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
- References:
- How to find the machine number nearest to an exact number -- N[] fails
- From: Alain Cochard <alain@geophysik.uni-muenchen.de>
- How to find the machine number nearest to an exact number -- N[] fails