NetWare update from Ananth Kesari.
[p5sagit/p5-mst-13.2.git] / lib / Math / BigInt.pm
index 3c142f2..591973e 100644 (file)
@@ -18,7 +18,7 @@ package Math::BigInt;
 my $class = "Math::BigInt";
 require 5.005;
 
-$VERSION = '1.55';
+$VERSION = '1.59';
 use Exporter;
 @ISA =       qw( Exporter );
 @EXPORT_OK = qw( objectify _swap bgcd blcm); 
@@ -524,7 +524,6 @@ sub bnan
     # otherwise do our own thing
     $self->{value} = $CALC->_zero();
     }
-  $self->{value} = $CALC->_zero();
   $self->{sign} = $nan;
   delete $self->{_a}; delete $self->{_p};      # rounding NaN is silly
   return $self;
@@ -903,12 +902,8 @@ sub badd
   my ($self,$x,$y,@r) = objectify(2,@_);
 
   return $x if $x->modify('badd');
-#  print "mbi badd ",join(' ',caller()),"\n";
-#  print "upgrade => ",$upgrade||'undef',
-#    " \$x (",ref($x),") \$y (",ref($y),")\n";
   return $upgrade->badd($x,$y,@r) if defined $upgrade &&
-    ((ref($x) eq $upgrade) || (ref($y) eq $upgrade));
-#  print "still badd\n";
+    ((!$x->isa($self)) || (!$y->isa($self)));
 
   $r[3] = $y;                          # no push!
   # inf and NaN handling
@@ -969,8 +964,10 @@ sub bsub
   my ($self,$x,$y,@r) = objectify(2,@_);
 
   return $x if $x->modify('bsub');
+
+# upgrade done by badd():
 #  return $upgrade->badd($x,$y,@r) if defined $upgrade &&
-#    ((ref($x) eq $upgrade) || (ref($y) eq $upgrade));
+#   ((!$x->isa($self)) || (!$y->isa($self)));
 
   if ($y->is_zero())
     { 
@@ -1295,8 +1292,9 @@ sub bdiv
   return $self->_div_inf($x,$y)
    if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/) || $y->is_zero());
 
+  #print "mbi bdiv $x $y\n";
   return $upgrade->bdiv($upgrade->new($x),$y,@r)
-   if defined $upgrade && $y->isa($upgrade);
+   if defined $upgrade && !$y->isa($self);
 
   $r[3] = $y;                                  # no push!
 
@@ -1358,6 +1356,9 @@ sub bdiv
   $x->round(@r); 
   }
 
