Upgrade to Math::BigInt 1.86

p4raw-id: //depot/perl@31159

17 files changed, 571 insertions, 132 deletions

diff --git a/lib/Math/BigFloat.pm b/lib/Math/BigFloat.pm
index 1242a37813..7c2794ced9 100644
--- a/lib/Math/BigFloat.pm
+++ b/lib/Math/BigFloat.pm
@@ -12,7 +12,7 @@ package Math::BigFloat;
 #   _a	: accuracy
 #   _p	: precision
 
-$VERSION = '1.54';
+$VERSION = '1.57';
 require 5.006002;
 
 require Exporter;
@@ -599,6 +599,8 @@ sub badd
     {
     ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
     }
+ 
+  return $x if $x->modify('badd');
 
   # inf and NaN handling
   if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
@@ -676,6 +678,8 @@ sub binc
   # increment arg by one
   my ($self,$x,@r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
 
+  return $x if $x->modify('binc');
+
   if ($x->{_es} eq '-')
     {
     return $x->badd($self->bone(),@r);	#  digits after dot
@@ -711,6 +715,8 @@ sub bdec
   # decrement arg by one
   my ($self,$x,@r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
 
+  return $x if $x->modify('bdec');
+
   if ($x->{_es} eq '-')
     {
     return $x->badd($self->bone('-'),@r);	#  digits after dot
@@ -748,6 +754,8 @@ sub blog
   {
   my ($self,$x,$base,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
 
+  return $x if $x->modify('blog');
+
   # $base > 0, $base != 1; if $base == undef default to $base == e
   # $x >= 0
 
@@ -777,7 +785,7 @@ sub blog
     }
 
   return $x->bzero(@params) if $x->is_one();
-  # base not defined => base == Euler's constant e
+  # base not defined => base == Euler's number e
   if (defined $base)
     {
     # make object, since we don't feed it through objectify() to still get the
@@ -878,11 +886,70 @@ sub blog
   $x;
   }
 
+sub _len_to_steps
+  {
+  # Given D (digits in decimal), compute N so that N! (N factorial) is
+  # at least D digits long. D should be at least 50.
+  my $d = shift;
+
+  # two constants for the Ramanujan estimate of ln(N!)
+  my $lg2 = log(2 * 3.14159265) / 2;
+  my $lg10 = log(10);
+
+  # D = 50 => N => 42, so L = 40 and R = 50
+  my $l = 40; my $r = $d;
+
+  # Otherwise this does not work under -Mbignum and we do not yet have "no bignum;" :(
+  $l = $l->numify if ref($l);
+  $r = $r->numify if ref($r);
+  $lg2 = $lg2->numify if ref($lg2);
+  $lg10 = $lg10->numify if ref($lg10);
+
+  # binary search for the right value (could this be written as the reverse of lg(n!)?)
+  while ($r - $l > 1)
+    {
+    my $n = int(($r - $l) / 2) + $l;
+    my $ramanujan = 
+      int(($n * log($n) - $n + log( $n * (1 + 4*$n*(1+2*$n)) ) / 6 + $lg2) / $lg10);
+    $ramanujan > $d ? $r = $n : $l = $n;
+    }
+  $l;
+  }
+
+sub bnok
+  {
+  # Calculate n over k (binomial coefficient or "choose" function) as integer.
+  # 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('bnok');
+
+  return $x->bnan() if $x->is_nan() || $y->is_nan();
+  return $x->binf() if $x->is_inf();
+
+  my $u = $x->as_int();
+  $u->bnok($y->as_int());
+
+  $x->{_m} = $u->{value};
+  $x->{_e} = $MBI->_zero();
+  $x->{_es} = '+';
+  $x->{sign} = '+';
+  $x->bnorm(@r);
+  }
+
 sub bexp
   {
-  # Calculate e ** X (Euler's constant to the power of X)
+  # Calculate e ** X (Euler's number to the power of X)
   my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
 
+  return $x if $x->modify('bexp');
+
   return $x->binf() if $x->{sign} eq '+inf';
   return $x->bzero() if $x->{sign} eq '-inf';
 
@@ -931,7 +998,6 @@ sub bexp
   local $Math::BigFloat::downgrade = undef;
 
   my $x_org = $x->copy();
-  delete $x_org->{_a}; delete $x_org->{_p};
 
   # We use the following Taylor series:
 
@@ -977,42 +1043,70 @@ sub bexp
   # -- + -- + -- + -- + -- + --- + --- + ---- = -----
   # 1     1    2    6   24   120   720   5040   5040
 
-  # Note that we cannot simple reduce 13700/5040 to 685/252.
+  # Note that we cannot simple reduce 13700/5040 to 685/252, but must keep A and B!
 
-  # So we start with A / B = 13490 / 5040 and F = 8
-  my $A = $MBI->_new(13700);
-  my $B = $MBI->_new(5040);
-  my $F = $MBI->_new(8);
-
-  # Code based on big rational math:
-  my $limit = $MBI->_new('1' . '0' x ($scale - 2));	# scale=5 => 100000
-
-  # while $B is not yet too big (making 1/$B too small)
-  while ($MBI->_acmp($B,$limit) < 0)
+  if ($scale <= 75)
     {
-    # calculate $a * $f + 1
-    $A = $MBI->_mul($A, $F);
-    $A = $MBI->_inc($A);
-    # calculate $b * $f
-    $B = $MBI->_mul($B, $F);
-    # increment f
-    $F = $MBI->_inc($F);
+    # set $x directly from a cached string form
+    $x->{_m} = $MBI->_new(
+    "27182818284590452353602874713526624977572470936999595749669676277240766303535476");
+    $x->{sign} = '+';
+    $x->{_es} = '-';
+    $x->{_e} = $MBI->_new(79);
     }
+  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 = _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 = $self->new($MBI->_str($A))->bdiv($MBI->_str($B), $scale);
+    # compute A/B with $scale digits in the result (truncate, not round)
+    $A = $MBI->_lsft( $A, $MBI->_new($scale), 10);
+    $A = $MBI->_div( $A, $B );
 
-  delete $x->{_a}; delete $x->{_p};
-  # raise $x to the wanted power
-  $x->bpow($x_org, $scale) unless $x_org->is_one();
+    $x->{_m} = $A;
+    $x->{sign} = '+';
+    $x->{_es} = '-';
+    $x->{_e} = $MBI->_new($scale);
+    }
 
-  # shortcut to not run through _find_round_parameters again
-  if (defined $params[0])
+  # $x contains now an estimate of e, with some surplus digits, so we can round
+  if (!$x_org->is_one())
     {
-    $x->bround($params[0],$params[2]);		# then round accordingly
+    # raise $x to the wanted power and round it in one step:
+    $x->bpow($x_org, @params);
     }
   else
     {
-    $x->bfround($params[1],$params[2]);		# then round accordingly
+    # 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)
     {
@@ -1466,6 +1560,8 @@ sub bmul
     ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
     }
 
+  return $x if $x->modify('bmul');
+
   return $x->bnan() if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
 
   # inf handling
@@ -1507,6 +1603,8 @@ sub bdiv
     ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
     }
 
+  return $x if $x->modify('bdiv');
+
   return $self->_div_inf($x,$y)
    if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/) || $y->is_zero());
 
@@ -1642,6 +1740,8 @@ sub bmod
     ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
     }
 
+  return $x if $x->modify('bmod');
+
   # handle NaN, inf, -inf
   if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
     {
@@ -1737,6 +1837,8 @@ sub broot
     ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
     }
 
+  return $x if $x->modify('broot');
+
   # NaN handling: $x ** 1/0, x or y NaN, or y inf/-inf or y == 0
   return $x->bnan() if $x->{sign} !~ /^\+/ || $y->is_zero() ||
          $y->{sign} !~ /^\+$/;
@@ -1861,6 +1963,8 @@ sub bsqrt
   # calculate square root
   my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
 
+  return $x if $x->modify('bsqrt');
+
   return $x->bnan() if $x->{sign} !~ /^[+]/;	# NaN, -inf or < 0
   return $x if $x->{sign} eq '+inf';		# sqrt(inf) == inf
   return $x->round($a,$p,$r) if $x->is_zero() || $x->is_one();
@@ -2030,7 +2134,9 @@ sub bfac
   # objectify is costly, so avoid it
   ($self,$x,@r) = objectify(1,@_) if !ref($x);
 
- return $x if $x->{sign} eq '+inf';	# inf => inf
+  # inf => inf
+  return $x if $x->modify('bfac') || $x->{sign} eq '+inf';	
+
   return $x->bnan() 
     if (($x->{sign} ne '+') ||		# inf, NaN, <0 etc => NaN
      ($x->{_es} ne '+'));		# digits after dot?
@@ -2161,6 +2267,8 @@ sub bpow
     ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
     }
 
+  return $x if $x->modify('bpow');
+
   return $x->bnan() if $x->{sign} eq $nan || $y->{sign} eq $nan;
   return $x if $x->{sign} =~ /^[+-]inf$/;
   
@@ -2539,6 +2647,7 @@ sub import
   my $self = shift;
   my $l = scalar @_;
   my $lib = ''; my @a;
+  my $lib_kind = 'try';
   $IMPORT=1;
   for ( my $i = 0; $i < $l ; $i++)
     {
@@ -2560,10 +2669,11 @@ sub import
       $downgrade = $_[$i+1];		# or undef to disable
       $i++;
       }
-    elsif ($_[$i] eq 'lib')
+    elsif ($_[$i] =~ /^(lib|try|only)\z/)
       {
       # alternative library
       $lib = $_[$i+1] || '';		# default Calc
+      $lib_kind = $1;			# lib, try or only
       $i++;
       }
     elsif ($_[$i] eq 'with')
@@ -2585,7 +2695,7 @@ sub import
   if ((defined $mbilib) && ($MBI eq 'Math::BigInt::Calc'))
     {
     # MBI already loaded
-    Math::BigInt->import('try',"$lib,$mbilib", 'objectify');
+    Math::BigInt->import( $lib_kind, "$lib,$mbilib", 'objectify');
     }
   else
     {
@@ -2598,7 +2708,7 @@ sub import
     # Perl < 5.6.0 dies with "out of memory!" when eval() and ':constant' is
     # used in the same script, or eval inside import(). So we require MBI:
     require Math::BigInt;
-    Math::BigInt->import( try => $lib, 'objectify' );
+    Math::BigInt->import( $lib_kind => $lib, 'objectify' );
     }
   if ($@)
     {
@@ -2722,6 +2832,8 @@ sub as_number
   # return copy as a bigint representation of this BigFloat number
   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
 
+  return $x if $x->modify('as_number');
+
   my $z = $MBI->_copy($x->{_m});
   if ($x->{_es} eq '-')			# < 0
     {
diff --git a/lib/Math/BigInt.pm b/lib/Math/BigInt.pm
index badf947ba6..f73af00c66 100644
--- a/lib/Math/BigInt.pm
+++ b/lib/Math/BigInt.pm
@@ -18,7 +18,7 @@ package Math::BigInt;
 my $class = "Math::BigInt";
 use 5.006002;
 
-$VERSION = '1.82';
+$VERSION = '1.86';
 
 @ISA = qw(Exporter);
 @EXPORT_OK = qw(objectify bgcd blcm); 
@@ -1265,6 +1265,57 @@ sub blog
   $x->round(@r);
   }
 
+sub bnok
+  {
+  # Calculate n over k (binomial coefficient or "choose" function) as integer.
+  # 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('bnok');
+  return $x->bnan() if $x->{sign} eq 'NaN' || $y->{sign} eq 'NaN';
+  return $x->binf() if $x->{sign} eq '+inf';
+
+  # k > n or k < 0 => 0
+  my $cmp = $x->bacmp($y);
+  return $x->bzero() if $cmp < 0 || $y->{sign} =~ /^-/;
+  # k == n => 1
+  return $x->bone(@r) if $cmp == 0;
+
+  if ($CALC->can('_nok'))
+    {
+    $x->{value} = $CALC->_nok($x->{value},$y->{value});
+    }
+  else
+    {
+    # ( 7 )    7!          7*6*5 * 4*3*2*1   7 * 6 * 5
+    # ( - ) = --------- =  --------------- = ---------
+    # ( 3 )   3! (7-3)!    3*2*1 * 4*3*2*1   3 * 2 * 1 
+
+    # compute n - k + 2 (so we start with 5 in the example above)
+    my $z = $x - $y;
+    if (!$z->is_one())
+      {
+      $z->binc();
+      my $r = $z->copy(); $z->binc();
+      my $d = $self->new(2);
+      while ($z->bacmp($x) <= 0)		# f < x ?
+        {
+        $r->bmul($z); $r->bdiv($d);
+        $z->binc(); $d->binc();
+        }
+      $x->{value} = $r->{value}; $x->{sign} = '+';
+      }
+    else { $x->bone(); }
+    }
+  $x->round(@r);
+  }
+
 sub bexp
   {
   # Calculate e ** $x (Euler's number to the power of X), truncated to
@@ -1281,6 +1332,7 @@ sub bexp
   my $u;
   {
     # run through Math::BigFloat unless told otherwise
+    require Math::BigFloat unless defined $upgrade;
     local $upgrade = 'Math::BigFloat' unless defined $upgrade;
     # calculate result, truncate it to integer
     $u = $upgrade->bexp($upgrade->new($x),@r);
@@ -2880,6 +2932,8 @@ Math::BigInt - Arbitrary size integer/float math package
   $x->broot($y);	   # $y'th root of $x (e.g. $y == 3 => cubic root)
   $x->bfac();		   # factorial of $x (1*2*3*4*..$x)
 
+  $x->bnok($y);		   # x over y (binomial coefficient n over k)
+
   $x->blog();		   # logarithm of $x to base e (Euler's number)
   $x->blog($base);	   # logarithm of $x to base $base (f.i. 2)
   $x->bexp();		   # calculate e ** $x where e is Euler's number
@@ -3418,6 +3472,19 @@ This method was added in v1.82 of Math::BigInt (April 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 v1.84 of Math::BigInt (April 2007).
+
 =head2 blsft()
 
 	$x->blsft($y);		# left shift in base 2
diff --git a/lib/Math/BigInt/Calc.pm b/lib/Math/BigInt/Calc.pm
index 6fb21b00d2..359736702e 100644
--- a/lib/Math/BigInt/Calc.pm
+++ b/lib/Math/BigInt/Calc.pm
@@ -4,9 +4,7 @@ use 5.006002;
 use strict;
 # use warnings;	# dont use warnings for older Perls
 
-use vars qw/$VERSION/;
-
-$VERSION = '0.48';
+our $VERSION = '0.50';
 
 # Package to store unsigned big integers in decimal and do math with them
 
@@ -33,7 +31,7 @@ $VERSION = '0.48';
 # global constants, flags and accessory
 
 # announce that we are compatible with MBI v1.70 and up
-sub api_version () { 1; }
+sub api_version () { 2; }
  
 # constants for easier life
 my ($BASE,$BASE_LEN,$MBASE,$RBASE,$MAX_VAL,$BASE_LEN_SMALL);
@@ -42,8 +40,9 @@ my ($AND_MASK,$XOR_MASK,$OR_MASK);
 
 sub _base_len 
   {
-  # set/get the BASE_LEN and assorted other, connected values
-  # used only be the testsuite, set is used only by the BEGIN block below
+  # Set/get the BASE_LEN and assorted other, connected values.
+  # Used only by the testsuite, the set variant is used only by the BEGIN
+  # block below:
   shift;
 
   my $b = shift;
@@ -76,10 +75,8 @@ sub _base_len
     undef &_mul;
     undef &_div;
 
-    # $caught & 1 != 0 => cannot use MUL
-    # $caught & 2 != 0 => cannot use DIV
-    # The parens around ($caught & 1) were important, indeed, if we would use
-    # & here.
+    # ($caught & 1) != 0 => cannot use MUL
+    # ($caught & 2) != 0 => cannot use DIV
     if ($caught == 2)				# 2
       {
       # must USE_MUL since we cannot use DIV
@@ -172,6 +169,9 @@ BEGIN
   $AND_MASK = __PACKAGE__->_new( ( 2 ** $AND_BITS ));
   $XOR_MASK = __PACKAGE__->_new( ( 2 ** $XOR_BITS ));
   $OR_MASK = __PACKAGE__->_new( ( 2 ** $OR_BITS ));
+
+  # We can compute the approximate lenght no faster than the real length:
+  *_alen = \&_len;
   }
 
 ###############################################################################
@@ -200,6 +200,16 @@ sub _ten
   [ 10 ];
   }
 
+sub _1ex
+  {
+  # create a 1Ex
+  my $rem = $_[1] % $BASE_LEN;		# remainder
+  my $parts = $_[1] / $BASE_LEN;	# parts
+
+  # 000000, 000000, 100 
+  [ (0) x $parts, '1' . ('0' x $rem) ];
+  }
+
 sub _copy
   {
   # make a true copy
@@ -425,7 +435,8 @@ sub _mul_use_mul
     $xi = shift @prod || 0;	# || 0 makes v5.005_3 happy
     }
   push @$xv, @prod;
-  __strip_zeros($xv);
+  # can't have leading zeros
+#  __strip_zeros($xv);
   $xv;
   }                                                                             
 
@@ -435,7 +446,7 @@ sub _mul_use_div
   # multiply two numbers in internal representation
   # modifies first arg, second need not be different from first
   my ($c,$xv,$yv) = @_;
- 
+
   if (@$yv == 1)
     {
     # shortcut for two small numbers, also handles $x == 0
@@ -460,7 +471,9 @@ sub _mul_use_div
     my $y = $yv->[0]; my $car = 0;
     foreach my $i (@$xv)
       {
+      # old, slower code (before use integer;)
       $i = $i * $y + $car; $car = int($i / $MBASE); $i -= $car * $MBASE;
+      #$i = $i * $y + $car; $i -= ($car = $i / $MBASE) * $MBASE;
       }
     push @$xv, $car if $car != 0;
     return $xv;
@@ -480,14 +493,14 @@ sub _mul_use_div
     for $yi (@$yv)
       {
       $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
-      $prod[$cty++] =
-       $prod - ($car = int($prod / $MBASE)) * $MBASE;
+      $prod[$cty++] = $prod - ($car = int($prod / $MBASE)) * $MBASE;
       }
     $prod[$cty] += $car if $car; # need really to check for 0?
     $xi = shift @prod || 0;	# || 0 makes v5.005_3 happy
     }
   push @$xv, @prod;
-  __strip_zeros($xv);
+  # can't have leading zeros
+#  __strip_zeros($xv);
   $xv;
   }                                                                             
 
@@ -907,7 +920,7 @@ sub _acmp
 
 sub _len
   {
-  # compute number of digits
+  # compute number of digits in base 10
 
   # int() because add/sub sometimes leaves strings (like '00005') instead of
   # '5' in this place, thus causing length() to report wrong length
@@ -1242,18 +1255,133 @@ sub _pow
   $cx;
   }
 
+sub _nok
+  {
+  # n over k
+  # ref to array, return ref to array
+  my ($c,$n,$k) = @_;
+
+  # ( 7 )    7!          7*6*5 * 4*3*2*1   7 * 6 * 5
+  # ( - ) = --------- =  --------------- = ---------
+  # ( 3 )   3! (7-3)!    3*2*1 * 4*3*2*1   3 * 2 * 1 
+
+  # compute n - k + 2 (so we start with 5 in the example above)
+  my $x = _copy($c,$n);
+
+  _sub($c,$n,$k);
+  if (!_is_one($c,$n))
+    {
+    _inc($c,$n);
+    my $f = _copy($c,$n); _inc($c,$f);		# n = 5, f = 6, d = 2
+    my $d = _two($c);
+    while (_acmp($c,$f,$x) <= 0)		# f < n ?
+      {
+      # n = (n * f / d) == 5 * 6 / 2 => n == 3
+      $n = _mul($c,$n,$f); $n = _div($c,$n,$d);
+      # f = 7, d = 3
+      _inc($c,$f); _inc($c,$d);
+      }
+    }
+  else 
+    {
+    # keep ref to $n and set it to 1
+    splice (@$n,1); $n->[0] = 1;
+    }
+  $n;
+  }
+
+my @factorials = (
+  1,
+  1,
+  2,
+  2*3,
+  2*3*4,
+  2*3*4*5,
+  2*3*4*5*6,
+  2*3*4*5*6*7,
+);
+
 sub _fac
   {
   # factorial of $x
   # ref to array, return ref to array
   my ($c,$cx) = @_;
 
-  if ((@$cx == 1) && ($cx->[0] <= 2))
+  if ((@$cx == 1) && ($cx->[0] <= 7))
     {
-    $cx->[0] ||= 1;		# 0 => 1, 1 => 1, 2 => 2
+    $cx->[0] = $factorials[$cx->[0]];		# 0 => 1, 1 => 1, 2 => 2 etc.
     return $cx;
     }
 
+  if ((@$cx == 1) && 		# we do this only if $x >= 12 and $x <= 7000
+      ($cx->[0] >= 12 && $cx->[0] < 7000))
+    {
+
+  # Calculate (k-j) * (k-j+1) ... k .. (k+j-1) * (k + j)
+  # See http://blogten.blogspot.com/2007/01/calculating-n.html
+  # The above series can be expressed as factors:
+  #   k * k - (j - i) * 2
+  # We cache k*k, and calculate (j * j) as the sum of the first j odd integers
+
+  # This will not work when N exceeds the storage of a Perl scalar, however,
+  # in this case the algorithm would be way to slow to terminate, anyway.
+
+  # As soon as the last element of $cx is 0, we split it up and remember
+  # how many zeors we got so far. The reason is that n! will accumulate
+  # zeros at the end rather fast.
+  my $zero_elements = 0;
+
+  # If n is even, set n = n -1
+  my $k = _num($c,$cx); my $even = 1;
+  if (($k & 1) == 0)
+    {
+    $even = $k; $k --;
+    }
+  # set k to the center point
+  $k = ($k + 1) / 2;
+#  print "k $k even: $even\n";
+  # now calculate k * k
+  my $k2 = $k * $k;
+  my $odd = 1; my $sum = 1;
+  my $i = $k - 1;
+  # keep reference to x
+  my $new_x = _new($c, $k * $even);
+  @$cx = @$new_x;
+  if ($cx->[0] == 0)
+    {
+    $zero_elements ++; shift @$cx;
+    }
+#  print STDERR "x = ", _str($c,$cx),"\n";
+  my $BASE2 = int(sqrt($BASE))-1;
+  my $j = 1; 
+  while ($j <= $i)
+    {
+    my $m = ($k2 - $sum); $odd += 2; $sum += $odd; $j++;
+    while ($j <= $i && ($m < $BASE2) && (($k2 - $sum) < $BASE2))
+      {
+      $m *= ($k2 - $sum);
+      $odd += 2; $sum += $odd; $j++;
+#      print STDERR "\n k2 $k2 m $m sum $sum odd $odd\n"; sleep(1);
+      }
+    if ($m < $BASE)
+      {
+      _mul($c,$cx,[$m]);
+      }
+    else
+      {
+      _mul($c,$cx,$c->_new($m));
+      }
+    if ($cx->[0] == 0)
+      {
+      $zero_elements ++; shift @$cx;
+      }
+#    print STDERR "Calculate $k2 - $sum = $m (x = ", _str($c,$cx),")\n";
+    }
+  # multiply in the zeros again
+  unshift @$cx, (0) x $zero_elements; 
+  return $cx;
+  }
+
   # go forward until $base is exceeded
   # limit is either $x steps (steps == 100 means a result always too high) or
   # $base.
@@ -1281,19 +1409,32 @@ sub _fac
     $n = _copy($c,$cx);
     }
 
-  $cx->[0] = $last; splice (@$cx,1);		# keep ref to $x
+  # Set $cx to the last result below $BASE (but keep ref to $x)
+  $cx->[0] = $last; splice (@$cx,1);
+  # As soon as the last element of $cx is 0, we split it up and remember
+  # how many zeors we got so far. The reason is that n! will accumulate
+  # zeros at the end rather fast.
   my $zero_elements = 0;
 
   # do left-over steps fit into a scalar?
   if (ref $n eq 'ARRAY')
     {
     # No, so use slower inc() & cmp()
+    # ($n is at least $BASE here)
+    my $base_2 = int(sqrt($BASE)) - 1;
+    #print STDERR "base_2: $base_2\n"; 
+    while ($step < $base_2)
+      {
+      if ($cx->[0] == 0)
+        {
+        $zero_elements ++; shift @$cx;
+        }
+      my $b = $step * ($step + 1); $step += 2;
+      _mul($c,$cx,[$b]);
+      }
     $step = [$step];
-    while (_acmp($step,$n) <= 0)
+    while (_acmp($c,$step,$n) <= 0)
       {
-      # as soon as the last element of $cx is 0, we split it up and remember
-      # how many zeors we got so far. The reason is that n! will accumulate
-      # zeros at the end rather fast.
       if ($cx->[0] == 0)
         {
         $zero_elements ++; shift @$cx;
@@ -1304,23 +1445,45 @@ sub _fac
   else
     {
     # Yes, so we can speed it up slightly
-    while ($step <= $n)
+  
+#    print "# left over steps $n\n";
+
+    my $base_4 = int(sqrt(sqrt($BASE))) - 2;
+    #print STDERR "base_4: $base_4\n";
+    my $n4 = $n - 4; 
+    while ($step < $n4 && $step < $base_4)
       {
-      # When the last element of $cx is 0, we split it up and remember
-      # how many we got so far. The reason is that n! will accumulate
-      # zeros at the end rather fast.
       if ($cx->[0] == 0)
         {
         $zero_elements ++; shift @$cx;
         }
+      my $b = $step * ($step + 1); $step += 2; $b *= $step * ($step + 1); $step += 2;
+      _mul($c,$cx,[$b]);
+      }
+    my $base_2 = int(sqrt($BASE)) - 1;
+    my $n2 = $n - 2; 
+    #print STDERR "base_2: $base_2\n"; 
+    while ($step < $n2 && $step < $base_2)
+      {
+      if ($cx->[0] == 0)
+        {
+        $zero_elements ++; shift @$cx;
+        }
+      my $b = $step * ($step + 1); $step += 2;
+      _mul($c,$cx,[$b]);
+      }
+    # do what's left over
+    while ($step <= $n)
+      {
       _mul($c,$cx,[$step]); $step++;
+      if ($cx->[0] == 0)
+        {
+        $zero_elements ++; shift @$cx;
+        }
       }
     }
   # multiply in the zeros again
-  while ($zero_elements-- > 0)
-    {
-    unshift @$cx, 0; 
-    }
+  unshift @$cx, (0) x $zero_elements;
   $cx;			# return result
   }
 
@@ -1430,7 +1593,7 @@ sub _sqrt
 
   if (scalar @$x == 1)
     {
-    # fit's into one Perl scalar, so result can be computed directly
+    # fits into one Perl scalar, so result can be computed directly
     $x->[0] = int(sqrt($x->[0]));
     return $x;
     } 
@@ -1513,7 +1676,7 @@ sub _root
       }
     else
       {
-      # fit's into one Perl scalar, so result can be computed directly
+      # fits into one Perl scalar, so result can be computed directly
       # cannot use int() here, because it rounds wrongly (try 
       # (81 ** 3) ** (1/3) to see what I mean)
       #$x->[0] = int( $x->[0] ** (1 / $n->[0]) );
@@ -1722,7 +1885,7 @@ sub _as_hex
   # convert a decimal number to hex (ref to array, return ref to string)
   my ($c,$x) = @_;
 
-  # fit's into one element (handle also 0x0 case)
+  # fits into one element (handle also 0x0 case)
   return sprintf("0x%x",$x->[0]) if @$x == 1;
 
   my $x1 = _copy($c,$x);
@@ -1752,7 +1915,7 @@ sub _as_bin
   # convert a decimal number to bin (ref to array, return ref to string)
   my ($c,$x) = @_;
 
-  # fit's into one element (and Perl recent enough), handle also 0b0 case
+  # fits into one element (and Perl recent enough), handle also 0b0 case
   # handle zero case for older Perls
   if ($] <= 5.005 && @$x == 1 && $x->[0] == 0)
     {
@@ -1790,7 +1953,7 @@ sub _as_oct
   # convert a decimal number to octal (ref to array, return ref to string)
   my ($c,$x) = @_;
 
-  # fit's into one element (handle also 0 case)
+  # fits into one element (handle also 0 case)
   return sprintf("0%o",$x->[0]) if @$x == 1;
 
   my $x1 = _copy($c,$x);
@@ -1810,7 +1973,7 @@ sub _as_oct
 
 sub _from_oct
   {
-  # convert a octal number to decimal (ref to string, return ref to array)
+  # convert a octal number to decimal (string, return ref to array)
   my ($c,$os) = @_;
 
   # for older Perls, play safe
@@ -1836,7 +1999,7 @@ sub _from_oct
 
 sub _from_hex
   {
-  # convert a hex number to decimal (ref to string, return ref to array)
+  # convert a hex number to decimal (string, return ref to array)
   my ($c,$hs) = @_;
 
   my $m = _new($c, 0x10000000);			# 28 bit at a time (<32 bit!)
@@ -1874,7 +2037,7 @@ sub _from_hex
 
 sub _from_bin
   {
-  # convert a hex number to decimal (ref to string, return ref to array)
+  # convert a hex number to decimal (string, return ref to array)
   my ($c,$bs) = @_;
 
   # instead of converting X (8) bit at a time, it is faster to "convert" the
@@ -2007,7 +2170,7 @@ version like 'Pari'.
 The following functions MUST be defined in order to support the use by
 Math::BigInt v1.70 or later:
 
-	api_version()	return API version, minimum 1 for v1.70
+	api_version()	return API version, 1 for v1.70, 2 for v1.83
 	_new(string)	return ref to new object from ref to decimal string
 	_zero()		return a new object with value 0
 	_one()		return a new object with value 1
@@ -2054,9 +2217,9 @@ Math::BigInt v1.70 or later:
 	_check(obj)	check whether internal representation is still intact
 			return 0 for ok, otherwise error message as string
 
-	_from_hex(str)	return ref to new object from ref to hexadecimal string
-	_from_bin(str)	return ref to new object from ref to binary string
-	_from_oct(str)	return ref to new object from ref to octal string
+	_from_hex(str)	return new object from a hexadecimal string
+	_from_bin(str)	return new object from a binary string
+	_from_oct(str)	return new object from an octal string
 	
 	_as_hex(str)	return string containing the value as
 			unsigned hex string, with the '0x' prepended.
@@ -2073,11 +2236,11 @@ Math::BigInt v1.70 or later:
 	_and(obj1,obj2)	AND (bit-wise) object 1 with object 2
 	_or(obj1,obj2)	OR (bit-wise) object 1 with object 2
 
-	_mod(obj,obj)	Return remainder of div of the 1st by the 2nd object
+	_mod(obj1,obj2)	Return remainder of div of the 1st by the 2nd object
 	_sqrt(obj)	return the square root of object (truncated to int)
 	_root(obj)	return the n'th (n >= 3) root of obj (truncated to int)
 	_fac(obj)	return factorial of object 1 (1*2*3*4..)
-	_pow(obj,obj)	return object 1 to the power of object 2
+	_pow(obj1,obj2)	return object 1 to the power of object 2
 			return undef for NaN
 	_zeros(obj)	return number of trailing decimal zeros
 	_modinv		return inverse modulus
@@ -2090,6 +2253,14 @@ Math::BigInt v1.70 or later:
 			 undef : unknown whether result is exactly RESULT
         _gcd(obj,obj)	return Greatest Common Divisor of two objects
 
+The following functions are REQUIRED for an api_version of 2 or greater:
+
+	_1ex($x)	create the number 1Ex where x >= 0
+	_alen(obj)	returns approximate count of the decimal digits of the
+			object. This estimate MUST always be greater or equal
+			to what _len() returns.
+        _nok(n,k)	calculate n over k (binomial coefficient)
+
 The following functions are optional, and can be defined if the underlying lib
 has a fast way to do them. If undefined, Math::BigInt will use pure Perl (hence
 slow) fallback routines to emulate these:
@@ -2098,7 +2269,6 @@ slow) fallback routines to emulate these:
 	_signed_and
 	_signed_xor
 
-
 Input strings come in as unsigned but with prefix (i.e. as '123', '0xabc'
 or '0b1101').
 
@@ -2109,7 +2279,7 @@ cases.
 
 The first parameter can be modified, that includes the possibility that you
 return a reference to a completely different object instead. Although keeping
-the reference and just changing it's contents is preferred over creating and
+the reference and just changing its contents is preferred over creating and
 returning a different reference.
 
 Return values are always references to objects, strings, or true/false for
@@ -2145,7 +2315,7 @@ Fixed, speed-up, streamlined and enhanced by Tels 2001 - 2007.
 
 =head1 SEE ALSO
 
-L<Math::BigInt>, L<Math::BigFloat>, L<Math::BigInt::BitVect>,
+L<Math::BigInt>, L<Math::BigFloat>,
 L<Math::BigInt::GMP>, L<Math::BigInt::FastCalc> and L<Math::BigInt::Pari>.
 
 =cut
diff --git a/lib/Math/BigInt/t/bare_mbf.t b/lib/Math/BigInt/t/bare_mbf.t
index 3a167e2021..4d471e2cf7 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 => 2042;
+  plan tests => 2064;
   }
 
 use Math::BigFloat lib => 'BareCalc';
diff --git a/lib/Math/BigInt/t/bare_mbi.t b/lib/Math/BigInt/t/bare_mbi.t
index 13512854ea..edb482e531 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 => 3055;
+  plan tests => 3093;
   }
 
 use Math::BigInt lib => 'BareCalc';
@@ -37,7 +37,7 @@ use vars qw ($class $try $x $y $f @args $ans $ans1 $ans1_str $setup $CL);
 $class = "Math::BigInt";
 $CL = "Math::BigInt::BareCalc";
 
-my $version = '1.61';	# for $VERSION tests, match current release (by hand!)
+my $version = '1.84';	# for $VERSION tests, match current release (by hand!)
 
 require 'bigintpm.inc';	# perform same tests as bigintpm
 
diff --git a/lib/Math/BigInt/t/bigfltpm.inc b/lib/Math/BigInt/t/bigfltpm.inc
index 45f48acc5b..9006afe1ca 100644
--- a/lib/Math/BigInt/t/bigfltpm.inc
+++ b/lib/Math/BigInt/t/bigfltpm.inc
@@ -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
diff --git a/lib/Math/BigInt/t/bigfltpm.t b/lib/Math/BigInt/t/bigfltpm.t
index 65f791d9fa..0c32908cb7 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 => 2042
+  plan tests => 2064
 	+ 3;		# own tests
   }
 
diff --git a/lib/Math/BigInt/t/bigintc.t b/lib/Math/BigInt/t/bigintc.t
index 89b515fcf2..44c43645df 100644
--- a/lib/Math/BigInt/t/bigintc.t
+++ b/lib/Math/BigInt/t/bigintc.t
@@ -13,7 +13,7 @@ BEGIN
     print "1..0\n";
     exit(0);
     }
-  plan tests => 322;
+  plan tests => 375;
   }
 
 use Math::BigInt::Calc;
@@ -97,31 +97,44 @@ ok ($C->_is_odd($C->_one()),1); ok ($C->_is_odd($C->_zero())||0,0);
 ok ($C->_is_even($C->_one()) || 0,0); ok ($C->_is_even($C->_zero()),1);
 
 # _len
-$x = $C->_new("1"); ok ($C->_len($x),1);
-$x = $C->_new("12"); ok ($C->_len($x),2);
-$x = $C->_new("123"); ok ($C->_len($x),3);
-$x = $C->_new("1234"); ok ($C->_len($x),4);
-$x = $C->_new("12345"); ok ($C->_len($x),5);
-$x = $C->_new("123456"); ok ($C->_len($x),6);
-$x = $C->_new("1234567"); ok ($C->_len($x),7);
-$x = $C->_new("12345678"); ok ($C->_len($x),8);
-$x = $C->_new("123456789"); ok ($C->_len($x),9);
-
-$x = $C->_new("8"); ok ($C->_len($x),1);
-$x = $C->_new("21"); ok ($C->_len($x),2);
-$x = $C->_new("321"); ok ($C->_len($x),3);
-$x = $C->_new("4321"); ok ($C->_len($x),4);
-$x = $C->_new("54321"); ok ($C->_len($x),5);
-$x = $C->_new("654321"); ok ($C->_len($x),6);
-$x = $C->_new("7654321"); ok ($C->_len($x),7);
-$x = $C->_new("87654321"); ok ($C->_len($x),8);
-$x = $C->_new("987654321"); ok ($C->_len($x),9);
-
-for (my $i = 1; $i < 9; $i++)
+for my $method (qw/_alen _len/)
   {
-  my $a = "$i" . '0' x ($i-1);
-  $x = $C->_new($a); 
-  print "# Tried len '$a'\n" unless ok ($C->_len($x),$i);
+  $x = $C->_new("1"); ok ($C->$method($x),1);
+  $x = $C->_new("12"); ok ($C->$method($x),2);
+  $x = $C->_new("123"); ok ($C->$method($x),3);
+  $x = $C->_new("1234"); ok ($C->$method($x),4);
+  $x = $C->_new("12345"); ok ($C->$method($x),5);
+  $x = $C->_new("123456"); ok ($C->$method($x),6);
+  $x = $C->_new("1234567"); ok ($C->$method($x),7);
+  $x = $C->_new("12345678"); ok ($C->$method($x),8);
+  $x = $C->_new("123456789"); ok ($C->$method($x),9);
+
+  $x = $C->_new("8"); ok ($C->$method($x),1);
+  $x = $C->_new("21"); ok ($C->$method($x),2);
+  $x = $C->_new("321"); ok ($C->$method($x),3);
+  $x = $C->_new("4321"); ok ($C->$method($x),4);
+  $x = $C->_new("54321"); ok ($C->$method($x),5);
+  $x = $C->_new("654321"); ok ($C->$method($x),6);
+  $x = $C->_new("7654321"); ok ($C->$method($x),7);
+  $x = $C->_new("87654321"); ok ($C->$method($x),8);
+  $x = $C->_new("987654321"); ok ($C->$method($x),9);
+
+  $x = $C->_new("0"); ok ($C->$method($x),1);
+  $x = $C->_new("20"); ok ($C->$method($x),2);
+  $x = $C->_new("320"); ok ($C->$method($x),3);
+  $x = $C->_new("4320"); ok ($C->$method($x),4);
+  $x = $C->_new("54320"); ok ($C->$method($x),5);
+  $x = $C->_new("654320"); ok ($C->$method($x),6);
+  $x = $C->_new("7654320"); ok ($C->$method($x),7);
+  $x = $C->_new("87654320"); ok ($C->$method($x),8);
+  $x = $C->_new("987654320"); ok ($C->$method($x),9);
+
+  for (my $i = 1; $i < 9; $i++)
+    {
+    my $a = "$i" . '0' x ($i-1);
+    $x = $C->_new($a); 
+    print "# Tried len '$a'\n" unless ok ($C->_len($x),$i);
+    }
   }
 
 # _digit
@@ -275,10 +288,12 @@ for my $i (2 .. 9)
   print "# _pow( ", '9' x $i, ", 2) \n" unless
    ok ($C->_str($C->_pow($x,$n)),$rc);
  
-  if ($i <= 7)
+  # if $i > $BASE_LEN, the test takes a really long time:
+  if ($i <= $BASE_LEN)
     {
     $x = '9' x $i; $x = $C->_new($x);
     $n = '9' x $i; $n = $C->_new($n);
+    print "# _root( ", '9' x $i, ", ", 9 x $i, ") \n";
     print "# _root( ", '9' x $i, ", ", 9 x $i, ") \n" unless
      ok ($C->_str($C->_root($x,$n)),'1');
 
@@ -287,6 +302,11 @@ for my $i (2 .. 9)
     print "# _root( ", '9' x $i, ", ", 9 x $i, ") \n" unless
      ok ($C->_str($C->_root($x,$n)), $res->[$i-2]);
     }
+  else
+    {
+    ok ("skipped $i", "skipped $i");
+    ok ("skipped $i", "skipped $i");
+    } 
   }
 
 ##############################################################################
@@ -395,6 +415,13 @@ ok ($C->_as_hex( $C->_new("12")), '0xc');
 ok ($C->_as_bin( $C->_new("12")), '0b1100');
 ok ($C->_as_oct( $C->_new("64")), '0100');
 
+# _1ex
+ok ($C->_str($C->_1ex(0)), "1");
+ok ($C->_str($C->_1ex(1)), "10");
+ok ($C->_str($C->_1ex(2)), "100");
+ok ($C->_str($C->_1ex(12)), "1000000000000");
+ok ($C->_str($C->_1ex(16)), "10000000000000000");
+
 # _check
 $x = $C->_new("123456789");
 ok ($C->_check($x),0);
diff --git a/lib/Math/BigInt/t/bigintpm.inc b/lib/Math/BigInt/t/bigintpm.inc
index c62d73e53d..49a68fc9be 100644
--- a/lib/Math/BigInt/t/bigintpm.inc
+++ b/lib/Math/BigInt/t/bigintpm.inc
@@ -89,6 +89,8 @@ while (<DATA>)
     $try .= '$m = $m->bstr(); $m = "NaN" if !defined $m;';
     $try .= '$e = $e->bstr(); $e = "NaN" if !defined $e;';
     $try .= '"$m,$e";';
+   }elsif ($f eq "bexp"){
+    $try .= "\$x->bexp();";
    } else {
     # binary ops
     $try .= "\$y = $class->new('$args[1]');";
@@ -154,6 +156,8 @@ while (<DATA>)
           {
           $try .= "\$x >> \$y;";
           }
+      }elsif ($f eq "bnok"){
+        $try .= "\$x->bnok(\$y);";
       }elsif ($f eq "broot"){
         $try .= "\$x->broot(\$y);";
       }elsif ($f eq "blog"){
@@ -2038,6 +2042,9 @@ NaNfac:NaN
 10:3628800
 11:39916800
 12:479001600
+20:2432902008176640000
+22:1124000727777607680000
+69:171122452428141311372468338881272839092270544893520369393648040923257279754140647424000000000000000
 &bpow
 abc:12:NaN
 12:abc:NaN
@@ -2189,7 +2196,7 @@ NaN:inf:NaN
 18446744073709551616:128:1
 # 213 ** 15
 84274086103068221283760416414557757:15:213
-# see t/bigroot for more tests
+# see t/bigroot.t for more tests
 &bsqrt
 145:12
 144:12
@@ -2219,6 +2226,29 @@ NaN:inf:NaN
 Nan:NaN
 +inf:inf
 -inf:NaN
+# see t/biglog.t for more tests
+&bexp
+NaN:NaN
+inf:inf
+1:2
+2:7
+# see t/bignok.t for more tests
+&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
+100:95:75287520
 &bround
 $round_mode('trunc')
 0:12:0
diff --git a/lib/Math/BigInt/t/bigintpm.t b/lib/Math/BigInt/t/bigintpm.t
index 5b7302c853..606143c7b2 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 => 3055;
+  plan tests => 3093;
   }
 
 use Math::BigInt lib => 'Calc';
