2

(I suspect this question is a duplicate, but I didn't find a sufficiently similar question with an answer to it.)

I'm having trouble with comparisons of symbolic Reals that are equal, but which Mathematica has trouble recognising them as equal, apparently because it uses N to compare these (apparently after some simple symbolic manipulation). Typically these issues creep up in conditions like one (greatly simplified) below:

x > y /. {x -> Sqrt[(2 - Sqrt[3])^2 + (-3 + 2*Sqrt[3])^2], y -> 4 - 2*Sqrt[3]}

N::meprec: "Internal precision limit $MaxExtraPrecision = 50.` reached while evaluating -4+2\ Sqrt[3]+Sqrt[(2-Sqrt[3])^2+(-3+2 Sqrt[3])^2]."

Sqrt[(2 - Sqrt[3])^2 + (-3 + 2*Sqrt[3])^2] > 4 - 2*Sqrt[3]

This is in no way undocumented feature; it is discussed under Possible Issues section of $MaxExtraPrecision.

What I really want is Mathematica to try a bit harder on solving these (in)equalities numerically in a block of code. How do I accomplish this? As a workaround for specific problem above, this does work (while Simplify on inequality part doesn't):

x > y /. Simplify @ {x -> Sqrt[(2 - Sqrt[3])^2 + (-3 + 2*Sqrt[3])^2], y -> 4 - 2*Sqrt[3]}

False

I would be happy to see a solution that could wrap up around a block of code, and which could possibly work also on non-algebraics that FullSimplify can handle.

EDIT:

To clarify my question: I want first example to evaluate like the expression below, and in general Greater to evaluate for NumericQ argument list similarly inside a code block where I want this feature to be used:

  Simplify[Sqrt[(2 - Sqrt[3])^2 + (-3 + 2*Sqrt[3])^2]] > Simplify[4 - 2*Sqrt[3]]

False

kirma
  • 19,056
  • 1
  • 51
  • 93
  • @Artes My specific question here is how to make comparisons work this way; typically this stuff creeps up in my case in complicated equation systems Mathematica has generated itself. – kirma Sep 07 '13 at 09:26
  • @Artes N appears only in the error message, not in my code. This is exactly the inconvenience I'm trying to get rid of. – kirma Sep 07 '13 at 09:31
  • 2
    I do not know if following will work for you or not, but myself, I like to keep everything in MachinePrecision only to be able to compare things with Matlab and Fortran when I want to. You can force all computation to be in MachinePrecision by putting this at the top of the notebook or in your init.m $MinPrecision = $MachinePrecision; $MaxPrecision = $MachinePrecision; Now your example runs with no warnings: Mathematica graphics again, do not know if this is what you want, or this will work for what you are doing. – Nasser Sep 07 '13 at 09:31
  • @Nasser Not a bad workaround for applications where minuscule differences are known not to matter; but I guess it turns comparison into machine-precision numeric operation? I'd wish to make Less try a bit more in symbolic manner, and return unevaluated if it's indecisive, instead of providing a possibly wrong result through limited-precision numeric evaluation. – kirma Sep 07 '13 at 09:37
  • 2
    @kirma You should use PossibleZeroQ, with Method -> "ExactAlgebraics". This post Most efficient way to determine conclusively whether an algebraic number is zero discusses the issue. – Artes Sep 07 '13 at 09:38
  • @Artes I'm aware of PossibleZeroQ and "ExactAlgebraics". :)

    What I'd wish here would be a similar "knob" for functions like Less; apparently it isn't present directly. Something that'd pre-evaluate Simplify for all arguments of these comparison functions before performing them.

     – kirma Sep 07 '13 at 09:42
  • 1
    @Nasser regardless of its applicability to this question, that sort of comparison may be inherently unreliable. MATLAB and Fortran follow IEEE754, while Mathematica does not, so one cannot expect to obtain the same results from different systems. – Oleksandr R. Sep 08 '13 at 14:50
  • @Oleksandr R. For elementary operations in machine arithmetic, Mathematica is almost certainly in compliance with IEEE754, for the simple reason that it only runs on compliant processors/OSs. The same is likely to hold for those function evaluations for which library code is used, since (I believe) many such libraries are in compliance. Not sure whether there is a compliance issue with some special functions insofar as I'm not sure the standard addresses that area in entirety... – Daniel Lichtblau Sep 09 '13 at 15:36
  • @Oleksandr R. ...If you had in mind bignum arithmetic rather than machine precision, then Mathematica is certainly different from IEEE754. For purposes of zero comparison it is (far) less likely (than IEEE-based numeric comparison codes) to claim that mathematically equivalent numeric expressions are not equal... – Daniel Lichtblau Sep 09 '13 at 15:40
  • @Oleksandr R. ... Of course there are other aspects of arithmetic besides zero testing, and IEEE standard for fixed precision bignums is not, at this time, quite attainable within Mathematica (the issue being, as you probably know, that various nondefault rounding modes are not supported). I am not sure to what extent it might be emulated-- this will give an idea of how close I was able to come, but maybe others have better approaches. – Daniel Lichtblau Sep 09 '13 at 15:42
  • @DanielLichtblau the main area of departure as I see it is the tolerance applied on the comparison operators and the fact that Infinity and Indeterminate behave differently to IEEE infs and NaNs. While disabling the tolerance is possible in our own routines, for users it's hard to know where it's applied internally, hence where and when an error (versus strict IEEE arithmetic) of up to 7 ULPs might be able to accumulate. In fact, knowing what is done with a library call or using (potentially mixed precision) bignum arithmetic is difficult in itself. By saying it is not fully IEEE-compliant – Oleksandr R. Sep 10 '13 at 14:46
  • @DanielLichtblau I don't mean to criticize Mathematica's approach to numerics on its own terms, but just to note that I don't think it's realistic to expect bit-for-bit reproducibility of non-trivial results obtained in dedicated numerics codes. – Oleksandr R. Sep 10 '13 at 14:53
  • @Oleksandr R. Good point about the comparison operators. So conditionals based thereon will likely not behave identically to, say, Fortran library code. That said, I suspect that pretty much any math-related code depending on comparisons will be dealing with granularity issues-- I don't recall offhand whether IEEE has much to say about that. What I'm getting at is bit reproducibility of any one program by another is likely to be difficult. But yes, if "the other" is Mathematica, then the difficulty maybe does increase. – Daniel Lichtblau Sep 10 '13 at 21:38

1 Answers1

3

This solution uses trick to redefine internal functions presented in this earlier answer by Leonid Shifrin.

(* Code run as an argument of WithFullySimplifiedComparisons tries
   harder to symbolically resolve equality and inequality tests. *)

ClearAll[WithFullySimplifiedComparisons];
WithFullySimplifiedComparisons =
  Module[
   {withReplacedFunctions, fullSimplifiedTest, testMappings, 
    splitInequality, wrapper},

   (* Runs code with functions in { old -> new, ... } mappings list
      executing new function in place of old, but old definitions
      working inside new implementations. *)

   SetAttributes[withReplacedFunctions, HoldRest];
   withReplacedFunctions[mappings_List, code_] :=
    Module[{outer = True},
     Internal`InheritedBlock[Evaluate@mappings[[All, 1]],
      (Unprotect[#1];
         #1[args___] :=
          Block[{outer = False}, #2[args]] /; outer;
         Protect[#1];) & @@@ mappings;
      code]];

   (* Returns a version of comparison function fun that attempts to
      simplify numeric comparisons in domain dom with FullSimplify. *)

   fullySimplifiedTest[fun_, dom_] :=
    Module[{test},
     test[] := True;
     test[_] := True;
     test[a_, b_, r___] :=
      With[{bn = FullSimplify[b]},
        If[fun[0, FullSimplify[bn - a]],
          test[bn, r], False, test[a, bn] && test[bn, r]]] /; 
       NumericQ[a] && NumericQ[b] && (a | b) \[Element] dom;
     test[a_, b_, r___] := fun[a, b] && test[b, r];
     test];

   (* Mappings for normal comparison functions. *)

   testMappings = 
    #1 -> fullySimplifiedTest[##] & @@@ 
     Join[{#, Complexes} & /@ {Equal, Unequal},
       {#, Reals} & /@ {Greater, GreaterEqual, Less, LessEqual}];

   (* Conversion function from Inequality to chain of simple tests. *)

   splitInequality[_] := True;
   splitInequality[a_, test_, b_, r___] := test[a, b] && splitInequality[b, r];

   (* Actual wrapper to be used.
      splitInequality uses new comparison functions. *)

   SetAttributes[wrapper, HoldFirst];
   wrapper[code_] :=
     withReplacedFunctions[testMappings,
       withReplacedFunctions[{Inequality -> splitInequality},
         code]];

   wrapper];

Now earlier problematic tests can run purely symbolically:

WithFullySimplifiedComparisons[
  x > y /. {x -> Sqrt[(2 - Sqrt[3])^2 + (-3 + 2*Sqrt[3])^2], y -> 4 - 2*Sqrt[3]}]

False

Community
  • 1
kirma
  • 19,056
  • 1
  • 51
  • 93
  • Accepted as answer as it solves the example, and no other attempts appeared... I'm a bit suspicious of reliability of withReplacedFunctions, but it's an attempt better than nothing. – kirma Sep 09 '13 at 09:21