pidigits benchmark and bpi() method in Math::BigFloat/Math::BigInt, take 7 [PATCH]

Message-Id: <200706081929.44888@bloodgate.com>

p4raw-id: //depot/perl@31364

11 files changed, 236 insertions, 28 deletions

diff --git a/lib/Math/BigFloat.pm b/lib/Math/BigFloat.pm
index 3498eb8698..dbc929d369 100644
--- a/lib/Math/BigFloat.pm
+++ b/lib/Math/BigFloat.pm
@@ -16,7 +16,8 @@ $VERSION = '1.58';
 require 5.006002;
 
 require Exporter;
-@ISA =       qw(Math::BigInt);
+@ISA		= qw/Math::BigInt/;
+@EXPORT_OK	= qw/bpi/;
 
 use strict;
 # $_trap_inf/$_trap_nan are internal and should never be accessed from outside
@@ -2324,6 +2325,137 @@ sub bpow
   }
 
 ###############################################################################
+# trigonometric functions
+
+# helper function for bpi()
+
+sub _signed_sub
+  {
+  my ($a, $s, $b) = @_;
+ 
+  if ($s == 0)
+    {
+    # $a and $b are negativ: -> add 
+    $MBI->_add($a, $b);
+    }
+  else
+    {
+    my $c = $MBI->_acmp($a,$b);
+    # $a positiv, $b negativ
+    if ($c >= 0)
+      {
+      $MBI->_sub($a,$b);
+      }
+    else
+      {
+      # make negativ
+      $a = $MBI->_sub( $MBI->_copy($b), $a);
+      $s = 0;
+      }
+    }
+
+  ($a,$s);
+  }
+
+sub _signed_add
+  {
+  my ($a, $s, $b) = @_;
+ 
+  if ($s == 1)
+    {
+    # $a and $b are positiv: -> add 
+    $MBI->_add($a, $b);
+    }
+  else
+    {
+    my $c = $MBI->_acmp($a,$b);
+    # $a positiv, $b negativ
+    if ($c >= 0)
+      {
+      $MBI->_sub($a,$b);
+      }
+    else
+      {
+      # make positiv
+      $a = $MBI->_sub( $MBI->_copy($b), $a);
+      $s = 1;
+      }
+    }
+
+  ($a,$s);
+  }
+
+sub bpi
+  {
+  # Calculate PI to N digits (including the 3 before the dot).
+
+  # The basic algorithm is the one implemented in:
+
+  # The Computer Language Shootout
+  #   http://shootout.alioth.debian.org/
+  #
+  #   contributed by Robert Bradshaw
+  #      modified by Ruud H.G.van Tol
+  #      modified by Emanuele Zeppieri
+ 
+  # We re-implement it here by using the low-level library directly. Also,
+  # the functions consume() and extract_digit() were inlined and some
+  # rendundand operations ( like *= 1 ) were removed.
+
+  my ($self,$n) = @_;
+  if (@_ == 1)
+    {
+    # called like Math::BigFloat::bpi(10);
+    $n = $self; $self = $class;
+    }
+  $self = ref($self) if ref($self);
+  $n = 40 if !defined $n || $n < 1;
+
+  my $z0 = $MBI->_one();
+  my $z1 = $MBI->_zero();
+  my $z2 = $MBI->_one();
+  my $ten = $MBI->_ten();
+  my $three = $MBI->_new(3);
+  my ($s, $d, $e, $r); my $k = 0; my $z1_sign = 0;
+
+  # main loop
+  for (1..$n)
+    {
+    while ( 1 < 3 )
+      {
+      if ($MBI->_acmp($z0,$z2) != 1)
+        {
+        # my $o = $z0 * 3 + $z1;
+        my $o = $MBI->_mul( $MBI->_copy($z0), $three);
+        $z1_sign == 0 ? $MBI->_sub( $o, $z1) : $MBI->_add( $o, $z1);
+        ($d,$r) = $MBI->_div( $MBI->_copy($o), $z2 );
+        $d = $MBI->_num($d);
+        $e = $MBI->_num( scalar $MBI->_div( $MBI->_add($o, $z0), $z2 ) );
+        last if $d == $e;
+        }
+      $k++;
+      my $k2 = $MBI->_new( 2*$k+1 );
+      # mul works regardless of the sign of $z1 since $k is always positive
+      $MBI->_mul( $z1, $k2 );
+      ($z1, $z1_sign) = _signed_add($z1, $z1_sign, $MBI->_mul( $MBI->_new(4*$k+2), $z0 ) );
+      $MBI->_mul( $z0, $MBI->_new($k) );
+      $MBI->_mul( $z2, $k2 );
+      }
+    $MBI->_mul( $z1, $ten );
+    ($z1, $z1_sign) = _signed_sub($z1, $z1_sign, $MBI->_mul( $MBI->_new(10*$d), $z2 ) );
+    $MBI->_mul( $z0, $ten );
+    $s .= $d;
+    }
+
+  my $x = $self->new(0);
+  $x->{_es} = '-';
+  $x->{_e} = $MBI->_new(length($s)-1);
+  $x->{_m} = $MBI->_new($s);
+
+  $x;
+  }
+
+###############################################################################
 # rounding functions
 
 sub bfround
