[Caffeine-Patch] Math::BigInt 1.87 take 11 (add batan2, fix batan, speedup bpi()

Message-Id: <200706231530.49865@bloodgate.com>

p4raw-id: //depot/perl@31459

1 files changed, 210 insertions, 103 deletions

diff --git a/lib/Math/BigFloat.pm b/lib/Math/BigFloat.pm
index 3c3cf523e7..3762574e9d 100644
--- a/lib/Math/BigFloat.pm
+++ b/lib/Math/BigFloat.pm
@@ -50,7 +50,8 @@ use overload
 # accessor methods instead.
 
 # class constants, use Class->constant_name() to access
-$round_mode = 'even'; # one of 'even', 'odd', '+inf', '-inf', 'zero' or 'trunc'
+# one of 'even', 'odd', '+inf', '-inf', 'zero', 'trunc' or 'common'
+$round_mode = 'even';
 $accuracy   = undef;
 $precision  = undef;
 $div_scale  = 40;
@@ -2438,81 +2439,142 @@ sub bmodpow
 ###############################################################################
 # trigonometric functions
 
-# helper function for bpi()
+# helper function for bpi() and batan2(), calculates arcus tanges (1/x)
 
-sub _signed_sub
+sub _atan_inv
   {
-  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;
-      }
-    }
+  # return a/b so that a/b approximates atan(1/x) to at least limit digits
+  my ($self, $x, $limit) = @_;
 
-  ($a,$s);
-  }
+  # Taylor:       x^3   x^5   x^7   x^9
+  #    atan = x - --- + --- - --- + --- - ...
+  #                3     5     7     9 
 
