use vars qw/@ISA $VERSION/;
@ISA = qw(Exporter);
-$VERSION = '0.20';
+$VERSION = '0.22';
# Package to store unsigned big integers in decimal and do math with them
# This routine clobbers up array x, but not y.
my ($c,$x,$y) = @_;
+
+ return $x if (@$y == 1) && $y->[0] == 0; # $x + 0 => $x
+ if ((@$x == 1) && $x->[0] == 0) # 0 + $y => $y->copy
+ {
+ # twice as slow as $x = [ @$y ], but necc. to retain $x as ref :(
+ @$x = @$y; return $x;
+ }
# for each in Y, add Y to X and carry. If after that, something is left in
# X, foreach in X add carry to X and then return X, carry
# modifies first arg, second need not be different from first
my ($c,$xv,$yv) = @_;
- # shortcut for two very short numbers
- # +0 since part maybe string '00001' from new()
+ # shortcut for two very short numbers (improved by Nathan Zook)
# works also if xv and yv are the same reference
- if ((@$xv == 1) && (@$yv == 1)
- && (length($xv->[0]+0) <= $BASE_LEN2)
- && (length($yv->[0]+0) <= $BASE_LEN2))
- {
- $xv->[0] *= $yv->[0];
- return $xv;
- }
-
+ if ((@$xv == 1) && (@$yv == 1))
+ {
+ if (($xv->[0] *= $yv->[0]) >= $MBASE)
+ {
+ $xv->[0] = $xv->[0] - ($xv->[1] = int($xv->[0] * $RBASE)) * $MBASE;
+ };
+ return $xv;
+ }
+ # shortcut for result == 0
+ if ( ((@$xv == 1) && ($xv->[0] == 0)) ||
+ ((@$yv == 1) && ($yv->[0] == 0)) )
+ {
+ @$xv = (0);
+ return $xv;
+ }
+
# since multiplying $x with $x fails, make copy in this case
$yv = [@$xv] if "$xv" eq "$yv"; # same references?
if ($LEN_CONVERT != 0)
# modifies first arg, second need not be different from first
my ($c,$xv,$yv) = @_;
- # shortcut for two very short numbers
- # +0 since part maybe string '00001' from new()
+ # shortcut for two very short numbers (improved by Nathan Zook)
# works also if xv and yv are the same reference
- if ((@$xv == 1) && (@$yv == 1)
- && (length($xv->[0]+0) <= $BASE_LEN2)
- && (length($yv->[0]+0) <= $BASE_LEN2))
- {
- $xv->[0] *= $yv->[0];
- return $xv;
- }
+ if ((@$xv == 1) && (@$yv == 1))
+ {
+ if (($xv->[0] *= $yv->[0]) >= $MBASE)
+ {
+ $xv->[0] =
+ $xv->[0] - ($xv->[1] = int($xv->[0] / $MBASE)) * $MBASE;
+ };
+ return $xv;
+ }
+ # shortcut for result == 0
+ if ( ((@$xv == 1) && ($xv->[0] == 0)) ||
+ ((@$yv == 1) && ($yv->[0] == 0)) )
+ {
+ @$xv = (0);
+ return $xv;
+ }
+
# since multiplying $x with $x fails, make copy in this case
$yv = [@$xv] if "$xv" eq "$yv"; # same references?
$cx;
}
-sub _sqrt1
+sub _fac
+ {
+ # factorial of $x
+ # ref to array, return ref to array
+ my ($c,$cx) = @_;
+
+ if ((@$cx == 1) && ($cx->[0] <= 2))
+ {
+ $cx->[0] = 1 * ($cx->[0]||1); # 0,1 => 1, 2 => 2
+ return $cx;
+ }
+
+ # go forward until $base is exceeded
+ # limit is either $x or $base (x == 100 means as result too high)
+ my $steps = 100; $steps = $cx->[0] if @$cx == 1;
+ my $r = 2; my $cf = 3; my $step = 1; my $last = $r;
+ while ($r < $BASE && $step < $steps)
+ {
+ $last = $r; $r *= $cf++; $step++;
+ }
+ if ((@$cx == 1) && ($step == $cx->[0]))
+ {
+ # completely done
+ $cx = [$last];
+ return $cx;
+ }
+ my $n = _copy($c,$cx);
+ $cx = [$last];
+
+ #$cx = _one();
+ while (!(@$n == 1 && $n->[0] == $step))
+ {
+ _mul($c,$cx,$n); _dec($c,$n);
+ }
+ $cx;
+ }
+
+use constant DEBUG => 0;
+
+my $steps = 0;
+
+sub steps { $steps };
+
+sub _sqrt
{
# square-root of $x
# ref to array, return ref to array
return $x;
}
my $y = _copy($c,$x);
- my $l = _len($c,$x) / 2; # hopefully _len/2 is < $BASE
- # my $l2 = [ _len($c,$x) / 2 ]; # old way: hopefully _len/2 is < $BASE
-
- splice @$x,0; $x->[0] = 1; # keep ref($x), but modify it
-
- # old way
- # _lsft($c,$x,$l2,10);
+ # hopefully _len/2 is < $BASE, the -1 is to always undershot the guess
+ # since our guess will "grow"
+ my $l = int((_len($c,$x)-1) / 2);
+
+ my $lastelem = $x->[-1]; # for guess
+ my $elems = scalar @$x - 1;
+ # not enough digits, but could have more?
+ if ((length($lastelem) <= 3) && ($elems > 1))
+ {
+ # right-align with zero pad
+ my $len = length($lastelem) & 1;
+ print "$lastelem => " if DEBUG;
+ $lastelem .= substr($x->[-2] . '0' x $BASE_LEN,0,$BASE_LEN);
+ # former odd => make odd again, or former even to even again
+ $lastelem = $lastelem / 10 if (length($lastelem) & 1) != $len;
+ print "$lastelem\n" if DEBUG;
+ }
# construct $x (instead of _lsft($c,$x,$l,10)
my $r = $l % $BASE_LEN; # 10000 00000 00000 00000 ($BASE_LEN=5)
$l = int($l / $BASE_LEN);
- $x->[$l--] = int('1' . '0' x $r);
- $x->[$l--] = 0 while ($l >= 0);
+ print "l = $l " if DEBUG;
+
+ splice @$x,$l; # keep ref($x), but modify it
+
+ # we make the first part of the guess not '1000...0' but int(sqrt($lastelem))
+ # that gives us:
+ # 14400 00000 => sqrt(14400) => 120
+ # 144000 000000 => sqrt(144000) => 379
+
+ # $x->[$l--] = int('1' . '0' x $r); # old way of guessing
+ print "$lastelem (elems $elems) => " if DEBUG;
+ $lastelem = $lastelem / 10 if ($elems & 1 == 1); # odd or even?
+ my $g = sqrt($lastelem); $g =~ s/\.//; # 2.345 => 2345
+ $r -= 1 if $elems & 1 == 0; # 70 => 7
+
+ # padd with zeros if result is too short
+ $x->[$l--] = int(substr($g . '0' x $r,0,$r+1));
+ print "now ",$x->[-1] if DEBUG;
+ print " would have been ", int('1' . '0' x $r),"\n" if DEBUG;
+
+ # If @$x > 1, we could compute the second elem of the guess, too, to create
+ # an even better guess. Not implemented yet.
+ $x->[$l--] = 0 while ($l >= 0); # all other digits of guess are zero
+ print "start x= ",${_str($c,$x)},"\n" if DEBUG;
my $two = _two();
my $last = _zero();
my $lastlast = _zero();
+ $steps = 0 if DEBUG;
while (_acmp($c,$last,$x) != 0 && _acmp($c,$lastlast,$x) != 0)
{
+ $steps++ if DEBUG;
$lastlast = _copy($c,$last);
$last = _copy($c,$x);
_add($c,$x, _div($c,_copy($c,$y),$x));
_div($c,$x, $two );
+ print " x= ",${_str($c,$x)},"\n" if DEBUG;
}
+ print "\nsteps in sqrt: $steps, " if DEBUG;
_dec($c,$x) if _acmp($c,$y,_mul($c,_copy($c,$x),$x)) < 0; # overshot?
+ print " final ",$x->[-1],"\n" if DEBUG;
$x;
}
_mod(obj,obj) Return remainder of div of the 1st by the 2nd object
_sqrt(obj) return the square root of object (truncate 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
_gcd(obj,obj) return Greatest Common Divisor of two objects
projects / p5sagit/p5-mst-13.2.git / blobdiff