Math::BigInt v1.87 take 10
Tels [Fri, 22 Jun 2007 19:02:22 +0000 (21:02 +0200)]
Message-Id: <200706221902.22487@bloodgate.com>

p4raw-id: //depot/perl@31449

lib/Math/BigFloat.pm
lib/Math/BigInt.pm
lib/Math/BigInt/t/bare_mbf.t
lib/Math/BigInt/t/bigfltpm.inc
lib/Math/BigInt/t/bigfltpm.t
lib/Math/BigInt/t/fallback.t
lib/Math/BigInt/t/sub_mbf.t
lib/Math/BigInt/t/with_sub.t

index bc4ca90..3c3cf52 100644 (file)
@@ -2244,13 +2244,13 @@ sub bfac
 
 sub _pow
   {
-  # Calculate a power where $y is a non-integer, like 2 ** 0.5
-  my ($x,$y,$a,$p,$r) = @_;
+  # Calculate a power where $y is a non-integer, like 2 ** 0.3
+  my ($x,$y,@r) = @_;
   my $self = ref($x);
 
   # if $y == 0.5, it is sqrt($x)
   $HALF = $self->new($HALF) unless ref($HALF);
-  return $x->bsqrt($a,$p,$r,$y) if $y->bcmp($HALF) == 0;
+  return $x->bsqrt(@r,$y) if $y->bcmp($HALF) == 0;
 
   # Using:
   # a ** x == e ** (x * ln a)
@@ -2264,7 +2264,7 @@ sub _pow
   # we need to limit the accuracy to protect against overflow
   my $fallback = 0;
   my ($scale,@params);
-  ($x,@params) = $x->_find_round_parameters($a,$p,$r);
+  ($x,@params) = $x->_find_round_parameters(@r);
     
   return $x if $x->is_nan();           # error in _find_round_parameters?
 
@@ -2275,7 +2275,7 @@ sub _pow
     $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;                   # round mode by caller or undef
+    $params[2] = $r[2];                        # round mode by caller or undef
     $fallback = 1;                     # to clear a/p afterwards
     }
   else
@@ -2380,7 +2380,7 @@ sub bpow
     }
   if ($x_is_zero)
     {
-    return $x->bone() if $y_is_zero;
+    #return $x->bone() if $y_is_zero;
     return $x if $y->{sign} eq '+';    # 0**y => 0 (if not y <= 0)
     # 0 ** -y => 1 / (0 ** y) => 1 / 0! (1 / 0 => +inf)
     return $x->binf();
@@ -2522,6 +2522,8 @@ sub bpi
     {
     # called like Math::BigFloat::bpi(10);
     $n = $self; $self = $class;
+    # called like Math::BigFloat->bpi();
+    $n = undef if $n eq 'Math::BigFloat';
     }
   $self = ref($self) if ref($self);
   $n = 40 if !defined $n || $n < 1;
@@ -2570,6 +2572,310 @@ sub bpi
   $x;
   }
 
