From: Tels <nospam-abuse@bloodgate.com>
Date: Mon, 9 Apr 2007 20:59:22 +0000 (+0000)
Subject: BigInt, FastCalc, BitRat, bignum released to CPAN [PATCH]
X-Git-Url: http://git.shadowcat.co.uk/gitweb/gitweb.cgi?a=commitdiff_plain;h=7d193e396ed9e1516565a568311b86ae5b3466a3;p=p5sagit%2Fp5-mst-13.2.git

BigInt, FastCalc, BitRat, bignum released to CPAN [PATCH]
Message-Id: <200704092059.24058@bloodgate.com>

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

diff --git a/ext/Math/BigInt/FastCalc/FastCalc.pm b/ext/Math/BigInt/FastCalc/FastCalc.pm
index 160a218..5b2ea2f 100644
--- a/ext/Math/BigInt/FastCalc/FastCalc.pm
+++ b/ext/Math/BigInt/FastCalc/FastCalc.pm
@@ -11,7 +11,7 @@ use vars qw/@ISA $VERSION $BASE $BASE_LEN/;
 
 @ISA = qw(DynaLoader);
 
-$VERSION = '0.12_01';
+$VERSION = '0.13';
 
 bootstrap Math::BigInt::FastCalc $VERSION;
 
diff --git a/ext/Math/BigInt/FastCalc/FastCalc.xs b/ext/Math/BigInt/FastCalc/FastCalc.xs
index 3e53876..b00ed05 100644
--- a/ext/Math/BigInt/FastCalc/FastCalc.xs
+++ b/ext/Math/BigInt/FastCalc/FastCalc.xs
@@ -18,6 +18,20 @@ MODULE = Math::BigInt::FastCalc		PACKAGE = Math::BigInt::FastCalc
  #  * added __strip_zeros(), _copy()
  # 2004-08-13 0.07 Tels
  #  * added _is_two(), _is_ten(), _ten()
+ # 2007-04-02 0.08 Tels
+ #  * plug leaks by creating mortals
+
+#define RETURN_MORTAL_INT(value)		\
+      ST(0) = sv_2mortal(newSViv(value));	\
+      XSRETURN(1);
+
+#define RETURN_MORTAL_BOOL(temp, comp)			\
+      ST(0) = sv_2mortal(boolSV( SvIV(temp) == comp));
+
+#define CONSTANT_OBJ(int)			\
+    RETVAL = newAV();				\
+    sv_2mortal((SV*)RETVAL);			\
+    av_push (RETVAL, newSViv( int ));
 
 void 
 _set_XS_BASE(BASE, BASE_LEN)
@@ -230,11 +244,6 @@ _num(class,x)
 
 ##############################################################################
 
-#define CONSTANT_OBJ(int)		\
-    RETVAL = newAV();			\
-    sv_2mortal((SV*)RETVAL);		\
-    av_push (RETVAL, newSViv( int ));
-
 AV *
 _zero(class)
   CODE:
@@ -281,7 +290,7 @@ _is_even(class, x)
   CODE:
     a = (AV*)SvRV(x);		/* ref to aray, don't check ref */
     temp = *av_fetch(a, 0, 0);	/* fetch first element */
-    ST(0) = boolSV((SvIV(temp) & 1) == 0);
+    ST(0) = sv_2mortal(boolSV((SvIV(temp) & 1) == 0));
 
 ##############################################################################
 
@@ -295,7 +304,7 @@ _is_odd(class, x)
   CODE:
     a = (AV*)SvRV(x);		/* ref to aray, don't check ref */
     temp = *av_fetch(a, 0, 0);	/* fetch first element */
-    ST(0) = boolSV((SvIV(temp) & 1) != 0);
+    ST(0) = sv_2mortal(boolSV((SvIV(temp) & 1) != 0));
 
 ##############################################################################
 
@@ -314,7 +323,7 @@ _is_one(class, x)
       XSRETURN(1);			/* len != 1, can't be '1' */
       }
     temp = *av_fetch(a, 0, 0);		/* fetch first element */
-    ST(0) = boolSV((SvIV(temp) == 1));
+    RETURN_MORTAL_BOOL(temp, 1);
 
 ##############################################################################
 
@@ -333,7 +342,7 @@ _is_two(class, x)
       XSRETURN(1);			/* len != 1, can't be '2' */
       }
     temp = *av_fetch(a, 0, 0);		/* fetch first element */
-    ST(0) = boolSV((SvIV(temp) == 2));
+    RETURN_MORTAL_BOOL(temp, 2);
 
 ##############################################################################
 
