Upgrade to Math-BigRat-0.20
authorSteve Peters <steve@fisharerojo.org>
Wed, 18 Jul 2007 18:52:12 +0000 (18:52 +0000)
committerSteve Peters <steve@fisharerojo.org>
Wed, 18 Jul 2007 18:52:12 +0000 (18:52 +0000)
p4raw-id: //depot/perl@31629

MANIFEST
lib/Math/BigRat.pm
lib/Math/BigRat/t/bigfltpm.inc
lib/Math/BigRat/t/bigfltrt.t
lib/Math/BigRat/t/biglog.t [new file with mode: 0644]
lib/Math/BigRat/t/bigroot.t [new file with mode: 0644]

index e7b7c0a..74494ae 100644 (file)
--- a/MANIFEST
+++ b/MANIFEST
@@ -2072,10 +2072,12 @@ lib/Math/BigRat.pm              Math::BigRat
 lib/Math/BigRat/t/big_ap.t             Math::BigRat test
 lib/Math/BigRat/t/bigfltpm.inc         Math::BigRat test
 lib/Math/BigRat/t/bigfltrt.t           Math::BigRat test
+lib/Math/BigRat/t/biglog.t             Math::BigRat test
 lib/Math/BigRat/t/bigratpm.inc         Math::BigRat test
 lib/Math/BigRat/t/bigratpm.t           Math::BigRat test
 lib/Math/BigRat/t/bigrat.t             Math::BigRat test
 lib/Math/BigRat/t/bigratup.t   test under $Math::BigInt::upgrade
+lib/Math/BigRat/t/bigroot.t            Math::BigRat test
 lib/Math/BigRat/t/requirer.t   see if require works properly
 lib/Math/BigRat/t/trap.t       see if trap_nan and trap_inf work
 lib/Math/Complex.pm            A Complex package
index 7732c36..a2f0796 100644 (file)
@@ -23,7 +23,7 @@ use vars qw($VERSION @ISA $upgrade $downgrade
 
 @ISA = qw(Math::BigFloat);
 
-$VERSION = '0.19';
+$VERSION = '0.20';
 
 use overload;                  # inherit overload from Math::BigFloat
 
@@ -338,6 +338,12 @@ sub config
   # return (later set?) configuration data as hash ref
   my $class = shift || 'Math::BigRat';
 
+  if (@_ == 1 && ref($_[0]) ne 'HASH')
+    {
+    my $cfg = $class->SUPER::config();
+    return $cfg->{$_[0]};
+    }
+
   my $cfg = $class->SUPER::config(@_);
 
   # now we need only to override the ones that are different from our parent
@@ -1011,6 +1017,130 @@ sub blog
   $x->_new_from_float( $x->_as_float()->blog(Math::BigFloat->new("$y"),@r) );
   }
 