@@ -2728,11 +2860,8 @@ sub import
 
   # register us with MBI to get notified of future lib changes
   Math::BigInt::_register_callback( $self, sub { $MBI = $_[0]; } );
-   
-  # any non :constant stuff is handled by our parent, Exporter
-  # even if @_ is empty, to give it a chance
-  $self->SUPER::import(@a);      	# for subclasses
-  $self->export_to_level(1,$self,@a);	# need this, too
+
+  $self->export_to_level(1,$self,@a);		# export wanted functions
   }
 
 sub bnorm
@@ -2887,13 +3016,16 @@ Math::BigFloat - Arbitrary size floating point math package
   use Math::BigFloat;
 
   # Number creation
-  $x = Math::BigFloat->new($str);	# defaults to 0
-  $nan  = Math::BigFloat->bnan();	# create a NotANumber
-  $zero = Math::BigFloat->bzero();	# create a +0
-  $inf = Math::BigFloat->binf();	# create a +inf
-  $inf = Math::BigFloat->binf('-');	# create a -inf
-  $one = Math::BigFloat->bone();	# create a +1
-  $one = Math::BigFloat->bone('-');	# create a -1
+  my $x = Math::BigFloat->new($str);	# defaults to 0
+  my $y = $x->copy();			# make a true copy
+  my $nan  = Math::BigFloat->bnan();	# create a NotANumber
+  my $zero = Math::BigFloat->bzero();	# create a +0
+  my $inf = Math::BigFloat->binf();	# create a +inf
+  my $inf = Math::BigFloat->binf('-');	# create a -inf
+  my $one = Math::BigFloat->bone();	# create a +1
+  my $mone = Math::BigFloat->bone('-');	# create a -1
+
+  my $pi = Math::BigFloat->bpi(100);	# PI to 100 digits
 
   # Testing
   $x->is_zero();		# true if arg is +0
@@ -3253,6 +3385,14 @@ function. The result is equivalent to:
 
 This method was added in v1.84 of Math::BigInt (April 2007).
 
+=head2 bpi()
+
+	print Math::BigFloat->bpi(100), "\n";
+
+Calculate PI to N digits (including the 3 before the dot).
+
+This method was added in v1.87 of Math::BigInt (June 2007).
+
 =head1 Autocreating constants
 
 After C<use Math::BigFloat ':constant'> all the floating point constants
@@ -3329,6 +3469,14 @@ request a different storage class for use with Math::BigFloat:
 However, this request is ignored, as the current code now uses the low-level
 math libary for directly storing the number parts.
 
+=head1 EXPORTS
+
+C<Math::BigFloat> exports nothing by default, but can export the C<bpi()> method:
+
+	use Math::BigFloat qw/bpi/;
+
+	print bpi(10), "\n";
+
 =head1 BUGS
 
 Please see the file BUGS in the CPAN distribution Math::BigInt for known bugs.
diff --git a/lib/Math/BigInt.pm b/lib/Math/BigInt.pm
index f99eea617e..10aeb1b025 100644
--- a/lib/Math/BigInt.pm
+++ b/lib/Math/BigInt.pm
@@ -2831,6 +2831,27 @@ sub __lcm
   }
 
 ###############################################################################
+# trigonometric functions
+
+sub bpi
+  {
+  # Calculate PI to N digits. Unless upgrading is in effect, returns the
+  # result truncated to an integer, that is, always returns '3'.
+  my ($self,$n) = @_;
+  if (@_ == 1)
+    {
+    # called like Math::BigInt::bpi(10);
+    $n = $self; $self = $class;
+    }
+  $self = ref($self) if ref($self);
+
+  return $upgrade->new($n) if defined $upgrade;
+
+  # hard-wired to "3"
+  $self->new(3);
+  }
+
+###############################################################################
 # this method returns 0 if the object can be modified, or 1 if not.
 # We use a fast constant sub() here, to avoid costly calls. Subclasses
 # may override it with special code (f.i. Math::BigInt::Constant does so)