diff --git a/lib/Math/BigInt/t/biglog.t b/lib/Math/BigInt/t/biglog.t
index a3f2e62cac..6c99cb5535 100644
--- a/lib/Math/BigInt/t/biglog.t
+++ b/lib/Math/BigInt/t/biglog.t
@@ -37,7 +37,7 @@ BEGIN
     }
   print "# INC = @INC\n";
 
-  plan tests => 66;
+  plan tests => 68;
   }
 
 use Math::BigFloat;
@@ -171,13 +171,23 @@ test_bpow ('9.86902225','0.5',undef, '3.1415');
 test_bpow ('0.2','0.41',10,   '0.5169187652');
 
 #############################################################################
-# test bexp()
+# test bexp() with cached results
 
 is ($cl->new(1)->bexp(), '2.718281828459045235360287471352662497757', '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'); 
+
 # all done
 1;
 
diff --git a/lib/Math/BigInt/t/bigroot.t b/lib/Math/BigInt/t/bigroot.t
index e48d659803..41fee89970 100644
--- a/lib/Math/BigInt/t/bigroot.t
+++ b/lib/Math/BigInt/t/bigroot.t
@@ -8,7 +8,7 @@
 # But it is better to test the numerical functionality, instead of not testing
 # it at all.
 
-use Test;
+use Test::More;
 use strict;
 
 BEGIN