+sub bexp
+  {
+  # set up parameters
+  my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
+
+  # objectify is costly, so avoid it
+  if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
+    {
+    ($self,$x,$y,$a,$p,$r) = objectify(2,$class,@_);
+    }
+
+  return $x->binf() if $x->{sign} eq '+inf';
+  return $x->bzero() if $x->{sign} eq '-inf';
+
+  # 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);
+
+  # also takes care of the "error in _find_round_parameters?" case
+  return $x if $x->{sign} eq '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;                 # P = undef
+    $scale = $params[0]+4;              # at least four more for proper round
+    $params[2] = $r;                    # 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's not enough...
+    $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
+    }
+
+  return $x->bone(@params) if $x->is_zero();
+
+  # See the comments in Math::BigFloat on how this algorithm works.
+  # Basically we calculate A and B (where B is faculty(N)) so that A/B = e
+
+  my $x_org = $x->copy();
+  if ($scale <= 75)
+    {
+    # set $x directly from a cached string form
+    $x->{_n} = $MBI->_new("90933395208605785401971970164779391644753259799242");
+    $x->{_d} = $MBI->_new("33452526613163807108170062053440751665152000000000");
+    $x->{sign} = '+';
+    }
+  else
+    {
+    # compute A and B so that e = A / B.
+
+    # After some terms we end up with this, so we use it as a starting point:
+    my $A = $MBI->_new("90933395208605785401971970164779391644753259799242");
+    my $F = $MBI->_new(42); my $step = 42;
+
+    # Compute how many steps we need to take to get $A and $B sufficiently big
+    my $steps = Math::BigFloat::_len_to_steps($scale - 4);
+#    print STDERR "# Doing $steps steps for ", $scale-4, " digits\n";
+    while ($step++ <= $steps)
+      {
+      # calculate $a * $f + 1
+      $A = $MBI->_mul($A, $F);
+      $A = $MBI->_inc($A);
+      # increment f
+      $F = $MBI->_inc($F);
+      }
+    # compute $B as factorial of $steps (this is faster than doing it manually)
+    my $B = $MBI->_fac($MBI->_new($steps));
+
+#  print "A ", $MBI->_str($A), "\nB ", $MBI->_str($B), "\n";
+
+    $x->{_n} = $A;
+    $x->{_d} = $B;
+    $x->{sign} = '+';
+    }
+
+  # $x contains now an estimate of e, with some surplus digits, so we can round
+  if (!$x_org->is_one())
+    {
+    # raise $x to the wanted power and round it in one step:
+    $x->bpow($x_org, @params);
+    }
+  else
+    {
+    # else just round the already computed result
+    delete $x->{_a}; delete $x->{_p};
+    # 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};
+    }
+
+  $x;
+  }
+
+sub bnok
+  {
+  # 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,$class,@_);
+    }
+
+  # do it with floats
+  $x->_new_from_float( $x->_as_float()->bnok(Math::BigFloat->new("$y"),@r) );
+  }
+
 sub _float_from_part
   {
   my $x = shift;
@@ -1602,10 +1732,6 @@ Calculates the factorial of $x. For instance:
 
 Works currently only for integers.
 
-=head2 blog()
-
-Is not yet implemented.
-
 =head2 bround()/round()/bfround()
 
 Are not yet implemented.
@@ -1739,6 +1865,30 @@ Compute $x ** $y.
 
 Please see the documentation in L<Math::BigInt> for further details.
 
+=head2 bexp()
+
+       $x->bexp($accuracy);            # calculate e ** X
+
+Calculates two integers A and B so that A/B is equal to C<e ** $x>, where C<e> is
+Euler's number.
+
+This method was added in v0.20 of Math::BigRat (May 2007).
+
+See also L<blog()>.
+
+=head2 bnok()
+
+       $x->bnok($y);              # x over y (binomial coefficient n over k)
+
+Calculates the binomial coefficient n over k, also called the "choose"
+function. The result is equivalent to:
+
+       ( n )      n!
+       | - |  = -------
+       ( k )    k!(n-k)!
+
+This method was added in v0.20 of Math::BigRat (May 2007).
+
 =head2 config()
 
         use Data::Dumper;
index 45f48ac..9006afe 100644 (file)
@@ -115,6 +115,8 @@ while (<DATA>)
         $try .= '$x->facmp($y);';
       } elsif ($f eq "fpow") {
         $try .= '$x ** $y;';
+      } elsif ($f eq "bnok") {
+        $try .= '$x->bnok($y);';
       } elsif ($f eq "froot") {
         $try .= "$setup; \$x->froot(\$y);";
       } elsif ($f eq "fadd") {
@@ -397,6 +399,21 @@ abc:+0:NaN
 +27:+90:270
 +1034:+804:415668
 $div_scale = 40;
+&bnok
++inf:10:inf
+NaN:NaN:NaN
+NaN:1:NaN
+1:NaN:NaN
+1:1:1
+# k > n
+1:2:0
+2:3:0
+# k < 0
+1:-2:0
+# 7 over 3 = 35
+7:3:35
+7:6:1
+100:90:17310309456440
 &flog
 0::NaN
 -1::NaN
index 3ad98db..ff911b9 100755 (executable)
@@ -29,7 +29,7 @@ BEGIN
   plan tests => 1;
   }
 
-use Math::BigRat::Test;                # test via this Subclass 
+use Math::BigRat::Test lib => 'Calc';  # test via this Subclass 
 
 use vars qw ($class $try $x $y $f @args $ans $ans1 $ans1_str $setup $CL);
 $class = "Math::BigRat::Test";
@@ -37,5 +37,5 @@ $CL = "Math::BigInt::Calc";
  
 ok (1,1);
 
-# fails stil 185 tests
+# fails still too many tests
 #require 'bigfltpm.inc';               # all tests here for sharing