@@ -2865,14 +2886,17 @@ Math::BigInt - Arbitrary size integer/float math package
   my $n = 1; my $sign = '-';
 
   # Number creation	
-  $x = Math::BigInt->new($str);		# defaults to 0
-  $y = $x->copy();			# make a true copy
-  $nan  = Math::BigInt->bnan(); 	# create a NotANumber
-  $zero = Math::BigInt->bzero();	# create a +0
-  $inf = Math::BigInt->binf();		# create a +inf
-  $inf = Math::BigInt->binf('-');	# create a -inf
-  $one = Math::BigInt->bone();		# create a +1
-  $one = Math::BigInt->bone('-');	# create a -1
+  my $x = Math::BigInt->new($str);	# defaults to 0
+  my $y = $x->copy();			# make a true copy
+  my $nan  = Math::BigInt->bnan(); 	# create a NotANumber
+  my $zero = Math::BigInt->bzero();	# create a +0
+  my $inf = Math::BigInt->binf();	# create a +inf
+  my $inf = Math::BigInt->binf('-');	# create a -inf
+  my $one = Math::BigInt->bone();	# create a +1
+  my $mone = Math::BigInt->bone('-');	# create a -1
+
+  my $pi = Math::BigInt->bpi();		# returns '3'
+					# see Math::BigFloat::bpi()
 
   $h = Math::BigInt->new('0x123');	# from hexadecimal
   $b = Math::BigInt->new('0b101');	# from binary
@@ -3497,6 +3521,23 @@ function. The result is equivalent to:
 
 This method was added in v1.84 of Math::BigInt (April 2007).
 
+=head2 bpi()
+
+	print Math::BigInt->bpi(100), "\n";		# 3
+
+Returns PI truncated to an integer, with the argument being ignored. that
+is it always returns C<3>.
+
+If upgrading is in effect, returns PI to N digits (including the "3"
+before the dot):
+
+	use Math::BigFloat;
+	use Math::BigInt upgrade => Math::BigFloat;
+	print Math::BigInt->bpi(3), "\n";		# 3.14
+	print Math::BigInt->bpi(100), "\n";		# 3.1415....
+
+This method was added in v1.87 of Math::BigInt (June 2007).
+
 =head2 blsft()
 
 	$x->blsft($y);		# left shift in base 2
@@ -4385,6 +4426,13 @@ All other methods upgrade themselves only when one (or all) of their
 arguments are of the class mentioned in $upgrade (This might change in later
 versions to a more sophisticated scheme):
 
+=head1 EXPORTS
+
+C<Math::BigInt> exports nothing by default, but can export the following methods:
+
+	bgcd
+	blcm
+
 =head1 BUGS
 
 =over 2
diff --git a/lib/Math/BigInt/t/bare_mbf.t b/lib/Math/BigInt/t/bare_mbf.t
index 4d471e2cf7..b2f1fd0a18 100644
--- a/lib/Math/BigInt/t/bare_mbf.t
+++ b/lib/Math/BigInt/t/bare_mbf.t
@@ -27,7 +27,7 @@ BEGIN
     }
   print "# INC = @INC\n";
 
-  plan tests => 2064;
+  plan tests => 2070;
   }
 
 use Math::BigFloat lib => 'BareCalc';
diff --git a/lib/Math/BigInt/t/bare_mbi.t b/lib/Math/BigInt/t/bare_mbi.t
index edb482e531..7c359c8948 100644
--- a/lib/Math/BigInt/t/bare_mbi.t
+++ b/lib/Math/BigInt/t/bare_mbi.t
@@ -26,7 +26,7 @@ BEGIN
     }
   print "# INC = @INC\n";
 
-  plan tests => 3093;
+  plan tests => 3099;
   }
 
 use Math::BigInt lib => 'BareCalc';
