[Maxima-discuss] e, pi, gamma, log(2)

[Volker v. N. <vol...@gm...>]

Greetings.

I attached new implementations of float.lisp/fpe1, fppi1, fpgamma1,
comp-log2 as a proposal.

I discovered a paper by Bruno Haible and Thomas Papanikolaou about the
technique of binary splitting to approximate infinite series.
See http://www.ginac.de/CLN/binsplit.pdf for details if you are interested.

With this technique I succeeded in implementing fast functions for
computing %e, %pi, %gamma and log(2).

[ I also implemented exp(x), log(x), sin(x) and friends for bigfloat x.
In contrast to fpe1, .. these implementions are still some kind of raw
material and need more error analysis. And of course it's a lot more
difficult to merge them into Maxima. My idea is to put them into
/share/contrib for now. ]

Below I want to describe the improvements for the 4 constants and I
would like to propose these changes and commit them if there are no
complains.

The values of %e, %pi, %gamma and log(2) computed by the described
functions are checked against reference values which were parsed and
rounded to the same precision. For all 0 < $fpprec < 10000 there are no
differences to the reference values.


1. Euler's number E

The simple unmodified Taylor-series fits best to the binary splitting
technique.

Timing results for a call to (fpe1) with GCL on 32 bit Ubuntu:
$fpprec =  500000
5.34.1:   2.039 secs
proposed: 0.370 secs
$fpprec = 1000000
5.34.1:   6.420 secs
proposed: 0.790 secs

2. PI

Best results with binary splitting and the Chudnovsky series.

Recently there are two versions of the square-root part for different
Lisp versions. I optimized the Heron/Newton version which is now the
fastest for all Lisps.

Timing results for a call to (fppi1):
$fpprec =  100000
5.34.1:   1.720 secs
proposed: 0.140 secs
$fpprec =  200000
5.34.1:   6.900 secs
proposed: 0.349 secs

3. GAMMA

In 2013 Brent and Johansson succeeded in proving the Brent-McMillan
algorithm which was a conjecture in 1980. They also proved an error
bound. See http://arxiv.org/pdf/1312.0039v1.pdf

The Brent-McMillan algorithm is a refinement of the Bessel function
approach which is used in Maxima 5.34.1.

Haible and Papanikolaou described how to compute the Bessel function
approach via binary splitting and it wasn't very hard to extend this to
Brent-McMillan.

Timing results for a call to (fpgamma1):
$fpprec =  2000
5.34.1:   1.049 secs
proposed: 0.029 secs
$fpprec =  4000
5.34.1:   6.590 secs
proposed: 0.070 secs

4. LOG(2)

Same algorithm as in 5.34.1. Just compute the Taylor series via binary
splitting.

Timing results for a call to (comp-log2):
$fpprec =  5000
5.34.1:   1.230 secs
proposed: 0.009 secs
$fpprec =  10000
5.34.1:   6.599 secs
proposed: 0.019 secs


Volker van Nek