+###############################################################################
+# modulus functions
+
 sub bmod 
   {
   # modulus (or remainder)
@@ -1380,7 +1381,12 @@ sub bmod
       {
       my $xsign = $x->{sign};
       $x->{sign} = $y->{sign};
-      $x = $y-$x if $xsign ne $y->{sign};      # one of them '-'
+      if ($xsign ne $y->{sign})
+        {
+        my $t = [ @{$x->{value}} ];                    # copy $x
+        $x->{value} = [ @{$y->{value}} ];              # copy $y to $x
+        $x->{value} = $CALC->_sub($y->{value},$t,1);   # $y-$x
+        }
       }
     else
       {
@@ -1397,6 +1403,117 @@ sub bmod
   $x;
   }
 
+sub bmodinv
+  {
+  # modular inverse.  given a number which is (hopefully) relatively
+  # prime to the modulus, calculate its inverse using Euclid's
+  # alogrithm.  if the number is not relatively prime to the modulus
+  # (i.e. their gcd is not one) then NaN is returned.
+
+  my ($self,$num,$mod,@r) = objectify(2,@_);
+
+  return $num if $num->modify('bmodinv');
+
+  return $num->bnan()
+       if ($mod->{sign} ne '+'                         # -, NaN, +inf, -inf
+         || $num->is_zero()                            # or num == 0
+        || $num->{sign} !~ /^[+-]$/                    # or num NaN, inf, -inf
+        );
+
+  # put least residue into $num if $num was negative, and thus make it positive
+  $num->bmod($mod) if $num->{sign} eq '-';
+
+  if ($CALC->can('_modinv'))
+    {
+    $num->{value} = $CALC->_modinv($num->{value},$mod->{value});
+    $num->bnan() if !defined $num->{value} ;            # in case there was no
+    return $num;
+    }
+
+  my ($u, $u1) = ($self->bzero(), $self->bone());
+  my ($a, $b) = ($mod->copy(), $num->copy());
+
+  # first step need always be done since $num (and thus $b) is never 0
+  # Note that the loop is aligned so that the check occurs between #2 and #1
+  # thus saving us one step #2 at the loop end. Typical loop count is 1. Even
+  # a case with 28 loops still gains about 3% with this layout.
+  my $q;
+  ($a, $q, $b) = ($b, $a->bdiv($b));                    # step #1
+  # Euclid's Algorithm
+  while (!$b->is_zero())
+    {
+    ($u, $u1) = ($u1, $u->bsub($u1->copy()->bmul($q))); # step #2
+    ($a, $q, $b) = ($b, $a->bdiv($b));                  # step #1 again
+    }
+
+  # if the gcd is not 1, then return NaN!  It would be pointless to
+  # have called bgcd to check this first, because we would then be performing
+  # the same Euclidean Algorithm *twice*
+  return $num->bnan() unless $a->is_one();
+
+  $u1->bmod($mod);
+  $num->{value} = $u1->{value};
+  $num->{sign} = $u1->{sign};
+  $num;
+  }
+
+sub bmodpow
+  {
+  # takes a very large number to a very large exponent in a given very
+  # large modulus, quickly, thanks to binary exponentation.  supports
+  # negative exponents.
+  my ($self,$num,$exp,$mod,@r) = objectify(3,@_);
+
+  return $num if $num->modify('bmodpow');
+
+  # check modulus for valid values
+  return $num->bnan() if ($mod->{sign} ne '+'          # NaN, - , -inf, +inf
+                       || $mod->is_zero());
+
+  # check exponent for valid values
+  if ($exp->{sign} =~ /\w/) 
+    {
+    # i.e., if it's NaN, +inf, or -inf...
+    return $num->bnan();
+    }
+
+  $num->bmodinv ($mod) if ($exp->{sign} eq '-');
+
+  # check num for valid values (also NaN if there was no inverse but $exp < 0)
+  return $num->bnan() if $num->{sign} !~ /^[+-]$/;
+
+  if ($CALC->can('_modpow'))
+    {
+    # $mod is positive, sign on $exp is ignored, result also positive
+    $num->{value} = $CALC->_modpow($num->{value},$exp->{value},$mod->{value});
+    return $num;
+    }
+
+  # in the trivial case,
+  return $num->bzero() if $mod->is_one();
+  return $num->bone() if $num->is_zero() or $num->is_one();
+
+  # $num->bmod($mod);           # if $x is large, make it smaller first
+  my $acc = $num->copy();      # but this is not really faster...
+
+  $num->bone(); # keep ref to $num
+
+  my $expbin = $exp->as_bin(); $expbin =~ s/^[-]?0b//; # ignore sign and prefix
+  my $len = length($expbin);
+  while (--$len >= 0)
+    {
+    if( substr($expbin,$len,1) eq '1')
+      {
+      $num->bmul($acc)->bmod($mod);
+      }
+    $acc->bmul($acc)->bmod($mod);
+    }
+
+  $num;
+  }
+
+###############################################################################
+
 sub bfac
   {
   # (BINT or num_str, BINT or num_str) return BINT
@@ -1436,7 +1553,7 @@ sub bpow
   return $x if $x->modify('bpow');
 
   return $upgrade->bpow($upgrade->new($x),$y,@r)
-   if defined $upgrade && $y->isa($upgrade);
+   if defined $upgrade && !$y->isa($self);
 
   $r[3] = $y;                                  # no push!
   return $x if $x->{sign} =~ /^[+-]inf$/;      # -inf/+inf ** x
@@ -1478,15 +1595,14 @@ sub bpow
 #    }
 
   my $pow2 = $self->__one();
-  my $y1 = $class->new($y);
-  my $two = $self->new(2);
-  while (!$y1->is_one())
+  my $y_bin = $y->as_bin(); $y_bin =~ s/^0b//;
+  my $len = length($y_bin);
+  while (--$len > 0)
     {
-    $pow2->bmul($x) if $y1->is_odd();
-    $y1->bdiv($two);
+    $pow2->bmul($x) if substr($y_bin,$len,1) eq '1';   # is odd?
     $x->bmul($x);
     }
-  $x->bmul($pow2) unless $pow2->is_one();
+  $x->bmul($pow2);
   $x->round(@r);
   }
 
@@ -1539,7 +1655,7 @@ sub brsft
       $bin =~ s/^-0b//;                        # strip '-0b' prefix
       $bin =~ tr/10/01/;               # flip bits
       # now shift
-      if (length($bin) <= $y)
+      if (CORE::length($bin) <= $y)
         {
        $bin = '0';                     # shifting to far right creates -1
                                        # 0, because later increment makes 
@@ -2074,7 +2190,6 @@ sub objectify
   # currently it tries 'Math::BigInt' + 1, which will not work.
 
   # some shortcut for the common cases
-
   # $x->unary_op();
   return (ref($_[1]),$_[1]) if (@_ == 2) && ($_[0]||0 == 1) && ref($_[1]);
 
@@ -2092,6 +2207,7 @@ sub objectify
     $a[0] = $class;
     $a[0] = shift if $_[0] =~ /^[A-Z].*::/;    # classname as first?
     }
+
   no strict 'refs';
   # disable downgrading, because Math::BigFLoat->foo('1.0','2.0') needs floats
   if (defined ${"$a[0]::downgrade"})
@@ -2100,6 +2216,7 @@ sub objectify
     ${"$a[0]::downgrade"} = undef;
     }
 