+sub bcos
+  {
+  # Calculate a cosinus of x.
+  my ($self,$x,@r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
+
+  # Taylor:      x^2   x^4   x^6   x^8
+  #    cos = 1 - --- + --- - --- + --- ...
+  #               2!    4!    6!    8!
+
+  # we need to limit the accuracy to protect against overflow
+  my $fallback = 0;
+  my ($scale,@params);
+  ($x,@params) = $x->_find_round_parameters(@r);
+    
+  #         constant object       or error in _find_round_parameters?
+  return $x if $x->modify('bcos') || $x->is_nan();
+
+  return $x->bone(@r) if $x->is_zero();
+
+  # 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
+    }
+
+  # when user set globals, they would interfere with our calculation, so
+  # disable them and later re-enable them
+  no strict 'refs';
+  my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
+  my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
+  # we also need to disable any set A or P on $x (_find_round_parameters took
+  # them already into account), since these would interfere, too
+  delete $x->{_a}; delete $x->{_p};
+  # need to disable $upgrade in BigInt, to avoid deep recursion
+  local $Math::BigInt::upgrade = undef;
+  my $last = 0;
+  my $over = $x * $x;                   # X ^ 2
+  my $x2 = $over->copy();               # X ^ 2; difference between terms
+  my $sign = 1;                         # start with -=
+  my $below = $self->new(2); my $factorial = $self->new(3);
+  $x->bone(); delete $x->{a}; delete $x->{p};
+
+  my $limit = $self->new("1E-". ($scale-1));
+  #my $steps = 0;
+  while (3 < 5)
+    {
+    # we calculate the next term, and add it to the last
+    # when the next term is below our limit, it won't affect the outcome
+    # anymore, so we stop:
+    my $next = $over->copy()->bdiv($below,$scale);
+    last if $next->bacmp($limit) <= 0;
+
+    if ($sign == 0)
+      {
+      $x->badd($next);
+      }
+    else
+      {
+      $x->bsub($next);
+      }
+    $sign = 1-$sign;                                   # alternate
+    # calculate things for the next term
+    $over->bmul($x2);                                  # $x*$x
+    $below->bmul($factorial); $factorial->binc();      # n*(n+1)
+    $below->bmul($factorial); $factorial->binc();      # n*(n+1)
+    }
+
+  # shortcut to not run through _find_round_parameters again
+  if (defined $params[0])
+    {
+    $x->bround($params[0],$params[2]);         # then round accordingly
+    }
+  else
+    {
+    $x->bfround($params[1],$params[2]);                # then round accordingly
+    }
+  if ($fallback)
+    {
+    # clear a/p after round, since user did not request it
+    delete $x->{_a}; delete $x->{_p};
+    }
+  # restore globals
+  $$abr = $ab; $$pbr = $pb;
+  $x;
+  }
+
+sub bsin
+  {
+  # Calculate a sinus of x.
+  my ($self,$x,@r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
+
+  # taylor:      x^3   x^5   x^7   x^9
+  #    sin = x - --- + --- - --- + --- ...
+  #               3!    5!    7!    9!
+
+  # we need to limit the accuracy to protect against overflow
+  my $fallback = 0;
+  my ($scale,@params);
+  ($x,@params) = $x->_find_round_parameters(@r);
+    
+  #         constant object       or error in _find_round_parameters?
+  return $x if $x->modify('bsin') || $x->is_nan();
+
+  return $x->bzero(@r) if $x->is_zero();
+
+  # 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
+    }
+
+  # when user set globals, they would interfere with our calculation, so
+  # disable them and later re-enable them
+  no strict 'refs';
+  my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
+  my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
+  # we also need to disable any set A or P on $x (_find_round_parameters took
+  # them already into account), since these would interfere, too
+  delete $x->{_a}; delete $x->{_p};
+  # need to disable $upgrade in BigInt, to avoid deep recursion
+  local $Math::BigInt::upgrade = undef;
+  my $last = 0;
+  my $over = $x * $x;                  # X ^ 2
+  my $x2 = $over->copy();              # X ^ 2; difference between terms
+  $over->bmul($x);                     # X ^ 3 as starting value
+  my $sign = 1;                                # start with -=
+  my $below = $self->new(6); my $factorial = $self->new(4);
+  delete $x->{a}; delete $x->{p};
+
+  my $limit = $self->new("1E-". ($scale-1));
+  #my $steps = 0;
+  while (3 < 5)
+    {
+    # we calculate the next term, and add it to the last
+    # when the next term is below our limit, it won't affect the outcome
+    # anymore, so we stop:
+    my $next = $over->copy()->bdiv($below,$scale);
+    last if $next->bacmp($limit) <= 0;
+
+    if ($sign == 0)
+      {
+      $x->badd($next);
+      }
+    else
+      {
+      $x->bsub($next);
+      }
+    $sign = 1-$sign;                                   # alternate
+    # calculate things for the next term
+    $over->bmul($x2);                                  # $x*$x
+    $below->bmul($factorial); $factorial->binc();      # n*(n+1)
+    $below->bmul($factorial); $factorial->binc();      # n*(n+1)
+    }
+
+  # shortcut to not run through _find_round_parameters again
+  if (defined $params[0])
+    {
+    $x->bround($params[0],$params[2]);         # then round accordingly
+    }
+  else
+    {
+    $x->bfround($params[1],$params[2]);                # then round accordingly
+    }
+  if ($fallback)
+    {
+    # clear a/p after round, since user did not request it
+    delete $x->{_a}; delete $x->{_p};
+    }
+  # restore globals
+  $$abr = $ab; $$pbr = $pb;
+  $x;
+  }
+
+sub batan
+  {
+  # Calculate a arcus tangens of x.
+  my ($x,@r) = @_;
+  my $self = ref($x);
+
+  # taylor:       x^3   x^5   x^7   x^9
+  #    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()!");
+    }
+
+  # we need to limit the accuracy to protect against overflow
+  my $fallback = 0;
+  my ($scale,@params);
+  ($x,@params) = $x->_find_round_parameters(@r);
+    
+  #         constant object       or error in _find_round_parameters?
+  return $x if $x->modify('batan') || $x->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
+    }
+
+  # when user set globals, they would interfere with our calculation, so
+  # disable them and later re-enable them
+  no strict 'refs';
+  my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
+  my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
+  # we also need to disable any set A or P on $x (_find_round_parameters took
+  # them already into account), since these would interfere, too
+  delete $x->{_a}; delete $x->{_p};
+  # need to disable $upgrade in BigInt, to avoid deep recursion
+  local $Math::BigInt::upgrade = undef;
+  my $last = 0;
+  my $over = $x * $x;                  # X ^ 2
+  my $x2 = $over->copy();              # X ^ 2; difference between terms
+  $over->bmul($x);                     # X ^ 3 as starting value
+  my $sign = 1;                                # start with -=
+  my $below = $self->new(3);
+  my $two = $self->new(2);
+  $x->bone(); delete $x->{a}; delete $x->{p};
+
+  my $limit = $self->new("1E-". ($scale-1));
+  #my $steps = 0;
+  while (3 < 5)
+    {
+    # we calculate the next term, and add it to the last
+    # when the next term is below our limit, it won't affect the outcome
+    # anymore, so we stop:
+    my $next = $over->copy()->bdiv($below,$scale);
+    last if $next->bacmp($limit) <= 0;
+
+    if ($sign == 0)
+      {
+      $x->badd($next);
+      }
+    else
+      {
+      $x->bsub($next);
+      }
+    $sign = 1-$sign;                                   # alternate
+    # calculate things for the next term
+    $over->bmul($x2);                                  # $x*$x
+    $below->badd($two);                                        # n += 2
+    }
+
+  # shortcut to not run through _find_round_parameters again
+  if (defined $params[0])
+    {
+    $x->bround($params[0],$params[2]);         # then round accordingly
+    }
+  else
+    {
+    $x->bfround($params[1],$params[2]);                # then round accordingly
+    }
+  if ($fallback)
+    {
+    # clear a/p after round, since user did not request it
+    delete $x->{_a}; delete $x->{_p};
+    }
+  # restore globals
+  $$abr = $ab; $$pbr = $pb;
+  $x;
+  }
+
 ###############################################################################
 # rounding functions
 
