X-Git-Url: http://git.shadowcat.co.uk/gitweb/gitweb.cgi?a=blobdiff_plain;f=lib%2FMath%2FBigInt%2FCalcEmu.pm;h=9f7fd16434a22870a946d47430139d8cac0cc8ed;hb=29489e7c741791873ea464cb7e13d2b5a19577a7;hp=d5d1734b220ffe6df0644023c3fe7313f495c5be;hpb=3a427a117ca05296a5c3d5d88415dc040792253b;p=p5sagit%2Fp5-mst-13.2.git diff --git a/lib/Math/BigInt/CalcEmu.pm b/lib/Math/BigInt/CalcEmu.pm index d5d1734..9f7fd16 100644 --- a/lib/Math/BigInt/CalcEmu.pm +++ b/lib/Math/BigInt/CalcEmu.pm @@ -1,12 +1,13 @@ -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. @@ -17,141 +18,6 @@ BEGIN $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) = @_; @@ -193,26 +59,29 @@ sub __emu_band $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) { @@ -223,20 +92,12 @@ sub __emu_band $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; @@ -316,16 +177,13 @@ sub __emu_bior $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; @@ -405,16 +263,9 @@ sub __emu_bxor $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} = '+'; @@ -425,170 +276,6 @@ sub __emu_bxor $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; - } - ############################################################################## ############################################################################## @@ -621,7 +308,7 @@ the same terms as Perl itself. =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