@@ -352,7 +361,7 @@ _is_ten(class, x)
       XSRETURN(1);			/* len != 1, can't be '10' */
       }
     temp = *av_fetch(a, 0, 0);		/* fetch first element */
-    ST(0) = boolSV((SvIV(temp) == 10));
+    RETURN_MORTAL_BOOL(temp, 10);
 
 ##############################################################################
 
@@ -371,7 +380,7 @@ _is_zero(class, x)
       XSRETURN(1);			/* len != 1, can't be '0' */
       }
     temp = *av_fetch(a, 0, 0);		/* fetch first element */
-    ST(0) = boolSV((SvIV(temp) == 0));
+    RETURN_MORTAL_BOOL(temp, 0);
 
 ##############################################################################
 
@@ -390,7 +399,7 @@ _len(class,x)
     temp = *av_fetch(a, elems, 0);	/* fetch last element */
     SvPV(temp, len);			/* convert to string & store length */
     len += (IV) XS_BASE_LEN * elems;
-    ST(0) = newSViv(len);
+    ST(0) = sv_2mortal(newSViv(len));
 
 ##############################################################################
 
@@ -418,13 +427,11 @@ _acmp(class, cx, cy);
 
     if (diff > 0)
       {
-      ST(0) = newSViv(1);		/* len differs: X > Y */
-      XSRETURN(1);
+      RETURN_MORTAL_INT(1);		/* len differs: X > Y */
       }
-    if (diff < 0)
+    else if (diff < 0)
       {
-      ST(0) = newSViv(-1);		/* len differs: X < Y */
-      XSRETURN(1);
+      RETURN_MORTAL_INT(-1);		/* len differs: X < Y */
       }
     /* both have same number of elements, so check length of last element
        and see if it differes */
@@ -435,13 +442,11 @@ _acmp(class, cx, cy);
     diff_str = (I32)lenx - (I32)leny;
     if (diff_str > 0)
       {
-      ST(0) = newSViv(1);		/* same len, but first elems differs in len */
-      XSRETURN(1);
+      RETURN_MORTAL_INT(1);		/* same len, but first elems differs in len */
       }
     if (diff_str < 0)
       {
-      ST(0) = newSViv(-1);		/* same len, but first elems differs in len */
-      XSRETURN(1);
+      RETURN_MORTAL_INT(-1);		/* same len, but first elems differs in len */
       }
     /* same number of digits, so need to make a full compare */
     diff_nv = 0;
@@ -458,13 +463,11 @@ _acmp(class, cx, cy);
       } 
     if (diff_nv > 0)
       {
-      ST(0) = newSViv(1);
-      XSRETURN(1);
+      RETURN_MORTAL_INT(1);
       }
     if (diff_nv < 0)
       {
-      ST(0) = newSViv(-1);
-      XSRETURN(1);
+      RETURN_MORTAL_INT(-1);
       }
-    ST(0) = newSViv(0);		/* equal */
+    ST(0) = sv_2mortal(newSViv(0));		/* X and Y are equal */
 
diff --git a/ext/Math/BigInt/FastCalc/t/leak.t b/ext/Math/BigInt/FastCalc/t/leak.t
index c7cae8b..b9cc596 100644
--- a/ext/Math/BigInt/FastCalc/t/leak.t
+++ b/ext/Math/BigInt/FastCalc/t/leak.t
@@ -1,6 +1,9 @@
 #!/usr/bin/perl -w
 
-# Test for memory leaks from _zero() and friends.
+# Test for memory leaks.
+
+# XXX TODO: This test file doesn't actually seem to work! If you remove
+# the sv_2mortal() in the XS file, it still happily passes all tests...
 
 use Test::More;
 use strict;
@@ -10,7 +13,7 @@ BEGIN
   $| = 1;
   chdir 't' if -d 't';
   unshift @INC, ('../lib', '../blib/arch');	# for running manually
-  plan tests => 4;
+  plan tests => 20;
   }
 
 #############################################################################
@@ -18,7 +21,7 @@ package Math::BigInt::FastCalc::LeakCheck;
 use base qw(Math::BigInt::FastCalc);
 
 my $destroyed = 0;
-sub DESTROY { $destroyed++ }
+sub DESTROY { $destroyed++; }
 
 #############################################################################
 package main;
@@ -32,3 +35,48 @@ for my $method (qw(_zero _one _two _ten))
     }
   is ($destroyed, 1, "$method does not leak memory");
   }