@@ -3141,6 +3447,11 @@ Math::BigFloat - Arbitrary size floating point math package
 
   my $pi = Math::BigFloat->bpi(100);   # PI to 100 digits
 
+  # the following examples compute their result to 100 digits accuracy:
+  my $cos  = Math::BigFloat->new(1)->bcos(100);                # cosinus(1)
+  my $sin  = Math::BigFloat->new(1)->bsin(100);                # sinus(1)
+  my $atan = Math::BigFloat->new(1)->batan(100);       # arcus tangens(1)
+
   # Testing
   $x->is_zero();               # true if arg is +0
   $x->is_nan();                        # true if arg is NaN
@@ -3508,6 +3819,33 @@ Calculate PI to N digits (including the 3 before the dot).
 
 This method was added in v1.87 of Math::BigInt (June 2007).
 
+=head2 bcos()
+
+       my $x = Math::BigFloat->new(1);
+       print $x->bcos(100), "\n";
+
+Calculate the cosinus of $x, modifying $x in place.
+
+This method was added in v1.87 of Math::BigInt (June 2007).
+
+=head2 bsin()
+
+       my $x = Math::BigFloat->new(1);
+       print $x->bsin(100), "\n";
+
+Calculate the sinus of $x, modifying $x in place.
+
+This method was added in v1.87 of Math::BigInt (June 2007).
+
+=head2 batan()
+
+       my $x = Math::BigFloat->new(1);
+       print $x->batan(100), "\n";
+
+Calculate the arcus tanges of $x, modifying $x in place.
+
+This method was added in v1.87 of Math::BigInt (June 2007).
+
 =head2 bmuladd()
 
        $x->bmuladd($y,$z);             
