From: Tels Date: Tue, 23 Dec 2003 01:09:23 +0000 (+0100) Subject: BigInt v1.68 - pre-release X-Git-Url: http://git.shadowcat.co.uk/gitweb/gitweb.cgi?a=commitdiff_plain;h=b282a5527464951004e354d07709b58fcb3bdad0;p=p5sagit%2Fp5-mst-13.2.git BigInt v1.68 - pre-release Message-Id: <200312230106.27661@bloodgate.com> p4raw-id: //depot/perl@21956 --- diff --git a/MANIFEST b/MANIFEST index acf2f20..e851d1d 100644 --- a/MANIFEST +++ b/MANIFEST @@ -1325,10 +1325,12 @@ lib/look.pl A "look" equivalent lib/Math/BigFloat.pm An arbitrary precision floating-point arithmetic package lib/Math/BigFloat/Trace.pm bignum tracing lib/Math/BigInt/Calc.pm Pure Perl module to support Math::BigInt +lib/Math/BigInt/CalcEmu.pm Pure Perl module to support Math::BigInt lib/Math/BigInt.pm An arbitrary precision integer arithmetic package lib/Math/BigInt/t/bare_mbf.t Test MBF under Math::BigInt::BareCalc lib/Math/BigInt/t/bare_mbi.t Test MBI under Math::BigInt::BareCalc lib/Math/BigInt/t/bare_mif.t Rounding tests under BareCalc +lib/Math/BigInt/t/alias.inc Support for BigInt tests lib/Math/BigInt/t/bigfltpm.inc Shared tests for bigfltpm.t and sub_mbf.t lib/Math/BigInt/t/bigfltpm.t See if BigFloat.pm works lib/Math/BigInt/t/bigintc.t See if BigInt/Calc.pm works @@ -1344,8 +1346,10 @@ lib/Math/BigInt/t/downgrade.t Test if use Math::BigInt(); under downgrade works lib/Math/BigInt/t/fallback.t Test Math::BigInt lib/Math/BigInt/t/inf_nan.t Special tests for inf and NaN handling lib/Math/BigInt/t/isa.t Test for Math::BigInt inheritance +lib/Math/BigInt/t/mbf_ali.t Tests for BigFloat lib/Math/BigInt/t/mbimbf.inc Actual BigInt/BigFloat accuracy, precision and fallback, round_mode tests lib/Math/BigInt/t/mbimbf.t BigInt/BigFloat accuracy, precision and fallback, round_mode +lib/Math/BigInt/t/mbi_ali.t Tests for BigInt lib/Math/BigInt/t/mbi_rand.t Test Math::BigInt randomly lib/Math/BigInt/Trace.pm bignum tracing lib/Math/BigInt/t/req_mbf0.t test: require Math::BigFloat; ->bzero(); @@ -1355,6 +1359,7 @@ lib/Math/BigInt/t/req_mbfi.t test: require Math::BigFloat; ->binf(); lib/Math/BigInt/t/req_mbfn.t test: require Math::BigFloat; ->new(); lib/Math/BigInt/t/req_mbfw.t require Math::BigFloat; import ( with => ); lib/Math/BigInt/t/require.t Test if require Math::BigInt works +lib/Math/BigInt/t/sub_ali.t Tests for aliases in BigInt subclasses lib/Math/BigInt/t/sub_mbf.t Empty subclass test of BigFloat lib/Math/BigInt/t/sub_mbi.t Empty subclass test of BigInt lib/Math/BigInt/t/sub_mif.t Test A & P with subclasses using mbimbf.inc diff --git a/lib/Math/BigFloat.pm b/lib/Math/BigFloat.pm index 3b8d5a6..9071648 100644 --- a/lib/Math/BigFloat.pm +++ b/lib/Math/BigFloat.pm @@ -12,16 +12,15 @@ package Math::BigFloat; # _p: precision # _f: flags, used to signal MBI not to touch our private parts -$VERSION = '1.41'; +$VERSION = '1.42'; require 5.005; use Exporter; @ISA = qw(Exporter Math::BigInt); use strict; -use vars qw/$AUTOLOAD $accuracy $precision $div_scale $round_mode $rnd_mode/; -use vars qw/$upgrade $downgrade/; -# the following are internal and should never be accessed from the outside -use vars qw/$_trap_nan $_trap_inf/; +# $_trap_inf and $_trap_nan are internal and should never be accessed from the outside +use vars qw/$AUTOLOAD $accuracy $precision $div_scale $round_mode $rnd_mode + $upgrade $downgrade $_trap_nan $_trap_inf/; my $class = "Math::BigFloat"; use overload @@ -50,7 +49,7 @@ my $MBI = 'Math::BigInt'; # the package we are using for our private parts # the following are private and not to be used from the outside: -use constant MB_NEVER_ROUND => 0x0001; +sub MB_NEVER_ROUND () { 0x0001; } # are NaNs ok? (otherwise it dies when encountering an NaN) set w/ config() $_trap_nan = 0; @@ -151,6 +150,7 @@ sub new return $self->bnorm(); } #print "new string '$wanted'\n"; + my ($mis,$miv,$mfv,$es,$ev) = Math::BigInt::_split(\$wanted); if (!ref $mis) { @@ -172,10 +172,33 @@ sub new # undef,undef to signal MBI that we don't need no bloody rounding $self->{_e} = $MBI->new("$$es$$ev",undef,undef); # exponent $self->{_m} = $MBI->new("$$miv$$mfv",undef,undef); # create mant. - # print $self->{_e}, " ", $self->{_m},"\n"; + + # this is to prevent automatically rounding when MBI's globals are set + $self->{_m}->{_f} = MB_NEVER_ROUND; + $self->{_e}->{_f} = MB_NEVER_ROUND; + # 3.123E0 = 3123E-3, and 3.123E-2 => 3123E-5 - $self->{_e} -= CORE::length($$mfv) if CORE::length($$mfv) != 0; + $self->{_e}->bsub( $MBI->new(CORE::length($$mfv),undef,undef)) + if CORE::length($$mfv) != 0; $self->{sign} = $$mis; + + #print "$$miv$$mfv $$es$$ev\n"; + + # we can only have trailing zeros on the mantissa of $$mfv eq '' + if (CORE::length($$mfv) == 0) + { + my $zeros = $self->{_m}->_trailing_zeros(); # correct for trailing zeros + if ($zeros != 0) + { + $self->{_m}->brsft($zeros,10); $self->{_e}->badd($MBI->new($zeros)); + } + } +# else +# { + # for something like 0Ey, set y to 1, and -0 => +0 + $self->{sign} = '+', $self->{_e}->bone() if $self->{_m}->is_zero(); +# } + return $self->round(@r) if !$downgrade; } # if downgrade, inf, NaN or integers go down @@ -352,8 +375,8 @@ sub bsstr sub numify { # Make a number from a BigFloat object - # simple return string and let Perl's atoi()/atof() handle the rest - my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_); + # simple return a string and let Perl's atoi()/atof() handle the rest + my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_); $x->bsstr(); } @@ -361,7 +384,7 @@ sub numify # public stuff (usually prefixed with "b") # tels 2001-08-04 -# todo: this must be overwritten and return NaN for non-integer values +# XXX TODO this must be overwritten and return NaN for non-integer values # band(), bior(), bxor(), too #sub bnot # { @@ -371,7 +394,6 @@ sub numify sub bcmp { # Compares 2 values. Returns one of undef, <0, =0, >0. (suitable for sort) - # (BFLOAT or num_str, BFLOAT or num_str) return cond_code # set up parameters my ($self,$x,$y) = (ref($_[0]),@_); @@ -437,7 +459,6 @@ sub bacmp { # Compares 2 values, ignoring their signs. # Returns one of undef, <0, =0, >0. (suitable for sort) - # (BFLOAT or num_str, BFLOAT or num_str) return cond_code # set up parameters my ($self,$x,$y) = (ref($_[0]),@_); @@ -573,7 +594,6 @@ sub bsub ($self,$x,$y,$a,$p,$r) = objectify(2,@_); } - # XXX TODO: remove? if ($y->is_zero()) # still round for not adding zero { return $x->round($a,$p,$r); @@ -589,42 +609,45 @@ sub bsub sub binc { # increment arg by one - my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_); + my ($self,$x,@r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_); if ($x->{_e}->sign() eq '-') { - return $x->badd($self->bone(),$a,$p,$r); # digits after dot + return $x->badd($self->bone(),@r); # digits after dot } - if (!$x->{_e}->is_zero()) + if (!$x->{_e}->is_zero()) # _e == 0 for NaN, inf, -inf { + # 1e2 => 100, so after the shift below _m has a '0' as last digit $x->{_m}->blsft($x->{_e},10); # 1e2 => 100 - $x->{_e}->bzero(); + $x->{_e}->bzero(); # normalize + # we know that the last digit of $x will be '1' or '9', depending on the + # sign } # now $x->{_e} == 0 if ($x->{sign} eq '+') { $x->{_m}->binc(); - return $x->bnorm()->bround($a,$p,$r); + return $x->bnorm()->bround(@r); } elsif ($x->{sign} eq '-') { $x->{_m}->bdec(); $x->{sign} = '+' if $x->{_m}->is_zero(); # -1 +1 => -0 => +0 - return $x->bnorm()->bround($a,$p,$r); + return $x->bnorm()->bround(@r); } # inf, nan handling etc - $x->badd($self->bone(),$a,$p,$r); # does round + $x->badd($self->bone(),@r); # badd() does round } sub bdec { # decrement arg by one - my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_); + my ($self,$x,@r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_); if ($x->{_e}->sign() eq '-') { - return $x->badd($self->bone('-'),$a,$p,$r); # digits after dot + return $x->badd($self->bone('-'),@r); # digits after dot } if (!$x->{_e}->is_zero()) @@ -640,16 +663,16 @@ sub bdec $x->{_m}->binc(); $x->{sign} = '-' if $zero; # 0 => 1 => -1 $x->{sign} = '+' if $x->{_m}->is_zero(); # -1 +1 => -0 => +0 - return $x->bnorm()->round($a,$p,$r); + return $x->bnorm()->round(@r); } # > 0 elsif ($x->{sign} eq '+') { $x->{_m}->bdec(); - return $x->bnorm()->round($a,$p,$r); + return $x->bnorm()->round(@r); } # inf, nan handling etc - $x->badd($self->bone('-'),$a,$p,$r); # does round + $x->badd($self->bone('-'),@r); # does round } sub DEBUG () { 0; } @@ -718,15 +741,40 @@ sub blog $x = Math::BigFloat->new($x); $self = ref($x); } - # first calculate the log to base e (using reduction by 10 (and probably 2)) - $self->_log_10($x,$scale); - - # and if a different base was requested, convert it - if (defined $base) + + my $done = 0; + + # If the base is defined and an integer, try to calculate integer result + # first. This is very fast, and in case the real result was found, we can + # stop right here. + if (defined $base && $base->is_int() && $x->is_int()) + { + my $int = $x->{_m}->copy(); + $int->blsft($x->{_e},10) unless $x->{_e}->is_zero(); + $int->blog($base->as_number()); + # if ($exact) + if ($base->copy()->bpow($int) == $x) + { + # found result, return it + $x->{_m} = $int; + $x->{_e} = $MBI->bzero(); + $x->bnorm(); + $done = 1; + } + } + + if ($done == 0) { - $base = Math::BigFloat->new($base) unless $base->isa('Math::BigFloat'); - # not ln, but some other base (don't modify $base) - $x->bdiv( $base->copy()->blog(undef,$scale), $scale ); + # first calculate the log to base e (using reduction by 10 (and probably 2)) + $self->_log_10($x,$scale); + + # and if a different base was requested, convert it + if (defined $base) + { + $base = Math::BigFloat->new($base) unless $base->isa('Math::BigFloat'); + # not ln, but some other base (don't modify $base) + $x->bdiv( $base->copy()->blog(undef,$scale), $scale ); + } } # shortcut to not run through _find_round_parameters again @@ -1541,20 +1589,23 @@ sub bfac { # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT # compute factorial number, modifies first argument - my ($self,$x,@r) = objectify(1,@_); + # set up parameters + my ($self,$x,@r) = (ref($_[0]),@_); + # objectify is costly, so avoid it + ($self,$x,@r) = objectify(1,@_) if !ref($x); + + return $x if $x->{sign} eq '+inf'; # inf => inf return $x->bnan() if (($x->{sign} ne '+') || # inf, NaN, <0 etc => NaN ($x->{_e}->{sign} ne '+')); # digits after dot? # use BigInt's bfac() for faster calc - if (! _is_zero_or_one($x->{_e})) + if (! $x->{_e}->is_zero()) { - $x->{_m}->blsft($x->{_e},10); # unnorm - $x->{_e}->bzero(); # norm again + $x->{_m}->blsft($x->{_e},10); # change 12e1 to 120e0 + $x->{_e}->bzero(); } - $x->{_m}->blsft($x->{_e},10); # un-norm m - $x->{_e}->bzero(); # norm again $x->{_m}->bfac(); # calculate factorial $x->bnorm()->round(@r); # norm again and round result } @@ -1948,7 +1999,7 @@ sub blsft sub DESTROY { - # going through AUTOLOAD for every DESTROY is costly, so avoid it by empty sub + # going through AUTOLOAD for every DESTROY is costly, avoid it by empty sub } sub AUTOLOAD @@ -2123,16 +2174,19 @@ sub bnorm return $x if $x->{sign} !~ /^[+-]$/; # inf, nan etc -# if (!$x->{_m}->is_odd()) -# { - my $zeros = $x->{_m}->_trailing_zeros(); # correct for trailing zeros - if ($zeros != 0) - { - $x->{_m}->brsft($zeros,10); $x->{_e}->badd($zeros); - } - # for something like 0Ey, set y to 1, and -0 => +0 + my $zeros = $x->{_m}->_trailing_zeros(); # correct for trailing zeros + if ($zeros != 0) + { + my $z = $MBI->new($zeros,undef,undef); + $x->{_m}->brsft($z,10); $x->{_e}->badd($z); + } + else + { + # $x can only be 0Ey if there are no trailing zeros ('0' has 0 trailing + # zeros). So, for something like 0Ey, set y to 1, and -0 => +0 $x->{sign} = '+', $x->{_e}->bone() if $x->{_m}->is_zero(); -# } + } + # this is to prevent automatically rounding when MBI's globals are set $x->{_m}->{_f} = MB_NEVER_ROUND; $x->{_e}->{_f} = MB_NEVER_ROUND; diff --git a/lib/Math/BigInt.pm b/lib/Math/BigInt.pm index 1b16600..aeee0e2 100644 --- a/lib/Math/BigInt.pm +++ b/lib/Math/BigInt.pm @@ -18,14 +18,14 @@ package Math::BigInt; my $class = "Math::BigInt"; require 5.005; -$VERSION = '1.67'; +$VERSION = '1.68'; use Exporter; @ISA = qw( Exporter ); @EXPORT_OK = qw( objectify bgcd blcm); -use vars qw/$round_mode $accuracy $precision $div_scale $rnd_mode/; -use vars qw/$upgrade $downgrade/; -# the following are internal and should never be accessed from the outside -use vars qw/$_trap_nan $_trap_inf/; +# _trap_inf and _trap_nan are internal and should never be accessed from the +# outside +use vars qw/$round_mode $accuracy $precision $div_scale $rnd_mode + $upgrade $downgrade $_trap_nan $_trap_inf/; use strict; # Inside overload, the first arg is always an object. If the original code had @@ -66,11 +66,6 @@ use overload "$_[1]" cmp $_[0]->bstr() : $_[0]->bstr() cmp "$_[1]" }, -#'cos' => sub { -# require Math::Big; -# return Math::Big::cos($_[0], ref($_[0])->accuracy()); -# }, - # make cos()/sin()/exp() "work" with BigInt's or subclasses 'cos' => sub { cos($_[0]->numify()) }, 'sin' => sub { sin($_[0]->numify()) }, @@ -152,7 +147,7 @@ $downgrade = undef; # default is no downgrade # these are internally, and not to be used from the outside -use constant MB_NEVER_ROUND => 0x0001; +sub MB_NEVER_ROUND () { 0x0001; } $_trap_nan = 0; # are NaNs ok? set w/ config() $_trap_inf = 0; # are infs ok? set w/ config() @@ -164,6 +159,9 @@ my %CAN; # cache for $CALC->can(...) my $IMPORT = 0; # was import() called yet? # used to make require work +my $EMU_LIB = 'Math/BigInt/CalcEmu.pm'; # emulate low-level math +my $EMU = 'Math::BigInt::CalcEmu'; # emulate low-level math + ############################################################################## # the old code had $rnd_mode, so we need to support it, too @@ -172,7 +170,16 @@ sub TIESCALAR { my ($class) = @_; bless \$round_mode, $class; } sub FETCH { return $round_mode; } sub STORE { $rnd_mode = $_[0]->round_mode($_[1]); } -BEGIN { tie $rnd_mode, 'Math::BigInt'; } +BEGIN + { + # tie to enable $rnd_mode to work transparently + tie $rnd_mode, 'Math::BigInt'; + + # set up some handy alias names + *as_int = \&as_number; + *is_pos = \&is_positive; + *is_neg = \&is_negative; + } ############################################################################## @@ -798,8 +805,9 @@ sub bsstr return 'inf'; # +inf } my ($m,$e) = $x->parts(); - my $sign = 'e+'; # e can only be positive - return $m->bstr().$sign.$e->bstr(); + #$m->bstr() . 'e+' . $e->bstr(); # e can only be positive in BigInt + # 'e+' because E can only be positive in BigInt + $m->bstr() . 'e+' . ${$CALC->_str($e->{value})}; } sub bstr @@ -814,7 +822,7 @@ sub bstr return 'inf'; # +inf } my $es = ''; $es = $x->{sign} if $x->{sign} eq '-'; - return $es.${$CALC->_str($x->{value})}; + $es.${$CALC->_str($x->{value})}; } sub numify @@ -834,7 +842,7 @@ sub numify sub sign { # return the sign of the number: +/-/-inf/+inf/NaN - my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_); + my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_); $x->{sign}; } @@ -960,7 +968,7 @@ sub round $r = ${"$c\::round_mode"} unless defined $r; if ($r !~ /^(even|odd|\+inf|\-inf|zero|trunc)$/) { - + require Carp; Carp::croak ("Unknown round mode '$r'"); } # now round, by calling either fround or ffround: @@ -979,7 +987,7 @@ sub bnorm { # (numstr or BINT) return BINT # Normalize number -- no-op here - my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_); + my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_); $x; } @@ -1050,7 +1058,7 @@ sub bcmp } # $x && $y both < 0 - $CALC->_acmp($y->{value},$x->{value}); # swaped (lib returns 0,1,-1) + $CALC->_acmp($y->{value},$x->{value}); # swaped acmp (lib returns 0,1,-1) } sub bacmp @@ -1116,12 +1124,11 @@ sub badd return $x; } - my ($sx, $sy) = ( $x->{sign}, $y->{sign} ); # get signs + my ($sx, $sy) = ( $x->{sign}, $y->{sign} ); # get signs if ($sx eq $sy) { $x->{value} = $CALC->_add($x->{value},$y->{value}); # same sign, abs add - $x->{sign} = $sx; } else { @@ -1140,7 +1147,6 @@ sub badd else # a < 0 { $x->{value} = $CALC->_sub($x->{value}, $y->{value}); # abs sub - $x->{sign} = $sx; } } $x->round(@r) if !exists $x->{_f} || $x->{_f} & MB_NEVER_ROUND == 0; @@ -1204,29 +1210,32 @@ sub binc sub bdec { # decrement arg by one - my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_); + my ($self,$x,@r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_); return $x if $x->modify('bdec'); - my $zero = $CALC->_is_zero($x->{value}) && $x->{sign} eq '+'; - # <= 0 - if (($x->{sign} eq '-') || $zero) + if ($x->{sign} eq '-') { + # < 0 $x->{value} = $CALC->_inc($x->{value}); - $x->{sign} = '-' if $zero; # 0 => 1 => -1 - $x->{sign} = '+' if $CALC->_is_zero($x->{value}); # -1 +1 => -0 => +0 - $x->round($a,$p,$r) if !exists $x->{_f} || $x->{_f} & MB_NEVER_ROUND == 0; - return $x; - } - # > 0 - elsif ($x->{sign} eq '+') + } + else { - $x->{value} = $CALC->_dec($x->{value}); - $x->round($a,$p,$r) if !exists $x->{_f} || $x->{_f} & MB_NEVER_ROUND == 0; - return $x; + return $x->badd($self->bone('-'),@r) unless $x->{sign} eq '+'; # inf/NaN + # >= 0 + if ($CALC->_is_zero($x->{value})) + { + # == 0 + $x->{value} = $CALC->_one(); $x->{sign} = '-'; # 0 => -1 + } + else + { + # > 0 + $x->{value} = $CALC->_dec($x->{value}); + } } - # inf, nan handling etc - $x->badd($self->bone('-'),$a,$p,$r); # badd does round - } + $x->round(@r) if !exists $x->{_f} || $x->{_f} & MB_NEVER_ROUND == 0; + $x; + } sub blog { @@ -1238,7 +1247,7 @@ sub blog # objectify is costly, so avoid it if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1]))) { - ($self,$x,$base,@r) = objectify(2,@_); + ($self,$x,$base,@r) = objectify(2,$class,@_); } # inf, -inf, NaN, <0 => NaN @@ -1256,26 +1265,8 @@ sub blog return $x->round(@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); + require $EMU_LIB; + __emu_blog($self,$x,$base,@r); } sub blcm @@ -1661,53 +1652,15 @@ sub bmodinv { my $sign; ($x->{value},$sign) = $CALC->_modinv($x->{value},$y->{value}); - $x->bnan() if !defined $x->{value}; # in case no GCD found - return $x if !defined $sign; # already real result - $x->{sign} = $sign; # flip/flop see below - $x->bmod($y); # calc real result + return $x->bnan() if !defined $x->{value}; # in case no GCD found + return $x if !defined $sign; # already real result + $x->{sign} = $sign; # flip/flop see below + $x->bmod($y); # calc real result return $x; } - 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*. - 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; + require $EMU_LIB; + __emu_bmodinv($self,$x,$y,@r); } sub bmodpow @@ -1742,24 +1695,8 @@ sub bmodpow return $num; } - # 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; + require $EMU_LIB; + __emu_bmodpow($self,$num,$exp,$mod,@r); } ############################################################################### @@ -1772,7 +1709,8 @@ sub bfac return $x if $x->modify('bfac'); - return $x->bnan() if $x->{sign} ne '+'; # inf, NnN, <0 etc => NaN + return $x if $x->{sign} eq '+inf'; # inf => inf + return $x->bnan() if $x->{sign} ne '+'; # NaN, <0 etc => NaN if ($CAN{fac}) { @@ -1780,17 +1718,8 @@ sub bfac return $x->round(@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 + require $EMU_LIB; + __emu_bfac($self,$x,@r); } sub bpow @@ -1815,8 +1744,9 @@ sub bpow $r[3] = $y; # no push! return $x if $x->{sign} =~ /^[+-]inf$/; # -inf/+inf ** x return $x->bnan() if $x->{sign} eq $nan || $y->{sign} eq $nan; - return $x->bone('+',@r) if $y->is_zero(); - return $x->round(@r) if $x->is_one() || $y->is_one(); + + # cases 0 ** Y, X ** 0, X ** 1, 1 ** Y are handled by Calc or Emu + if ($x->{sign} eq '-' && $CALC->_is_one($x->{value})) { # if $x == -1 and odd/even y => +1/-1 @@ -1825,44 +1755,18 @@ sub bpow } # 1 ** -y => 1 / (1 ** |y|) # so do test for negative $y after above's clause - return $x->bnan() if $y->{sign} eq '-'; - return $x->round(@r) if $x->is_zero(); # 0**y => 0 (if not y <= 0) + return $x->bnan() if $y->{sign} eq '-' && !$x->is_one(); if ($CAN{pow}) { $x->{value} = $CALC->_pow($x->{value},$y->{value}); + $x->{sign} = '+' if $CALC->_is_zero($y->{value}); $x->round(@r) if !exists $x->{_f} || $x->{_f} & MB_NEVER_ROUND == 0; return $x; } -# based on the assumption that shifting in base 10 is fast, and that mul -# works faster if numbers are small: we count trailing zeros (this step is -# O(1)..O(N), but in case of O(N) we save much more time due to this), -# stripping them out of the multiplication, and add $count * $y zeros -# afterwards like this: -# 300 ** 3 == 300*300*300 == 3*3*3 . '0' x 2 * 3 == 27 . '0' x 6 -# creates deep recursion since brsft/blsft use bpow sometimes. -# my $zeros = $x->_trailing_zeros(); -# if ($zeros > 0) -# { -# $x->brsft($zeros,10); # remove zeros -# $x->bpow($y); # recursion (will not branch into here again) -# $zeros = $y * $zeros; # real number of zeros to add -# $x->blsft($zeros,10); -# return $x->round(@r); -# } - - 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; + require $EMU_LIB; + __emu_bpow($self,$x,$y,@r); } sub blsft @@ -1890,7 +1794,7 @@ sub blsft $x->{value} = $t; return $x->round(@r); } # fallback - return $x->bmul( $self->bpow($n, $y, @r), @r ); + $x->bmul( $self->bpow($n, $y, @r), @r ); } sub brsft @@ -1946,6 +1850,7 @@ sub brsft $x->{value} = $res->{value}; # take over value return $x->round(@r); # we are done now, magic, isn't? } + # x < 0, n == 2, y == 1 $x->bdec(); # n == 2, but $y == 1: this fixes it } @@ -1976,12 +1881,11 @@ sub band return $x if $x->modify('band'); $r[3] = $y; # no push! - local $Math::BigInt::upgrade = undef; return $x->bnan() if ($x->{sign} !~ /^[+-]$/ || $y->{sign} !~ /^[+-]$/); - my $sx = 1; $sx = -1 if $x->{sign} eq '-'; - my $sy = 1; $sy = -1 if $y->{sign} eq '-'; + my $sx = $x->{sign} eq '+' ? 1 : -1; + my $sy = $y->{sign} eq '+' ? 1 : -1; if ($CAN{and} && $sx == 1 && $sy == 1) { @@ -1994,92 +1898,9 @@ sub band $x->{value} = $CALC->_signed_and($x->{value},$y->{value},$sx,$sy); return $x->round(@r); } - - return $x->bzero(@r) if $y->is_zero() || $x->is_zero(); - - my $sign = 0; # sign of result - $sign = 1 if ($x->{sign} eq '-') && ($y->{sign} eq '-'); - - my ($bx,$by); - - if ($sx == -1) # if x is negative - { - # two's complement: inc and flip all "bits" in $bx - $bx = $x->binc()->as_hex(); # -1 => 0, -2 => 1, -3 => 2 etc - $bx =~ s/-?0x//; - $bx =~ tr/0123456789abcdef/\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00/; - } - else - { - $bx = $x->as_hex(); # get binary representation - $bx =~ s/-?0x//; - $bx =~ tr/fedcba9876543210/\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00/; - } - if ($sy == -1) # if y is negative - { - # two's complement: inc and flip all "bits" in $by - $by = $y->copy()->binc()->as_hex(); # -1 => 0, -2 => 1, -3 => 2 etc - $by =~ s/-?0x//; - $by =~ tr/0123456789abcdef/\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00/; - } - else - { - $by = $y->as_hex(); # get binary representation - $by =~ s/-?0x//; - $by =~ tr/fedcba9876543210/\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00/; - } - # now we have bit-strings from X and Y, reverse them for padding - $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) - my $diff = CORE::length($bx) - CORE::length($by); - if ($diff > 0) - { - $bx = substr($bx,0,CORE::length($by)); - } - elsif ($diff < 0) - { - $by = substr($by,0,CORE::length($bx)); - } - - # 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 - # 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/0123456789abcdef/; - } - else - { - $bx =~ tr/\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00/fedcba9876543210/; - } - - $bx = '0x' . $bx; - if ($CAN{from_hex}) - { - $x->{value} = $CALC->_from_hex( \$bx ); - } - else - { - $r = $self->new($bx); - $x->{value} = $r->{value}; - } - - # calculate sign of result - $x->{sign} = '+'; - $x->{sign} = '-' if $sx == $sy && $sx == -1 && !$x->is_zero(); - - $x->bdec() if $sign == 1; - - $x->round(@r); + + require $EMU_LIB; + __emu_band($self,$x,$y,$sx,$sy,@r); } sub bior @@ -2102,8 +1923,8 @@ sub bior return $x->bnan() if ($x->{sign} !~ /^[+-]$/ || $y->{sign} !~ /^[+-]$/); - my $sx = 1; $sx = -1 if $x->{sign} eq '-'; - my $sy = 1; $sy = -1 if $y->{sign} eq '-'; + my $sx = $x->{sign} eq '+' ? 1 : -1; + my $sy = $y->{sign} eq '+' ? 1 : -1; # the sign of X follows the sign of X, e.g. sign of Y irrelevant for bior() @@ -2114,96 +1935,15 @@ sub bior return $x->round(@r); } - # if lib can do negatvie values, so use it + # if lib can do negative values, let it handle this if ($CAN{signed_or}) { $x->{value} = $CALC->_signed_or($x->{value},$y->{value},$sx,$sy); return $x->round(@r); } - return $x->round(@r) if $y->is_zero(); - - my $sign = 0; # sign of result - $sign = 1 if ($x->{sign} eq '-') || ($y->{sign} eq '-'); - - my ($bx,$by); - - if ($sx == -1) # if x is negative - { - # two's complement: inc and flip all "bits" in $bx - $bx = $x->binc()->as_hex(); # -1 => 0, -2 => 1, -3 => 2 etc - $bx =~ s/-?0x//; - $bx =~ tr/0123456789abcdef/\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00/; - } - else - { - $bx = $x->as_hex(); # get binary representation - $bx =~ s/-?0x//; - $bx =~ tr/fedcba9876543210/\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00/; - } - if ($sy == -1) # if y is negative - { - # two's complement: inc and flip all "bits" in $by - $by = $y->copy()->binc()->as_hex(); # -1 => 0, -2 => 1, -3 => 2 etc - $by =~ s/-?0x//; - $by =~ tr/0123456789abcdef/\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00/; - } - else - { - $by = $y->as_hex(); # get binary representation - $by =~ s/-?0x//; - $by =~ tr/fedcba9876543210/\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00/; - } - # now we have bit-strings from X and Y, reverse them for padding - $bx = reverse $bx; - $by = reverse $by; - - # 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) - { - $by .= $yy x $diff; - } - elsif ($diff < 0) - { - $bx .= $xx x abs($diff); - } - - # or 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 - # 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/0123456789abcdef/; - } - else - { - $bx =~ tr/\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00/fedcba9876543210/; - } - - $bx = '0x' . $bx; - if ($CAN{from_hex}) - { - $x->{value} = $CALC->_from_hex( \$bx ); - } - else - { - $r = $self->new($bx); - $x->{value} = $r->{value}; - } - - # if one of X or Y was negative, we need to decrement result - $x->bdec() if $sign == 1; - - $x->round(@r); + require $EMU_LIB; + __emu_bior($self,$x,$y,$sx,$sy,@r); } sub bxor @@ -2222,12 +1962,10 @@ sub bxor return $x if $x->modify('bxor'); $r[3] = $y; # no push! - local $Math::BigInt::upgrade = undef; - return $x->bnan() if ($x->{sign} !~ /^[+-]$/ || $y->{sign} !~ /^[+-]$/); - my $sx = 1; $sx = -1 if $x->{sign} eq '-'; - my $sy = 1; $sy = -1 if $y->{sign} eq '-'; + my $sx = $x->{sign} eq '+' ? 1 : -1; + my $sy = $y->{sign} eq '+' ? 1 : -1; # don't use lib for negative values if ($CAN{xor} && $sx == 1 && $sy == 1) @@ -2236,104 +1974,20 @@ sub bxor return $x->round(@r); } - # if lib can do negatvie values, so use it + # if lib can do negative values, let it handle this if ($CAN{signed_xor}) { $x->{value} = $CALC->_signed_xor($x->{value},$y->{value},$sx,$sy); return $x->round(@r); } - return $x->round(@r) if $y->is_zero(); - - my $sign = 0; # sign of result - $sign = 1 if $x->{sign} ne $y->{sign}; - - my ($bx,$by); - - if ($sx == -1) # if x is negative - { - # two's complement: inc and flip all "bits" in $bx - $bx = $x->binc()->as_hex(); # -1 => 0, -2 => 1, -3 => 2 etc - $bx =~ s/-?0x//; - $bx =~ tr/0123456789abcdef/\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00/; - } - else - { - $bx = $x->as_hex(); # get binary representation - $bx =~ s/-?0x//; - $bx =~ tr/fedcba9876543210/\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00/; - } - if ($sy == -1) # if y is negative - { - # two's complement: inc and flip all "bits" in $by - $by = $y->copy()->binc()->as_hex(); # -1 => 0, -2 => 1, -3 => 2 etc - $by =~ s/-?0x//; - $by =~ tr/0123456789abcdef/\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00/; - } - else - { - $by = $y->as_hex(); # get binary representation - $by =~ s/-?0x//; - $by =~ tr/fedcba9876543210/\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00/; - } - # now we have bit-strings from X and Y, reverse them for padding - $bx = reverse $bx; - $by = reverse $by; - - # 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) - { - $by .= $yy x $diff; - } - elsif ($diff < 0) - { - $bx .= $xx x abs($diff); - } - - # xor 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 - # 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/0123456789abcdef/; - } - else - { - $bx =~ tr/\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00/fedcba9876543210/; - } - - $bx = '0x' . $bx; - if ($CAN{from_hex}) - { - $x->{value} = $CALC->_from_hex( \$bx ); - } - else - { - $r = $self->new($bx); - $x->{value} = $r->{value}; - } - - # calculate sign of result - $x->{sign} = '+'; - $x->{sign} = '-' if $sx != $sy && !$x->is_zero(); - - $x->bdec() if $sign == 1; - - $x->round(@r); + require $EMU_LIB; + __emu_bxor($self,$x,$y,$sx,$sy,@r); } sub length { - my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_); + my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_); my $e = $CALC->_len($x->{value}); wantarray ? ($e,0) : $e; @@ -2349,7 +2003,7 @@ sub digit sub _trailing_zeros { - # return the amount of trailing zeros in $x + # return the amount of trailing zeros in $x (as scalar) my $x = shift; $x = $class->new($x) unless ref $x; @@ -2381,28 +2035,8 @@ sub bsqrt return $x->round(@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); + require $EMU_LIB; + __emu_bsqrt($self,$x,@r); } sub broot @@ -2437,48 +2071,8 @@ sub broot return $x->round(@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($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 - { - # Should compute a guess of the result (by rule of thumb), then improve it - # via Newton's method or something similiar. - # XXX TODO - warn ('broot() not fully implemented in BigInt.'); - } - return $x->round(@r); + require $EMU_LIB; + __emu_broot($self,$x,$y,@r); } sub exponent @@ -2488,13 +2082,12 @@ sub exponent if ($x->{sign} !~ /^[+-]$/) { - my $s = $x->{sign}; $s =~ s/^[+-]//; - return $self->new($s); # -inf,+inf => inf + my $s = $x->{sign}; $s =~ s/^[+-]//; # NaN, -inf,+inf => NaN or inf + return $self->new($s); } - my $e = $class->bzero(); - return $e->binc() if $x->is_zero(); - $e += $x->_trailing_zeros(); - $e; + return $self->bone() if $x->is_zero(); + + $self->new($x->_trailing_zeros()); } sub mantissa @@ -2504,10 +2097,11 @@ sub mantissa if ($x->{sign} !~ /^[+-]$/) { - return $self->new($x->{sign}); # keep + or - sign + # for NaN, +inf, -inf: keep the sign + return $self->new($x->{sign}); } - my $m = $x->copy(); - # that's inefficient + my $m = $x->copy(); delete $m->{_p}; delete $m->{_a}; + # that's a bit inefficient: my $zeros = $m->_trailing_zeros(); $m->brsft($zeros,10) if $zeros != 0; $m; @@ -2529,18 +2123,14 @@ sub bfround # precision: round to the $Nth digit left (+$n) or right (-$n) from the '.' # $n == 0 || $n == 1 => round to integer my $x = shift; $x = $class->new($x) unless ref $x; + my ($scale,$mode) = $x->_scale_p($x->precision(),$x->round_mode(),@_); - return $x if !defined $scale; # no-op - return $x if $x->modify('bfround'); + + return $x if !defined $scale || $x->modify('bfround'); # no-op # no-op for BigInts if $n <= 0 - if ($scale <= 0) - { - $x->{_a} = undef; # clear an eventual set A - $x->{_p} = $scale; return $x; - } + $x->bround( $x->length()-$scale, $mode) if $scale > 0; - $x->bround( $x->length()-$scale, $mode); $x->{_a} = undef; # bround sets {_a} $x->{_p} = $scale; # so correct it $x; @@ -2549,9 +2139,7 @@ sub bfround sub _scan_for_nonzero { # internal, used by bround() - my $x = shift; - my $pad = shift; - my $xs = shift; + my ($x,$pad,$xs) = @_; my $len = $x->length(); return 0 if $len == 1; # '5' is trailed by invisible zeros @@ -2708,35 +2296,15 @@ sub as_hex return $x->bstr() if $x->{sign} !~ /^[+-]$/; # inf, nan etc - my $es = ''; my $s = ''; + my $s = ''; $s = $x->{sign} if $x->{sign} eq '-'; if ($CAN{as_hex}) { - $es = ${$CALC->_as_hex($x->{value})}; + return $s . ${$CALC->_as_hex($x->{value})}; } - else - { - return '0x0' if $x->is_zero(); - 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($h,pack('v',$xr->numify())); - } - $es = reverse $es; - $es =~ s/^[0]+//; # strip leading zeros - $s .= '0x'; - } - $s . $es; + require $EMU_LIB; + __emu_as_hex(ref($x),$x,$s); } sub as_bin @@ -2746,34 +2314,15 @@ sub as_bin return $x->bstr() if $x->{sign} !~ /^[+-]$/; # inf, nan etc - my $es = ''; my $s = ''; - $s = $x->{sign} if $x->{sign} eq '-'; + my $s = ''; $s = $x->{sign} if $x->{sign} eq '-'; if ($CAN{as_bin}) { - $es = ${$CALC->_as_bin($x->{value})}; - } - else - { - return '0b0' if $x->is_zero(); - 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($b,pack('v',$xr->numify())); - } - $es = reverse $es; - $es =~ s/^[0]+//; # strip leading zeros - $s .= '0b'; + return $s . ${$CALC->_as_bin($x->{value})}; } - $s . $es; + + require $EMU_LIB; + __emu_as_bin(ref($x),$x,$s); + } ############################################################################## @@ -2977,17 +2526,17 @@ sub __from_hex my $sign = '+'; $sign = '-' if ($$hs =~ /^-/); $$hs =~ s/^[+-]//; # strip sign - if ($CAN{'_from_hex'}) + if ($CAN{'from_hex'}) { $x->{value} = $CALC->_from_hex($hs); } else { # fallback to pure perl - my $mul = Math::BigInt->bzero(); $mul++; + my $mul = Math::BigInt->bone(); my $x65536 = Math::BigInt->new(65536); - my $len = CORE::length($$hs)-2; - $len = int($len/4); # 4-digit parts, w/o '0x' + my $len = CORE::length($$hs)-2; # minus 2 for 0x + $len = int($len/4); # 4-digit parts, w/o '0x' my $val; my $i = -4; while ($len >= 0) { @@ -3016,15 +2565,15 @@ sub __from_bin my $sign = '+'; $sign = '-' if ($$bs =~ /^\-/); $$bs =~ s/^[+-]//; # strip sign - if ($CAN{'_from_bin'}) + if ($CAN{'from_bin'}) { $x->{value} = $CALC->_from_bin($bs); } else { - my $mul = Math::BigInt->bzero(); $mul++; + my $mul = Math::BigInt->bone(); my $x256 = Math::BigInt->new(256); - my $len = CORE::length($$bs)-2; + my $len = CORE::length($$bs)-2; # minus 2 for 0b $len = int($len/8); # 8-digit parts, w/o '0b' my $val; my $i = -8; while ($len >= 0) @@ -3142,7 +2691,7 @@ sub __gcd ############################################################################### # this method return 0 if the object can be modified, or 1 for not -# We use a fast use constant statement here, to avoid costly calls. Subclasses +# We use a fast constant sub() here, to avoid costly calls. Subclasses # may override it with special code (f.i. Math::BigInt::Constant does so) sub modify () { 0; } @@ -3182,8 +2731,8 @@ Math::BigInt - Arbitrary size integer math package $x->is_one('-'); # if $x is -1 $x->is_odd(); # if $x is odd $x->is_even(); # if $x is even - $x->is_positive(); # if $x >= 0 - $x->is_negative(); # if $x < 0 + $x->is_pos(); # if $x >= 0 + $x->is_neg(); # if $x < 0 $x->is_inf(sign); # if $x is +inf, or -inf (sign is default '+') $x->is_int(); # if $x is an integer (not a float) @@ -3260,14 +2809,15 @@ Math::BigInt - Arbitrary size integer math package $x->mantissa(); # return (signed) mantissa as BigInt $x->parts(); # return (mantissa,exponent) as BigInt $x->copy(); # make a true copy of $x (unlike $y = $x;) - $x->as_number(); # return as BigInt (in BigInt: same as copy()) + $x->as_int(); # return as BigInt (in BigInt: same as copy()) + $x->numify(); # return as scalar (might overflow!) # conversation to string (do not modify their argument) $x->bstr(); # normalized string $x->bsstr(); # normalized string in scientific notation $x->as_hex(); # as signed hexadecimal string with prefixed 0x $x->as_bin(); # as signed binary string with prefixed 0b - + # precision and accuracy (see section about rounding for more) $x->precision(); # return P of $x (or global, if P of $x undef) @@ -3546,10 +3096,10 @@ like: if ($x == 0) -=head2 is_positive()/is_negative() +=head2 is_pos()/is_neg() - $x->is_positive(); # true if >= 0 - $x->is_negative(); # true if < 0 + $x->is_pos(); # true if >= 0 + $x->is_neg(); # true if < 0 The methods return true if the argument is positive or negative, respectively. C is neither positive nor negative, while C<+inf> counts as positive, and @@ -3557,6 +3107,11 @@ C<-inf> is negative. A C is positive. These methods are only testing the sign, and not the value. +C and C are aliase to C and +C, respectively. C and C were +introduced in v1.36, while C and C were only introduced +in v1.68. + =head2 is_odd()/is_even()/is_int() $x->is_odd(); # true if odd, false for even @@ -3785,13 +3340,21 @@ Return the signed mantissa of $x as BigInt. $x->copy(); # make a true copy of $x (unlike $y = $x;) -=head2 as_number +=head2 as_int - $x->as_number(); # return as BigInt (in BigInt: same as copy()) + $x->as_int(); + +Returns $x as a BigInt (truncated towards zero). In BigInt this is the same as +C. + +C is an alias to this method. C was introduced in +v1.22, while C was only introduced in v1.68. -=head2 bsrt +=head2 bstr + + $x->bstr(); - $x->bstr(); # return normalized string +Returns a normalized string represantation of C<$x>. =head2 bsstr diff --git a/lib/Math/BigInt/Calc.pm b/lib/Math/BigInt/Calc.pm index 02770e2..e1cae77 100644 --- a/lib/Math/BigInt/Calc.pm +++ b/lib/Math/BigInt/Calc.pm @@ -4,11 +4,9 @@ use 5.005; use strict; # use warnings; # dont use warnings for older Perls -require Exporter; -use vars qw/@ISA $VERSION/; -@ISA = qw(Exporter); +use vars qw/$VERSION/; -$VERSION = '0.37'; +$VERSION = '0.38'; # Package to store unsigned big integers in decimal and do math with them @@ -194,6 +192,10 @@ sub _new # 1ex format. Assumes normalized value as input. my $d = $_[1]; my $il = length($$d)-1; + + # < BASE_LEN due len-1 above + return [ int($$d) ] if $il < $BASE_LEN; # shortcut for short numbers + # this leaves '00000' instead of int 0 and will be corrected after any op [ reverse(unpack("a" . ($il % $BASE_LEN+1) . ("a$BASE_LEN" x ($il / $BASE_LEN)), $$d)) ]; @@ -1240,6 +1242,22 @@ sub _pow # ref to array, ref to array, return ref to array my ($c,$cx,$cy) = @_; + if (scalar @$cy == 1 && $cy->[0] == 0) + { + splice (@$cx,1); $cx->[0] = 1; # y == 0 => x => 1 + return $cx; + } + if ((scalar @$cx == 1 && $cx->[0] == 1) || # x == 1 + (scalar @$cy == 1 && $cy->[0] == 1)) # or y == 1 + { + return $cx; + } + if (scalar @$cx == 1 && $cx->[0] == 0) + { + splice (@$cx,1); $cx->[0] = 0; # 0 ** y => 0 (if not y <= 0) + return $cx; + } + my $pow2 = _one(); my $y_bin = ${_as_bin($c,$cy)}; $y_bin =~ s/^0b//; @@ -1346,8 +1364,7 @@ sub _log_int return if (scalar @$x == 1 && $x->[0] == 0); # BASE 0 or 1 => NaN return if (scalar @$base == 1 && $base->[0] < 2); - my $cmp = _acmp($c,$x,$base); - # X == BASE => 1 + my $cmp = _acmp($c,$x,$base); # X == BASE => 1 if ($cmp == 0) { splice (@$x,1); $x->[0] = 1; @@ -1366,11 +1383,43 @@ sub _log_int my $x_org = _copy($c,$x); # preserve x splice(@$x,1); $x->[0] = 1; # keep ref to $x + my $trial = _copy($c,$base); + + # XXX TODO this only works if $base has only one element + if (scalar @$base == 1) + { + # compute int ( length_in_base_10(X) / ( log(base) / log(10) ) ) + my $len = _len($c,$x_org); + my $res = int($len / (log($base->[0]) / log(10))) || 1; # avoid $res == 0 + + $x->[0] = $res; + $trial = _pow ($c, _copy($c, $base), $x); + my $a = _acmp($x,$trial,$x_org); + return ($x,1) if $a == 0; + # we now that $res is too small + if ($res < 0) + { + _mul($c,$trial,$base); _add($c, $x, [1]); + } + else + { + # or too big + _div($c,$trial,$base); _sub($c, $x, [1]); + } + # did we now get the right result? + $a = _acmp($x,$trial,$x_org); + return ($x,1) if $a == 0; # yes, exactly + # still too big + if ($a > 0) + { + _div($c,$trial,$base); _sub($c, $x, [1]); + } + } + # simple loop that increments $x by two in each step, possible overstepping # the real result by one - # use a loop that keeps $x as scalar as long as possible (this is faster) - my $trial = _copy($c,$base); my $a; + my $a; my $base_mul = _mul($c, _copy($c,$base), $base); while (($a = _acmp($x,$trial,$x_org)) < 0) @@ -1981,6 +2030,7 @@ slow) fallback routines to emulate these: _root(obj) return the n'th (n >= 3) root of obj (truncated to int) _fac(obj) return factorial of object 1 (1*2*3*4..) _pow(obj,obj) return object 1 to the power of object 2 + return undef for NaN _gcd(obj,obj) return Greatest Common Divisor of two objects _zeros(obj) return number of trailing decimal zeros diff --git a/lib/Math/BigInt/CalcEmu.pm b/lib/Math/BigInt/CalcEmu.pm new file mode 100644 index 0000000..4ec244e --- /dev/null +++ b/lib/Math/BigInt/CalcEmu.pm @@ -0,0 +1,594 @@ +package Math::BigInt; + +use 5.005; +use strict; +# use warnings; # dont use warnings for older Perls + +use vars qw/$VERSION/; + +$VERSION = '0.01'; + +# See SYNOPSIS below. + +my $CALC_EMU; + +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) = @_; + + return $x->bzero(@r) if $y->is_zero() || $x->is_zero(); + + my $sign = 0; # sign of result + $sign = 1 if $sx == -1 && $sy == -1; + + my ($bx,$by); + + if ($sx == -1) # if x is negative + { + # two's complement: inc and flip all "bits" in $bx + $bx = $x->binc()->as_hex(); # -1 => 0, -2 => 1, -3 => 2 etc + $bx =~ s/-?0x//; + $bx =~ tr/0123456789abcdef/\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00/; + } + else + { + $bx = $x->as_hex(); # get binary representation + $bx =~ s/-?0x//; + $bx =~ tr/fedcba9876543210/\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00/; + } + if ($sy == -1) # if y is negative + { + # two's complement: inc and flip all "bits" in $by + $by = $y->copy()->binc()->as_hex(); # -1 => 0, -2 => 1, -3 => 2 etc + $by =~ s/-?0x//; + $by =~ tr/0123456789abcdef/\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00/; + } + else + { + $by = $y->as_hex(); # get binary representation + $by =~ s/-?0x//; + $by =~ tr/fedcba9876543210/\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00/; + } + # now we have bit-strings from X and Y, reverse them for padding + $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) + my $diff = CORE::length($bx) - CORE::length($by); + if ($diff > 0) + { + $bx = substr($bx,0,CORE::length($by)); + } + elsif ($diff < 0) + { + $by = substr($by,0,CORE::length($bx)); + } + + # 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 + # 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/0123456789abcdef/; + } + else + { + $bx =~ tr/\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00/fedcba9876543210/; + } + + $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}; + } + + # 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; + + $x->round(@r); + } + +sub __emu_bior + { + my ($self,$x,$y,$sx,$sy,@r) = @_; + + return $x->round(@r) if $y->is_zero(); + + my $sign = 0; # sign of result + $sign = 1 if ($sx == -1) || ($sy == -1); + + my ($bx,$by); + + if ($sx == -1) # if x is negative + { + # two's complement: inc and flip all "bits" in $bx + $bx = $x->binc()->as_hex(); # -1 => 0, -2 => 1, -3 => 2 etc + $bx =~ s/-?0x//; + $bx =~ tr/0123456789abcdef/\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00/; + } + else + { + $bx = $x->as_hex(); # get binary representation + $bx =~ s/-?0x//; + $bx =~ tr/fedcba9876543210/\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00/; + } + if ($sy == -1) # if y is negative + { + # two's complement: inc and flip all "bits" in $by + $by = $y->copy()->binc()->as_hex(); # -1 => 0, -2 => 1, -3 => 2 etc + $by =~ s/-?0x//; + $by =~ tr/0123456789abcdef/\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00/; + } + else + { + $by = $y->as_hex(); # get binary representation + $by =~ s/-?0x//; + $by =~ tr/fedcba9876543210/\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00/; + } + # now we have bit-strings from X and Y, reverse them for padding + $bx = reverse $bx; + $by = reverse $by; + + # 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) + { + $by .= $yy x $diff; + } + elsif ($diff < 0) + { + $bx .= $xx x abs($diff); + } + + # or 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 + # 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/0123456789abcdef/; + } + else + { + $bx =~ tr/\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00/fedcba9876543210/; + } + + $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}; + } + + # if one of X or Y was negative, we need to decrement result + $x->bdec() if $sign == 1; + + $x->round(@r); + } + +sub __emu_bxor + { + my ($self,$x,$y,$sx,$sy,@r) = @_; + + return $x->round(@r) if $y->is_zero(); + + my $sign = 0; # sign of result + $sign = 1 if $x->{sign} ne $y->{sign}; + + my ($bx,$by); + + if ($sx == -1) # if x is negative + { + # two's complement: inc and flip all "bits" in $bx + $bx = $x->binc()->as_hex(); # -1 => 0, -2 => 1, -3 => 2 etc + $bx =~ s/-?0x//; + $bx =~ tr/0123456789abcdef/\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00/; + } + else + { + $bx = $x->as_hex(); # get binary representation + $bx =~ s/-?0x//; + $bx =~ tr/fedcba9876543210/\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00/; + } + if ($sy == -1) # if y is negative + { + # two's complement: inc and flip all "bits" in $by + $by = $y->copy()->binc()->as_hex(); # -1 => 0, -2 => 1, -3 => 2 etc + $by =~ s/-?0x//; + $by =~ tr/0123456789abcdef/\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00/; + } + else + { + $by = $y->as_hex(); # get binary representation + $by =~ s/-?0x//; + $by =~ tr/fedcba9876543210/\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00/; + } + # now we have bit-strings from X and Y, reverse them for padding + $bx = reverse $bx; + $by = reverse $by; + + # 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) + { + $by .= $yy x $diff; + } + elsif ($diff < 0) + { + $bx .= $xx x abs($diff); + } + + # xor 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 + # 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/0123456789abcdef/; + } + else + { + $bx =~ tr/\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00/fedcba9876543210/; + } + + $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}; + } + + # calculate sign of result + $x->{sign} = '+'; + $x->{sign} = '-' if $sx != $sy && !$x->is_zero(); + + $x->bdec() if $sign == 1; + + $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 + { + # Should compute a guess of the result (by rule of thumb), then improve it + # via Newton's method or something similiar. + # XXX TODO + warn ('broot() not fully implemented in BigInt.'); + } + $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; + } + +############################################################################## +############################################################################## + +1; +__END__ + +=head1 NAME + +Math::BigInt::CalcEmu - Emulate low-level math with BigInt code + +=head1 SYNOPSIS + +Contains routines that emulate low-level math functions in BigInt, e.g. +optional routines the low-level math package does not provide on it's own. + +Will be loaded on demand and automatically by BigInt. + +Stuff here is really low-priority to optimize, +since it is far better to implement the operation in the low-level math +libary directly, possible even using a call to the native lib. + +=head1 DESCRIPTION + +=head1 METHODS + +=head1 LICENSE + +This program is free software; you may redistribute it and/or modify it under +the same terms as Perl itself. + +=head1 AUTHORS + +(c) Tels http://bloodgate.com 2003 - based on BigInt code by +Tels from 2001-2003. + +=head1 SEE ALSO + +L, L, L, +L and L. + +=cut diff --git a/lib/Math/BigInt/t/alias.inc b/lib/Math/BigInt/t/alias.inc new file mode 100644 index 0000000..84310fc --- /dev/null +++ b/lib/Math/BigInt/t/alias.inc @@ -0,0 +1,12 @@ + +# alias subroutine testing, included by sub_ali.t and mbi_ali.t + +my $x = $CL->new(123); + +is ($x->is_pos(), 1, '123 is positive'); +is ($x->is_neg(), 0, '123 is not negative'); +is ($x->as_int(), 123, '123 is 123 as int'); +is (ref($x->as_int()), $CL, '123 is scalar as int'); +$x->bneg(); +is ($x->is_pos(), 0, '-123 is not positive'); +is ($x->is_neg(), 1, '-123 is negative'); diff --git a/lib/Math/BigInt/t/bigfltpm.inc b/lib/Math/BigInt/t/bigfltpm.inc index 60a8f08..1a05a66 100644 --- a/lib/Math/BigInt/t/bigfltpm.inc +++ b/lib/Math/BigInt/t/bigfltpm.inc @@ -1201,7 +1201,7 @@ abc:1:abc:NaN &ffac Nanfac:NaN -1:NaN -+inf:NaN ++inf:inf -inf:NaN 0:1 1:1 diff --git a/lib/Math/BigInt/t/bigintc.t b/lib/Math/BigInt/t/bigintc.t index 8d352eb..1f0804c 100644 --- a/lib/Math/BigInt/t/bigintc.t +++ b/lib/Math/BigInt/t/bigintc.t @@ -14,7 +14,7 @@ use Math::BigInt::Calc; BEGIN { - plan tests => 296; + plan tests => 300; } my ($BASE_LEN, $AND_BITS, $XOR_BITS, $OR_BITS, $BASE_LEN_SMALL, $MAX_VAL) = @@ -224,6 +224,15 @@ $x = $C->_new(\"81"); $n = $C->_new(\"4"); # 3*3*3*3 == 81 ok (${$C->_str($C->_root($x,$n))},'3'); # _pow (and _root) +$x = $C->_new(\"0"); $n = $C->_new(\"3"); # 0 ** y => 0 +ok (${$C->_str($C->_pow($x,$n))}, 0); +$x = $C->_new(\"3"); $n = $C->_new(\"0"); # x ** 0 => 1 +ok (${$C->_str($C->_pow($x,$n))}, 1); +$x = $C->_new(\"1"); $n = $C->_new(\"3"); # 1 ** y => 1 +ok (${$C->_str($C->_pow($x,$n))}, 1); +$x = $C->_new(\"5"); $n = $C->_new(\"1"); # x ** 1 => x +ok (${$C->_str($C->_pow($x,$n))}, 5); + $x = $C->_new(\"81"); $n = $C->_new(\"3"); # 81 ** 3 == 531441 ok (${$C->_str($C->_pow($x,$n))},81 ** 3); diff --git a/lib/Math/BigInt/t/bigintpm.inc b/lib/Math/BigInt/t/bigintpm.inc index db52553..3cbb993 100644 --- a/lib/Math/BigInt/t/bigintpm.inc +++ b/lib/Math/BigInt/t/bigintpm.inc @@ -1891,7 +1891,7 @@ abc:NaN,NaN &bfac -1:NaN NaNfac:NaN -+inf:NaN ++inf:inf -inf:NaN 0:1 1:1 diff --git a/lib/Math/BigInt/t/mbf_ali.t b/lib/Math/BigInt/t/mbf_ali.t new file mode 100644 index 0000000..1ca4315 --- /dev/null +++ b/lib/Math/BigInt/t/mbf_ali.t @@ -0,0 +1,42 @@ +#!/usr/bin/perl -w + +# test that the new alias names work + +use Test::More; +use strict; + +BEGIN + { + $| = 1; + # to locate the testing files + my $location = $0; $location =~ s/mbf_ali.t//i; + if ($ENV{PERL_CORE}) + { + # testing with the core distribution + @INC = qw(../t/lib); + } + unshift @INC, qw(../lib); + if (-d 't') + { + chdir 't'; + require File::Spec; + unshift @INC, File::Spec->catdir(File::Spec->updir, $location); + } + else + { + unshift @INC, $location; + } + print "# INC = @INC\n"; + + plan tests => 6; + } + +use Math::BigFloat; + +use vars qw/$x $CL/; + +$CL = 'Math::BigFloat'; + +require 'alias.inc'; + + diff --git a/lib/Math/BigInt/t/mbi_ali.t b/lib/Math/BigInt/t/mbi_ali.t new file mode 100644 index 0000000..4028017 --- /dev/null +++ b/lib/Math/BigInt/t/mbi_ali.t @@ -0,0 +1,42 @@ +#!/usr/bin/perl -w + +# test that the new alias names work + +use Test::More; +use strict; + +BEGIN + { + $| = 1; + # to locate the testing files + my $location = $0; $location =~ s/mbi_ali.t//i; + if ($ENV{PERL_CORE}) + { + # testing with the core distribution + @INC = qw(../t/lib); + } + unshift @INC, qw(../lib); + if (-d 't') + { + chdir 't'; + require File::Spec; + unshift @INC, File::Spec->catdir(File::Spec->updir, $location); + } + else + { + unshift @INC, $location; + } + print "# INC = @INC\n"; + + plan tests => 6; + } + +use Math::BigInt; + +use vars qw/$x $CL/; + +$CL = 'Math::BigInt'; + +require 'alias.inc'; + + diff --git a/lib/Math/BigInt/t/sub_ali.t b/lib/Math/BigInt/t/sub_ali.t new file mode 100644 index 0000000..93620a9 --- /dev/null +++ b/lib/Math/BigInt/t/sub_ali.t @@ -0,0 +1,40 @@ +#!/usr/bin/perl -w + +# test that the new alias names work + +use Test::More; +use strict; + +BEGIN + { + $| = 1; + # to locate the testing files + my $location = $0; $location =~ s/sub_ali.t//i; + if ($ENV{PERL_CORE}) + { + # testing with the core distribution + @INC = qw(../t/lib); + } + unshift @INC, qw(../lib); + if (-d 't') + { + chdir 't'; + require File::Spec; + unshift @INC, File::Spec->catdir(File::Spec->updir, $location); + } + else + { + unshift @INC, $location; + } + print "# INC = @INC\n"; + + plan tests => 6; + } + +use Math::BigInt::Subclass; + +use vars qw/$CL $x/; +$CL = 'Math::BigInt::Subclass'; + +require 'alias.inc'; + diff --git a/lib/Math/BigInt/t/upgrade.inc b/lib/Math/BigInt/t/upgrade.inc index 0b66640..49dbf91 100644 --- a/lib/Math/BigInt/t/upgrade.inc +++ b/lib/Math/BigInt/t/upgrade.inc @@ -1245,7 +1245,7 @@ abc:NaN,NaN &bfac -1:NaN NaNfac:NaN -+inf:NaN ++inf:inf -inf:NaN 0:1 1:1