MathGroup Archive 2010

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

Search the Archive

Re: Re: Numerical Problem

  • To: mathgroup at
  • Subject: [mg107115] Re: [mg107044] Re: Numerical Problem
  • From: Daniel Lichtblau <danl at>
  • Date: Wed, 3 Feb 2010 06:09:49 -0500 (EST)
  • References: <> <> <hk17nr$oja$> <>

Richard Fateman wrote:
> Thanks, Tony, for providing a nice example of how arithmetic matters.
> Your  3-line program got pretty plausible results in version 5.0
> Clear[epcm,fepcm]
> epcm[h_,y_List,f_List]:=y+0.5 h (f.y+f.(y+h f.y))
> fepcm[h_,t0_,tmax_,y0_List,f_List]:=NestList[epcm[h,#,f]&,y0,10]
> app1=fepcm[.1,0.,1,{2.,-1.},{{998,1998},{-999,-1999}}]
> But apparently not in version 7, where DanL suggests that Intel's
> numerical library may have sent it off in the wrong direction.
> Let's try using Mathematica's significance arithmetic, and do everything
> "exactly" except for one number,  0.5.  Let us use N[1/2,40].  That is,
> 40 decimal digits or so.  Observe that 1/2 is exactly representable in 
> binary, so that SetPrecision[ [N[1/2,40], 100]  is
> 0.50000000000000000000000000000000000000000000000000000000000000000000\
> 00000000000000000000000000000000
> We change the program then to this, below, using software high-precision 
> and therefore significance arithmetic.
> Clear[epcms,fepcms] (*software floats *)
> epcms[h_,y_List,f_List]:= y+N[1/2,40]*h (f.y+f.(y+h f.y))
> fepcms[h_,t0_,tmax_,y0_List,f_List]:=NestList[epcms[h,#,f]&,y0,10]
> app1=fepcms[1/10,0,1,{2,-1},{{998,1998},{-999,-1999}}]
> If you look at the trailing few numbers using InputFormat you find that 
> for some of them, the numbers that approximate the exactly correct 
> answer for about 18 decimal digits are believe wrong after 3 or 4 
> decimal digits. And the last two numbers, which print as 0.x10^1 are,
> how can we put this delicately... wrong?
> This is using 40 digits -- an extra 22 decimal digits compared to 
> machine precision. But machine precision on Mathematica 5.2 apparently 
> gets the numbers pretty much right.  I am curious to see how version 7 
> can be using Intel's library to get different answers.
> NOte that if you carry even MORE digits, the problem eventually goes 
> away,  but for that to work you would have to be looking for this 
> problem to appear, and you might not notice it if the results of this 
> program were fed in to yet another program.

Somebody's library has a big effect on it. As for why a library code 
would be invoked at all, I think this is how Dot with machine numbers is 
handled in the Mathematica kernel (that is, by using e.g. MKL BLAS).

Anyway, I can now emulate the "good" behavior in a way that explains 
what is happening. It is, as RJF suggests, the use of extended precision 
arithmetic in intermediate computations that makes the difference.

In brief, some platforms support, and some machine double libraries use, 
higher precision register arithmetic. Specifically, machine doubles 
support mantissas of 53 bits, and extended precision registers use 64, 
via 80 bit registers (the other 16 bits are for sign and exponent). But 
the results are eventually stored in 64 bit memory, hence the mantissas 
get altered (and on't ask what happens when the exponents don't fit 
their 11 bit allotments...).

It took me a while to get the emulation correct because I was truncating 
  the 64 bit results, and zero-padding the emulated "true" doubles (the 
things with 53 bits). The latter was appropriate, but the former really 
requires rounding, so I was getting junk in the last bits of the 
doubles, and the ill conditioning magnified it, and...

Anyway, here it is.

SetAttributes[emulateFloat, Listable];
SetAttributes[emulateExtended, Listable];

emulateFloat[rr_Real] :=
  Module[{rf = RealDigits[SetPrecision[rr, $MachinePrecision], 2]},
   rf[[1]] = PadRight[rf[[1]], 53]; (*probably no longer necessary*)
   Sign[rr]*N[FromDigits[rf, 2], $MachinePrecision]]

emulateExtended[rr_Real] :=
  Module[{rf = RealDigits[rr, 2]}, rf[[1]] = PadRight[rf[[1]], 64];
   Sign[rr]*N[FromDigits[rf, 2], 64*Log[10., 2]]]

The iteration, with extended precision emulated, involves bumping floats 
to 64-bit-mantissa bignums, doing the operation, and forcing 
$MachinePrecision digits (converted to bits, that is, 53) to be retained.

epcm[h_, y_List, f_List] := y + 1/2*h*(f.y + f.(y + h*f.y))
app[0] = {2., -1.};
mat = {{998, 1998}, {-999, -1999}};

    app[j] =
     emulateFloat[epcm[1/10, emulateExtended[app[j - 1]], mat]]]], {j,
   1, 16}]

Do[Print[app[j] =
    emulateFloat[epcm[1/10, emulateExtended[app[j - 1]], mat]]], {j, 1,


One might ask whether this "magic" is really worth the bother. Maybe it 
is but I do not know that it is all that valuable for this particular 
problem. As Steve Schochet pointed out in a different response, changing 
the input by even an ULP will cause the instability to manifest in the 
iterations, and that holds true regardless of whether or not extended 
precision registers are used. Numeric ill-conditioning trumps clever 
arithmetic every time.

Daniel Lichtblau
Wolfram Research

  • Prev by Date: Re: position of sequence of numbers in list
  • Next by Date: Re: Re: Numerical Problem
  • Previous by thread: Re: Re: Numerical Problem
  • Next by thread: Re: Re: Numerical Problem