index ea4876c..134a4b7 100644 (file)
@@ -81,9 +81,8 @@ use overload
                "$_[1]" cmp $_[0]->bstr() :
                $_[0]->bstr() cmp "$_[1]" },
 
-# make cos()/sin()/atan2() "work" with BigInt's or subclasses
-'cos'  =>      sub { cos($_[0]->numify()) }, 
-'sin'  =>      sub { sin($_[0]->numify()) }, 
+'cos'  =>      sub { $_[0]->copy->bcos(); }, 
+'sin'  =>      sub { $_[0]->copy->bsin(); }, 
 'atan2'        =>      sub { $_[2] ?
                        atan2($_[1],$_[0]->numify()) :
                        atan2($_[0]->numify(),$_[1]) },
@@ -2927,6 +2926,65 @@ sub bpi
   $self->new(3);
   }
 
+sub bcos
+  {
+  # Calculate cosinus(x) to N digits. Unless upgrading is in effect, returns the
+  # result truncated to an integer.
+  my ($self,$x,@r) = ref($_[0]) ? (undef,@_) : objectify(1,@_);
+
+  return $x if $x->modify('bcos');
+
+  return $x->bnan() if $x->{sign} !~ /^[+-]\z/;        # -inf +inf or NaN => NaN
+
+  return $upgrade->new($x)->bcos(@r) if defined $upgrade;
+
+  # calculate the result and truncate it to integer
+  my $t = Math::BigFloat->new($x)->bcos(@r)->as_int();
+
+  $x->bone() if $t->is_one();
+  $x->bzero() if $t->is_zero();
+  $x->round(@r);
+  }
+
+sub bsin
+  {
+  # Calculate sinus(x) to N digits. Unless upgrading is in effect, returns the
+  # result truncated to an integer.
+  my ($self,$x,@r) = ref($_[0]) ? (undef,@_) : objectify(1,@_);
+
+  return $x if $x->modify('bsin');
+
+  return $x->bnan() if $x->{sign} !~ /^[+-]\z/;        # -inf +inf or NaN => NaN
+
+  return $upgrade->new($x)->bsin(@r) if defined $upgrade;
+
+  # calculate the result and truncate it to integer
+  my $t = Math::BigFloat->new($x)->bsin(@r)->as_int();
+
+  $x->bone() if $t->is_one();
+  $x->bzero() if $t->is_zero();
+  $x->round(@r);
+  }
+
+sub batan
+  {
+  # Calculate arcus tangens of x to N digits. Unless upgrading is in effect, returns the
+  # result truncated to an integer.
+  my ($self,$x,@r) = ref($_[0]) ? (undef,@_) : objectify(1,@_);
+
+  return $x if $x->modify('batan');
+
+  return $x->bnan() if $x->{sign} !~ /^[+-]\z/;        # -inf +inf or NaN => NaN
+
+  return $upgrade->new($x)->batan(@r) if defined $upgrade;
+
+  # calculate the result and truncate it to integer
+  my $t = Math::BigFloat->new($x)->batan(@r);
+
+  $x->{value} = $CALC->_new( $x->as_int()->bstr() );
+  $x->round(@r);
+  }
+
 ###############################################################################
 # 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
@@ -3624,6 +3682,33 @@ before the dot):
 
 This method was added in v1.87 of Math::BigInt (June 2007).
 
+=head2 bcos()
+
+       my $x = Math::BigFloat->new(1);
+       print $x->bcos(100), "\n";
+
+Calculate the cosinus of $x, modifying $x in place.
+
+This method was added in v1.87 of Math::BigInt (June 2007).
+
+=head2 bsin()
+
+       my $x = Math::BigFloat->new(1);
+       print $x->bsin(100), "\n";
+
+Calculate the sinus of $x, modifying $x in place.
+
+This method was added in v1.87 of Math::BigInt (June 2007).
+
+=head2 batan()
+
+       my $x = Math::BigFloat->new(0.5);
+       print $x->batan(100), "\n";
+
+Calculate the arcus tanges of $x, modifying $x in place.
+
+This method was added in v1.87 of Math::BigInt (June 2007).
+
 =head2 blsft()
 
        $x->blsft($y);          # left shift in base 2
