Services & ResourcesWolfram ForumsMathGroup Archive

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

Re: Floating-Point Computing

• To: mathgroup at smc.vnet.net
• Subject: [mg93927] Re: Floating-Point Computing
• From: Szabolcs Horvát <szhorvat at gmail.com>
• Date: Fri, 28 Nov 2008 07:11:18 -0500 (EST)
• Organization: University of Bergen
• References: <ggofnk$rsh$1@smc.vnet.net>

Antonio wrote:
> Dear List,
> 
> I came across this Sun article on floating point:
> http://developers.sun.com/solaris/articles/fp_errors.html
> 
> which calculates the following code:
> 
> main()
> {
>  register i;
>  double f=0.0;
>  double sqrt();
> 
>  for(i=0;i<=1000000000;i++)
>  {
>     f+=sqrt( (double) i);
>  }
>  printf("Finished %20.14f\n",f);
> }
> 
> whose correct answer using Euler-Maclaurin formula is:
> f=21081851083600.37596259382529338
> 
> Tried to implement in mathematica, just for fun, so that later I could
> play with precision, but the kernel has no more memory and shuts down.
> 
> N[Total[Sqrt[Range[1000000000]]]]
> 

If all those numbers are stored as 32-bit integers in a packed array, 
they would take up about 4*10^9/1024^3 ~= 3.7 GB space.  That is already 
too much for a 32-bit machine, and more than most present-day personal 
computers can handle.  Mapping Sqrt to the list unpacks it, eating up 
even more memory.  Creating a symbolic sum with Total and only 
numericising at the end doesn't help either.

The only reasonable approach here is to write code that is equivalent to 
the C version, i.e. doesn't pre-compute the list of numbers to be 
summed.  But that will run too slowly in Mathematica.  So Mathematica is 
not really the right tool for summing 10^9 numbers.

Two possible versions:

fun1[n_Integer] := Module[{f = 0.},
   Do[f += N@Sqrt[i], {i, n}];
   f
   ]

fun2 = Compile[{{n, _Integer}}, Module[{f = 0.},
    Do[f += N@Sqrt[i], {i, n}];
    f
    ]
   ]

fun2[10^9] can finish in a reasonable time, but it only works with 
machine precision numbers (because of Compile).

If you sum fewer numbers, use Total@Sqrt@N@Range[n] instead of 
N@Total@Sqrt@Range[n].

• Follow-Ups:
  • Re: Re: Floating-Point Computing
    • From: "Hobbs, Sylvia (DPH)" <Sylvia.Hobbs@state.ma.us>