@@ -62,10 +62,8 @@ sub test_broot
   my ($x,$n,$y,$scale,$result) = @_;
 
   my $s = $scale || 'undef';
-  print "# Try: $cl $x->bpow($n)->broot($y,$s) == $result:\n"; 
-   ok ($cl->new($x)->bpow($n)->broot($y,$scale),$result);
+  is ($cl->new($x)->bpow($n)->broot($y,$scale),$result, "Try: $cl $x->bpow($n)->broot($y,$s) == $result");
   $result =~ s/\..*//;
-  print "# Try: $c $x->bpow($n)->broot($y,$s) == $result:\n"; 
-   ok ($c->new($x)->bpow($n)->broot($y,$scale),$result);
+  is ($c->new($x)->bpow($n)->broot($y,$scale),$result, "Try: $c $x->bpow($n)->broot($y,$s) == $result");
   }
 
diff --git a/lib/Math/BigInt/t/lib_load.t b/lib/Math/BigInt/t/lib_load.t
index 3aff7c4037..ff3972effa 100644
--- a/lib/Math/BigInt/t/lib_load.t
+++ b/lib/Math/BigInt/t/lib_load.t
@@ -26,7 +26,7 @@ BEGIN
     }
   print "# INC = @INC\n";
 
-  plan tests => 2;
+  plan tests => 4;
   }
 
 # first load BigInt with Calc
