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:
# _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;
# 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} !~ /^[+-]$/))
{
$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++;
{
$x->bsqrt(); $x->bmul($xc); $d->bdec(); # 1
}
- print "at $x\n";
}
# assume fraction ends in 1
$x->bsqrt(); # 1
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)
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);
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);
$x->round(@r);
}
+###############################################################################
+# modulus functions
+
sub bmod
{
# modulus (or remainder)
$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
${"$a[0]::downgrade"} = undef;
}
+ my $up = ${"$a[0]::upgrade"};
# print "Now in objectify, my class is today $a[0]\n";
if ($count == 0)
{
{
$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);
{
$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);
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
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
# 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
$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<NaN>' is
+returned unless C<$num> is relatively prime to C<$mod>, i.e. unless
+C<bgcd($num, $mod)==1>.
+
+=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<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)
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
}
# 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)
{
# 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)
{
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);
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);
$x;
}
+sub _modinv
+ {
+ # inverse modulus
+ }
+
+sub _modpow
+ {
+ # modulus of power ($x ** $y) % $z
+ }
+
##############################################################################
##############################################################################
_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').
$| = 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
while (<DATA>)
{
- chop;
+ chomp;
$_ =~ s/#.*$//; # remove comments
$_ =~ s/\s+$//; # trailing spaces
next if /^$/; # skip empty lines & comments
while (<DATA>)
{
- chop;
+ chomp;
next if /^#/; # skip comments
if (s/^&//)
{
$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;
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
my ($func,@args,$ans,$rc,$class,$try);
while (<DATA>)
{
- chop;
+ chomp;
next if /^#/; # skip comments
if (s/^&//)
{
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())
{
# $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);
}
my $CALC = Math::BigInt->config()->{lib};
while (<DATA>)
{
- chop;
+ chomp;
next if /^\s*(#|$)/; # skip comments and empty lines
if (s/^&//)
{
while (<DATA>)
{
- chop;
+ chomp;
next if /^#/; # skip comments
if (s/^&//)
{
projects / p5sagit/p5-mst-13.2.git / commitdiff