+
+my $num = Math::BigInt::FastCalc->_zero();
+for my $method (qw(_is_zero _is_one _is_two _is_ten _num))
+  {
+  $destroyed = 0;
+    {
+    my $rc = Math::BigInt::FastCalc->$method($num);
+    bless \$rc, "Math::BigInt::FastCalc::LeakCheck";
+    }
+  is ($destroyed, 1, "$method does not leak memory");
+  }
+
+my $num_10 = Math::BigInt::FastCalc->_ten();
+my $num_2 = Math::BigInt::FastCalc->_two();
+
+my $num_long   = Math::BigInt::FastCalc->_new("1234567890");
+my $num_long_2 = Math::BigInt::FastCalc->_new("12345678900987654321");
+
+# to hit all possible code branches
+_test_acmp($num, $num);
+_test_acmp($num_10, $num_10);
+_test_acmp($num, $num_10);
+_test_acmp($num_10, $num);
+_test_acmp($num, $num_2);
+_test_acmp($num_2, $num);
+_test_acmp($num_long, $num);
+_test_acmp($num, $num_long);
+_test_acmp($num_long, $num_long);
+_test_acmp($num_long, $num_long_2);
+_test_acmp($num_long_2, $num_long);
+
+sub _test_acmp
+  {
+  my ($n1,$n2) = @_;
+
+  $destroyed = 0;
+    {
+    my $rc = Math::BigInt::FastCalc->_acmp($n1,$n2);
+    bless \$rc, "Math::BigInt::FastCalc::LeakCheck";
+    }
+  my $n_1 = Math::BigInt::FastCalc->_str($n1);
+  my $n_2 = Math::BigInt::FastCalc->_str($n2);
+  is ($destroyed, 1, "_acmp($n_1,$n_2) does not leak memory");
+  }
+
diff --git a/lib/Math/BigFloat.pm b/lib/Math/BigFloat.pm
index f569036..1242a37 100644
--- a/lib/Math/BigFloat.pm
+++ b/lib/Math/BigFloat.pm
@@ -12,7 +12,7 @@ package Math::BigFloat;
 #   _a	: accuracy
 #   _p	: precision
 
-$VERSION = '1.53';
+$VERSION = '1.54';
 require 5.006002;
 
 require Exporter;
