meanwhile there are some timing results by Raymond Toy on the list, but I want to add my
results, which convinced me to commit the new functions.
And, of course, the justification is needed.
I did the computations on a Windows 2000 notebook with the Maxima standard installation,
which means GCL is used. The CPU was running at 600MHz and the computations used
about 128 MB of my 512 MB RAM.
(%i1) fpprintprec:10$
(%i2) fpprec:200000$
(%i3) pi_ref: bfloat(%pi)$
(%i4) time(%o3);
(%o4) [308.94]
this is the result with the original 5.13.0 code
(%i5) :lisp (setq bigfloat%pi '((bigfloat simp 56.) 56593902016227522. 2))
((BIGFLOAT SIMP 56) 56593902016227522 2)
(%i5) load("E:\\home\\maxima\\chudnovsky_pi.o")$
this sets %pi back to the default value (there is a storing mechanism in fppi) and loads the
new functions
(%i6) pi: bfloat(%pi)$
(%i7) time(%o6);
(%o7) [23.94]
to justify the result I took a higher precision and compared
(%i8) fpprec:200010$
(%i9) pi_10: bfloat(%pi)$
(%i10) pi_10-pi;
(%o10) -6.417479737b-200001
(%i11) pi_10-pi_ref;
(%o11) 5.416844384b-199997
this proves, that the new code is even more accurate
I also checked the first 5000 digits against a printed table.
here is the same for Eulers number:
(%i1) fpprintprec:10$
(%i2) fpprec:50000$
(%i3) e_ref: bfloat(%e)$
(%i4) time(%o3);
(%o4) [570.74]
(%i5) :lisp (setq bigfloat%e '((bigfloat simp 56.) 48968212118944587. 2))
((BIGFLOAT SIMP 56) 48968212118944587 2)
(%i5) load("E:\\home\\maxima\\taylor_e.o");
(%o5) E:\home\maxima\taylor_e.o
(%i6) e: bfloat(%e)$
(%i7) time(%o6);
(%o7) [0.21]
(%i8) fpprec:50010$
(%i9) e_10: bfloat(%e)$
(%i10) e_10-e;
(%o10) -2.989818907b-50001
(%i11) e_10-e_ref;
(%o11) -4.270544836b-50000
for this computation I basically use the Taylor series
E = 1 + 1/1 + 1/(2*1) + 1/(3*2*1) + 1/(4*3*2*1) + 1/(5*4*3*2*1) + etc..
the computation takes remarkably less time, if you summarize a definite number of
summands, e.g.
1/1 + 1/(2*1) + 1/(3*2*1) + 1/(4*3*2*1)+ 1/(5*4*3*2*1)
= [ 5*4*3*2 + 5*4*3 + 5*4 + 5 + 1 ] / (5*4*3*2*1)
In this case 5 - 1 time consuming divisions are avoided. That's the trick.
The number of 1001, that I use in the code, isn't something special. It just says, take a lot
for great precisions, but don't take too much for small precisions. Perhaps a precision
dependend number would be better. Maybe (*quo (isqrt fpprec) 10)
the PI-computations published by Bailey et al. (e.g. Brent-Salamin, Borwein-Borwein) using
recursive algorithms aren't so fast in Maxima because the bigfloat computations (lots of nth-
roots) with full precision take too much time.
the Chudnovsky algorithm only needs two single bigfloat computations. the rest is bignum
arithmetic.
beside Chudnovsky and Ramanujan I got good results with atan-series (published by
Haenel)
Regards
Volker van Nek
Am 4 Nov 2007 um 16:27 hat Richard Fateman geschrieben:
> Could you tell us about your timing results? On what computer and which
> lisp?
> Comppi does not use bigfloats at all.
>
> I think it uses only integer calculations, though there is a "quotient"
> which is like truncate, I think.
>
> I think that there are fast results from Chudnovskys, but I think there is a
> more recent survey of such results in
> "Experimental Mathematics" by Bailey et al.
>
> Similarly, fpexp1 uses a taylor series, so your comment that you changed it
> seems to require some justification (timing, accuracy, etc.)
>
> I think that changes to code that works already should be held to extremely
> high standards.
>
> Regards
> RJF
>
>
> > -----Original Message-----
> > From: maxima-bounces at math.utexas.edu
> > [mailto:maxima-bounces at math.utexas.edu] On Behalf Of van Nek
> > Sent: Sunday, November 04, 2007 1:32 PM
> > To: fateman at cs.berkeley.edu
> > Cc: Maxima at math.utexas.edu
> > Subject: Re: [Maxima] comppi
> >
> > Am 13 Oct 2007 um 22:01 hat Richard Fateman geschrieben:
> >
> > > I haven't looked at it, but I think I copied it from some
> > continued fraction
> > > approximation in Hakmem, probably written out by Bill Gosper.
> > > There are likely faster-converging formulas known now.
> >
> > It took me some time to find out, what series it was.
> > I found, that it was the Newton series based on
> > pi = 6*asin(1/2)
> >
> > Meanwhile I was looking for a faster way to compute PI in Maxima.
> > The best I found is a some kind of Ramanujan series by
> > Chudnovsky & Chudnovsky.
> > The computation mainly uses bignum calculations and is
> > therefor faster than other
> > algorithms which are based on bigfloat calculations (e.g.
> > Borwein's recursive algorithms).
> >
> > I changed float.lisp/comppi and also wrote a new function
> > float.lisp/compe which computes
> > E by some modified Taylor series.
> >
> > Volker van Nek
> > _______________________________________________
> > Maxima mailing list
> > Maxima at math.utexas.edu
> > http://www.math.utexas.edu/mailman/listinfo/maxima
> >
> _______________________________________________
> Maxima mailing list
> Maxima at math.utexas.edu
> http://www.math.utexas.edu/mailman/listinfo/maxima