@@ -43,3 +43,11 @@ Math::BigInt->import( lib => 'BareCalc' );
 
 is (Math::BigFloat::config()->{lib}, 'Math::BigInt::BareCalc', 'BigFloat was notified');
 
+# See that Math::BigFloat supports "only"
+eval "Math::BigFloat->import('only' => 'Calc')";
+is (Math::BigFloat::config()->{lib}, 'Math::BigInt::Calc', '"only" worked');
+
+# See that Math::BigFloat supports "try"
+eval "Math::BigFloat->import('try' => 'BareCalc')";
+is (Math::BigFloat::config()->{lib}, 'Math::BigInt::BareCalc', '"try" worked');
+
diff --git a/lib/Math/BigInt/t/mbi_rand.t b/lib/Math/BigInt/t/mbi_rand.t
index d24920f68b..e2bf6637de 100644
--- a/lib/Math/BigInt/t/mbi_rand.t
+++ b/lib/Math/BigInt/t/mbi_rand.t
@@ -1,6 +1,6 @@
 #!/usr/bin/perl -w
 
-use Test;
+use Test::More;
 use strict;
 
 my $count;
@@ -57,7 +57,7 @@ for (my $i = 0; $i < $count; $i++)
   print "# A $A\n# B $B\n";
   if ($A->is_zero() || $B->is_zero())
     {
-    for (1..4) { ok (1,1); } next;
+    for (1..4) { is (1,1, 'skipped this test'); } next;
     }
 
   # check that int(A/B)*B + A % B == A holds for all inputs