index 7428535..1ff4b52 100644 (file)
@@ -27,7 +27,7 @@ BEGIN
     }
   print "# INC = @INC\n";
 
-  plan tests => 2200;
+  plan tests => 2226;
   }
 
 use Math::BigFloat lib => 'BareCalc';
index c99cc70..873a4ca 100644 (file)
@@ -119,6 +119,12 @@ while (<DATA>)
         $try .= '$x ** $y;';
       } elsif ($f eq "bnok") {
         $try .= '$x->bnok($y);';
+      } elsif ($f eq "bcos") {
+        $try .= '$x->bcos($y);';
+      } elsif ($f eq "bsin") {
+        $try .= '$x->bsin($y);';
+      } elsif ($f eq "batan") {
+        $try .= '$x->batan($y);';
       } elsif ($f eq "froot") {
         $try .= "$setup; \$x->froot(\$y);";
       } elsif ($f eq "fadd") {
@@ -142,9 +148,9 @@ while (<DATA>)
        $try .= "\$z = $class->new(\"$args[2]\");";
 
        if( $f eq "bmodpow") {
-         $try .= "\$x->bmodpow(\$y,\$z);";
+          $try .= '$x->bmodpow($y,$z);';
          } elsif ($f eq "bmuladd"){
-          $try .= "\$x->bmuladd(\$y,\$z);";
+          $try .= '$x->bmuladd($y,$z);';
          } else { warn "Unknown op '$f'"; }
        }
     }
@@ -410,6 +416,21 @@ abc:+0:NaN
 +27:+90:270
 +1034:+804:415668
 $div_scale = 40;
+&bcos
+1.2:10:0.3623577545
+2.4:12:-0.737393715541
+0:10:1
+0:20:1
+1:10:0.5403023059
+1:12:0.540302305868
+&bsin
+1:10:0.8414709848
+0:10:0
+0:20:0
+2.1:12:0.863209366649
+1.2:13:0.9320390859672
+0.2:13:0.1986693307951
+3.2:12:-0.0583741434276
 &bpi
 77:3.1415926535897932384626433832795028841971693993751058209749445923078164062862
 +0:3.141592653589793238462643383279502884197
index 1b9127e..bcaf6e1 100755 (executable)
@@ -26,7 +26,7 @@ BEGIN
     }
   print "# INC = @INC\n";
 
-  plan tests => 2200
+  plan tests => 2226
        + 5;            # own tests
   }
 
index 010ff8a..c69a3ff 100644 (file)
@@ -28,7 +28,7 @@ BEGIN
     }
   print "# INC = @INC\n";
 
-  plan tests => 9;
+  plan tests => 7;
   }
 
 # The tests below test that cos(BigInt) = cos(Scalar) which is DWIM, but not
@@ -41,18 +41,14 @@ use Math::BigFloat;
 
 my $bi = Math::BigInt->new(1);
 
-ok (cos($bi), cos(1));
-ok (sin($bi), sin(1));
+ok (cos($bi), int(cos(1)));
+ok (sin($bi), int(sin(1)));
 ok (atan2($bi,$bi), atan2(1,1));
 
 my $bf = Math::BigInt->new(0);
 
-ok (cos($bf), cos(0));
-ok (sin($bf), sin(0));
+ok (cos($bf), int(cos(0)));
+ok (sin($bf), int(sin(0)));
 ok (atan2($bi,$bf), atan2(1,0));
 ok (atan2($bf,$bi), atan2(0,1));
 
-my $bone = Math::BigInt->new(1);
-ok (cos($bone), cos(1));
-ok (sin($bone), sin(1));
-
index 08b1ae4..8c15747 100755 (executable)
@@ -26,7 +26,7 @@ BEGIN
     }
   print "# INC = @INC\n"; 
   
-  plan tests => 2200
+  plan tests => 2226
     + 6;       # + our own tests
   }
 
index a8d3ec8..695fb08 100644 (file)
@@ -28,7 +28,7 @@ BEGIN
     }
   print "# INC = @INC\n";
 
-  plan tests => 2200
+  plan tests => 2226
        + 1;
   }