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

Re: Numerical accuracy/precision - this is a bug or a feature?

• To: mathgroup at smc.vnet.net
• Subject: [mg120083] Re: Numerical accuracy/precision - this is a bug or a feature?
• From: Richard Fateman <fateman at cs.berkeley.edu>
• Date: Fri, 8 Jul 2011 04:50:39 -0400 (EDT)
• References: <ius5op$2g7$1@smc.vnet.net> <ius7hm$358$1@smc.vnet.net> <iuuko3$eqi$1@smc.vnet.net> <iv45mo$f3n$1@smc.vnet.net>

On 7/7/2011 4:35 AM, Daniel Lichtblau wrote:
> On Jul 5, 4:15 am, Richard Fateman<fate... at cs.berkeley.edu>  wrote:
>> In brief:  yes to both.
>>
>> On 7/4/2011 4:18 AM, Kevin J. McCann wrote:
>> ...
>>
>>
>>
>>>> 1. N[2.0,20] should give 2 with accuracy/precision/whatsever about 20
>>>> decimal digits, i.e. 2.00000000000000000000
>>
>> WRI will defend this as a feature.
>>
>> You thought that the semantics of N[] were the same as SetPrecision
>> e.g.
>>
>> N[SetPrecision[2.0, 20]*N[Sqrt[2], 20], 20]
>>
>> works as you expected.
>>
>> So from your perspective, and from the perspective of anyone else who
>> thinks along the same lines, it is a bug.  I would prefer to call it
>> a design error.
>>
>> 2.0  (indeed, any floating point number) is a perfectly respectable way
>> of denoting a number that can be expressed in higher precision by adding
>> binary zeros. WRI doesn't agree.
>>
>> RJF
>
> I don't think this holds up for all floats in the sense you seem to
> indicate. Take 2.1 for example.

My contention is that 2.1, written that way, is a perfectly respectable 
way of denoting a number that can be expressed in higher precision by 
adding binary zeros.  2.1 is a double-float in most systems today.
What is that number?  In decimal, to higher precision, it looks like
2.10000000000000008881784197001....

A program that decodes a float into its component fraction, exponent, 
and sign gives this information for the fraction:
4728779608739021  or in binary
10000110011001100110011001100110011001100110011001101

That is the number to which we add binary zeros.

Let us call this number P.





>We can (and will) use SetPrecision to
> pad with binary zeroes, just as you use it above on 2.0.
>
> In[606]:= N[SetPrecision[2.1, 20]*N[Sqrt[2], 20], 20]
>
> Out[606]= 2.9698484809834997281
>
> But this is not an accurate representation of 21/10 * Sqrt[2] to 20
> places.


Just because you are doing arithmetic on numbers represented to 20 
digits does not mean that the result of the arithmetic will be right to 
20 digits, however, in the single operation of multiplication, you 
should not lose any digits.  What's going on here?

You simply computed two different quantities.  Let S= sqrt(2) to 20 
decimal places.  In line 606 you multiplied P*S and rounded it to 20 
digits  [ actually you did something to the binary reps, but it doesn't
matter here].
In line  607, below, you multiplied 21/10 by S.   Since P is not equal 
to 21/10, why should the results be the same?  Oh, Mathematica thinks 
that P ==21/10,  but that is merely the terribly broken notion of 
equality that Mathematica uses. We used to call that "close enough for 
government work".   But it is terribly broken because == fails to 
satisfy the axioms of "equivalence relations" such as a==b  and b==c 
implies a==c.


>
> In[607]:= N[21/10*Sqrt[2], 20]
> Out[607]= 2.9698484809834996025
>
> They are off by arould a machine epsilon. No surprise here, I think.
>
> In[608]:= % - %%
> Out[608]= -1.256*10^-16
>
> I also think that this sort of numeric discrepancy is inevitable* in
> any program that uses decimal input and binary representation.

Absent any interior calculations, I would tend to agree.  But there is 
something else going on. I think that Mathematica has embedded in it 
several design decisions that elevate the (typical) 
one-half-unit-in-last-place error in binary-to-decimal and 
decimal-to-binary conversion that happens in input or output,  into a 
global slush bucket that infects and propagates into all interior 
computations with "software floats".
>
> Daniel Lichtblau
> Wolfram Research
>
> *Perhaps more inevitable than these threads. Now that's pretty
> inevitable.
>
maybe so.