From: Tels <nospam-abuse@bloodgate.com>
Date: Tue, 26 Jun 2007 21:00:53 +0000 (+0200)
Subject: Math::BigInt take 12 [PATCH]
X-Git-Url: http://git.shadowcat.co.uk/gitweb/gitweb.cgi?a=commitdiff_plain;h=30afc38d44f16307c535215ac91b9b6259f70b14;p=p5sagit%2Fp5-mst-13.2.git

Math::BigInt take 12 [PATCH]
Message-Id: <200706262100.54138@bloodgate.com>

p4raw-id: //depot/perl@31478
---

diff --git a/lib/Math/BigFloat.pm b/lib/Math/BigFloat.pm
index 3762574..7fb9863 100644
--- a/lib/Math/BigFloat.pm
+++ b/lib/Math/BigFloat.pm
@@ -2488,11 +2488,11 @@ sub _atan_inv
   #        sign = 1 - sign
 
   my $a = $MBI->_one();
-  my $b = $MBI->_new($x);
+  my $b = $MBI->_copy($x);
  
-  my $x2  = $MBI->_mul( $MBI->_new($x), $b);		# x2 = x * x
+  my $x2  = $MBI->_mul( $MBI->_copy($x), $b);		# x2 = x * x
   my $d   = $MBI->_new( 3 );				# d = 3
-  my $c   = $MBI->_mul( $MBI->_new($x), $x2);		# c = x ^ 3
+  my $c   = $MBI->_mul( $MBI->_copy($x), $x2);		# c = x ^ 3
   my $two = $MBI->_new( 2 );
 
   # run the first step unconditionally
@@ -2567,7 +2567,7 @@ sub _atan_inv
 
     }
 
-#  print "Took $step steps for $x\n";
+#  print "Took $step steps for ", $MBI->_str($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);
@@ -2597,12 +2597,12 @@ sub bpi
   # 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);
+  my ($a,$b) = $self->_atan_inv( $MBI->_new(239),$n);
+  my ($c,$d) = $self->_atan_inv( $MBI->_new(1023),$n);
+  my ($e,$f) = $self->_atan_inv( $MBI->_new(5832),$n);
+  my ($g,$h) = $self->_atan_inv( $MBI->_new(110443),$n);
+  my ($i,$j) = $self->_atan_inv( $MBI->_new(4841182),$n);
+  my ($k,$l) = $self->_atan_inv( $MBI->_new(6826318),$n);
 
   $MBI->_mul($a, $MBI->_new(732));
   $MBI->_mul($c, $MBI->_new(128));
