24
$\begingroup$

In this answer, Oleksandr R. mentioned an undocumented function Internal`CompareNumeric and briefly explained its usage as follows:

Internal`CompareNumeric[prec, a, b] returns -1, 0, or 1 according to whether a is less, equal, or greater than b when compared at the precision of a or b (whichever is less) minus prec decimal digits of "tolerance". It is the fundamental operation underlying Less, Equal, Greater, LessEqual etc. for finite-precision numeric types.

…To be honest, I failed to understand this description. In my eyes, it seems to suggest that, Internal`CompareNumeric will compare

With[{minPre = Precision /@ {a, b} // Min}, N[{a, b}, minPre - prec]]

but it's not true, because

Internal`CompareNumeric[1, 1.1`2, 1.2`2]
(* -1 *)

while the output should be 0 in my understanding.

Can someone explain the usage of Internal`CompareNumeric in a more detailed way?

$\endgroup$
  • $\begingroup$ my guess the comparison is actually made directly on the binary representation and the "decimal digits" count only approximately translates to a specific number of base 2 digits. $\endgroup$ – george2079 Mar 20 '17 at 22:03
  • $\begingroup$ @george2079 Hmm… but Internal`CompareNumeric[1, 11`2, 12`2] still gives -1. $\endgroup$ – xzczd Mar 21 '17 at 3:13
  • 4
    $\begingroup$ FYI, the split point of -1 and 0 is around7.3456427360042 $\endgroup$ – happy fish Mar 26 '17 at 9:10
  • $\begingroup$ I was discussing this with yode the other day. Internal`CompareNumeric[n,a,b] basically answers the question, "how do a and b related ignoring n points of precision". So if you have two numbers and drop $MachinePrecision points of precision it's basically comparing them at the integer level. e.g. Internal`CompareNumeric[IntegerPart@$MachinePrecision, 2.1, 2] (because my $MachinePrecision isn't an int) gives 0, but Internal`CompareNumeric[IntegerPart@$MachinePrecision - 1, 2.1, 2] gives 1. $\endgroup$ – b3m2a1 Mar 29 '17 at 6:13
  • 1
    $\begingroup$ Additional question: Why does Internal`CompareNumeric give an error message when prec is set smaller than $MinPrecision? $\endgroup$ – mmeent Mar 29 '17 at 14:21
7
$\begingroup$

This is just to highlight another aspect of the other answer given. I would have expected a linear progression from the definition. I.e., if Precision[n]=p then we'd expect the switch point for CompareNumeric[x, n, n+epsilon] to be p. And this holds once you get past $MachinePrecision. Probably happens when Mathematica has to go from 64 bit numbers to whatever its internal infinite-precision numeric is. We can see this here:

pr[x_] :=
  SelectFirst[
   Range[1, 100, .1],
   Internal`CompareNumeric[#, SetPrecision[1.1`, x], 1.2`100] == 0 &
   ];

Show[
 Plot[x, {x, 0, 25}, PlotStyle -> Gray],
 {#, pr@#} & /@ Range[1, 25, .01] // ListPlot
 ]

blah

I looked at the part for n < $MachinePrecision and thought: "oh, a log-type function (probably log2)". But no simple scaling on Log2 worked for me.

Show[ListLinePlot@Table[Table[x*N@Log2[n], {n, 15}], {x, 5}], 
 ListPlot[{#, pr@#} & /@ Range[1, 15, .01]]]

log2

So maybe someone else can tackle what sort of pre-form that is. Once we have that we can probably piece back how Mathematica is interpreting things for its sub-machine precision numbers.

Update

As found by happy fish, $MachinePrecision - $MachinePrecision/x fits this really well:

Plot[{
  pr[x],
  $MachinePrecision - $MachinePrecision/x
  },
 {x, 1, $MachinePrecision}
 ]

fitted

Noting that $MachinePrecision is 53Log10[2] by simple log rearrangements we have Log10[2^(53*(1 - 1/x))]. Since $MachinePrecision comes from how a 64 bit number is stored (there are 53 bits of precision), this could be saying it's using 53*(1 - 1/x) bits of precision for these numbers with Precision specified below $MachinePrecision. Why that would be, I have no idea. But it's a possibility.

Update 2

Shadowray points out this is even odder when you have fewer points of precision available before they would be equivalent

For example, make pr multivariate:

pr[x_, y_: 1] :=
  SelectFirst[
   Range[1, 100, .1],
   With[{
     a = SetPrecision[1. + (10^(-y)), x],
     b = SetPrecision[1. + (2*10^(-y)), 100]
     },
    Internal`CompareNumeric[#, a, b] == 0 &
    ]
   ];

Now plotting this for different amounts of reduction:

Show[
 Plot[x, {x, 0, 25}, PlotStyle -> Gray],
 Table[
   {#, pr[#, n]} & /@ Range[1, 25, .05],
   {n, 17}
   ] // ListPlot
 ]

changingprec

Note that the first 6 curves actually continue. They just jump up to y=x after $MachinePrecision (as I would naively have predicted).

If we use a near-continuous change:

Show[
 Plot[x, {x, 0, 25}, PlotStyle -> Gray],
 Table[
   {#, pr[#, n]} & /@ Range[1, 25, .05],
   {n, Range[1, 17, .1]}
   ] // ListPlot
 ]

continuous

The pre-discontinuity part can still be easily formulated:

Show[
 Plot[x, {x, 0, 25}, PlotStyle -> Gray],
 Plot[
  Evaluate@Table[$MachinePrecision (1 - n/x),
    {n, 16}],
  {x, 1, $MachinePrecision}],
 Table[
   {#, pr[#, n]} & /@ Range[1, 25, .1],
   {n, 17}
   ] // ListPlot
 ]

shifting

What's happening at the discontinuities is somewhat less clear. Note that the oddity at 17 is from the fact that 10^-17 is smaller than 10^-$MachinePrecision and so the system just treats the addition as doing nothing. The jumps need more explanation though.

In any case, if we return to the storage argument, now we have 53 n/x bits being dropped, assuming the interpreting the pre-jump part of the curves as showing Log10[2^(53 (1 - n/x))] is valid.

$\endgroup$
  • $\begingroup$ @xzczd made the pr much more sensitive. Shows the real discontinuity at $MachinePrecision now. $\endgroup$ – b3m2a1 Mar 29 '17 at 8:06
  • $\begingroup$ I got similar plots. Although I used FindFormula to get very close fits, but those forms and coefficients do not make sense... $\endgroup$ – happy fish Mar 29 '17 at 8:13
  • $\begingroup$ @happyfish It honestly looks like a -1/r thing (but shifted vertically). Simple (16+-1/r) is way too steep though. $\endgroup$ – b3m2a1 Mar 29 '17 at 8:21
  • $\begingroup$ $MachinePrecision -17.217894/x no idea $\endgroup$ – happy fish Mar 29 '17 at 8:27
  • $\begingroup$ @happyfish Wow. That fits perfectly. Who knows if it's right. Maybe someone brighter than myself will see that and be able to explain why it works. Maybe it's supposed to be $MachinePrecision - $MachinePrecision / x? That's just as good a fit, I think. $\endgroup$ – b3m2a1 Mar 29 '17 at 8:29
5
$\begingroup$

After some experiments I found some of the "rules" with which Internal`CompareNumeric[c, a, b] operates. If you find any exceptions please notify me. Lets denote minpr = Min[Precision[a],Precision[b]]. If c is greater than pr[minpr] then Internal`CompareNumeric[c, a, b] kind of throw away anything after the point. The function pr depends on a and b, the only symmetry I was able to detect is that if minpr>15, pr[minpr]=minpr (maybe connected to MachinePrecision?). The following figure is a possible pr[].

enter image description here

$\endgroup$
  • 1
    $\begingroup$ What is pr[]? $\endgroup$ – Michael E2 Mar 29 '17 at 11:28
  • $\begingroup$ @MichaelE2 An assumed function relationship between minpr and the threshold? $\endgroup$ – xzczd Mar 30 '17 at 9:25
  • $\begingroup$ @xzczd There's no definition of pr[] in your OP, nor code for it. I notice MB1965 has a code for it (but no definition other than the code). Is it the same as used here? Since the ordering of answers changes over time (this was the top answer when I read it), it would be helpful to point to the definition of pr[] or give the code, if it's different. In the answer, "The following figure is a possible pr[]" suggests to me that there might be various ways to define pr[] and the plot is a possible one. $\endgroup$ – Michael E2 Mar 30 '17 at 11:03
  • $\begingroup$ @MichaelE2 In my understanding, pr[] is more of a guess in this answer,. (BTW this is indeed the earlest answer. ) Anyway let's wait for Most Perfect Username's clarification :) $\endgroup$ – xzczd Mar 30 '17 at 11:10
  • $\begingroup$ @xzczd A guess at what? There's no def. But yes, I think the OP can easily clarify. $\endgroup$ – Michael E2 Mar 30 '17 at 11:13
5
+500
$\begingroup$

In this answer I will present an overall simulation of Internal`CompareNumeric's behavior. Throughout this answer, there are four variables: Internal`CompareNumeric[tol, a, b], first=Min[Abs[a],Abs[b]] diff=Abs[a-b] and prec is the minimum precision of the two specified. Here's my findings so far

  • For exact arithmetic, Internal`CompareNumeric gives correct answer regardless of tol.
  • For machine precision integers (1., 2.), the split point is always $MachinePrecision.
  • The order and sign of a, b does not matter, the minimum precision is used for calculation(this can be implied from my definition above).
  • The smallest prec that makes the two numbers comparable is $-\log_{10}(\frac{diff}{first+diff})$ (this number comes from the definition of Precision), i.e. when $prec>-\log_{10}(\frac{diff}{first+diff})$, there exist a tol, such that the function returns -1 or 1; when prec is less than this number, the function returns 0 for all tol settings.
  • When $prec\in [-\log_{10}(\frac{diff}{first+diff}),\$MachinePrecision]$, there exist a split point $y$, such that when $tol<y$, the function returns -1 or 1, when $tol>y$, the function returns 0. The analytical formula with first, diff, prec as parameters is $$y(first,diff,prec)=\$MachinePrecision (1+\frac{1}{prec}\log_{10}(\frac{diff}{first+diff}))$$
  • For prec>$MachinePrecision, the split point of tol is exactly the same as prec.
  • For prec>$MachinePrecision, the formula above is incorrect. Some experiments show that when the ratio of first and diff is big enough, the comparison is converted to exact arithmetic(independent of tol). There should be two such points when diff/first is small and big, I haven't worked on it yet. A quick fixed of the formula is y[prec,first,diff]=prec-If[#>8,#,0]&[-Log10[diff/(first+diff)]], I will make careful research later when I have time.

The following is the code version of the findings above:

fishCompareNumeric[tol_, a_, b_] := 
 With[{first = Min[Abs@a, Abs@b], diff = Abs[a - b], prec = Precision@{a, b}},
  Module[{split},
   Catch[
    split = Piecewise[{
       {Throw@Sign[a - b], prec === Infinity},
       {$MachinePrecision, prec === MachinePrecision},
       {Throw@0, prec < -Quiet@Log10[diff/(first + diff)]},
       {$MachinePrecision (1 + Log10[diff/(first + diff)]/prec), -Log10[diff/(
           first + diff)] <= prec <= $MachinePrecision},
       {prec, prec > $MachinePrecision}
       }];
    Piecewise[{
      {Sign[a - b], tol < split}
      }]]]]

Here is a plot of my binary search for approximate values, with parameters(OP's example) first=1.1, diff=0.1, horizontal axis being prec and vertical axis being split point (of tol).

enter image description here

You may notice the two vertical lines at around 1.07 and 16, they are introduced by the sudden value changes, at the two points in my formula.

Here is a plot of my exact formula with same parameters enter image description here

Code for binary search approximation:

splitPointBinary[{{val1_, val2_}, {low_, up_}}] := 
 With[{mid = Internal`CompareNumeric[(low + up)/2, val1, val2]}, 
  If[Internal`CompareNumeric[low, val1, val2] == 
    mid, {{val1, val2}, {(low + up)/2, up}}, {{val1, 
     val2}, {low, (low + up)/2}}]]
splitPointApprox[first_, diff_, prec_] := 
 FixedPoint[
   splitPointBinary, {{SetPrecision[first, prec], 
     SetPrecision[first + diff, prec]}, {0.1, 1000}}][[2, 1]]

Code for my formula:

splitPoint[first_, diff_, prec_] := 
 Piecewise[{{$MachinePrecision (1 + Log10[diff/(first + diff)]/prec), 
    -Log10[diff/(first + diff)] < prec < $MachinePrecision}, {prec, 
    prec > $MachinePrecision}}]

For OP's example,

splitPointApprox[1.1, 0.1, 2]
(*7.34564*)
splitPoint[1.1, 0.1, 2]
(*7.34564*)

Another one,

splitPointApprox[2.2, 0.5, 6]
(*14.0071*)
splitPoint[2.2, 0.5, 6]
(*14.0071*)
$\endgroup$
  • $\begingroup$ This looks strange for splitPointApprox[1.1, 0.0000000000001, 2] $\endgroup$ – Shadowray Mar 29 '17 at 11:14
  • $\begingroup$ @xzczd answer updated. $\endgroup$ – happy fish Mar 31 '17 at 6:26
  • $\begingroup$ I took the liberty to add a function interpreted from your findings to your answer, feel free to edit or roll back if you don't like it :) . $\endgroup$ – xzczd Apr 2 '17 at 14:03
  • $\begingroup$ @xzczd I like it :) $\endgroup$ – happy fish Apr 2 '17 at 14:41
  • $\begingroup$ Since I haven't think out a systematic enough way to judge the answer, I'll leave it unaccepted for a while, but this is undoubtfully the most plausible answer I received so far, +51 :) . $\endgroup$ – xzczd Apr 3 '17 at 3:21
2
$\begingroup$

Let's consider a slightly different approach than above. For Internal`CompareNumeric[prec, a, a+ε], lets fix prec and a and find the threshold value for ε for which this switches from 0 to -1. For this purpose I define the following function

tf[p_?NumericQ, y_?NumericQ, z_?NumericQ] :=
   Internal`CompareNumeric[
     p,
     SetPrecision[1, y],
     SetPrecision[1, y] + 10^-SetPrecision[z, 100]
]

That is we take a = SetPrecision[1, y] and ε = 10^-SetPrecision[p, 100]. (The SetPrecision in the definition of ε is there to avoid that the precision of z in the input can affect the outcome in anyway.)

We can then find the threshold value using simple bisection. For example abusing FindRoot:

tf2[p_?NumericQ, y_?NumericQ] := (z /. 
  FindRoot[tf[p, y, z] == -1/2, {z, -10, p + y + 1}, 
  WorkingPrecision -> 100, MaxIterations -> 20])

Plotting the threshold value as a function of y for values of p we get enter image description here

The top curve (p=0) is simply given by z = y. More generally the behavior of tf2 is modelled precisely by

Max[
 Piecewise[
  {
   {y - p, y > MP},
   {((MP - p)/MP) y, y < MP}
  }
 ],
 MP/2 UnitStep[y - Max[MP, p]]
]

where

MP=$MachinePrecision

What does this tell us about the behavior of Internal`CompareNumeric[prec,a,b]? Well, if the precision of a and b is set higher than machine precision and prec is not too large, the behavior is exactly as described by Oleksandr R.

It becomes a bit weirder when prec is large (compared to the precision of a and b), in which case there is an apparent upperbound of on the threshold difference of ε = 10^(-MP/2).

I don't quite understand the behavior of the threshold when the precision of a and b is smaller than machine precision. Maybe somebody else can shed some light.

Update:

Based on the above we can venture to propose a mock version of Internal`CompareNumeric[prec, a, b]:

MockCompareNumeric[p_,a_,b_]:=With[
 {ip = Min[Precision[a], Precision[b]]},
  With[
   {ep = Which[
    ip === MachinePrecision, MP-p
    ip >= Max[MP, p], Max[ip - p, MP/2],
    True, Max[(1 - p/MP) ip, 0]
    ]
   },
   If[N[ Abs[a - b]/Max[Abs[a], Abs[b]]] < N[10^-ep], 0, Sign[a-b]] 
 ]
];

I have tested the behavior of this mock function with the actual Internal`CompareNumeric over a large range of prec, a b, and they match. (edit 3.4.2017: adjusted the mock function to get the correct behavior for MachinePrecision input.)

Addendum

Ofcourse, the actual Internal`CompareNumeric would not work like this. Instead it probably uses some low level functions working directly on the internal representation of the arbitrary precision numbers. (It certainly is a lot faster) However, the mock function gives us a fairly good idea of the behavior of Internal`CompareNumeric[prec,a,b].

In short it determines some effective precision (ep in the mock function) from prec and the minimal precision of the inputs a and b and then determines whether a > b within a precision tolerance set by ep.

The difference between the behavior of ep dependent on whether the input precision is larger or smaller than machine precision probably arises due to differences in implementation leveraging the fact that for precision smaller than machine precision, Mathematica represents its arbitrary precision numbers simply as a pair of machine numbers representing the value and the precision of the number, in which case it can leverage various machine operations to do the comparison. I suspect that similarly the lower bounds on ep arise some checks that can be down efficiently using low level functions. The scaling of ep when ip < MP presumably is some design choice relating to the desired functionality, possibly related with wanting the effective precision to be larger than zero for reasonable values of prec.

$\endgroup$
  • $\begingroup$ Have you seen this, by any chance? $\endgroup$ – J. M. is away Mar 29 '17 at 13:25
  • $\begingroup$ @ J.M. No I hadn't. I am also not sure why you are pointing to it. Could you explain? $\endgroup$ – mmeent Mar 29 '17 at 13:55
  • $\begingroup$ You mentioned wanting to find the threshold value where two numbers are no longer the same, so I think a function that would give the "next" floating point number might suit your needs. $\endgroup$ – J. M. is away Mar 29 '17 at 13:58
  • $\begingroup$ @ J.M. Ah, I see. Well, in this case we are not interested in finding the threshold where the two numbers are no longer considered the "same" by Internal`CompareNumeric which is a slightly different question, than find the threshold where the internal representation becomes the same. (Arbitrary precision numbers in Mathematica can have many more digits in their internal representation than their actual precision value.) $\endgroup$ – mmeent Mar 29 '17 at 14:06
  • $\begingroup$ What if one evaluates MockCompareNumeric[5, -0.9252872333493522561`2., 0.9289648686205458361`3.] ? $\endgroup$ – Shadowray Mar 30 '17 at 9:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.