my $class = "Math::BigInt";
require 5.005;
-$VERSION = '1.57';
+$VERSION = '1.59';
use Exporter;
@ISA = qw( Exporter );
@EXPORT_OK = qw( objectify _swap bgcd blcm);
{
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
{
$x;
}
-sub bmodinv_not_yet_implemented
+sub bmodinv
{
# modular inverse. given a number which is (hopefully) relatively
# prime to the modulus, calculate its inverse using Euclid's
|| $num->is_zero() # or num == 0
|| $num->{sign} !~ /^[+-]$/ # or num NaN, inf, -inf
);
-# return $num # i.e., NaN or some kind of infinity,
-# if ($num->{sign} =~ /\w/);
- # the remaining case, nonpositive case, $num < 0, is addressed below.
+ # 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());
- # put least residue into $b if $num was negative
- $b %= $mod if $b->{sign} eq '-';
-
- # Euclid's Algorithm
- while( ! $b->is_zero()) {
- ($a, my $q, $b) = ($b, $self->bdiv( $a->copy(), $b));
- ($u, $u1) = ($u1, $u - $u1 * $q);
+ # 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 first, because we would then be performing the
- # same Euclidean Algorithm *twice*
- return $self->bnan() unless $a->is_one();
+ # 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();
- $u %= $mod;
- return $u;
+ $u1->bmod($mod);
+ $num->{value} = $u1->{value};
+ $num->{sign} = $u1->{sign};
+ $num;
}
-sub bmodpow_not_yet_implemented
+sub bmodpow
{
# takes a very large number to a very large exponent in a given very
# large modulus, quickly, thanks to binary exponentation. supports
# i.e., if it's NaN, +inf, or -inf...
return $num->bnan();
}
- elsif ($exp->{sign} eq '-')
+
+ $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'))
{
- $exp->babs();
- $num->bmodinv ($mod);
- return $num if $num->{sign} !~ /^[+-]/; # i.e. if there was no inverse
+ # $mod is positive, sign on $exp is ignored, result also positive
+ $num->{value} = $CALC->_modpow($num->{value},$exp->{value},$mod->{value});
+ return $num;
}
- # check num for valid values
- return $num->bnan() if $num->{sign} !~ /^[+-]$/;
+ # in the trivial case,
+ return $num->bzero() if $mod->is_one();
+ return $num->bone() if $num->is_zero() or $num->is_one();
- # 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...
- my $acc = $num->copy(); $num->bone(); # keep ref to $num
+ $num->bone(); # keep ref to $num
- print "$num $acc $exp\n";
- while( !$exp->is_zero() ) {
- if( $exp->is_odd() ) {
- $num->bmul($acc)->bmod($mod);
+ 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);
- $exp->brsft( 1, 2); # remove last (binary) digit
- print "$num $acc $exp\n";
+ $acc->bmul($acc)->bmod($mod);
}
- return $num;
+
+ $num;
}
###############################################################################
# }
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);
}
$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
{
$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 '';
}
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
}
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
# 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
=head2 bmodinv
-Not yet implemented.
-
- bmodinv($num,$mod); # modular inverse (no OO style)
+ $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
=head2 bmodpow
-Not yet implemented.
-
- bmodpow($num,$exp,$mod); # modular exponentation ($num**$exp % $mod)
+ $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
projects / p5sagit/p5-mst-13.2.git / blobdiff