@@ -89,7 +89,7 @@ sub STORE       { $rnd_mode = $_[0]->round_mode($_[1]); }
 
 BEGIN
   {
-  # when someone set's $rnd_mode, we catch this and check the value to see
+  # when someone sets $rnd_mode, we catch this and check the value to see
   # whether it is valid or not. 
   $rnd_mode   = 'even'; tie $rnd_mode, 'Math::BigFloat';
 
@@ -103,8 +103,8 @@ BEGIN
   # valid method aliases for AUTOLOAD
   my %methods = map { $_ => 1 }  
    qw / fadd fsub fmul fdiv fround ffround fsqrt fmod fstr fsstr fpow fnorm
-        fint facmp fcmp fzero fnan finf finc fdec flog ffac fneg
-	fceil ffloor frsft flsft fone flog froot
+        fint facmp fcmp fzero fnan finf finc fdec ffac fneg
+	fceil ffloor frsft flsft fone flog froot fexp
       /;
   # valid methods that can be handed up (for AUTOLOAD)
   my %hand_ups = map { $_ => 1 }  
@@ -812,6 +812,7 @@ sub blog
   local $Math::BigFloat::downgrade = undef;
 
   # upgrade $x if $x is not a BigFloat (handle BigInt input)
+  # XXX TODO: rebless!
   if (!$x->isa('Math::BigFloat'))
     {
     $x = Math::BigFloat->new($x);
@@ -844,7 +845,8 @@ sub blog
 
   if ($done == 0)
     {
-    # first calculate the log to base e (using reduction by 10 (and probably 2))
+    # base is undef, so base should be e (Euler's number), so first calculate the
+    # log to base e (using reduction by 10 (and probably 2)):
     $self->_log_10($x,$scale);
 
     # and if a different base was requested, convert it
@@ -876,6 +878,153 @@ sub blog
   $x;
   }
 
+sub bexp
+  {
+  # Calculate e ** X (Euler's constant to the power of X)
+  my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
+
+  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();
+
+  if (!$x->isa('Math::BigFloat'))
+    {
+    $x = Math::BigFloat->new($x);
+    $self = ref($x);
+    }
+  
+  # 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;
+  local $Math::BigFloat::downgrade = undef;
+
+  my $x_org = $x->copy();
+  delete $x_org->{_a}; delete $x_org->{_p};
+
+  # We use the following Taylor series:
+
+  #           x    x^2   x^3   x^4
+  #  e = 1 + --- + --- + --- + --- ...
+  #           1!    2!    3!    4!
+
+  # The difference for each term is X and N, which would result in:
+  # 2 copy, 2 mul, 2 add, 1 inc, 1 div operations per term
+
+  # But it is faster to compute exp(1) and then raising it to the
+  # given power, esp. if $x is really big and an integer because:
+
+  #  * The numerator is always 1, making the computation faster
+  #  * the series converges faster in the case of x == 1
+  #  * We can also easily check when we have reached our limit: when the
+  #    term to be added is smaller than "1E$scale", we can stop - f.i.
+  #    scale == 5, and we have 1/40320, then we stop since 1/40320 < 1E-5.
+  #  * we can compute the *exact* result by simulating bigrat math:
+
+  #  1   1    gcd(3,4) = 1    1*24 + 1*6    5
+  #  - + -                  = ---------- =  --                 
+  #  6   24                      6*24       24
+
+  # We do not compute the gcd() here, but simple do:
+  #  1   1    1*24 + 1*6   30
+  #  - + -  = --------- =  --                 
+  #  6   24       6*24     144
+
+  # In general:
+  #  a   c    a*d + c*b 	and note that c is always 1 and d = (b*f)
+  #  - + -  = ---------
+  #  b   d       b*d
+
+  # This leads to:         which can be reduced by b to:
+  #  a   1     a*b*f + b    a*f + 1
+  #  - + -   = --------- =  -------
+  #  b   b*f     b*b*f        b*f
+
+  # The first terms in the series are:
+
+  # 1     1    1    1    1    1     1     1     13700
+  # -- + -- + -- + -- + -- + --- + --- + ---- = -----
+  # 1     1    2    6   24   120   720   5040   5040
+
+  # Note that we cannot simple reduce 13700/5040 to 685/252.
+
+  # 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)
+    {
+    # 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);
+    }
+
+  $x = $self->new($MBI->_str($A))->bdiv($MBI->_str($B), $scale);
+
+  delete $x->{_a}; delete $x->{_p};
+  # raise $x to the wanted power
+  $x->bpow($x_org, $scale) unless $x_org->is_one();
+
+  # 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;						# return modified $x
+  }
+
 sub _log
   {
   # internal log function to calculate ln() based on Taylor series.
@@ -885,6 +1034,8 @@ sub _log
   # in case of $x == 1, result is 0
   return $x->bzero() if $x->is_one();
 
+  # XXX TODO: rewrite this in a similiar manner to bexp()
+
   # http://www.efunda.com/math/taylor_series/logarithmic.cfm?search_string=log
 
   # u = x-1, v = x+1
@@ -931,7 +1082,7 @@ sub _log
     # if we truncated $over and $below we might get 0.12345. Does this matter
     # for the end result? So we give $over and $below 4 more digits to be
     # on the safe side (unscientific error handling as usual... :+D
-    
+
     $next = $over->copy->bround($scale+4)->bdiv(
       $below->copy->bmul($factor)->bround($scale+4), 
       $scale);
@@ -950,8 +1101,8 @@ sub _log
       $steps++; print "step $steps = $x\n" if $steps % 10 == 0;
       }
     }
-  $x->bmul($f);					# $x *= 2
   print "took $steps steps\n" if DEBUG;
+  $x->bmul($f);					# $x *= 2
   }
 
 sub _log_10
@@ -960,19 +1111,27 @@ sub _log_10
   # and then "correcting" the result to the proper one. Modifies $x in place.
   my ($self,$x,$scale) = @_;
 
-  # taking blog() from numbers greater than 10 takes a *very long* time, so we
+  # Taking blog() from numbers greater than 10 takes a *very long* time, so we
   # break the computation down into parts based on the observation that:
-  #  blog(x*y) = blog(x) + blog(y)
-  # We set $y here to multiples of 10 so that $x is below 1 (the smaller $x is
-  # the faster it get's, especially because 2*$x takes about 10 times as long,
-  # so by dividing $x by 10 we make it at least factor 100 faster...)
-
-  # The same observation is valid for numbers smaller than 0.1 (e.g. computing
-  # log(1) is fastest, and the farther away we get from 1, the longer it takes)
-  # so we also 'break' this down by multiplying $x with 10 and subtract the
+  #  blog(X*Y) = blog(X) + blog(Y)
+  # We set Y here to multiples of 10 so that $x becomes below 1 - the smaller
+  # $x is the faster it gets. Since 2*$x takes about 10 times as
+  # long, we make it faster by about a factor of 100 by dividing $x by 10.
+
+  # The same observation is valid for numbers smaller than 0.1, e.g. computing
+  # log(1) is fastest, and the further away we get from 1, the longer it takes.
+  # So we also 'break' this down by multiplying $x with 10 and subtract the
   # log(10) afterwards to get the correct result.
 
