MathGroup Archive 2009

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: How to find the machine number nearest to an exact

  • To: mathgroup at smc.vnet.net
  • Subject: [mg103898] Re: [mg103869] How to find the machine number nearest to an exact
  • From: danl at wolfram.com
  • Date: Sun, 11 Oct 2009 08:06:33 -0400 (EDT)
  • References: <200910101107.HAA12596@smc.vnet.net>

> Hello.
>
>   6.0 for Linux x86 (32-bit) (February 7, 2008)
>
> I had always thought that N[exact number] gives the machine number
> closest to that exact number.  But consider the following exact number
>
>   ex=14539079527638557142952619724749887685325578333570\
>      9865657255467349279813823794598008431302162071157\
>     /                                                 \
>      8000000000000000000000000000000000000000000000000\
>     0000000000000000000000000000000000000000000000000;
>
> find its machine number representation:
>
>   num=N[ex];
>
> and determine the corresponding rounding error:
>
>   err=N[Abs[ex-  SetPrecision[num,Infinity]]  ]
>
> For me this gives 3.59393 10^-16, which is larger than Eps =
> $MachineEpsilon = 2^-52 ~ 2.22 10^-16.  But, between 1 and 2, the
> spacing between machine numbers is Eps.  So the rounding error should
> be smaller than or equal to Eps/2 ~ 1.11*10^-16 if the rounding is to
> the nearest machine number.
>
> To investigate further, consider the binary expansions of ex:
>
>   binex=RealDigits[ex,2,55]
>
> which gives:
>
>   Out[4]= {{1, 1, 1, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0,
> 0,
>
>   >     0, 1, 1, 1, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 0, 0, 1, 0, 0, 1, 1, 0,
> 1, 0,
>
>   >     0, 1, 0, 1, 1, 0, 0, 1, 0}, 1}
>                           ^
>                           ^-53rd digit
>
> Thus I would expect that the nearest binary number should end up
> either like (rounding below):
>
>   >     0, 1, 0, 1, 1, 0, 0}, 1}
>
> or like (rounding above):
>
>   >     0, 1, 0, 1, 1, 0, 1}, 1}
>
> whereas the actual binary expansion is given by:
>
>   binum=RealDigits[num,2]
>
> and ends as
>
>   >     0, 1, 0, 1, 0, 1, 1}, 1}
>
> In order to check, let us artificially construct those 2 other
> rounding numbers:
>
>   binum2=binum;
>   binum2[[1,51]]=1;
>   binum2[[1,52]]=0;
>   binum2[[1,53]]=0;
>
>   binum3=binum2;
>   binum3[[1,53]]=1;
>
> and compute again the rounding errors:
>
>   err2= N[ Abs[ ex-FromDigits[binum2,2] ]]
>   err3= N[ Abs[ ex-FromDigits[binum3,2] ]]
>
> This gives 1.37348 10^-16 for err2, which is this time less than Eps
> and 8.46964 10^-17 for err3, which is less then Eps/2, as expected.
>
> So, in this case, not only does N[] not give the nearest machine
> number, but it "jumps" a machine number.
>
> Is this a normal behavior?  Is there a flaw in my reasoning above? Is
> there an easy way to automatically generate the nearest machine
> number?
>
> Many thanks in advance.
> Alain
>
> PS: all outputs are identical if num is evaluated as
> SetPrecision[ex,MachinePrecision] instead of by N[ex].

The behavior strikes me as correct. It follows from rounding error. First
recognize that N applies itself to arguments working from the inside out.
So it is effectively doing N[numerator]/N[denominator].

N[numerator] turns out to be low, though I rather suspect it indeed gives
the closest representable double. It turns out that 1/N[denominator] is
also slightly low. I do not know offhand how machine double inversion is
done, but I suspect it is giving a result in accordance with IEEE specs.
The product has an error maybe 3 times what you expected. But notice that
representing operands as n+d_n (as in numerator plus delta_numerator) and
d+d_d (where I use d actually for 1/denominator), the product is n*d +
n*d_d + d*d_n up to first order error terms. If d_n and d_d are each
ballpark of half machine epsilon, the multipliers (n and d that is) can
bring this total to around 1.5*epsilon.

Here is an interval approach that shows how large we might fairly stray.

ex = 14539079527638557142952619724749887685325578333570986565725546734\
9279813823794598008431302162071157/
   8000000000000000000000000000000000000000000000000000000000000000000\
0000000000000000000000000000000;

nex = Numerator[ex];
dex = Denominator[ex];

In[146]:= InputForm[in =
  Interval[{N[nex, $MachinePrecision], N[nex, $MachinePrecision]}]]

Out[146]//InputForm=
Interval[
 {1.4539079527638554107532475622048214373667`15.954589770191005*^98,
  1.4539079527638560178372763827451560996985`15.954589770191005*^98}]

In[159]:= InputForm[id =
  Interval[{1/N[dex, $MachinePrecision], 1/N[dex, $MachinePrecision]}]]

Out[159]//InputForm=
Interval[
 {1.2499999999999998375717224117984497104`15.954589770191005*^-98,
  1.2500000000000001624282775882015502895`15.954589770191005*^-98}]

In[148]:= InputForm[in*id]

Out[148]//InputForm=
Interval[{1.81738494095481858319658506347207658741`15.653559774527023,
  1.81738494095482070254156986771549394137`15.653559774527023}]

As for obtaining the actual nearest double, you could do something like

In[62]:= InputForm[i2 = N[N[ex, 60]]]

Out[62]//InputForm=
1.8173849409548197


Daniel Lichtblau
Wolfram Research






  • Prev by Date: Re: Re: How to find which variable caused the trigger in Manipulate[]
  • Next by Date: Re: Plot with axis with the same scaling factor
  • Previous by thread: How to find the machine number nearest to an exact number -- N[] fails
  • Next by thread: update of Synergetics coordinates Notebook