-package Math::BigInt;
+package Math::BigInt::CalcEmu;
use 5.005;
use strict;
# use warnings; # dont use warnings for older Perls
-
use vars qw/$VERSION/;
-$VERSION = '0.02';
+$VERSION = '0.04';
+
+package Math::BigInt;
# See SYNOPSIS below.
$CALC_EMU = Math::BigInt->config()->{'lib'};
}
-sub __emu_blog
- {
- my ($self,$x,$base,@r) = @_;
-
- return $x->bnan() if $x->is_zero() || $base->is_zero() || $base->is_one();
-
- my $acmp = $x->bacmp($base);
- return $x->bone('+',@r) if $acmp == 0;
- return $x->bzero(@r) if $acmp < 0 || $x->is_one();
-
- # blog($x,$base) ** $base + $y = $x
-
- # this trial multiplication is very fast, even for large counts (like for
- # 2 ** 1024, since this still requires only 1024 very fast steps
- # (multiplication of a large number by a very small number is very fast))
- # See Calc for an even faster algorightmn
- my $x_org = $x->copy(); # preserve orgx
- $x->bzero(); # keep ref to $x
- my $trial = $base->copy();
- while ($trial->bacmp($x_org) <= 0)
- {
- $trial->bmul($base); $x->binc();
- }
- $x->round(@r);
- }
-
-sub __emu_bmodinv
- {
- my ($self,$x,$y,@r) = @_;
-
- my ($u, $u1) = ($self->bzero(), $self->bone());
- my ($a, $b) = ($y->copy(), $x->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 (calculate GCD of ($a,$b) in $a and also calculate
- # two values in $u and $u1, we use only $u1 afterwards)
- my $sign = 1; # flip-flop
- while (!$b->is_zero()) # found GCD if $b == 0
- {
- # the original algorithm had:
- # ($u, $u1) = ($u1, $u->bsub($u1->copy()->bmul($q))); # step #2
- # The following creates exact the same sequence of numbers in $u1,
- # except for the sign ($u1 is now always positive). Since formerly
- # the sign of $u1 was alternating between '-' and '+', the $sign
- # flip-flop will take care of that, so that at the end of the loop
- # we have the real sign of $u1. Keeping numbers positive gains us
- # speed since badd() is faster than bsub() and makes it possible
- # to have the algorithmn in Calc for even more speed.
-
- ($u, $u1) = ($u1, $u->badd($u1->copy()->bmul($q))); # step #2
- $sign = - $sign; # flip sign
-
- ($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* in case the gcd is 1.
- return $x->bnan() unless $a->is_one();
-
- $u1->bneg() if $sign != 1; # need to flip?
-
- $u1->bmod($y); # calc result
- $x->{value} = $u1->{value}; # and copy over to $x
- $x->{sign} = $u1->{sign}; # to modify in place
- $x->round(@r);
- }
-
-sub __emu_bmodpow
- {
- my ($self,$num,$exp,$mod,@r) = @_;
-
- # in the trivial case,
- return $num->bzero(@r) if $mod->is_one();
- return $num->bone('+',@r) 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 = CORE::length($expbin);
- while (--$len >= 0)
- {
- $num->bmul($acc)->bmod($mod) if substr($expbin,$len,1) eq '1';
- $acc->bmul($acc)->bmod($mod);
- }
-
- $num->round(@r);
- }
-
-sub __emu_bfac
- {
- my ($self,$x,@r) = @_;
-
- return $x->bone('+',@r) if $x->is_zero() || $x->is_one(); # 0 or 1 => 1
-
- my $n = $x->copy();
- $x->bone();
- # seems we need not to temp. clear A/P of $x since the result is the same
- my $f = $self->new(2);
- while ($f->bacmp($n) < 0)
- {
- $x->bmul($f); $f->binc();
- }
- $x->bmul($f,@r); # last step and also round result
- }
-
-sub __emu_bpow
- {
- my ($self,$x,$y,@r) = @_;
-
- return $x->bone('+',@r) if $y->is_zero();
- return $x->round(@r) if $x->is_one() || $y->is_one();
- return $x->round(@r) if $x->is_zero(); # 0**y => 0 (if not y <= 0)
-
- my $pow2 = $self->bone();
- my $y_bin = $y->as_bin(); $y_bin =~ s/^0b//;
- my $len = CORE::length($y_bin);
- while (--$len > 0)
- {
- $pow2->bmul($x) if substr($y_bin,$len,1) eq '1'; # is odd?
- $x->bmul($x);
- }
- $x->bmul($pow2);
- $x->round(@r) if !exists $x->{_f} || $x->{_f} & MB_NEVER_ROUND == 0;
- $x;
- }
-
sub __emu_band
{
my ($self,$x,$y,$sx,$sy,@r) = @_;
$bx = reverse $bx;
$by = reverse $by;
- # cut the longer string to the length of the shorter one (the result would
- # be 0 due to AND anyway)
+ # padd the shorter string
+ my $xx = "\x00"; $xx = "\x0f" if $sx == -1;
+ my $yy = "\x00"; $yy = "\x0f" if $sy == -1;
my $diff = CORE::length($bx) - CORE::length($by);
if ($diff > 0)
{
- $bx = substr($bx,0,CORE::length($by));
+ # if $yy eq "\x00", we can cut $bx, otherwise we need to padd $by
+ $by .= $yy x $diff;
}
elsif ($diff < 0)
{
- $by = substr($by,0,CORE::length($bx));
+ # if $xx eq "\x00", we can cut $by, otherwise we need to padd $bx
+ $bx .= $xx x abs($diff);
}
-
+
# and the strings together
my $r = $bx & $by;
# and reverse the result again
$bx = reverse $r;
- # one of $x or $y was negative, so need to flip bits in the result
- # in both cases (one or two of them negative, or both positive) we need
+ # One of $x or $y was negative, so need to flip bits in the result.
+ # In both cases (one or two of them negative, or both positive) we need
# to get the characters back.
if ($sign == 1)
{
$bx =~ tr/\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00/fedcba9876543210/;
}
+ # leading zeros will be stripped by _from_hex()
$bx = '0x' . $bx;
- if ($CALC_EMU->can('_from_hex'))
- {
- $x->{value} = $CALC_EMU->_from_hex( \$bx );
- }
- else
- {
- $r = $self->new($bx);
- $x->{value} = $r->{value};
- }
+ $x->{value} = $CALC_EMU->_from_hex( $bx );
# calculate sign of result
$x->{sign} = '+';
- #$x->{sign} = '-' if $sx == $sy && $sx == -1 && !$x->is_zero();
$x->{sign} = '-' if $sign == 1 && !$x->is_zero();
$x->bdec() if $sign == 1;
$bx =~ tr/\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00/fedcba9876543210/;
}
+ # leading zeros will be stripped by _from_hex()
$bx = '0x' . $bx;
- if ($CALC_EMU->can('_from_hex'))
- {
- $x->{value} = $CALC_EMU->_from_hex( \$bx );
- }
- else
- {
- $r = $self->new($bx);
- $x->{value} = $r->{value};
- }
+ $x->{value} = $CALC_EMU->_from_hex( $bx );
+
+ # calculate sign of result
+ $x->{sign} = '+';
+ $x->{sign} = '-' if $sign == 1 && !$x->is_zero();
# if one of X or Y was negative, we need to decrement result
$x->bdec() if $sign == 1;
$bx =~ tr/\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00/fedcba9876543210/;
}
+ # leading zeros will be stripped by _from_hex()
$bx = '0x' . $bx;
- if ($CALC_EMU->can('_from_hex'))
- {
- $x->{value} = $CALC_EMU->_from_hex( \$bx );
- }
- else
- {
- $r = $self->new($bx);
- $x->{value} = $r->{value};
- }
+ $x->{value} = $CALC_EMU->_from_hex( $bx );
# calculate sign of result
$x->{sign} = '+';
$x->round(@r);
}
-sub __emu_bsqrt
- {
- my ($self,$x,@r) = @_;
-
- # this is slow:
- return $x->round(@r) if $x->is_zero(); # 0,1 => 0,1
-
- return $x->bone('+',@r) if $x < 4; # 1,2,3 => 1
- my $y = $x->copy();
- my $l = int($x->length()/2);
-
- $x->bone(); # keep ref($x), but modify it
- $x->blsft($l,10) if $l != 0; # first guess: 1.('0' x (l/2))
-
- my $last = $self->bzero();
- my $two = $self->new(2);
- my $lastlast = $self->bzero();
- #my $lastlast = $x+$two;
- while ($last != $x && $lastlast != $x)
- {
- $lastlast = $last; $last = $x->copy();
- $x->badd($y / $x);
- $x->bdiv($two);
- }
- $x->bdec() if $x * $x > $y; # overshot?
- $x->round(@r);
- }
-
-sub __emu_broot
- {
- my ($self,$x,$y,@r) = @_;
-
- return $x->bsqrt() if $y->bacmp(2) == 0; # 2 => square root
-
- # since we take at least a cubic root, and only 8 ** 1/3 >= 2 (==2):
- return $x->bone('+',@r) if $x < 8; # $x=2..7 => 1
-
- my $num = $x->numify();
-
- if ($num <= 1000000)
- {
- $x = $self->new( int ( sprintf ("%.8f", $num ** (1 / $y->numify() ))));
- return $x->round(@r);
- }
-
- # if $n is a power of two, we can repeatedly take sqrt($X) and find the
- # proper result, because sqrt(sqrt($x)) == root($x,4)
- # See Calc.pm for more details
- my $b = $y->as_bin();
- if ($b =~ /0b1(0+)$/)
- {
- my $count = CORE::length($1); # 0b100 => len('00') => 2
- my $cnt = $count; # counter for loop
- my $shift = $self->new(6);
- $x->blsft($shift); # add some zeros (even amount)
- while ($cnt-- > 0)
- {
- # 'inflate' $X by adding more zeros
- $x->blsft($shift);
- # calculate sqrt($x), $x is now a bit too big, again. In the next
- # round we make even bigger, again.
- $x->bsqrt($x);
- }
- # $x is still to big, so truncate result
- $x->brsft($shift);
- }
- else
- {
- # trial computation by starting with 2,4,6,8,10 etc until we overstep
- my $step;
- my $trial = $self->new(2);
- my $two = $self->new(2);
- my $s_128 = $self->new(128);
-
- local undef $Math::BigInt::accuracy;
- local undef $Math::BigInt::precision;
-
- # while still to do more than X steps
- do
- {
- $step = $self->new(2);
- while ( $trial->copy->bpow($y)->bacmp($x) < 0)
- {
- $step->bmul($two);
- $trial->badd($step);
- }
-
- # hit exactly?
- if ( $trial->copy->bpow($y)->bacmp($x) == 0)
- {
- $x->{value} = $trial->{value}; # make copy while preserving ref to $x
- return $x->round(@r);
- }
- # overstepped, so go back on step
- $trial->bsub($step);
- } while ($step > $s_128);
-
- $step = $two->copy();
- while ( $trial->copy->bpow($y)->bacmp($x) < 0)
- {
- $trial->badd($step);
- }
-
- # not hit exactly?
- if ( $x->bacmp( $trial->copy()->bpow($y) ) < 0)
- {
- $trial->bdec();
- }
- # copy result into $x (preserve ref)
- $x->{value} = $trial->{value};
- }
- $x->round(@r);
- }
-
-sub __emu_as_hex
- {
- my ($self,$x,$s) = @_;
-
- return '0x0' if $x->is_zero();
-
- my $x1 = $x->copy()->babs(); my ($xr,$x10000,$h,$es);
- if ($] >= 5.006)
- {
- $x10000 = $self->new (0x10000); $h = 'h4';
- }
- else
- {
- $x10000 = $self->new (0x1000); $h = 'h3';
- }
- while (!$x1->is_zero())
- {
- ($x1, $xr) = bdiv($x1,$x10000);
- $es .= unpack($h,pack('v',$xr->numify()));
- }
- $es = reverse $es;
- $es =~ s/^[0]+//; # strip leading zeros
- $s . '0x' . $es;
- }
-
-sub __emu_as_bin
- {
- my ($self,$x,$s) = @_;
-
- return '0b0' if $x->is_zero();
-
- my $x1 = $x->copy()->babs(); my ($xr,$x10000,$b,$es);
- if ($] >= 5.006)
- {
- $x10000 = $self->new (0x10000); $b = 'b16';
- }
- else
- {
- $x10000 = $self->new (0x1000); $b = 'b12';
- }
- while (!$x1->is_zero())
- {
- ($x1, $xr) = bdiv($x1,$x10000);
- $es .= unpack($b,pack('v',$xr->numify()));
- }
- $es = reverse $es;
- $es =~ s/^[0]+//; # strip leading zeros
- $s . '0b' . $es;
- }
-
##############################################################################
##############################################################################
=head1 AUTHORS
-(c) Tels http://bloodgate.com 2003 - based on BigInt code by
+(c) Tels http://bloodgate.com 2003, 2004 - based on BigInt code by
Tels from 2001-2003.
=head1 SEE ALSO
projects / p5sagit/p5-mst-13.2.git / blobdiff