[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].
Prev by Date:
**Re: Floating-Point Computing**
Next by Date:
**Re: Package Problem**
Previous by thread:
**Re: Floating-Point Computing**
Next by thread:
**Re: Re: Floating-Point Computing**
| |