Re: Floating-Point Computing
- To: mathgroup at smc.vnet.net
- Subject: [mg94006] Re: [mg93885] Floating-Point Computing
- From: Sseziwa Mukasa <mukasa at jeol.com>
- Date: Tue, 2 Dec 2008 00:39:22 -0500 (EST)
- References: <200811281004.FAA28552@smc.vnet.net>
On Nov 28, 2008, at 5:04 AM, 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]]]] > > Any ideas? Here are a few: (Debug) In[56]:= Timing[N[Sum[Sqrt[i],{i,0,1000000000}],31]] Timing[Sum[Sqrt[N[i]],{i,0,1000000000}]] Timing[NSum[Sqrt[i],{i,0,1000000000},WorkingPrecision->31]] Timing[NSum[Sqrt[i],{i,0,1000000000},WorkingPrecision->40]] Timing[N[Zeta[-(1/2)]-Zeta[-(1/2),1000000001],31]] (Debug) Out[56]= {0.175015,2.108185108360037596259382529338*10^13} (Debug) Out[57]= {0.175027,Zeta[-(1/2)]-Zeta[-(1/2),1000000001]} (Debug) Out[58]= {0.030142,2.1081851083600375962594*10^13} (Debug) Out[59]= {0.093878,2.108185108360037596259382529338*10^13} (Debug) Out[60]= {0.000041,2.108185108360037596259382529338*10^13} Regards, Sseziwa Mukasa