@@ -2830,38 +2830,139 @@ sub bsin
 
 sub batan2
   { 
-  # calculate arcus tangens of ($x/$y)
+  # calculate arcus tangens of ($y/$x)
   
   # set up parameters
-  my ($self,$x,$y,@r) = (ref($_[0]),@_);
+  my ($self,$y,$x,@r) = (ref($_[0]),@_);
   # objectify is costly, so avoid it
   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
     {
-    ($self,$x,$y,@r) = objectify(2,@_);
+    ($self,$y,$x,@r) = objectify(2,@_);
     }
 
-  return $x if $x->modify('batan2');
+  return $y if $y->modify('batan2');
 
-  return $x->bnan() if (($x->{sign} eq $nan) ||
-			($y->{sign} eq $nan) ||
-			($x->is_zero() && $y->is_zero()));
+  return $y->bnan() if ($y->{sign} eq $nan) || ($x->{sign} eq $nan);
 
-  # inf handling
-  if (($x->{sign} =~ /^[+-]inf$/) || ($y->{sign} =~ /^[+-]inf$/))
+  return $upgrade->new($y)->batan2($upgrade->new($x),@r) if defined $upgrade;
+
+  # Y X
+  # 0 0 result is 0
+  # 0 +x result is 0
+  return $y->bzero(@r) if $y->is_zero() && $x->{sign} eq '+';
+
+  # Y X
+  # 0 -x result is PI
+  if ($y->is_zero())
     {
-    # XXX TODO:
-    return $x->bnan();
+    # calculate PI
+    my $pi = $self->bpi(@r);
+    # modify $x in place
+    $y->{_m} = $pi->{_m};
+    $y->{_e} = $pi->{_e};
+    $y->{_es} = $pi->{_es};
+    $y->{sign} = '+';
+    return $y;
+    }
+
+  # Y X
+  # +y 0 result is PI/2
+  # -y 0 result is -PI/2
+  if ($y->is_inf() || $x->is_zero())
+    {
+    # calculate PI/2
+    my $pi = $self->bpi(@r);
+    # modify $x in place
+    $y->{_m} = $pi->{_m};
+    $y->{_e} = $pi->{_e};
+    $y->{_es} = $pi->{_es};
+    # -y => -PI/2, +y => PI/2
+    $y->{sign} = substr($y->{sign},0,1);		# +inf => +
+    $MBI->_div($y->{_m}, $MBI->_new(2));
+    return $y;
     }
 
-  return $upgrade->new($x)->batan2($upgrade->new($y),@r) if defined $upgrade;
+  # we need to limit the accuracy to protect against overflow
+  my $fallback = 0;
+  my ($scale,@params);
+  ($y,@params) = $y->_find_round_parameters(@r);
+    
+  # error in _find_round_parameters?
+  return $y if $y->is_nan();
+
+  # no rounding at all, so must use fallback
+  if (scalar @params == 0)
+    {
+    # simulate old behaviour
+    $params[0] = $self->div_scale();	# and round to it as accuracy
+    $params[1] = undef;			# disable P
+    $scale = $params[0]+4; 		# at least four more for proper round
+    $params[2] = $r[2];			# round mode by caller or undef
+    $fallback = 1;			# to clear a/p afterwards
+    }
+  else
+    {
+    # the 4 below is empirical, and there might be cases where it is not
+    # enough...
+    $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
+    }
+
+  # inlined is_one() && is_one('-')
+  if ($MBI->_is_one($y->{_m}) && $MBI->_is_zero($y->{_e}))
+    {
+    # shortcut: 1 1 result is PI/4
+    # inlined is_one() && is_one('-')
+    if ($MBI->_is_one($x->{_m}) && $MBI->_is_zero($x->{_e}))
+      {
+      # 1,1 => PI/4
+      my $pi_4 = $self->bpi( $scale - 3);
+      # modify $x in place
+      $y->{_m} = $pi_4->{_m};
+      $y->{_e} = $pi_4->{_e};
+      $y->{_es} = $pi_4->{_es};
+      # 1 1 => +
+      # -1 1 => -
+      # 1 -1 => -
+      # -1 -1 => +
+      $y->{sign} = $x->{sign} eq $y->{sign} ? '+' : '-';
+      $MBI->_div($y->{_m}, $MBI->_new(4));
+      return $y;
+      }
+    # shortcut: 1 int(X) result is _atan_inv(X)
+
+    # is integer
+    if ($x->{_es} eq '+')
+      {
+      my $x1 = $MBI->_copy($x->{_m});
+      $MBI->_lsft($x1, $x->{_e},10) unless $MBI->_is_zero($x->{_e});
+
+      my ($a,$b) = $self->_atan_inv($x1, $scale);
+      my $y_sign = $y->{sign};
+      # calculate A/B
+      $y->bone(); $y->{_m} = $a; my $y_d = $self->bone(); $y_d->{_m} = $b;
+      $y->bdiv($y_d, @r);
+      $y->{sign} = $y_sign;
+      return $y;
+      }
+    }
+
+  # handle all other cases
+  #  X  Y
+  # +x +y 0 to PI/2
+  # -x +y PI/2 to PI
+  # +x -y 0 to -PI/2
+  # -x -y -PI/2 to -PI 
+
+  my $y_sign = $y->{sign};
 
   # divide $x by $y
-  $x->bdiv($y)->batan(@r);
+  $y->bdiv($x, $scale) unless $x->is_one();
+  $y->batan(@r);
 
-  # set the sign of $x depending on $y
-  $x->{sign} = '-' if $y->{sign} eq '-';
+  # restore sign
+  $y->{sign} = $y_sign;
 
-  $x;
+  $y;
   }
 
 sub batan
@@ -2874,17 +2975,6 @@ sub batan
   #    atan = x - --- + --- - --- + --- ...
   #                3     5     7     9 
 
-  # XXX TODO:
-  # This series is only valid if -1 < x < 1, so for other x we need to
-  # find a different way.
-
-  if ($x < -1 || $x > 1)
-    {
-    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);
@@ -2893,6 +2983,24 @@ sub batan
   #         constant object       or error in _find_round_parameters?
   return $x if $x->modify('batan') || $x->is_nan();
 