+  my $up = ${"$a[0]::upgrade"};
   # print "Now in objectify, my class is today $a[0]\n";
   if ($count == 0)
     {
@@ -2110,7 +2227,7 @@ sub objectify
         {
         $k = $a[0]->new($k);
         }
-      elsif (ref($k) ne $a[0])
+      elsif (!defined $up && ref($k) ne $a[0])
        {
        # foreign object, try to convert to integer
         $k->can('as_number') ?  $k = $k->as_number() : $k = $a[0]->new($k);
@@ -2128,7 +2245,7 @@ sub objectify
         {
         $k = $a[0]->new($k);
         }
-      elsif (ref($k) ne $a[0])
+      elsif (!defined $up && ref($k) ne $a[0])
        {
        # foreign object, try to convert to integer
         $k->can('as_number') ?  $k = $k->as_number() : $k = $a[0]->new($k);
@@ -2147,29 +2264,30 @@ sub import
   my $self = shift;
 
   $IMPORT++;
-  my @a = @_; my $l = scalar @_; my $j = 0;
-  for ( my $i = 0; $i < $l ; $i++,$j++ )
+  my @a; my $l = scalar @_;
+  for ( my $i = 0; $i < $l ; $i++ )
     {
     if ($_[$i] eq ':constant')
       {
       # this causes overlord er load to step in
       overload::constant integer => sub { $self->new(shift) };
       overload::constant binary => sub { $self->new(shift) };
-      splice @a, $j, 1; $j --;
       }
     elsif ($_[$i] eq 'upgrade')
       {
       # this causes upgrading
       $upgrade = $_[$i+1];             # or undef to disable
-      my $s = 2; $s = 1 if @a-$j < 2;  # avoid "can not modify non-existant..."
-      splice @a, $j, $s; $j -= $s;
+      $i++;
       }
     elsif ($_[$i] =~ /^lib$/i)
       {
       # this causes a different low lib to take care...
       $CALC = $_[$i+1] || '';
-      my $s = 2; $s = 1 if @a-$j < 2;  # avoid "can not modify non-existant..."
-      splice @a, $j, $s; $j -= $s;
+      $i++;
+      }
+    else
+      {
+      push @a, $_[$i];
       }
     }
   # any non :constant stuff is handled by our parent, Exporter
@@ -2183,15 +2301,18 @@ sub import
   $CALC = '';                                  # signal error
   foreach my $lib (@c)
     {
+    next if ($lib || '') eq '';
     $lib = 'Math::BigInt::'.$lib if $lib !~ /^Math::BigInt/i;
     $lib =~ s/\.pm$//;
     if ($] < 5.006)
       {
       # Perl < 5.6.0 dies with "out of memory!" when eval() and ':constant' is
       # used in the same script, or eval inside import().
-      (my $mod = $lib . '.pm') =~ s!::!/!g;
-      # require does not automatically :: => /, so portability problems arise
-      eval { require $mod; $lib->import( @c ); }
+      my @parts = split /::/, $lib;             # Math::BigInt => Math BigInt
+      my $file = pop @parts; $file .= '.pm';    # BigInt => BigInt.pm
+      require File::Spec;
+      $file = File::Spec->catfile (@parts, $file);
+      eval { require "$file"; $lib->import( @c ); }
       }
     else
       {
@@ -2331,7 +2452,8 @@ sub _split
     $es = $1; $ev = $2;
     # valid mantissa?
     return if $m eq '.' || $m eq '';
-    my ($mi,$mf) = split /\./,$m;
+    my ($mi,$mf,$last) = split /\./,$m;
+    return if defined $last;           # last defined => 1.2.3 or others
     $mi = '0' if !defined $mi;
     $mi .= '0' if $mi =~ /^[\-\+]?$/;
     $mf = '0' if !defined $mf || $mf eq '';
@@ -2372,12 +2494,19 @@ sub as_hex
     }
   else
     {
-    my $x1 = $x->copy()->babs(); my $xr;
-    my $x10000 = Math::BigInt->new (0x10000);
+    my $x1 = $x->copy()->babs(); my ($xr,$x10000,$h);
+    if ($] >= 5.006)
+      {
+      $x10000 = Math::BigInt->new (0x10000); $h = 'h4';
+      }
+    else
+      {
+      $x10000 = Math::BigInt->new (0x1000); $h = 'h3';
+      }
     while (!$x1->is_zero())
       {
       ($x1, $xr) = bdiv($x1,$x10000);
-      $es .= unpack('h4',pack('v',$xr->numify()));
+      $es .= unpack($h,pack('v',$xr->numify()));
       }
     $es = reverse $es;
     $es =~ s/^[0]+//;  # strip leading zeros
@@ -2402,12 +2531,19 @@ sub as_bin
     }
   else
     {
-    my $x1 = $x->copy()->babs(); my $xr;
-    my $x10000 = Math::BigInt->new (0x10000);
+    my $x1 = $x->copy()->babs(); my ($xr,$x10000,$b);
+    if ($] >= 5.006)
+      {
+      $x10000 = Math::BigInt->new (0x10000); $b = 'b16';
+      }
+    else
+      {
+      $x10000 = Math::BigInt->new (0x1000); $b = 'b12';
+      }
     while (!$x1->is_zero())
       {
       ($x1, $xr) = bdiv($x1,$x10000);
-      $es .= unpack('b16',pack('v',$xr->numify()));
+      $es .= unpack($b,pack('v',$xr->numify()));
       }
     $es = reverse $es; 
     $es =~ s/^[0]+//;  # strip leading zeros
@@ -2515,6 +2651,9 @@ Math::BigInt - Arbitrary size integer math package
                                # return (quo,rem) or quo if scalar
 
   $x->bmod($y);                        # modulus (x % y)
+  $x->bmodpow($exp,$mod);      # modular exponentation (($num**$exp) % $mod))
+  $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
   $x->brsft($y);               # right shift 
@@ -2557,6 +2696,8 @@ Math::BigInt - Arbitrary size integer math package
   $x->bsstr();                 # normalized string in scientific notation
   $x->as_hex();                        # as signed hexadecimal string with prefixed 0x
   $x->as_bin();                        # as signed binary string with prefixed 0b
+  
+  Math::BigInt->config();      # return hash containing configuration/version
 
 =head1 DESCRIPTION
 
@@ -2612,6 +2753,15 @@ Each of the methods below accepts three additional parameters. These arguments
 $A, $P and $R are accuracy, precision and round_mode. Please see more in the
 section about ACCURACY and ROUNDIND.
 
+=head2 config
+
+       use Data::Dumper;
+
+       print Dumper ( Math::BigInt->config() );
+
+Returns a hash containing the configuration, e.g. the version number, lib
+loaded etc.
+
 =head2 accuracy
 
        $x->accuracy(5);                # local for $x
@@ -2824,6 +2974,35 @@ numbers.
 
   $x->bmod($y);                        # modulus (x % y)
 
+=head2 bmodinv
+
+  $num->bmodinv($mod);         # modular inverse
+
+Returns the inverse of C<$num> in the given modulus C<$mod>.  'C<NaN>' is
+returned unless C<$num> is relatively prime to C<$mod>, i.e. unless
+C<bgcd($num, $mod)==1>.
+
+=head2 bmodpow
+
+  $num->bmodpow($exp,$mod);    # modular exponentation ($num**$exp % $mod)
+
+Returns the value of C<$num> taken to the power C<$exp> in the modulus
+C<$mod> using binary exponentation.  C<bmodpow> is far superior to
+writing
+
+  $num ** $exp % $mod
+
+because C<bmodpow> is much faster--it reduces internal variables into
+the modulus whenever possible, so it operates on smaller numbers.
+
+C<bmodpow> also supports negative exponents.
+
+  bmodpow($num, -1, $mod)
+
+is exactly equivalent to
+
+  bmodinv($num, $mod)
+
 =head2 bpow
 
   $x->bpow($y);                        # power of arguments (x ** y)