diff --git a/lib/Math/BigInt/t/bigfltpm.inc b/lib/Math/BigInt/t/bigfltpm.inc
index 9006afe1ca..7e71d7a7c6 100644
--- a/lib/Math/BigInt/t/bigfltpm.inc
+++ b/lib/Math/BigInt/t/bigfltpm.inc
@@ -63,6 +63,8 @@ while (<DATA>)
       # some is_xxx test function	
       } elsif ($f =~ /^is_(zero|one|negative|positive|odd|even|nan|int)$/) {
         $try .= "\$x->$f();";
+      } elsif ($f eq "bpi") {
+        $try .= '$class->bpi($x);';
       } elsif ($f eq "finc") {
         $try .= '++$x;';
       } elsif ($f eq "fdec") {
@@ -399,6 +401,10 @@ abc:+0:NaN
 +27:+90:270
 +1034:+804:415668
 $div_scale = 40;
+&bpi
+77:3.1415926535897932384626433832795028841971693993751058209749445923078164062862
++0:3.141592653589793238462643383279502884197
+11:3.1415926535
 &bnok
 +inf:10:inf
 NaN:NaN:NaN
diff --git a/lib/Math/BigInt/t/bigfltpm.t b/lib/Math/BigInt/t/bigfltpm.t
index 69686d009d..1645663f1c 100755
--- a/lib/Math/BigInt/t/bigfltpm.t
+++ b/lib/Math/BigInt/t/bigfltpm.t
@@ -26,7 +26,7 @@ BEGIN
     }
   print "# INC = @INC\n";
 
-  plan tests => 2064
+  plan tests => 2070
 	+ 5;		# own tests
   }
 
diff --git a/lib/Math/BigInt/t/bigintpm.inc b/lib/Math/BigInt/t/bigintpm.inc
index 49a68fc9be..b2065c6d65 100644
--- a/lib/Math/BigInt/t/bigintpm.inc
+++ b/lib/Math/BigInt/t/bigintpm.inc
@@ -91,6 +91,8 @@ while (<DATA>)
     $try .= '"$m,$e";';
    }elsif ($f eq "bexp"){
     $try .= "\$x->bexp();";
+   } elsif ($f eq "bpi"){
+    $try .= "$class\->bpi(\$x);";
    } else {
     # binary ops
     $try .= "\$y = $class->new('$args[1]');";
@@ -2232,6 +2234,10 @@ NaN:NaN
 inf:inf
 1:2
 2:7
+&bpi
+77:3
++0:3
+11:3
 # see t/bignok.t for more tests
 &bnok
 +inf:10:inf
diff --git a/lib/Math/BigInt/t/bigintpm.t b/lib/Math/BigInt/t/bigintpm.t
index 606143c7b2..f8e2edaa07 100755
--- a/lib/Math/BigInt/t/bigintpm.t
+++ b/lib/Math/BigInt/t/bigintpm.t
@@ -10,7 +10,7 @@ BEGIN
   my $location = $0; $location =~ s/bigintpm.t//;
   unshift @INC, $location; # to locate the testing files
   chdir 't' if -d 't';
-  plan tests => 3093;
+  plan tests => 3099;
   }
 
 use Math::BigInt lib => 'Calc';
diff --git a/lib/Math/BigInt/t/sub_mbf.t b/lib/Math/BigInt/t/sub_mbf.t
index f0761f89a2..5f25de80f4 100755
--- a/lib/Math/BigInt/t/sub_mbf.t
+++ b/lib/Math/BigInt/t/sub_mbf.t
@@ -26,7 +26,7 @@ BEGIN
     }
   print "# INC = @INC\n"; 
   
-  plan tests => 2064
+  plan tests => 2070
     + 6;	# + our own tests
   }
 
diff --git a/lib/Math/BigInt/t/sub_mbi.t b/lib/Math/BigInt/t/sub_mbi.t
index 798bfa2413..f35b819f55 100755
--- a/lib/Math/BigInt/t/sub_mbi.t
+++ b/lib/Math/BigInt/t/sub_mbi.t
@@ -26,7 +26,7 @@ BEGIN
     }
   print "# INC = @INC\n";
 
-  plan tests => 3093
+  plan tests => 3099
     + 5;	# +5 own tests
   }
 
diff --git a/lib/Math/BigInt/t/with_sub.t b/lib/Math/BigInt/t/with_sub.t
index a66b6e166c..c8c9fc993e 100644
--- a/lib/Math/BigInt/t/with_sub.t
+++ b/lib/Math/BigInt/t/with_sub.t
@@ -28,7 +28,7 @@ BEGIN
     }
   print "# INC = @INC\n";
 
-  plan tests => 2064
+  plan tests => 2070
 	+ 1;
   }