-  # calculate nr of digits before dot
+  # To get $x even closer to 1, we also divide by 2 and then use log(2) to
+  # correct for this. For instance if $x is 2.4, we use the formula:
+  #  blog(2.4 * 2) == blog (1.2) + blog(2)
+  # and thus calculate only blog(1.2) and blog(2), which is faster in total
+  # than calculating blog(2.4).
+
+  # In addition, the values for blog(2) and blog(10) are cached.
+
+  # Calculate nr of digits before dot:
   my $dbd = $MBI->_num($x->{_e});
   $dbd = -$dbd if $x->{_es} eq '-';
   $dbd += $MBI->_len($x->{_m});
@@ -990,9 +1149,10 @@ sub _log_10
     # we can use the cached value in these cases
     if ($scale <= $LOG_10_A)
       {
-      $x->bzero(); $x->badd($LOG_10);
+      $x->bzero(); $x->badd($LOG_10);		# modify $x in place
       $calc = 0; 				# no need to calc, but round
       }
+    # if we can't use the shortcut, we continue normally
     }
   else
     {
@@ -1003,9 +1163,10 @@ sub _log_10
       # we can use the cached value in these cases
       if ($scale <= $LOG_2_A)
         {
-        $x->bzero(); $x->badd($LOG_2);
+        $x->bzero(); $x->badd($LOG_2);		# modify $x in place
         $calc = 0; 				# no need to calc, but round
         }
+      # if we can't use the shortcut, we continue normally
       }
     }
 
@@ -1028,6 +1189,8 @@ sub _log_10
   my $l_10;				# value of ln(10) to A of $scale
   my $l_2;				# value of ln(2) to A of $scale
 
+  my $two = $self->new(2);
+
   # $x == 2 => 1, $x == 13 => 2, $x == 0.1 => 0, $x == 0.01 => -1
   # so don't do this shortcut for 1 or 0
   if (($dbd > 1) || ($dbd < 0))
@@ -1049,10 +1212,40 @@ sub _log_10
       }
     else
       {
-      # else: slower, compute it (but don't cache it, because it could be big)
+      # else: slower, compute and cache result
       # also disable downgrade for this code path
       local $Math::BigFloat::downgrade = undef;
-      $l_10 = $self->new(10)->blog(undef,$scale);	# scale+4, actually
+
+      # shorten the time to calculate log(10) based on the following:
+      # log(1.25 * 8) = log(1.25) + log(8)
+      #               = log(1.25) + log(2) + log(2) + log(2)
+
+      # first get $l_2 (and possible compute and cache log(2))
+      $LOG_2 = $self->new($LOG_2,undef,undef) unless ref $LOG_2;
+      if ($scale <= $LOG_2_A)
+        {
+        # use cached value
+        $l_2 = $LOG_2->copy();			# copy() for the mul below
+        }
+      else
+        {
+        # else: slower, compute and cache result
+        $l_2 = $two->copy(); $self->_log($l_2, $scale); # scale+4, actually
+        $LOG_2 = $l_2->copy();			# cache the result for later
+						# the copy() is for mul below
+        $LOG_2_A = $scale;
+        }
+
+      # now calculate log(1.25):
+      $l_10 = $self->new('1.25'); $self->_log($l_10, $scale); # scale+4, actually
+
+      # log(1.25) + log(2) + log(2) + log(2):
+      $l_10->badd($l_2);
+      $l_10->badd($l_2);
+      $l_10->badd($l_2);
+      $LOG_10 = $l_10->copy();		# cache the result for later
+					# the copy() is for mul below
+      $LOG_10_A = $scale;
       }
     $dbd-- if ($dbd > 1); 		# 20 => dbd=2, so make it dbd=1	
     $l_10->bmul( $self->new($dbd));	# log(10) * (digits_before_dot-1)
@@ -1075,31 +1268,33 @@ sub _log_10
   $HALF = $self->new($HALF) unless ref($HALF);
 
   my $twos = 0;				# default: none (0 times)	
-  my $two = $self->new(2);
-  while ($x->bacmp($HALF) <= 0)
+  while ($x->bacmp($HALF) <= 0)		# X <= 0.5
     {
     $twos--; $x->bmul($two);
     }
-  while ($x->bacmp($two) >= 0)
+  while ($x->bacmp($two) >= 0)		# X >= 2
     {
     $twos++; $x->bdiv($two,$scale+4);		# keep all digits
     }