+  if ($x->{sign} =~ /^[+-]inf\z/)
+    {
+    # +inf result is PI/2
+    # -inf result is -PI/2
+    # calculate PI/2
+    my $pi = $self->bpi(@r);
+    # modify $x in place
+    $x->{_m} = $pi->{_m};
+    $x->{_e} = $pi->{_e};
+    $x->{_es} = $pi->{_es};
+    # -y => -PI/2, +y => PI/2
+    $x->{sign} = substr($x->{sign},0,1);		# +inf => +
+    $MBI->_div($x->{_m}, $MBI->_new(2));
+    return $x;
+    }
+
+  return $x->bzero(@r) if $x->is_zero();
+
   # no rounding at all, so must use fallback
   if (scalar @params == 0)
     {
@@ -2910,6 +3018,35 @@ sub batan
     $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
     }
 
+  # 1 or -1 => PI/4
+  # inlined is_one() && is_one('-')
+  if ($MBI->_is_one($x->{_m}) && $MBI->_is_zero($x->{_e}))
+    {
+    my $pi = $self->bpi($scale - 3);
+    # modify $x in place
+    $x->{_m} = $pi->{_m};
+    $x->{_e} = $pi->{_e};
+    $x->{_es} = $pi->{_es};
+    # leave the sign of $x alone (+1 => +PI/4, -1 => -PI/4)
+    $MBI->_div($x->{_m}, $MBI->_new(4));
+    return $x;
+    }
+  
+  # This series is only valid if -1 < x < 1, so for other x we need to
+  # to calculate PI/ - atan(1/x):
+  my $one = $MBI->_new(1);
+  my $pi = undef;
+  if ($x->{_es} eq '+' && ($MBI->_acmp($x->{_m},$one) >= 0))
+    {
+    # calculate PI/2
+    $pi = $self->bpi($scale - 3);
+    $MBI->_div($pi->{_m}, $MBI->_new(2));
+    # calculate 1/$x:
+    my $x_copy = $x->copy();
+    # modify $x in place
+    $x->bone(); $x->bdiv($x_copy,$scale);
+    }
+
   # when user set globals, they would interfere with our calculation, so
   # disable them and later re-enable them
   no strict 'refs';
@@ -2954,6 +3091,17 @@ sub batan
     $below->badd($two);					# n += 2
     }
 
