Mapping algebraic numbers to machine precision numbers in Mathematica

The following email exchange deals with a problem in the N[] function of mathematica. The computer algebra system does not compute N[ Sqrt[ n]] correctly for large integers n. I argued with Wolfram research support about a in my opinion severe platform independent bug (or limitation) of the N[] function of Mathematica. The function N[] has problems with transforming exact algebraic numbers to numerical numbers. Here is an example
N[Sqrt[2^117 + 1]] - N[Sqrt[2^117 + 1], 100]
The result is -64. The behavior of N[] is quite random:
Table[N[Sqrt[2^k + 1]] - N[Sqrt[2^k + 1], 1000], {k, 1, 120}]
It is easy to bypass the problem by just avoiding N[] and giving a second argument. In comparison, Maple projects algebraic numbers correctly into machine precision numbers:
evalf(sqrt(2^117 + 1))
--------------------------------------------------------------------------
From: Oliver Knill 
Date: Wed, 17 Sep 2008 11:02:22 -0400 (EDT)
Subject: bug with N[] implementation in Mathematica
To: support@wolfram.com

Hi,
the following seems to point to a bug in the N[]  implementation

N[-10238998959279348 + 217141 Sqrt[2223470289078996563761]]

The error can be eliminated by giving a second argument to N[] like

N[-10238998959279348 + 217141 Sqrt[2223470289078996563761],20]

It seems to be platform independent (tested on Linux,Solaris and OS X)
appears in Mathematica 5 and Mathematica 6. Its not just a display
problem because also graphs display the wrong result
(ListPlot for example must use N[] internally and displays the wrong
values).

Oliver Knill
--------------------------------------------------------------------------

Date: Wed, 17 Sep 2008 23:30:03 UT
From: support@wolfram.com
To: knill@math.harvard.edu

Hello,

Thank you for the email.  Mathematica uses several different forms of
arithmetic. There is exact, reserved for whole numbers.  There is also
Machine Precision which is a little less than 16 digits of precision.  Then
there high precision which can be set in various ways including the N
function with a second argument.  Without setting the second argument
however it defaults to Machine Precision. I've attached a notebook
illustrating this idea and how to deal with it.

Sincerely,

Jason Martinez
Technical Support
Wolfram Research, Inc.
support@wolfram.com

Notebook

--------------------------------------------------------------------------
From: knill@math.harvard.edu (Oliver Knill)
Date: Wed, 17 Sep 2008 19:59:59 -0400 (EDT)
To: support@wolfram.com

Dear Jason,

I still think this is not supposed to be. Even with a default
Machine precision of 16 digits. Mathematica computes

N[Sqrt[111111111111111111111111111111111111111111111111111111]]

correctly. And I would expect that the projection from number
fields to the real numbers is done with the machine precision.
Any other computer algebra system does that correctly. And
Mathematica does this correctly also most of the time, even
with less precision:

217141 N[Sqrt[2223470289078996563761],3]
N[217141 Sqrt[2223470289078996563761],3]

give both the right result up to 3 digits.  Also

N[-10238998959279348 + 217141 Sqrt[2223470289078996563761],3]

is correct up to three digits and

N[-10238998959279348 + 217141 Sqrt[2223470289078996563761],16]

is correct up to 16 digits but

N[-10238998959279348 + 217141 Sqrt[2223470289078996563761]]

is completely off? I would expect the last two results to be the same up
to 15 digits.

Oliver Knill


P.S. the default N[] implementation seems really doing strange things:

(* Case where things go wrong *)
f[k_]:=N[k Sqrt[2223470289078996563761]-Floor[k Sqrt[2223470289078996563761]]]
ListPlot[Table[f[200*k],{k,1,1001}]]

(* Case where things are right, with 17 digits *)
f[k_]:=N[k Sqrt[222347028907899656]-Floor[k Sqrt[222347028907899656]]]
ListPlot[Table[f[200*k],{k,1,1001}]]

(* Case where things are even worse and x-floor(x) gets 8 *)
f[k_]:=N[k Sqrt[112223470289078996563761]-Floor[k Sqrt[112223470289078996563761]]]
ListPlot[Table[f[200*k],{k,1,1001}]]

(* Case where things are even worse and x-floor(x) gets -250 *)
s=112223470289078996563761888;
f[k_]:=N[k Sqrt[s]-Floor[k Sqrt[s]]]
ListPlot[Table[f[200*k],{k,1,1001}]]


-------------------------------------------------------------------------

From support@wolfram.com  Mon Sep 22 16:00:16 2008
Date: Mon, 22 Sep 2008 20:00:00 UT


Hello,

Thank you for the emails.  These all boil down to rounding errors, either
from using MachinePrecision or from converting from Exact numbers to
Machine Precision.  I've gone over this in detail in the attached
notebook.  Please let me know if you still have any questions about this.

Sincerely,

Jason Martinez
Technical Support
Wolfram Research, Inc.
support@wolfram.com

 Notebook 

-------------------------------------------------------------------------
From: knill@math.harvard.edu (Oliver Knill)
Date: Wed, 24 Sep 2008 17:16:12 -0400 (EDT)
To: support@wolfram.com

Hi Jason,

I understand the difference and how to fix the issues.
But you don't answer the question. My alert is about the
way mathematica does the conversion from exact to
machine precision.

I often make use of the exact arithmetic Mathematica use
and this is one of the strength of this computer algebra
system.

It is about the projection N[] from exact to machine precision
on algebraic numbers which is inconsistent and as I believe flawed.
The user expects the result of N[sqrt(b)] to be correct
up the machine precision for arbitrary integers.
This is obviously not the case. I think this is a bug which
needs to be fixed.

The example show that N[sqrt(b)]  is all over the place.
I don't see any mathematical reason why N[] should
sometimes be an integer, sometimes be a real and sometimes be 0.

Again, I don't have the problem to find a fix to

Table[N[Sqrt[2^k + 1]] - N[Sqrt[2^k + 1], 1000], {k, 1, 120}]

You mentioned the fix several times. My problem is that this situation
appears in Mathematica at all. And I'm even more bothered why there
is no explanation for the strange fact that the difference
can become -64 in

N[Sqrt[2^117 + 1]] - N[Sqrt[2^117 + 1], 100]

It would be nice (and not affect any exact calculation issues), if
N[sqrt(a)] would work properly for arbitrary large integers in
a future version of mathematica.

Oliver Knill

P.S. The bug or rather (limitation of Mathematica) leads to
interesting questions although, like what the following numbers are:

u=Union[Table[k (-1/64) (N[Sqrt[2^117 + k]]-N[Sqrt[2^117+k],100]),{k,1000}]]

As discussed previously, the sequence is independent of the second
argument of N[].

---------------------------------------------------------------------------

From support@wolfram.com  Fri Sep 26 16:30:06 2008
To: knill@math.harvard.edu

Hello,

Thank you for the email.

Thank you for taking the time to send us this suggestion for automatically
going to higher precision when dealing with large exact numbers. I have
forwarded your suggestion to the appropriate people in our development
group.

I have also included your contact information so that you can be notified
when this gets implemented.

Sincerely,

Jason Martinez
Technical Support
Wolfram Research, Inc.
support@wolfram.com

---------------------------------------------------------------------------