-  # $twos > 0 => did mul 2, < 0 => did div 2 (never both)
-  # calculate correction factor based on ln(2)
+  # $twos > 0 => did mul 2, < 0 => did div 2 (but we never did both)
+  # So calculate correction factor based on ln(2):
   if ($twos != 0)
     {
     $LOG_2 = $self->new($LOG_2,undef,undef) unless ref $LOG_2;
     if ($scale <= $LOG_2_A)
       {
       # use cached value
-      $l_2 = $LOG_2->copy();			# copy for mul
+      $l_2 = $LOG_2->copy();			# copy() for the mul below
       }
     else
       {
-      # else: slower, compute it (but don't cache it, because it could be big)
+      # else: slower, compute and cache result
       # also disable downgrade for this code path
       local $Math::BigFloat::downgrade = undef;
-      $l_2 = $two->blog(undef,$scale);	# scale+4, actually
+      $l_2 = $two->copy(); $self->_log($l_2, $scale); # scale+4, actually
+      $LOG_2 = $l_2->copy();			# cache the result for later
+						# the copy() is for mul below
+      $LOG_2_A = $scale;
       }
     $l_2->bmul($twos);		# * -2 => subtract, * 2 => add
     }
@@ -1107,7 +1302,9 @@ sub _log_10
   $self->_log($x,$scale);			# need to do the "normal" way
   $x->badd($l_10) if defined $l_10; 		# correct it by ln(10)
   $x->badd($l_2) if defined $l_2;		# and maybe by ln(2)
+
   # all done, $x contains now the result
+  $x;
   }
 
 sub blcm 
@@ -1919,7 +2116,7 @@ sub _pow
     {
     # 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
+    # anymore, so we stop:
     $next = $over->copy()->bdiv($below,$scale);
     last if $next->bacmp($limit) <= 0;
     $x->badd($next);
@@ -2627,6 +2824,7 @@ Math::BigFloat - Arbitrary size floating point math package
   
   $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
   
   $x->band($y);			# bit-wise and
   $x->bior($y);			# bit-wise inclusive or
diff --git a/lib/Math/BigInt.pm b/lib/Math/BigInt.pm
index 600970f..badf947 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.80';
+$VERSION = '1.82';
 
 @ISA = qw(Exporter);
 @EXPORT_OK = qw(objectify bgcd blcm); 
@@ -81,10 +81,9 @@ use overload
                "$_[1]" cmp $_[0]->bstr() :
                $_[0]->bstr() cmp "$_[1]" },
 
-# make cos()/sin()/exp() "work" with BigInt's or subclasses
+# make cos()/sin()/atan2() "work" with BigInt's or subclasses
 'cos'	=>	sub { cos($_[0]->numify()) }, 
 'sin'	=>	sub { sin($_[0]->numify()) }, 
-'exp'	=>	sub { exp($_[0]->numify()) }, 
 'atan2'	=>	sub { $_[2] ?
 			atan2($_[1],$_[0]->numify()) :
 			atan2($_[0]->numify(),$_[1]) },
@@ -95,6 +94,7 @@ use overload
 
 # log(N) is log(N, e), where e is Euler's number
 'log'	=>	sub { $_[0]->copy()->blog($_[1], undef); }, 
+'exp'	=>	sub { $_[0]->copy()->bexp($_[1]); }, 
 'int'	=>	sub { $_[0]->copy(); }, 
 'neg'	=>	sub { $_[0]->copy()->bneg(); }, 
 'abs'	=>	sub { $_[0]->copy()->babs(); },
@@ -1145,6 +1145,7 @@ sub bsub
   
   # set up parameters
   my ($self,$x,$y,@r) = (ref($_[0]),@_);
+
   # objectify is costly, so avoid it
   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
     {
@@ -1264,6 +1265,37 @@ sub blog
   $x->round(@r);
   }
 
+sub bexp
+  {
+  # Calculate e ** $x (Euler's number to the power of X), truncated to
+  # an integer value.
+  my ($self,$x,@r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
+  return $x if $x->modify('bexp');
+
+  # inf, -inf, NaN, <0 => NaN
+  return $x->bnan() if $x->{sign} eq 'NaN';
+  return $x->bone() if $x->is_zero();
+  return $x if $x->{sign} eq '+inf';
+  return $x->bzero() if $x->{sign} eq '-inf';
+
+  my $u;
+  {
+    # run through Math::BigFloat unless told otherwise
+    local $upgrade = 'Math::BigFloat' unless defined $upgrade;
+    # calculate result, truncate it to integer
+    $u = $upgrade->bexp($upgrade->new($x),@r);
+  }
+
+  if (!defined $upgrade)
+    {
+    $u = $u->as_int();
+    # modify $x in place
+    $x->{value} = $u->{value};
+    $x->round(@r);
+    }
+  else { $x = $u; }
+  }
+
 sub blcm 
   { 
   # (BINT or num_str, BINT or num_str) return BINT
@@ -2055,7 +2087,8 @@ sub exponent
     }
   return $self->bone() if $x->is_zero();
 
-  $self->new($x->_trailing_zeros());
+  # 12300 => 2 trailing zeros => exponent is 2
+  $self->new( $CALC->_zeros($x->{value}) );
   }
 
 sub mantissa
