Subject: | Math::BigFloat->bpi is very slow with many digits |
At one point I decided I wanted a million or so digits of Pi, and tried using Math::BigFloat. It's really slow. While the following patch certainly won't give y-cruncher or similar programs any competition, it does speed up large inputs by a great deal. If we really cared, allowing a _bpi method to the backend would allow C implementations (Pari would just get Pari's Pi constant).
I tried Wetherfield's newer arctan formula, which speeds things up about 10%. Not really worth upgrading.
I briefly looked at spigot formulas, but in pure Perl they were a lot slower than the existing arctan solution for me.
I tried the Chudnovsky's formula, which is just a little faster than the included patch using the Calc backend, but was substantially slower with the Pari or GMP backends. Someone may have a more efficient implementation.
I tried Arndt's quartic method, which runs about half the speed of the AGM method with my implementation. It's possible one of the Borwein-type formulas may be faster.
Timing for the included AGM method. We see a much better looking growth rate.
digits Time-orig (s) Time-AGM (s)
1000 0.007 0.008
10000 0.52 0.13
50000 26.0 1.7
100000 130.9 4.7
200000 631.8 12.6
400000 31.7
800000 89.1
I've used a large string from y-cruncher and compared the output at various digits to the rounded string. The only issues I found up to 2M or so are similar to the ones impacting the current code: no matter how many guard digits we use, eventually we will find a x499999... or x049999... that will go a different direction. Look, for example, at 761, 776, 16686, 820160 digits with different rounding modes. If we really, really cared we could use one of the digit formulas to add more digits as needed for proper rounding.
Subject: | bigfloat-bpi-performance.patch |
--- perl-5.19.5/lib/Math/BigFloat.pm 2013-08-26 03:14:41.000000000 -0700
+++ perl-5.19.5/lib/5.19.5/Math/BigFloat.pm 2013-11-04 00:27:23.983550954 -0800
@@ -2679,6 +2679,27 @@
my $fallback = defined $n ? 0 : 1;
$n = 40 if !defined $n || $n < 1;
+ # For large accuracy, the arctan formulas become very inefficient with
+ # Math::BigFloat. Switch to Brent-Salamin (aka AGM or Gauss-Legendre).
+ if ($n >= 1000) {
+ $n += 8;
+ $HALF = $self->new($HALF) unless ref($HALF);
+ my ($an, $bn, $tn, $pn) = ($self->bone, $HALF->copy->bsqrt($n),
+ $HALF->copy->bmul($HALF), $self->bone);
+ while ($pn < $n) {
+ my $prev_an = $an->copy;
+ $an->badd($bn)->bmul($HALF, $n);
+ $bn->bmul($prev_an)->bsqrt($n);
+ $prev_an->bsub($an);
+ $tn->bsub($pn * $prev_an * $prev_an);
+ $pn->badd($pn);
+ }
+ $an->badd($bn);
+ $an->bmul($an,$n)->bdiv(4*$tn, $n-8);
+ delete $an->{_a} if $fallback == 1;
+ return $an;
+ }
+
# after é»è¦å© (Hwang Chien-Lih) (1997)
# pi/4 = 183 * atan(1/239) + 32 * atan(1/1023) â 68 * atan(1/5832)
# + 12 * atan(1/110443) - 12 * atan(1/4841182) - 100 * atan(1/6826318)