@@ -69,8 +69,8 @@ for (my $i = 0; $i < $count; $i++)
 
   print "# seed $seed, ". join(' ',Math::BigInt::Calc->_base_len()),"\n".
         "# tried $ADB * $B + $two*$AMB - $AMB\n"
-   unless ok ($ADB*$B+$two*$AMB-$AMB,$As);
-  if (ok ($ADB*$B/$B,$ADB))
+   unless is ($ADB*$B+$two*$AMB-$AMB,$As, "ADB * B + 2 * AMB - AMB == A");
+  if (is ($ADB*$B/$B,$ADB, "ADB * B / B == ADB"))
     {
     print "# seed: $seed, \$ADB * \$B / \$B = ", $ADB * $B / $B, " != $ADB (\$B=$B)\n";
     if (Math::BigInt->config()->{lib} =~ /::Calc/)
@@ -84,10 +84,10 @@ for (my $i = 0; $i < $count; $i++)
   # print "check: $ADB $AMB";
   print "# seed $seed, ". join(' ',Math::BigInt::Calc->_base_len()),"\n".
         "# tried $ADB * $A + $two*$AMB - $AMB\n"
-   unless ok ($ADB*$A+$two*$AMB-$AMB,$Bs);
+   unless is ($ADB*$A+$two*$AMB-$AMB,$Bs, "ADB * A + 2 * AMB - AMB == B");
   print "# +$two * $AMB = ",$ADB * $A + $two * $AMB,"\n";
   print "# -$AMB = ",$ADB * $A + $two * $AMB - $AMB,"\n";
   print "# seed $seed, \$ADB * \$A / \$A = ", $ADB * $A / $A, " != $ADB (\$A=$A)\n"
-   unless ok ($ADB*$A/$A,$ADB);
+   unless is ($ADB*$A/$A,$ADB, "ADB * A/A == ADB");
   }
 
diff --git a/lib/Math/BigInt/t/sub_mbf.t b/lib/Math/BigInt/t/sub_mbf.t
index 53a3fcb737..f0761f89a2 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 => 2042
+  plan tests => 2064
     + 6;	# + our own tests
   }
 
diff --git a/lib/Math/BigInt/t/sub_mbi.t b/lib/Math/BigInt/t/sub_mbi.t
index 4cd9ff4d9c..798bfa2413 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 => 3055
+  plan tests => 3093
     + 5;	# +5 own tests
   }
 
diff --git a/lib/Math/BigInt/t/with_sub.t b/lib/Math/BigInt/t/with_sub.t
index 5601a207b6..a66b6e166c 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 => 2042
+  plan tests => 2064
 	+ 1;
   }