@@ -2069,8 +2102,9 @@ sub mantissa
     return $self->new($x->{sign});
     }
   my $m = $x->copy(); delete $m->{_p}; delete $m->{_a};
+
   # that's a bit inefficient:
-  my $zeros = $m->_trailing_zeros();
+  my $zeros = $CALC->_zeros($m->{value});
   $m->brsft($zeros,10) if $zeros != 0;
   $m;
   }
@@ -2343,7 +2377,7 @@ sub objectify
     }
 
   my $up = ${"$a[0]::upgrade"};
-  #print "Now in objectify, my class is today $a[0], count = $count\n";
+  # print STDERR "# Now in objectify, my class is today $a[0], count = $count\n";
   if ($count == 0)
     {
     while (@_)
@@ -2366,7 +2400,7 @@ sub objectify
     while ($count > 0)
       {
       $count--; 
-      $k = shift; 
+      $k = shift;
       if (!ref($k))
         {
         $k = $a[0]->new($k);
@@ -2830,8 +2864,8 @@ Math::BigInt - Arbitrary size integer/float math package
   $x->bmodinv($mod);	   # the inverse of $x in the given modulus $mod
 
   $x->bpow($y);		   # power of arguments (x ** y)
-  $x->blsft($y);	   # left shift in base 10
-  $x->brsft($y);	   # right shift in base 10
+  $x->blsft($y);	   # left shift in base 2
+  $x->brsft($y);	   # right shift in base 2
 			   # returns (quo,rem) or quo if in scalar context
   $x->blsft($y,$n);	   # left shift by $y places in base $n
   $x->brsft($y,$n);	   # right shift by $y places in base $n
@@ -2846,6 +2880,10 @@ 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->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
+  
   $x->round($A,$P,$mode);  # round to accuracy or precision using mode $mode
   $x->bround($n);	   # accuracy: preserve $n digits
   $x->bfround($n);	   # round to $nth digit, no-op for BigInts
@@ -3362,14 +3400,32 @@ is exactly equivalent to
 
 	$x->bpow($y);			# power of arguments (x ** y)
 
+=head2 blog()
+
+	$x->blog($base, $accuracy);	# logarithm of x to the base $base
+
+If C<$base> is not defined, Euler's number (e) is used:
+
+	print $x->blog(undef, 100);	# log(x) to 100 digits
+
+=head2 bexp()
+
+	$x->bexp($accuracy);		# calculate e ** X
+
+Calculates the expression C<e ** $x> where C<e> is Euler's number.
+
+This method was added in v1.82 of Math::BigInt (April 2007).
+
+See also L<blog()>.
+
 =head2 blsft()
 
-	$x->blsft($y);		# left shift
+	$x->blsft($y);		# left shift in base 2
 	$x->blsft($y,$n);	# left shift, in base $n (like 10)
 
 =head2 brsft()
 
-	$x->brsft($y);		# right shift 
+	$x->brsft($y);		# right shift in base 2
 	$x->brsft($y,$n);	# right shift, in base $n (like 10)
 
 =head2 band()
@@ -4240,6 +4296,8 @@ is in effect, they will always hand up their work:
 
 =item blog()
 
+=item bexp()
+
 =back
 
 Beware: This list is not complete.
diff --git a/lib/Math/BigInt/t/biglog.t b/lib/Math/BigInt/t/biglog.t
index 0958ddc..a3f2e62 100644
--- a/lib/Math/BigInt/t/biglog.t
+++ b/lib/Math/BigInt/t/biglog.t
@@ -1,6 +1,6 @@
 #!/usr/bin/perl -w
 
-# Test blog function (and bpow, since it uses blog).
+# Test blog function (and bpow, since it uses blog), as well as bexp().
 
 # It is too slow to be simple included in bigfltpm.inc, where it would get
 # executed 3 times. One time would be under BareCalc, which shouldn't make any
@@ -11,7 +11,7 @@
 # it at all (which did lead to wrong answers for 0 < $x < 1 in blog() in
 # versions up to v1.63, and for bsqrt($x) when $x << 1 for instance).
 
-use Test;
+use Test::More;
 use strict;
 
 BEGIN
@@ -37,7 +37,7 @@ BEGIN
     }
   print "# INC = @INC\n";
 
-  plan tests => 56;
+  plan tests => 66;
   }
 
 use Math::BigFloat;
@@ -45,16 +45,48 @@ use Math::BigInt;
 
 my $cl = "Math::BigInt";
 
+#############################################################################
 # test log($n) in BigInt (broken until 1.80)
 
-ok ($cl->new(2)->blog(), '0');
-ok ($cl->new(288)->blog(), '5');
-ok ($cl->new(2000)->blog(), '7');
+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) in BigInt
+
+is ($cl->new(1)->bexp(), '2', "bexp(1)");
+is ($cl->new(2)->bexp(), '7',"bexp(2)");
+is ($cl->new(3)->bexp(), '20', "bexp(3)");
 
 #############################################################################
