From: Jarkko Hietaniemi Date: Thu, 30 May 2002 20:39:30 +0000 (+0000) Subject: Upgrade to Math::BigInt 1.57. X-Git-Url: http://git.shadowcat.co.uk/gitweb/gitweb.cgi?a=commitdiff_plain;h=d614cd8b2519c84f1ee8ae0c9c71fba2ed16cfb3;p=p5sagit%2Fp5-mst-13.2.git Upgrade to Math::BigInt 1.57. p4raw-id: //depot/perl@16906 --- diff --git a/lib/Math/BigFloat.pm b/lib/Math/BigFloat.pm index 33cf3d1..ea78da5 100644 --- a/lib/Math/BigFloat.pm +++ b/lib/Math/BigFloat.pm @@ -1,7 +1,7 @@ package Math::BigFloat; # -# Mike grinned. 'Two down, infinity to go' - Mike Nostrus in Before and After +# Mike grinned. 'Two down, infinity to go' - Mike Nostrus in 'Before and After' # # The following hash values are internally used: @@ -12,7 +12,7 @@ package Math::BigFloat; # _p: precision # _f: flags, used to signal MBI not to touch our private parts -$VERSION = '1.31'; +$VERSION = '1.32'; require 5.005; use Exporter; use File::Spec; @@ -440,7 +440,6 @@ sub badd # return result as BFLOAT my ($self,$x,$y,$a,$p,$r) = objectify(2,@_); - #print "mbf badd $x $y\n"; # inf and NaN handling if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/)) { @@ -678,7 +677,7 @@ sub blog $next = $over->copy()->bdiv($below->copy()->bmul($factor),$scale); last if $next->bcmp($limit) <= 0; $x->badd($next); - # print "step $steps $x\n"; + # print "step $x\n"; # calculate things for the next term $over *= $u; $below *= $v; $factor->badd($f); #$steps++; @@ -1199,7 +1198,6 @@ sub _pow2 { $x->bsqrt(); $x->bmul($xc); $d->bdec(); # 1 } - print "at $x\n"; } # assume fraction ends in 1 $x->bsqrt(); # 1 @@ -1337,7 +1335,7 @@ sub bpow return $x->bone() if $y->is_zero(); return $x if $x->is_one() || $y->is_one(); - return $x->_pow2($y,$a,$p,$r) if !$y->is_int(); # non-integer power + return $x->_pow($y,$a,$p,$r) if !$y->is_int(); # non-integer power my $y1 = $y->as_number(); # make bigint # if ($x == -1) diff --git a/lib/Math/BigInt.pm b/lib/Math/BigInt.pm index dd6521e..a85a474 100644 --- a/lib/Math/BigInt.pm +++ b/lib/Math/BigInt.pm @@ -18,7 +18,7 @@ package Math::BigInt; my $class = "Math::BigInt"; require 5.005; -$VERSION = '1.56'; +$VERSION = '1.57'; use Exporter; @ISA = qw( Exporter ); @EXPORT_OK = qw( objectify _swap bgcd blcm); @@ -1292,6 +1292,7 @@ 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($self); @@ -1355,6 +1356,9 @@ sub bdiv $x->round(@r); } +############################################################################### +# modulus functions + sub bmod { # modulus (or remainder) @@ -1394,6 +1398,97 @@ sub bmod $x; } +sub bmodinv_not_yet_implemented + { + # 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 + ); +# return $num # i.e., NaN or some kind of infinity, +# if ($num->{sign} =~ /\w/); + + # the remaining case, nonpositive case, $num < 0, is addressed below. + + 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); + } + + # 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(); + + $u %= $mod; + return $u; + } + +sub bmodpow_not_yet_implemented + { + # 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(); + } + elsif ($exp->{sign} eq '-') + { + $exp->babs(); + $num->bmodinv ($mod); + return $num if $num->{sign} !~ /^[+-]/; # i.e. if there was no inverse + } + + # 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(); + + my $acc = $num->copy(); $num->bone(); # keep ref to $num + + print "$num $acc $exp\n"; + while( !$exp->is_zero() ) { + if( $exp->is_odd() ) { + $num->bmul($acc)->bmod($mod); + } + $acc->bmul($acc)->bmod($mod); + $exp->brsft( 1, 2); # remove last (binary) digit + print "$num $acc $exp\n"; + } + return $num; + } + +############################################################################### + sub bfac { # (BINT or num_str, BINT or num_str) return BINT @@ -2097,6 +2192,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) { @@ -2107,7 +2203,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); @@ -2125,7 +2221,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,7 +2243,6 @@ sub import my @a; my $l = scalar @_; for ( my $i = 0; $i < $l ; $i++ ) { -# print "at $_[$i]\n"; if ($_[$i] eq ':constant') { # this causes overlord er load to step in @@ -2171,7 +2266,6 @@ sub import push @a, $_[$i]; } } -# print "a @a\n"; # any non :constant stuff is handled by our parent, Exporter # even if @_ is empty, to give it a chance $self->SUPER::import(@a); # need it for subclasses @@ -2515,6 +2609,7 @@ Math::BigInt - Arbitrary size integer math package # return (quo,rem) or quo if scalar $x->bmod($y); # modulus (x % y) + $x->bpow($y); # power of arguments (x ** y) $x->blsft($y); # left shift $x->brsft($y); # right shift @@ -2835,6 +2930,39 @@ numbers. $x->bmod($y); # modulus (x % y) +=head2 bmodinv + +Not yet implemented. + + bmodinv($num,$mod); # modular inverse (no OO style) + +Returns the inverse of C<$num> in the given modulus C<$mod>. 'C' is +returned unless C<$num> is relatively prime to C<$mod>, i.e. unless +C. + +=head2 bmodpow + +Not yet implemented. + + bmodpow($num,$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 is far superior to +writing + + $num ** $exp % $mod + +because C is much faster--it reduces internal variables into +the modulus whenever possible, so it operates on smaller numbers. + +C 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) diff --git a/lib/Math/BigInt/Calc.pm b/lib/Math/BigInt/Calc.pm index f1ade92..a2fe812 100644 --- a/lib/Math/BigInt/Calc.pm +++ b/lib/Math/BigInt/Calc.pm @@ -8,7 +8,7 @@ require Exporter; use vars qw/@ISA $VERSION/; @ISA = qw(Exporter); -$VERSION = '0.27'; +$VERSION = '0.28'; # Package to store unsigned big integers in decimal and do math with them @@ -485,9 +485,11 @@ sub _mul_use_mul } # since multiplying $x with $x fails, make copy in this case - $yv = [@$xv] if "$xv" eq "$yv"; # same references? + $yv = [@$xv] if $xv == $yv; # same references? +# $yv = [@$xv] if "$xv" eq "$yv"; # same references? + # since multiplying $x with $x would fail here, use the faster squaring -# return _square($c,$xv) if "$xv" eq "$yv"; # same reference? +# return _square($c,$xv) if $xv == $yv; # same reference? if ($LEN_CONVERT != 0) { @@ -565,9 +567,10 @@ sub _mul_use_div # since multiplying $x with $x fails, make copy in this case - $yv = [@$xv] if "$xv" eq "$yv"; # same references? + $yv = [@$xv] if $xv == $yv; # same references? +# $yv = [@$xv] if "$xv" eq "$yv"; # same references? # since multiplying $x with $x would fail here, use the faster squaring -# return _square($c,$xv) if "$xv" eq "$yv"; # same reference? +# return _square($c,$xv) if $xv == $yv; # same reference? if ($LEN_CONVERT != 0) { @@ -642,7 +645,7 @@ sub _div_use_mul return $x; } - my $y = [ @$yorg ]; + my $y = [ @$yorg ]; # always make copy to preserve if ($LEN_CONVERT != 0) { $c->_to_small($x); $c->_to_small($y); @@ -782,7 +785,7 @@ sub _div_use_div return $x; } - my $y = [ @$yorg ]; + my $y = [ @$yorg ]; # always make copy to preserve if ($LEN_CONVERT != 0) { $c->_to_small($x); $c->_to_small($y); @@ -1571,6 +1574,16 @@ sub _from_bin $x; } +sub _modinv + { + # inverse modulus + } + +sub _modpow + { + # modulus of power ($x ** $y) % $z + } + ############################################################################## ############################################################################## @@ -1674,6 +1687,8 @@ slow) fallback routines to emulate these: _gcd(obj,obj) return Greatest Common Divisor of two objects _zeros(obj) return number of trailing decimal zeros + _modinv return inverse modulus + _modpow return modulus of power ($x ** $y) % $z Input strings come in as unsigned but with prefix (i.e. as '123', '0xabc' or '0b1101'). diff --git a/lib/Math/BigInt/t/bare_mbf.t b/lib/Math/BigInt/t/bare_mbf.t index a160c7c..8b682f9 100644 --- a/lib/Math/BigInt/t/bare_mbf.t +++ b/lib/Math/BigInt/t/bare_mbf.t @@ -8,6 +8,7 @@ BEGIN $| = 1; # to locate the testing files my $location = $0; $location =~ s/bare_mbf.t//i; + print "$0\n"; if ($ENV{PERL_CORE}) { # testing with the core distribution diff --git a/lib/Math/BigInt/t/bigfltpm.inc b/lib/Math/BigInt/t/bigfltpm.inc index 734b935..d455165 100644 --- a/lib/Math/BigInt/t/bigfltpm.inc +++ b/lib/Math/BigInt/t/bigfltpm.inc @@ -4,7 +4,7 @@ ok ($class->config()->{lib},$CL); while () { - chop; + chomp; $_ =~ s/#.*$//; # remove comments $_ =~ s/\s+$//; # trailing spaces next if /^$/; # skip empty lines & comments diff --git a/lib/Math/BigInt/t/bigintpm.inc b/lib/Math/BigInt/t/bigintpm.inc index 2bcf346..4d74a41 100644 --- a/lib/Math/BigInt/t/bigintpm.inc +++ b/lib/Math/BigInt/t/bigintpm.inc @@ -40,7 +40,7 @@ my ($f,$z,$a,$exp,@a,$m,$e,$round_mode,$expected_class); while () { - chop; + chomp; next if /^#/; # skip comments if (s/^&//) { @@ -163,9 +163,18 @@ while () $try .= "\$x ^ \$y;"; }elsif ($f eq "bpow"){ $try .= "\$x ** \$y;"; + } elsif( $f eq "bmodinv") { + $try .= "\$x->bmodinv(\$y);"; }elsif ($f eq "digit"){ $try .= "\$x->digit(\$y);"; - } else { warn "Unknown op '$f'"; } + } else { + $try .= "\$z = $class->new(\"$args[2]\");"; + + # Functions with three arguments + if( $f eq "bmodpow") { + $try .= "\$x->bmodpow(\$y,\$z);"; + } else { warn "Unknown op '$f'"; } + } } # end else all other ops $ans1 = eval $try; @@ -1311,6 +1320,42 @@ inf:0:inf 14:3:4 # bug in Calc with '99999' vs $BASE-1 10000000000000000000000000000000000000000000000000000000000000000000000000000000000:10000000375084540248994272022843165711074:999999962491547381984643365663244474111576 +#&bmodinv +## format: number:modulus:result +## bmodinv Data errors +#abc:abc:NaN +#abc:5:NaN +#5:abc:NaN +## bmodinv Expected Results from normal use +#1:5:1 +#3:5:2 +#-2:5:2 +#324958749843759385732954874325984357439658735983745:2348249874968739:1741662881064902 +## bmodinv Error cases / useless use of function +#3:-5:NaN +#inf:5:NaN +#&bmodpow +## format: number:exponent:modulus:result +## bmodpow Data errors +#abc:abc:abc:NaN +#5:abc:abc:NaN +#abc:5:abc:NaN +#abc:abc:5:NaN +#5:5:abc:NaN +#5:abc:5:NaN +#abc:5:5:NaN +## bmodpow Expected results +#0:0:2:1 +#1:0:2:1 +#0:0:1:0 +#8:7:5032:3840 +#8:-1:5033:4404 +#98436739867439843769485798542749827593285729587325:43698764986460981048259837659386739857456983759328457:6943857329857295827698367:3104744730915914415259518 +## bmodpow Error cases +#8:8:-5:NaN +#8:-1:16:NaN +#inf:5:13:NaN +#5:inf:13:NaN &bmod # inf handling, see table in doc 0:inf:0 diff --git a/lib/Math/BigInt/t/calling.t b/lib/Math/BigInt/t/calling.t index ca78fe9..8a0d2bb 100644 --- a/lib/Math/BigInt/t/calling.t +++ b/lib/Math/BigInt/t/calling.t @@ -63,7 +63,7 @@ my $version = '1.46'; # adjust manually to match latest release my ($func,@args,$ans,$rc,$class,$try); while () { - chop; + chomp; next if /^#/; # skip comments if (s/^&//) { diff --git a/lib/Math/BigInt/t/mbi_rand.t b/lib/Math/BigInt/t/mbi_rand.t index bbe73c8..69be2d4 100644 --- a/lib/Math/BigInt/t/mbi_rand.t +++ b/lib/Math/BigInt/t/mbi_rand.t @@ -13,33 +13,35 @@ BEGIN my $location = $0; $location =~ s/mbi_rand.t//; unshift @INC, $location; # to locate the testing files chdir 't' if -d 't'; - $count = 500; + $count = 128; plan tests => $count*2; } use Math::BigInt; my $c = 'Math::BigInt'; -my $length = 200; +my $length = 128; # If you get a failure here, please re-run the test with the printed seed # value as input: perl t/mbi_rand.t seed my $seed = int(rand(65537)); print "# seed: $seed\n"; srand($seed); -my ($A,$B,$ADB,$AMB,$la,$lb); +my ($A,$B,$As,$Bs,$ADB,$AMB,$la,$lb); +my $two = Math::BigInt->new(2); for (my $i = 0; $i < $count; $i++) { # length of A and B $la = int(rand($length)+1); $lb = int(rand($length)+1); - $A = ''; $B = ''; + $As = ''; $Bs = ''; # we create the numbers from "patterns", e.g. get a random number and a # random count and string them together. This means things like # "100000999999999999911122222222" are much more likely. If we just strung # together digits, we would end up with "1272398823211223" etc. - while (length($A) < $la) { $A .= int(rand(100)) x int(rand(16)); } - while (length($B) < $lb) { $B .= int(rand(100)) x int(rand(16)); } - $A = $c->new($A); $B = $c->new($B); + while (length($As) < $la) { $As .= int(rand(100)) x int(rand(16)); } + while (length($Bs) < $lb) { $Bs .= int(rand(100)) x int(rand(16)); } + $As =~ s/^0+//; $Bs =~ s/^0+//; + $A = $c->new($As); $B = $c->new($Bs); # print "# A $A\n# B $B\n"; if ($A->is_zero() || $B->is_zero()) { @@ -49,11 +51,11 @@ for (my $i = 0; $i < $count; $i++) # $X = ($A/$B)*$B + 2 * ($A % $B) - ($A % $B); ($ADB,$AMB) = $A->copy()->bdiv($B); print "# ". join(' ',Math::BigInt::Calc->_base_len()),"\n" - unless ok ($ADB*$B+2*$AMB-$AMB,$A); + unless ok ($ADB*$B+$two*$AMB-$AMB,$As); # swap 'em and try this, too # $X = ($B/$A)*$A + $B % $A; ($ADB,$AMB) = $B->copy()->bdiv($A); print "# ". join(' ',Math::BigInt::Calc->_base_len()),"\n" - unless ok ($ADB*$A+2*$AMB-$AMB,$B); + unless ok ($ADB*$A+$two*$AMB-$AMB,$Bs); } diff --git a/lib/Math/BigInt/t/mbimbf.inc b/lib/Math/BigInt/t/mbimbf.inc index 968d830..1460161 100644 --- a/lib/Math/BigInt/t/mbimbf.inc +++ b/lib/Math/BigInt/t/mbimbf.inc @@ -596,7 +596,7 @@ my ($ans1,$f,$a,$p,$xp,$yp,$xa,$ya,$try,$ans,@args); my $CALC = Math::BigInt->config()->{lib}; while () { - chop; + chomp; next if /^\s*(#|$)/; # skip comments and empty lines if (s/^&//) { diff --git a/lib/Math/BigInt/t/upgrade.inc b/lib/Math/BigInt/t/upgrade.inc index bf35261..fc70873 100644 --- a/lib/Math/BigInt/t/upgrade.inc +++ b/lib/Math/BigInt/t/upgrade.inc @@ -52,7 +52,7 @@ my ($f,$z,$a,$exp,@a,$m,$e,$round_mode,$expected_class); while () { - chop; + chomp; next if /^#/; # skip comments if (s/^&//) {