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>