-sub _signed_add
-  {
-  my ($a, $s, $b) = @_;
+  #               1      1         1        1
+  #    atan 1/x = - - ------- + ------- - ------- + ...
+  #               x   x^3 * 3   x^5 * 5   x^7 * 7 
+
+  #               1      1         1            1
+  #    atan 1/x = - - --------- + ---------- - ----------- + ... 
+  #               5    3 * 125     5 * 3125     7 * 78125
+
+  # Subtraction/addition of a rational:
+
+  #  5    7    5*3 +- 7*4
+  #  - +- -  = ----------
+  #  4    3       4*3
+
+  # Term:  N        N+1
+  #
+  #        a             1                  a * d * c +- b
+  #        ----- +- ------------------  =  ----------------
+  #        b           d * c                b * d * c
+
+  #  since b1 = b0 * (d-2) * c
+
+  #        a             1                  a * d +- b / c
+  #        ----- +- ------------------  =  ----------------
+  #        b           d * c                b * d 
+
+  # and  d = d + 2
+  # and  c = c * x * x
+
+  #        u = d * c
+  #        stop if length($u) > limit 
+  #        a = a * u +- b
+  #        b = b * u
+  #        d = d + 2
+  #        c = c * x * x
+  #        sign = 1 - sign
+
+  my $a = $MBI->_one();
+  my $b = $MBI->_new($x);
  
-  if ($s == 1)
-    {
-    # $a and $b are positiv: -> add 
-    $MBI->_add($a, $b);
-    }
-  else
+  my $x2  = $MBI->_mul( $MBI->_new($x), $b);		# x2 = x * x
+  my $d   = $MBI->_new( 3 );				# d = 3
+  my $c   = $MBI->_mul( $MBI->_new($x), $x2);		# c = x ^ 3
+  my $two = $MBI->_new( 2 );
+
+  # run the first step unconditionally
+  my $u = $MBI->_mul( $MBI->_copy($d), $c);
+  $a = $MBI->_mul($a, $u);
+  $a = $MBI->_sub($a, $b);
+  $b = $MBI->_mul($b, $u);
+  $d = $MBI->_add($d, $two);
+  $c = $MBI->_mul($c, $x2);
+
+  # a is now a * (d-3) * c
+  # b is now b * (d-2) * c
+
+  # run the second step unconditionally
+  $u = $MBI->_mul( $MBI->_copy($d), $c);
+  $a = $MBI->_mul($a, $u);
+  $a = $MBI->_add($a, $b);
+  $b = $MBI->_mul($b, $u);
+  $d = $MBI->_add($d, $two);
+  $c = $MBI->_mul($c, $x2);
+
+  # a is now a * (d-3) * (d-5) * c * c  
+  # b is now b * (d-2) * (d-4) * c * c
+
+  # so we can remove c * c from both a and b to shorten the numbers involved:
+  $a = $MBI->_div($a, $x2);
+  $b = $MBI->_div($b, $x2);
+  $a = $MBI->_div($a, $x2);
+  $b = $MBI->_div($b, $x2);
+
+#  my $step = 0; 
+  my $sign = 0;						# 0 => -, 1 => +
+  while (3 < 5)
     {
-    my $c = $MBI->_acmp($a,$b);
-    # $a positiv, $b negativ
-    if ($c >= 0)
+#    $step++;
+#    if (($i++ % 100) == 0)
+#      {
+#    print "a=",$MBI->_str($a),"\n";
+#    print "b=",$MBI->_str($b),"\n";
+#      }
+#    print "d=",$MBI->_str($d),"\n";
+#    print "x2=",$MBI->_str($x2),"\n";
+#    print "c=",$MBI->_str($c),"\n";
+
+    my $u = $MBI->_mul( $MBI->_copy($d), $c);
+    # use _alen() for libs like GMP where _len() would be O(N^2)
+    last if $MBI->_alen($u) > $limit;
+    my ($bc,$r) = $MBI->_div( $MBI->_copy($b), $c);
+    if ($MBI->_is_zero($r))
       {
-      $MBI->_sub($a,$b);
+      # b / c is an integer, so we can remove c from all terms
+      # this happens almost every time:
+      $a = $MBI->_mul($a, $d);
+      $a = $MBI->_sub($a, $bc) if $sign == 0;
+      $a = $MBI->_add($a, $bc) if $sign == 1;
+      $b = $MBI->_mul($b, $d);
       }
     else
       {
-      # make positiv
-      $a = $MBI->_sub( $MBI->_copy($b), $a);
-      $s = 1;
+      # b / c is not an integer, so we keep c in the terms
+      # this happens very rarely, for instance for x = 5, this happens only
+      # at the following steps:
+      # 1, 5, 14, 32, 72, 157, 340, ...
+      $a = $MBI->_mul($a, $u);
+      $a = $MBI->_sub($a, $b) if $sign == 0;
+      $a = $MBI->_add($a, $b) if $sign == 1;
+      $b = $MBI->_mul($b, $u);
       }
+    $d = $MBI->_add($d, $two);
+    $c = $MBI->_mul($c, $x2);
+    $sign = 1 - $sign;
+
     }
 
-  ($a,$s);
+#  print "Took $step steps for $x\n";
+#  print "a=",$MBI->_str($a),"\n"; print "b=",$MBI->_str($b),"\n";
+  # return a/b so that a/b approximates atan(1/x)
+  ($a,$b);
   }
 
 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 (@_ == 0)
     {
@@ -2528,48 +2590,45 @@ sub bpi
   $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;
+  # 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)
+
+  # a few more to prevent rounding errors
+  $n += 4;
+
+  my ($a,$b) = $self->_atan_inv(239,$n);
+  my ($c,$d) = $self->_atan_inv(1023,$n);
+  my ($e,$f) = $self->_atan_inv(5832,$n);
+  my ($g,$h) = $self->_atan_inv(110443,$n);
+  my ($i,$j) = $self->_atan_inv(4841182,$n);
+  my ($k,$l) = $self->_atan_inv(6826318,$n);
+
+  $MBI->_mul($a, $MBI->_new(732));
+  $MBI->_mul($c, $MBI->_new(128));
+  $MBI->_mul($e, $MBI->_new(272));
+  $MBI->_mul($g, $MBI->_new(48));
+  $MBI->_mul($i, $MBI->_new(48));
+  $MBI->_mul($k, $MBI->_new(400));
+
+  my $x = $self->bone(); $x->{_m} = $a; my $x_d = $self->bone(); $x_d->{_m} = $b;
+  my $y = $self->bone(); $y->{_m} = $c; my $y_d = $self->bone(); $y_d->{_m} = $d;
+  my $z = $self->bone(); $z->{_m} = $e; my $z_d = $self->bone(); $z_d->{_m} = $f;
+  my $u = $self->bone(); $u->{_m} = $g; my $u_d = $self->bone(); $u_d->{_m} = $h;
+  my $v = $self->bone(); $v->{_m} = $i; my $v_d = $self->bone(); $v_d->{_m} = $j;
+  my $w = $self->bone(); $w->{_m} = $k; my $w_d = $self->bone(); $w_d->{_m} = $l;
+  $x->bdiv($x_d, $n);
+  $y->bdiv($y_d, $n);
+  $z->bdiv($z_d, $n);
+  $u->bdiv($u_d, $n);
+  $v->bdiv($v_d, $n);
+  $w->bdiv($w_d, $n);
+
+  delete $x->{a}; delete $y->{a}; delete $z->{a};
+  delete $u->{a}; delete $v->{a}; delete $w->{a};
+  $x->badd($y)->bsub($z)->badd($u)->bsub($v)->bsub($w);
+
+  $x->round($n-4);
   }
 
 sub bcos