+#############################################################################
+# BigFloat tests
+
+#############################################################################
+# test log(2, N) where N > 67 (broken until 1.82)
 
 $cl = "Math::BigFloat";
 
+# These tests can take quite a while, but are nec. Maybe protect them with
+# some alarm()?
+
+# this triggers the calculation and caching of ln(2):
+ok ($cl->new(5)->blog(undef,71), 
+'1.6094379124341003746007593332261876395256013542685177219126478914741790');
+
+# if the cache was correct, we should get this result, fast:
+ok ($cl->new(2)->blog(undef,71), 
+'0.69314718055994530941723212145817656807550013436025525412068000949339362');
+
+ok ($cl->new(10)->blog(undef,71), 
+'2.3025850929940456840179914546843642076011014886287729760333279009675726');
+
+ok ($cl->new(21)->blog(undef,71), 
+'3.0445224377234229965005979803657054342845752874046106401940844835750742');
+
+#############################################################################
+
 # These tests are now really fast, since they collapse to blog(10), basically
 # Don't attempt to run them with older versions. You are warned.
 
@@ -107,12 +139,12 @@ ok ($cl->new('1.2')->bpow('0.3',10),  '1.056219968');
 ok ($cl->new('10')->bpow('0.6',10),   '3.981071706');
 
 # blog should handle bigint input
-ok (Math::BigFloat::blog(Math::BigInt->new(100),10), 2);
+is (Math::BigFloat::blog(Math::BigInt->new(100),10), 2, "blog(100)");
 
 # some integer results
-ok ($cl->new(2)->bpow(32)->blog(2),  '32');	# 2 ** 32
-ok ($cl->new(3)->bpow(32)->blog(3),  '32');	# 3 ** 32
-ok ($cl->new(2)->bpow(65)->blog(2),  '65');	# 2 ** 65
+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");
 
 # test for bug in bsqrt() not taking negative _e into account
 test_bpow ('200','0.5',10,      '14.14213562');
@@ -138,6 +170,18 @@ test_bpow ('9.86902225','0.5',undef, '3.1415');
 
 test_bpow ('0.2','0.41',10,   '0.5169187652');
 
+#############################################################################
+# test bexp()
+
+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)'); 
+
+# all done
+1;
+
+#############################################################################
 sub test_bpow
   {
   my ($x,$y,$scale,$result) = @_;
@@ -146,3 +190,4 @@ sub test_bpow
    unless ok ($cl->new($x)->bpow($y,$scale),$result);
   }
 
+
diff --git a/lib/Math/BigInt/t/fallback.t b/lib/Math/BigInt/t/fallback.t
index 00f1dfd..010ff8a 100644
--- a/lib/Math/BigInt/t/fallback.t
+++ b/lib/Math/BigInt/t/fallback.t
@@ -1,6 +1,6 @@
 #!/usr/bin/perl -w
 
-# test 'fallback' for overload cos/sin/atan2/exp
+# test 'fallback' for overload cos/sin/atan2
 
 use Test;
 use strict;
@@ -28,7 +28,7 @@ BEGIN
     }
   print "# INC = @INC\n";
 
-  plan tests => 12;
+  plan tests => 9;
   }
 
 # The tests below test that cos(BigInt) = cos(Scalar) which is DWIM, but not
@@ -43,19 +43,16 @@ my $bi = Math::BigInt->new(1);
 
 ok (cos($bi), cos(1));
 ok (sin($bi), sin(1));
-ok (exp($bi), exp(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 (exp($bf), exp(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));
-ok (exp($bone), exp(1));