diff --git a/lib/Math/BigRat/t/biglog.t b/lib/Math/BigRat/t/biglog.t
new file mode 100644 (file)
index 0000000..d201c41
--- /dev/null
@@ -0,0 +1,99 @@
+#!/usr/bin/perl -w
+
+# Test blog function (and bpow, since it uses blog), as well as bexp().
+
+use Test::More;
+use strict;
+
+BEGIN
+  {
+  $| = 1;
+  # to locate the testing files
+  my $location = $0; $location =~ s/biglog.t//i;
+  if ($ENV{PERL_CORE})
+    {
+    # testing with the core distribution
+    @INC = qw(../lib);
+    }
+  unshift @INC, '../lib';
+  if (-d 't')
+    {
+    chdir 't';
+    require File::Spec;
+    unshift @INC, File::Spec->catdir(File::Spec->updir, $location);
+    }
+  else
+    {
+    unshift @INC, $location;
+    }
+  print "# INC = @INC\n";
+
+  plan tests => 14;
+  }
+
+use Math::BigRat;
+
+my $cl = "Math::BigRat";
+
+#############################################################################
+# test log($n)
+
+# does not work yet
+#is ($cl->new(2)->blog(), '0', "blog(2)");
+#is ($cl->new(288)->blog(), '5',"blog(288)");
+#is ($cl->new(2000)->blog(), '7', "blog(2000)");
+
+#############################################################################
+# test exp($n)
+
+is ($cl->new(1)->bexp()->as_int(), '2', "bexp(1)");
+is ($cl->new(2)->bexp()->as_int(), '7',"bexp(2)");
+is ($cl->new(3)->bexp()->as_int(), '20', "bexp(3)");
+
+# rounding not implemented yet
+#is ($cl->new(3)->bexp(10), '20', "bexp(3,10)");
+
+# $x < 0 => NaN
+ok ($cl->new(-2)->blog(), 'NaN');
+ok ($cl->new(-1)->blog(), 'NaN');
+ok ($cl->new(-10)->blog(), 'NaN');
+ok ($cl->new(-2,2)->blog(), 'NaN');
+
+#############################################################################
+# test bexp() with cached results
+
+is ($cl->new(1)->bexp(), 
+  '90933395208605785401971970164779391644753259799242' . '/' .
+  '33452526613163807108170062053440751665152000000000',
+  'bexp(1)');
+#is ($cl->new(2)->bexp(40), $cl->new(1)->bexp(45)->bpow(2,40), 'bexp(2)'); 
+
+#is ($cl->new("12.5")->bexp(61), $cl->new(1)->bexp(65)->bpow(12.5,61), 'bexp(12.5)'); 
+
+#############################################################################
+# test bexp() with big values (non-cached)
+
+#is ($cl->new(1)->bexp(100), 
+#  '2.718281828459045235360287471352662497757247093699959574966967627724076630353547594571382178525166427',
+# 'bexp(100)');
+
+is ($cl->new("12.5")->bexp(91), $cl->new(1)->bexp(95)->bpow(12.5,91), 
+  'bexp(12.5) to 91 digits'); 
+
+#############################################################################
+# some integer results
+is ($cl->new(2)->bpow(32)->blog(2),  '32', "2 ** 32");
+is ($cl->new(3)->bpow(32)->blog(3),  '32', "3 ** 32");
+is ($cl->new(2)->bpow(65)->blog(2),  '65', "2 ** 65");
+
+my $x = Math::BigInt->new( '777' ) ** 256;
+my $base = Math::BigInt->new( '12345678901234' );
+is ($x->copy()->blog($base), 56, 'blog(777**256, 12345678901234)');
+
+$x = Math::BigInt->new( '777' ) ** 777;
+$base = Math::BigInt->new( '777' );
+is ($x->copy()->blog($base), 777, 'blog(777**777, 777)');
+
+# all done
+1;
+
diff --git a/lib/Math/BigRat/t/bigroot.t b/lib/Math/BigRat/t/bigroot.t
new file mode 100644 (file)
index 0000000..41fee89
--- /dev/null
@@ -0,0 +1,69 @@
+#!/usr/bin/perl -w
+
+# Test broot function (and bsqrt() function, since it is used by broot()).
+
+# It is too slow to be simple included in bigfltpm.inc, where it would get
+# executed 3 times.
+
+# But it is better to test the numerical functionality, instead of not testing
+# it at all.
+
+use Test::More;
+use strict;
+
+BEGIN
+  {
+  $| = 1;
+  # to locate the testing files
+  my $location = $0; $location =~ s/bigroot.t//i;
+  if ($ENV{PERL_CORE})
+    {
+    # testing with the core distribution
+    @INC = qw(../lib);
+    }
+  unshift @INC, '../lib';
+  if (-d 't')
+    {
+    chdir 't';
+    require File::Spec;
+    unshift @INC, File::Spec->catdir(File::Spec->updir, $location);
+    }
+  else
+    {
+    unshift @INC, $location;
+    }
+  print "# INC = @INC\n";
+
+  plan tests => 4 * 2;
+  }
+
+use Math::BigFloat;
+use Math::BigInt;
+
+my $cl = "Math::BigFloat";
+my $c = "Math::BigInt";
+
+# 2 ** 240 = 
+# 1766847064778384329583297500742918515827483896875618958121606201292619776
+
+# takes way too long
+#test_broot ('2','240', 8, undef,   '1073741824');
+#test_broot ('2','240', 9, undef,   '106528681.3099908308759836475139583940127');
+#test_broot ('2','120', 9, undef,   '10321.27324073880096577298929482324664787');
+#test_broot ('2','120', 17, undef,   '133.3268493632747279600707813049418888729');
+
+test_broot ('2','120', 8, undef,   '32768');
+test_broot ('2','60', 8, undef,   '181.0193359837561662466161566988413540569');
+test_broot ('2','60', 9, undef,   '101.5936673259647663841091609134277286651');
+test_broot ('2','60', 17, undef,   '11.54672461623965153271017217302844672562');
+
+sub test_broot
+  {
+  my ($x,$n,$y,$scale,$result) = @_;
+
+  my $s = $scale || 'undef';
+  is ($cl->new($x)->bpow($n)->broot($y,$scale),$result, "Try: $cl $x->bpow($n)->broot($y,$s) == $result");
+  $result =~ s/\..*//;
+  is ($c->new($x)->bpow($n)->broot($y,$scale),$result, "Try: $c $x->bpow($n)->broot($y,$s) == $result");
+  }
+