@@ -2769,6 +2828,42 @@ sub bsin
   $x;
   }
 
+sub batan2
+  { 
+  # calculate arcus tangens of ($x/$y)
+  
+  # set up parameters
+  my ($self,$x,$y,@r) = (ref($_[0]),@_);
+  # objectify is costly, so avoid it
+  if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
+    {
+    ($self,$x,$y,@r) = objectify(2,@_);
+    }
+
+  return $x if $x->modify('batan2');
+
+  return $x->bnan() if (($x->{sign} eq $nan) ||
+			($y->{sign} eq $nan) ||
+			($x->is_zero() && $y->is_zero()));
+
+  # inf handling
+  if (($x->{sign} =~ /^[+-]inf$/) || ($y->{sign} =~ /^[+-]inf$/))
+    {
+    # XXX TODO:
+    return $x->bnan();
+    }
+
+  return $upgrade->new($x)->batan2($upgrade->new($y),@r) if defined $upgrade;
+
+  # divide $x by $y
+  $x->bdiv($y)->batan(@r);
+
+  # set the sign of $x depending on $y
+  $x->{sign} = '-' if $y->{sign} eq '-';
+
+  $x;
+  }
+
 sub batan
   {
   # Calculate a arcus tangens of x.
@@ -2788,6 +2883,8 @@ sub batan
     die("$x is out of range for batan()!");
     }
 
+  return $x->bzero(@r) if $x->is_zero();
+
   # we need to limit the accuracy to protect against overflow
   my $fallback = 0;
   my ($scale,@params);
@@ -2831,7 +2928,7 @@ sub batan
   my $sign = 1;				# start with -=
   my $below = $self->new(3);
   my $two = $self->new(2);
-  $x->bone(); delete $x->{a}; delete $x->{p};
+  delete $x->{a}; delete $x->{p};
 
   my $limit = $self->new("1E-". ($scale-1));
   #my $steps = 0;
@@ -3815,7 +3912,8 @@ This method was added in v1.84 of Math::BigInt (April 2007).
 
 	print Math::BigFloat->bpi(100), "\n";
 
-Calculate PI to N digits (including the 3 before the dot).
+Calculate PI to N digits (including the 3 before the dot). The result is
+rounded according to the current rounding mode, which defaults to "even".
 
 This method was added in v1.87 of Math::BigInt (June 2007).
 
@@ -3846,6 +3944,15 @@ Calculate the arcus tanges of $x, modifying $x in place.
 
 This method was added in v1.87 of Math::BigInt (June 2007).
 
+=head2 batan2()
+
+	my $x = Math::BigFloat->new(0.5);
+	print $x->batan2(0.5), "\n";
+
+Calculate the arcus tanges of $x / $y, modifying $x in place.
+
+This method was added in v1.87 of Math::BigInt (June 2007).
+
 =head2 bmuladd()
 
 	$x->bmuladd($y,$z);