+  if (defined $pi)
+    {
+    my $x_copy = $x->copy();
+    # modify $x in place
+    $x->{_m} = $pi->{_m};
+    $x->{_e} = $pi->{_e};
+    $x->{_es} = $pi->{_es};
+    # PI/2 - $x
+    $x->bsub($x_copy);
+    }
+
   # shortcut to not run through _find_round_parameters again
   if (defined $params[0])
     {
@@ -3549,6 +3697,10 @@ Math::BigFloat - Arbitrary size floating point math package
   my $sin  = Math::BigFloat->new(1)->bsin(100);		# sinus(1)
   my $atan = Math::BigFloat->new(1)->batan(100);	# arcus tangens(1)
 
+  my $atan2 = Math::BigFloat->new(  1 )->batan2( 1 ,100); # batan(1)
+  my $atan2 = Math::BigFloat->new(  1 )->batan2( 8 ,100); # batan(1/8)
+  my $atan2 = Math::BigFloat->new( -2 )->batan2( 1 ,100); # batan(-2)
+
   # Testing
   $x->is_zero();		# true if arg is +0
   $x->is_nan();			# true if arg is NaN
@@ -3935,21 +4087,23 @@ Calculate the sinus of $x, modifying $x in place.
 
 This method was added in v1.87 of Math::BigInt (June 2007).
 
-=head2 batan()
+=head2 batan2()
 
-	my $x = Math::BigFloat->new(1);
-	print $x->batan(100), "\n";
+	my $y = Math::BigFloat->new(2);
+	my $x = Math::BigFloat->new(3);
+	print $y->batan2($x), "\n";
 
-Calculate the arcus tanges of $x, modifying $x in place.
+Calculate the arcus tanges of C<$y> divided by C<$x>, modifying $y in place.
+See also L<batan()>.
 
 This method was added in v1.87 of Math::BigInt (June 2007).
 
-=head2 batan2()
+=head2 batan()
 
-	my $x = Math::BigFloat->new(0.5);
-	print $x->batan2(0.5), "\n";
+	my $x = Math::BigFloat->new(1);
+	print $x->batan(100), "\n";
 
-Calculate the arcus tanges of $x / $y, modifying $x in place.
+Calculate the arcus tanges of $x, modifying $x in place. See also L<batan2()>.
 
 This method was added in v1.87 of Math::BigInt (June 2007).
 
diff --git a/lib/Math/BigInt.pm b/lib/Math/BigInt.pm
index fcb6a78..a9bc6c7 100644
--- a/lib/Math/BigInt.pm
+++ b/lib/Math/BigInt.pm
@@ -2970,33 +2970,30 @@ sub bsin
 
 sub batan2
   { 
-  # calculate arcus tangens of ($x/$y)
+  # calculate arcus tangens of ($y/$x)
  
   # set up parameters
-  my ($self,$x,$y,@r) = (ref($_[0]),@_);
+  my ($self,$y,$x,@r) = (ref($_[0]),@_);
   # objectify is costly, so avoid it
   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
     {
-    ($self,$x,$y,@r) = objectify(2,@_);
+    ($self,$y,$x,@r) = objectify(2,@_);
     }
 
-  return $x if $x->modify('batan2');
+  return $y if $y->modify('batan2');
 
-  return $x->bnan() if (($x->{sign} eq $nan) ||
-			($y->{sign} eq $nan) ||
-			($x->is_zero() && $y->is_zero()));
+  return $y->bnan() if ($y->{sign} eq $nan) || ($x->{sign} eq $nan);
+
+  return $y->bzero() if	$y->is_zero() && $x->{sign} eq '+';		# x >= 0
 
   # inf handling
-  if (($x->{sign} =~ /^[+-]inf$/) || ($y->{sign} =~ /^[+-]inf$/))
-    {
-    # XXX TODO:
-    return $x->bnan();
-    }
+  # +-inf => --PI/2 => +-1
+  return $y->bone( substr($y->{sign},0,1) ) if $y->{sign} =~ /^[+-]inf$/;
 
-  return $upgrade->new($x)->batan2($upgrade->new($y),@r) if defined $upgrade;
+  return $upgrade->new($y)->batan2($upgrade->new($x),@r) if defined $upgrade;
 
   require Math::BigFloat;
-  my $r = Math::BigFloat->new($x)->batan2(Math::BigFloat->new($y),@r)->as_int();
+  my $r = Math::BigFloat->new($y)->batan2(Math::BigFloat->new($x),@r)->as_int();
 
   $x->{value} = $r->{value};
   $x->{sign} = $r->{sign};
@@ -3744,24 +3741,25 @@ integer.
 
 This method was added in v1.87 of Math::BigInt (June 2007).
 
-=head2 batan()
+=head2 batan2()
 
 	my $x = Math::BigInt->new(1);
-	print $x->batan(100), "\n";
+	my $y = Math::BigInt->new(1);
+	print $y->batan2($x), "\n";
 
-Calculate the arcus tanges of $x, modifying $x in place.
+Calculate the arcus tangens of C<$y> divided by C<$x>, modifying $y in place.
 
 In BigInt, unless upgrading is in effect, the result is truncated to an
 integer.
 
 This method was added in v1.87 of Math::BigInt (June 2007).
 
-=head2 batan2()
+=head2 batan()
 
-	my $x = Math::BigInt->new(1);
-	print $x->batan2(5), "\n";
+	my $x = Math::BigFloat->new(0.5);
+	print $x->batan(100), "\n";
 
-Calculate the arcus tanges of $x / $y, modifying $x in place.
+Calculate the arcus tangens of $x, modifying $x in place.
 
 In BigInt, unless upgrading is in effect, the result is truncated to an
 integer.
diff --git a/lib/Math/BigInt/t/bare_mbf.t b/lib/Math/BigInt/t/bare_mbf.t
index b501695..53b9fb4 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 => 2236;
+  plan tests => 2292;
   }
 
 use Math::BigFloat lib => 'BareCalc';
diff --git a/lib/Math/BigInt/t/bare_mbi.t b/lib/Math/BigInt/t/bare_mbi.t
index 793fd69..6b572ee 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 => 3217;
+  plan tests => 3257;
   }
 
 use Math::BigInt lib => 'BareCalc';
