Services & ResourcesWolfram ForumsMathGroup Archive

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

Re: Re: Floating-Point Computing

• To: mathgroup at smc.vnet.net
• Subject: [mg93998] Re: [mg93927] Re: Floating-Point Computing
• From: danl at wolfram.com
• Date: Mon, 1 Dec 2008 07:02:19 -0500 (EST)
• 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].

These two approaches can be hybridized to give quite good speed. The one
caveat is that, when totalling a list, the semantics are not defined as to
order of summing. Though in actual practice I believe Total sums
first-to-last, which is in accord with the semantics of the original
question.

Anyway, here is code that sums square roots to n, using blocks of size m.

cfun = Compile[{{n, _Integer}, {m, _Integer}},
   Module[{tot = 0., flr = Floor[n/m]},
    Do[
     tot += Total[Sqrt[j*m + Range[m]]];
     , {j, 0, flr - 1}];
    tot += Total[Sqrt[Range[flr*m + 1, n]]];
    tot]];

It's about four times faster than what I had for summing one at a time
(essentially the same as fun2 above).

In[129]:= Timing[cfun[10^7, 10^5]]
Out[129]= {0.841, 2.10819\[Times]10^10}

In[131]:= Timing[cfun[10^8, 10^7]]
Out[131]= {10.005, 6.66667\[Times]10^11}

In[130]:= Timing[cfun[10^9, 10^7]]
Out[130]= {104.44, 2.10819\[Times]10^13}

My general experience with Compile suggest that this speed is perhaps
within a factor of 2-4 of the C code above. The test below also suggest
this, since it spends most of its time in C code (though it may be a
factor of 2 slower than need be, since it makes two passes over the list).

In[135]:= Timing[Total[Sqrt[Range[1., 10^7]]]]

Out[135]= {0.431, 2.10819\[Times]10^10}

Daniel Lichtblau
Wolfram Research

• Follow-Ups:
  • Re: Re: Re: Floating-Point Computing
    • From: Daniel Lichtblau <danl@wolfram.com>