diff --git a/lib/Math/BigInt/t/bigfltpm.inc b/lib/Math/BigInt/t/bigfltpm.inc
index 6225146..60b5cd2 100644
--- a/lib/Math/BigInt/t/bigfltpm.inc
+++ b/lib/Math/BigInt/t/bigfltpm.inc
@@ -434,11 +434,41 @@ $div_scale = 40;
 0.2:13:0.1986693307951
 3.2:12:-0.0583741434276
 &batan
-0.2:14:0.19739555984988
+NaN:10:NaN
+inf:14:1.5707963267949
+-inf:14:-1.5707963267949
 0.2:13:0.1973955598499
+0.2:14:0.19739555984988
+0:10:0
+1:14:0.78539816339744
+-1:14:-0.78539816339744
+# test an argument X > 1
+2:14:1.1071487177941
 &batan2
+NaN:1:10:NaN
+NaN:NaN:10:NaN
+1:NaN:10:NaN
+inf:1:14:1.5707963267949
+-inf:1:14:-1.5707963267949
 1:5:13:0.1973955598499
 1:5:14:0.19739555984988
+0:0:10:0
+0:1:14:0
+0:2:14:0
+1:0:14:1.5707963267949
+5:0:14:1.5707963267949
+-1:0:11:-1.5707963268
+-2:0:77:-1.5707963267948966192313216916397514420985846996875529104874722961539082031431
+2:0:77:1.5707963267948966192313216916397514420985846996875529104874722961539082031431
+-1:5:14:-0.19739555984988
+1:5:14:0.19739555984988
+-1:8:14:-0.12435499454676
+1:8:14:0.12435499454676
+-1:1:14:-0.78539816339744
+# test an argument X > 1 and one X < 1
+1:2:24:0.463647609000806116214256
+2:1:14:1.1071487177941
+-2:1:14:-1.1071487177941
 &bpi
 150:3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664709384460955058223172535940813
 77:3.1415926535897932384626433832795028841971693993751058209749445923078164062862
diff --git a/lib/Math/BigInt/t/bigfltpm.t b/lib/Math/BigInt/t/bigfltpm.t
index 83a5da3..0956883 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 => 2236
+  plan tests => 2292
 	+ 5;		# own tests
   }
 
diff --git a/lib/Math/BigInt/t/bigintpm.inc b/lib/Math/BigInt/t/bigintpm.inc
index dbcb469..f948156 100644
--- a/lib/Math/BigInt/t/bigintpm.inc
+++ b/lib/Math/BigInt/t/bigintpm.inc
@@ -176,6 +176,8 @@ while (<DATA>)
        $try .= "\$x->bmodinv(\$y);";
       }elsif ($f eq "digit"){
         $try .= "\$x->digit(\$y);";
+      }elsif ($f eq "batan2"){
+        $try .= "\$x->batan2(\$y);";
       } else {
        # Functions with three arguments
        $try .= "\$z = $class->new(\"$args[2]\");";
@@ -2296,6 +2298,27 @@ NaN:NaN
 inf:inf
 1:2
 2:7
+&batan2
+NaN:1:10:NaN
+NaN:NaN:10:NaN
+1:NaN:10:NaN
+inf:1:14:1
+-inf:1:14:-1
+1:5:13:0
+1:5:14:0
+0:0:10:0
+0:1:14:0
+0:2:14:0
+1:0:14:1
+5:0:14:1
+-1:0:11:-1
+-2:0:77:-1
+2:0:77:1
+-1:5:14:0
+1:5:14:0
+-1:8:14:0
+1:8:14:0
+-1:1:14:0
 &bpi
 77:3
 +0:3
diff --git a/lib/Math/BigInt/t/bigintpm.t b/lib/Math/BigInt/t/bigintpm.t
index 86a41c5..948f05b 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 => 3217;
+  plan tests => 3257;
   }
 
 use Math::BigInt lib => 'Calc';
diff --git a/lib/Math/BigInt/t/sub_mbf.t b/lib/Math/BigInt/t/sub_mbf.t
index 09a0b29..6d64cf7 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 => 2236
+  plan tests => 2292
     + 6;	# + our own tests
   }
 
diff --git a/lib/Math/BigInt/t/sub_mbi.t b/lib/Math/BigInt/t/sub_mbi.t
index b9c598a..ea38bff 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 => 3217
+  plan tests => 3257
     + 5;	# +5 own tests
   }
 
diff --git a/lib/Math/BigInt/t/with_sub.t b/lib/Math/BigInt/t/with_sub.t
index e6e4815..3e829b1 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 => 2236
+  plan tests => 2292
 	+ 1;
   }