| 1 | package Math::BigInt::Calc;
|
|---|
| 2 |
|
|---|
| 3 | use 5.005;
|
|---|
| 4 | use strict;
|
|---|
| 5 | # use warnings; # dont use warnings for older Perls
|
|---|
| 6 |
|
|---|
| 7 | use vars qw/$VERSION/;
|
|---|
| 8 |
|
|---|
| 9 | $VERSION = '0.47';
|
|---|
| 10 |
|
|---|
| 11 | # Package to store unsigned big integers in decimal and do math with them
|
|---|
| 12 |
|
|---|
| 13 | # Internally the numbers are stored in an array with at least 1 element, no
|
|---|
| 14 | # leading zero parts (except the first) and in base 1eX where X is determined
|
|---|
| 15 | # automatically at loading time to be the maximum possible value
|
|---|
| 16 |
|
|---|
| 17 | # todo:
|
|---|
| 18 | # - fully remove funky $# stuff in div() (maybe - that code scares me...)
|
|---|
| 19 |
|
|---|
| 20 | # USE_MUL: due to problems on certain os (os390, posix-bc) "* 1e-5" is used
|
|---|
| 21 | # instead of "/ 1e5" at some places, (marked with USE_MUL). Other platforms
|
|---|
| 22 | # BS2000, some Crays need USE_DIV instead.
|
|---|
| 23 | # The BEGIN block is used to determine which of the two variants gives the
|
|---|
| 24 | # correct result.
|
|---|
| 25 |
|
|---|
| 26 | # Beware of things like:
|
|---|
| 27 | # $i = $i * $y + $car; $car = int($i / $MBASE); $i = $i % $MBASE;
|
|---|
| 28 | # This works on x86, but fails on ARM (SA1100, iPAQ) due to whoknows what
|
|---|
| 29 | # reasons. So, use this instead (slower, but correct):
|
|---|
| 30 | # $i = $i * $y + $car; $car = int($i / $MBASE); $i -= $MBASE * $car;
|
|---|
| 31 |
|
|---|
| 32 | ##############################################################################
|
|---|
| 33 | # global constants, flags and accessory
|
|---|
| 34 |
|
|---|
| 35 | # announce that we are compatible with MBI v1.70 and up
|
|---|
| 36 | sub api_version () { 1; }
|
|---|
| 37 |
|
|---|
| 38 | # constants for easier life
|
|---|
| 39 | my ($BASE,$BASE_LEN,$MBASE,$RBASE,$MAX_VAL,$BASE_LEN_SMALL);
|
|---|
| 40 | my ($AND_BITS,$XOR_BITS,$OR_BITS);
|
|---|
| 41 | my ($AND_MASK,$XOR_MASK,$OR_MASK);
|
|---|
| 42 |
|
|---|
| 43 | sub _base_len
|
|---|
| 44 | {
|
|---|
| 45 | # set/get the BASE_LEN and assorted other, connected values
|
|---|
| 46 | # used only be the testsuite, set is used only by the BEGIN block below
|
|---|
| 47 | shift;
|
|---|
| 48 |
|
|---|
| 49 | my $b = shift;
|
|---|
| 50 | if (defined $b)
|
|---|
| 51 | {
|
|---|
| 52 | # find whether we can use mul or div or none in mul()/div()
|
|---|
| 53 | # (in last case reduce BASE_LEN_SMALL)
|
|---|
| 54 | $BASE_LEN_SMALL = $b+1;
|
|---|
| 55 | my $caught = 0;
|
|---|
| 56 | while (--$BASE_LEN_SMALL > 5)
|
|---|
| 57 | {
|
|---|
| 58 | $MBASE = int("1e".$BASE_LEN_SMALL);
|
|---|
| 59 | $RBASE = abs('1e-'.$BASE_LEN_SMALL); # see USE_MUL
|
|---|
| 60 | $caught = 0;
|
|---|
| 61 | $caught += 1 if (int($MBASE * $RBASE) != 1); # should be 1
|
|---|
| 62 | $caught += 2 if (int($MBASE / $MBASE) != 1); # should be 1
|
|---|
| 63 | last if $caught != 3;
|
|---|
| 64 | }
|
|---|
| 65 | # BASE_LEN is used for anything else than mul()/div()
|
|---|
| 66 | $BASE_LEN = $BASE_LEN_SMALL;
|
|---|
| 67 | $BASE_LEN = shift if (defined $_[0]); # one more arg?
|
|---|
| 68 | $BASE = int("1e".$BASE_LEN);
|
|---|
| 69 |
|
|---|
| 70 | $MBASE = int("1e".$BASE_LEN_SMALL);
|
|---|
| 71 | $RBASE = abs('1e-'.$BASE_LEN_SMALL); # see USE_MUL
|
|---|
| 72 | $MAX_VAL = $MBASE-1;
|
|---|
| 73 |
|
|---|
| 74 | # avoid redefinitions
|
|---|
| 75 |
|
|---|
| 76 | undef &_mul;
|
|---|
| 77 | undef &_div;
|
|---|
| 78 |
|
|---|
| 79 | # $caught & 1 != 0 => cannot use MUL
|
|---|
| 80 | # $caught & 2 != 0 => cannot use DIV
|
|---|
| 81 | # The parens around ($caught & 1) were important, indeed, if we would use
|
|---|
| 82 | # & here.
|
|---|
| 83 | if ($caught == 2) # 2
|
|---|
| 84 | {
|
|---|
| 85 | # must USE_MUL since we cannot use DIV
|
|---|
| 86 | *{_mul} = \&_mul_use_mul;
|
|---|
| 87 | *{_div} = \&_div_use_mul;
|
|---|
| 88 | }
|
|---|
| 89 | else # 0 or 1
|
|---|
| 90 | {
|
|---|
| 91 | # can USE_DIV instead
|
|---|
| 92 | *{_mul} = \&_mul_use_div;
|
|---|
| 93 | *{_div} = \&_div_use_div;
|
|---|
| 94 | }
|
|---|
| 95 | }
|
|---|
| 96 | return $BASE_LEN unless wantarray;
|
|---|
| 97 | return ($BASE_LEN, $AND_BITS, $XOR_BITS, $OR_BITS, $BASE_LEN_SMALL, $MAX_VAL, $BASE);
|
|---|
| 98 | }
|
|---|
| 99 |
|
|---|
| 100 | sub _new
|
|---|
| 101 | {
|
|---|
| 102 | # (ref to string) return ref to num_array
|
|---|
| 103 | # Convert a number from string format (without sign) to internal base
|
|---|
| 104 | # 1ex format. Assumes normalized value as input.
|
|---|
| 105 | my $il = length($_[1])-1;
|
|---|
| 106 |
|
|---|
| 107 | # < BASE_LEN due len-1 above
|
|---|
| 108 | return [ int($_[1]) ] if $il < $BASE_LEN; # shortcut for short numbers
|
|---|
| 109 |
|
|---|
| 110 | # this leaves '00000' instead of int 0 and will be corrected after any op
|
|---|
| 111 | [ reverse(unpack("a" . ($il % $BASE_LEN+1)
|
|---|
| 112 | . ("a$BASE_LEN" x ($il / $BASE_LEN)), $_[1])) ];
|
|---|
| 113 | }
|
|---|
| 114 |
|
|---|
| 115 | BEGIN
|
|---|
| 116 | {
|
|---|
| 117 | # from Daniel Pfeiffer: determine largest group of digits that is precisely
|
|---|
| 118 | # multipliable with itself plus carry
|
|---|
| 119 | # Test now changed to expect the proper pattern, not a result off by 1 or 2
|
|---|
| 120 | my ($e, $num) = 3; # lowest value we will use is 3+1-1 = 3
|
|---|
| 121 | do
|
|---|
| 122 | {
|
|---|
| 123 | $num = ('9' x ++$e) + 0;
|
|---|
| 124 | $num *= $num + 1.0;
|
|---|
| 125 | } while ("$num" =~ /9{$e}0{$e}/); # must be a certain pattern
|
|---|
| 126 | $e--; # last test failed, so retract one step
|
|---|
| 127 | # the limits below brush the problems with the test above under the rug:
|
|---|
| 128 | # the test should be able to find the proper $e automatically
|
|---|
| 129 | $e = 5 if $^O =~ /^uts/; # UTS get's some special treatment
|
|---|
| 130 | $e = 5 if $^O =~ /^unicos/; # unicos is also problematic (6 seems to work
|
|---|
| 131 | # there, but we play safe)
|
|---|
| 132 | $e = 5 if $] < 5.006; # cap, for older Perls
|
|---|
| 133 | $e = 7 if $e > 7; # cap, for VMS, OS/390 and other 64 bit systems
|
|---|
| 134 | # 8 fails inside random testsuite, so take 7
|
|---|
| 135 |
|
|---|
| 136 | __PACKAGE__->_base_len($e); # set and store
|
|---|
| 137 |
|
|---|
| 138 | use integer;
|
|---|
| 139 | # find out how many bits _and, _or and _xor can take (old default = 16)
|
|---|
| 140 | # I don't think anybody has yet 128 bit scalars, so let's play safe.
|
|---|
| 141 | local $^W = 0; # don't warn about 'nonportable number'
|
|---|
| 142 | $AND_BITS = 15; $XOR_BITS = 15; $OR_BITS = 15;
|
|---|
| 143 |
|
|---|
| 144 | # find max bits, we will not go higher than numberofbits that fit into $BASE
|
|---|
| 145 | # to make _and etc simpler (and faster for smaller, slower for large numbers)
|
|---|
| 146 | my $max = 16;
|
|---|
| 147 | while (2 ** $max < $BASE) { $max++; }
|
|---|
| 148 | {
|
|---|
| 149 | no integer;
|
|---|
| 150 | $max = 16 if $] < 5.006; # older Perls might not take >16 too well
|
|---|
| 151 | }
|
|---|
| 152 | my ($x,$y,$z);
|
|---|
| 153 | do {
|
|---|
| 154 | $AND_BITS++;
|
|---|
| 155 | $x = oct('0b' . '1' x $AND_BITS); $y = $x & $x;
|
|---|
| 156 | $z = (2 ** $AND_BITS) - 1;
|
|---|
| 157 | } while ($AND_BITS < $max && $x == $z && $y == $x);
|
|---|
| 158 | $AND_BITS --; # retreat one step
|
|---|
| 159 | do {
|
|---|
| 160 | $XOR_BITS++;
|
|---|
| 161 | $x = oct('0b' . '1' x $XOR_BITS); $y = $x ^ 0;
|
|---|
| 162 | $z = (2 ** $XOR_BITS) - 1;
|
|---|
| 163 | } while ($XOR_BITS < $max && $x == $z && $y == $x);
|
|---|
| 164 | $XOR_BITS --; # retreat one step
|
|---|
| 165 | do {
|
|---|
| 166 | $OR_BITS++;
|
|---|
| 167 | $x = oct('0b' . '1' x $OR_BITS); $y = $x | $x;
|
|---|
| 168 | $z = (2 ** $OR_BITS) - 1;
|
|---|
| 169 | } while ($OR_BITS < $max && $x == $z && $y == $x);
|
|---|
| 170 | $OR_BITS --; # retreat one step
|
|---|
| 171 |
|
|---|
| 172 | $AND_MASK = __PACKAGE__->_new( ( 2 ** $AND_BITS ));
|
|---|
| 173 | $XOR_MASK = __PACKAGE__->_new( ( 2 ** $XOR_BITS ));
|
|---|
| 174 | $OR_MASK = __PACKAGE__->_new( ( 2 ** $OR_BITS ));
|
|---|
| 175 | }
|
|---|
| 176 |
|
|---|
| 177 | ###############################################################################
|
|---|
| 178 |
|
|---|
| 179 | sub _zero
|
|---|
| 180 | {
|
|---|
| 181 | # create a zero
|
|---|
| 182 | [ 0 ];
|
|---|
| 183 | }
|
|---|
| 184 |
|
|---|
| 185 | sub _one
|
|---|
| 186 | {
|
|---|
| 187 | # create a one
|
|---|
| 188 | [ 1 ];
|
|---|
| 189 | }
|
|---|
| 190 |
|
|---|
| 191 | sub _two
|
|---|
| 192 | {
|
|---|
| 193 | # create a two (used internally for shifting)
|
|---|
| 194 | [ 2 ];
|
|---|
| 195 | }
|
|---|
| 196 |
|
|---|
| 197 | sub _ten
|
|---|
| 198 | {
|
|---|
| 199 | # create a 10 (used internally for shifting)
|
|---|
| 200 | [ 10 ];
|
|---|
| 201 | }
|
|---|
| 202 |
|
|---|
| 203 | sub _copy
|
|---|
| 204 | {
|
|---|
| 205 | # make a true copy
|
|---|
| 206 | [ @{$_[1]} ];
|
|---|
| 207 | }
|
|---|
| 208 |
|
|---|
| 209 | # catch and throw away
|
|---|
| 210 | sub import { }
|
|---|
| 211 |
|
|---|
| 212 | ##############################################################################
|
|---|
| 213 | # convert back to string and number
|
|---|
| 214 |
|
|---|
| 215 | sub _str
|
|---|
| 216 | {
|
|---|
| 217 | # (ref to BINT) return num_str
|
|---|
| 218 | # Convert number from internal base 100000 format to string format.
|
|---|
| 219 | # internal format is always normalized (no leading zeros, "-0" => "+0")
|
|---|
| 220 | my $ar = $_[1];
|
|---|
| 221 |
|
|---|
| 222 | my $l = scalar @$ar; # number of parts
|
|---|
| 223 | if ($l < 1) # should not happen
|
|---|
| 224 | {
|
|---|
| 225 | require Carp;
|
|---|
| 226 | Carp::croak("$_[1] has no elements");
|
|---|
| 227 | }
|
|---|
| 228 |
|
|---|
| 229 | my $ret = "";
|
|---|
| 230 | # handle first one different to strip leading zeros from it (there are no
|
|---|
| 231 | # leading zero parts in internal representation)
|
|---|
| 232 | $l --; $ret .= int($ar->[$l]); $l--;
|
|---|
| 233 | # Interestingly, the pre-padd method uses more time
|
|---|
| 234 | # the old grep variant takes longer (14 vs. 10 sec)
|
|---|
| 235 | my $z = '0' x ($BASE_LEN-1);
|
|---|
| 236 | while ($l >= 0)
|
|---|
| 237 | {
|
|---|
| 238 | $ret .= substr($z.$ar->[$l],-$BASE_LEN); # fastest way I could think of
|
|---|
| 239 | $l--;
|
|---|
| 240 | }
|
|---|
| 241 | $ret;
|
|---|
| 242 | }
|
|---|
| 243 |
|
|---|
| 244 | sub _num
|
|---|
| 245 | {
|
|---|
| 246 | # Make a number (scalar int/float) from a BigInt object
|
|---|
| 247 | my $x = $_[1];
|
|---|
| 248 |
|
|---|
| 249 | return 0+$x->[0] if scalar @$x == 1; # below $BASE
|
|---|
| 250 | my $fac = 1;
|
|---|
| 251 | my $num = 0;
|
|---|
| 252 | foreach (@$x)
|
|---|
| 253 | {
|
|---|
| 254 | $num += $fac*$_; $fac *= $BASE;
|
|---|
| 255 | }
|
|---|
| 256 | $num;
|
|---|
| 257 | }
|
|---|
| 258 |
|
|---|
| 259 | ##############################################################################
|
|---|
| 260 | # actual math code
|
|---|
| 261 |
|
|---|
| 262 | sub _add
|
|---|
| 263 | {
|
|---|
| 264 | # (ref to int_num_array, ref to int_num_array)
|
|---|
| 265 | # routine to add two base 1eX numbers
|
|---|
| 266 | # stolen from Knuth Vol 2 Algorithm A pg 231
|
|---|
| 267 | # there are separate routines to add and sub as per Knuth pg 233
|
|---|
| 268 | # This routine clobbers up array x, but not y.
|
|---|
| 269 |
|
|---|
| 270 | my ($c,$x,$y) = @_;
|
|---|
| 271 |
|
|---|
| 272 | return $x if (@$y == 1) && $y->[0] == 0; # $x + 0 => $x
|
|---|
| 273 | if ((@$x == 1) && $x->[0] == 0) # 0 + $y => $y->copy
|
|---|
| 274 | {
|
|---|
| 275 | # twice as slow as $x = [ @$y ], but necc. to retain $x as ref :(
|
|---|
| 276 | @$x = @$y; return $x;
|
|---|
| 277 | }
|
|---|
| 278 |
|
|---|
| 279 | # for each in Y, add Y to X and carry. If after that, something is left in
|
|---|
| 280 | # X, foreach in X add carry to X and then return X, carry
|
|---|
| 281 | # Trades one "$j++" for having to shift arrays
|
|---|
| 282 | my $i; my $car = 0; my $j = 0;
|
|---|
| 283 | for $i (@$y)
|
|---|
| 284 | {
|
|---|
| 285 | $x->[$j] -= $BASE if $car = (($x->[$j] += $i + $car) >= $BASE) ? 1 : 0;
|
|---|
| 286 | $j++;
|
|---|
| 287 | }
|
|---|
| 288 | while ($car != 0)
|
|---|
| 289 | {
|
|---|
| 290 | $x->[$j] -= $BASE if $car = (($x->[$j] += $car) >= $BASE) ? 1 : 0; $j++;
|
|---|
| 291 | }
|
|---|
| 292 | $x;
|
|---|
| 293 | }
|
|---|
| 294 |
|
|---|
| 295 | sub _inc
|
|---|
| 296 | {
|
|---|
| 297 | # (ref to int_num_array, ref to int_num_array)
|
|---|
| 298 | # Add 1 to $x, modify $x in place
|
|---|
| 299 | my ($c,$x) = @_;
|
|---|
| 300 |
|
|---|
| 301 | for my $i (@$x)
|
|---|
| 302 | {
|
|---|
| 303 | return $x if (($i += 1) < $BASE); # early out
|
|---|
| 304 | $i = 0; # overflow, next
|
|---|
| 305 | }
|
|---|
| 306 | push @$x,1 if (($x->[-1] || 0) == 0); # last overflowed, so extend
|
|---|
| 307 | $x;
|
|---|
| 308 | }
|
|---|
| 309 |
|
|---|
| 310 | sub _dec
|
|---|
| 311 | {
|
|---|
| 312 | # (ref to int_num_array, ref to int_num_array)
|
|---|
| 313 | # Sub 1 from $x, modify $x in place
|
|---|
| 314 | my ($c,$x) = @_;
|
|---|
| 315 |
|
|---|
| 316 | my $MAX = $BASE-1; # since MAX_VAL based on MBASE
|
|---|
| 317 | for my $i (@$x)
|
|---|
| 318 | {
|
|---|
| 319 | last if (($i -= 1) >= 0); # early out
|
|---|
| 320 | $i = $MAX; # underflow, next
|
|---|
| 321 | }
|
|---|
| 322 | pop @$x if $x->[-1] == 0 && @$x > 1; # last underflowed (but leave 0)
|
|---|
| 323 | $x;
|
|---|
| 324 | }
|
|---|
| 325 |
|
|---|
| 326 | sub _sub
|
|---|
| 327 | {
|
|---|
| 328 | # (ref to int_num_array, ref to int_num_array, swap)
|
|---|
| 329 | # subtract base 1eX numbers -- stolen from Knuth Vol 2 pg 232, $x > $y
|
|---|
| 330 | # subtract Y from X by modifying x in place
|
|---|
| 331 | my ($c,$sx,$sy,$s) = @_;
|
|---|
| 332 |
|
|---|
| 333 | my $car = 0; my $i; my $j = 0;
|
|---|
| 334 | if (!$s)
|
|---|
| 335 | {
|
|---|
| 336 | for $i (@$sx)
|
|---|
| 337 | {
|
|---|
| 338 | last unless defined $sy->[$j] || $car;
|
|---|
| 339 | $i += $BASE if $car = (($i -= ($sy->[$j] || 0) + $car) < 0); $j++;
|
|---|
| 340 | }
|
|---|
| 341 | # might leave leading zeros, so fix that
|
|---|
| 342 | return __strip_zeros($sx);
|
|---|
| 343 | }
|
|---|
| 344 | for $i (@$sx)
|
|---|
| 345 | {
|
|---|
| 346 | # we can't do an early out if $x is < than $y, since we
|
|---|
| 347 | # need to copy the high chunks from $y. Found by Bob Mathews.
|
|---|
| 348 | #last unless defined $sy->[$j] || $car;
|
|---|
| 349 | $sy->[$j] += $BASE
|
|---|
| 350 | if $car = (($sy->[$j] = $i-($sy->[$j]||0) - $car) < 0);
|
|---|
| 351 | $j++;
|
|---|
| 352 | }
|
|---|
| 353 | # might leave leading zeros, so fix that
|
|---|
| 354 | __strip_zeros($sy);
|
|---|
| 355 | }
|
|---|
| 356 |
|
|---|
| 357 | sub _mul_use_mul
|
|---|
| 358 | {
|
|---|
| 359 | # (ref to int_num_array, ref to int_num_array)
|
|---|
| 360 | # multiply two numbers in internal representation
|
|---|
| 361 | # modifies first arg, second need not be different from first
|
|---|
| 362 | my ($c,$xv,$yv) = @_;
|
|---|
| 363 |
|
|---|
| 364 | if (@$yv == 1)
|
|---|
| 365 | {
|
|---|
| 366 | # shortcut for two very short numbers (improved by Nathan Zook)
|
|---|
| 367 | # works also if xv and yv are the same reference, and handles also $x == 0
|
|---|
| 368 | if (@$xv == 1)
|
|---|
| 369 | {
|
|---|
| 370 | if (($xv->[0] *= $yv->[0]) >= $MBASE)
|
|---|
| 371 | {
|
|---|
| 372 | $xv->[0] = $xv->[0] - ($xv->[1] = int($xv->[0] * $RBASE)) * $MBASE;
|
|---|
| 373 | };
|
|---|
| 374 | return $xv;
|
|---|
| 375 | }
|
|---|
| 376 | # $x * 0 => 0
|
|---|
| 377 | if ($yv->[0] == 0)
|
|---|
| 378 | {
|
|---|
| 379 | @$xv = (0);
|
|---|
| 380 | return $xv;
|
|---|
| 381 | }
|
|---|
| 382 | # multiply a large number a by a single element one, so speed up
|
|---|
| 383 | my $y = $yv->[0]; my $car = 0;
|
|---|
| 384 | foreach my $i (@$xv)
|
|---|
| 385 | {
|
|---|
| 386 | $i = $i * $y + $car; $car = int($i * $RBASE); $i -= $car * $MBASE;
|
|---|
| 387 | }
|
|---|
| 388 | push @$xv, $car if $car != 0;
|
|---|
| 389 | return $xv;
|
|---|
| 390 | }
|
|---|
| 391 | # shortcut for result $x == 0 => result = 0
|
|---|
| 392 | return $xv if ( ((@$xv == 1) && ($xv->[0] == 0)) );
|
|---|
| 393 |
|
|---|
| 394 | # since multiplying $x with $x fails, make copy in this case
|
|---|
| 395 | $yv = [@$xv] if $xv == $yv; # same references?
|
|---|
| 396 |
|
|---|
| 397 | my @prod = (); my ($prod,$car,$cty,$xi,$yi);
|
|---|
| 398 |
|
|---|
| 399 | for $xi (@$xv)
|
|---|
| 400 | {
|
|---|
| 401 | $car = 0; $cty = 0;
|
|---|
| 402 |
|
|---|
| 403 | # slow variant
|
|---|
| 404 | # for $yi (@$yv)
|
|---|
| 405 | # {
|
|---|
| 406 | # $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
|
|---|
| 407 | # $prod[$cty++] =
|
|---|
| 408 | # $prod - ($car = int($prod * RBASE)) * $MBASE; # see USE_MUL
|
|---|
| 409 | # }
|
|---|
| 410 | # $prod[$cty] += $car if $car; # need really to check for 0?
|
|---|
| 411 | # $xi = shift @prod;
|
|---|
| 412 |
|
|---|
| 413 | # faster variant
|
|---|
| 414 | # looping through this if $xi == 0 is silly - so optimize it away!
|
|---|
| 415 | $xi = (shift @prod || 0), next if $xi == 0;
|
|---|
| 416 | for $yi (@$yv)
|
|---|
| 417 | {
|
|---|
| 418 | $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
|
|---|
| 419 | ## this is actually a tad slower
|
|---|
| 420 | ## $prod = $prod[$cty]; $prod += ($car + $xi * $yi); # no ||0 here
|
|---|
| 421 | $prod[$cty++] =
|
|---|
| 422 | $prod - ($car = int($prod * $RBASE)) * $MBASE; # see USE_MUL
|
|---|
| 423 | }
|
|---|
| 424 | $prod[$cty] += $car if $car; # need really to check for 0?
|
|---|
| 425 | $xi = shift @prod || 0; # || 0 makes v5.005_3 happy
|
|---|
| 426 | }
|
|---|
| 427 | push @$xv, @prod;
|
|---|
| 428 | __strip_zeros($xv);
|
|---|
| 429 | $xv;
|
|---|
| 430 | }
|
|---|
| 431 |
|
|---|
| 432 | sub _mul_use_div
|
|---|
| 433 | {
|
|---|
| 434 | # (ref to int_num_array, ref to int_num_array)
|
|---|
| 435 | # multiply two numbers in internal representation
|
|---|
| 436 | # modifies first arg, second need not be different from first
|
|---|
| 437 | my ($c,$xv,$yv) = @_;
|
|---|
| 438 |
|
|---|
| 439 | if (@$yv == 1)
|
|---|
| 440 | {
|
|---|
| 441 | # shortcut for two small numbers, also handles $x == 0
|
|---|
| 442 | if (@$xv == 1)
|
|---|
| 443 | {
|
|---|
| 444 | # shortcut for two very short numbers (improved by Nathan Zook)
|
|---|
| 445 | # works also if xv and yv are the same reference, and handles also $x == 0
|
|---|
| 446 | if (($xv->[0] *= $yv->[0]) >= $MBASE)
|
|---|
| 447 | {
|
|---|
| 448 | $xv->[0] =
|
|---|
| 449 | $xv->[0] - ($xv->[1] = int($xv->[0] / $MBASE)) * $MBASE;
|
|---|
| 450 | };
|
|---|
| 451 | return $xv;
|
|---|
| 452 | }
|
|---|
| 453 | # $x * 0 => 0
|
|---|
| 454 | if ($yv->[0] == 0)
|
|---|
| 455 | {
|
|---|
| 456 | @$xv = (0);
|
|---|
| 457 | return $xv;
|
|---|
| 458 | }
|
|---|
| 459 | # multiply a large number a by a single element one, so speed up
|
|---|
| 460 | my $y = $yv->[0]; my $car = 0;
|
|---|
| 461 | foreach my $i (@$xv)
|
|---|
| 462 | {
|
|---|
| 463 | $i = $i * $y + $car; $car = int($i / $MBASE); $i -= $car * $MBASE;
|
|---|
| 464 | }
|
|---|
| 465 | push @$xv, $car if $car != 0;
|
|---|
| 466 | return $xv;
|
|---|
| 467 | }
|
|---|
| 468 | # shortcut for result $x == 0 => result = 0
|
|---|
| 469 | return $xv if ( ((@$xv == 1) && ($xv->[0] == 0)) );
|
|---|
| 470 |
|
|---|
| 471 | # since multiplying $x with $x fails, make copy in this case
|
|---|
| 472 | $yv = [@$xv] if $xv == $yv; # same references?
|
|---|
| 473 |
|
|---|
| 474 | my @prod = (); my ($prod,$car,$cty,$xi,$yi);
|
|---|
| 475 | for $xi (@$xv)
|
|---|
| 476 | {
|
|---|
| 477 | $car = 0; $cty = 0;
|
|---|
| 478 | # looping through this if $xi == 0 is silly - so optimize it away!
|
|---|
| 479 | $xi = (shift @prod || 0), next if $xi == 0;
|
|---|
| 480 | for $yi (@$yv)
|
|---|
| 481 | {
|
|---|
| 482 | $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
|
|---|
| 483 | $prod[$cty++] =
|
|---|
| 484 | $prod - ($car = int($prod / $MBASE)) * $MBASE;
|
|---|
| 485 | }
|
|---|
| 486 | $prod[$cty] += $car if $car; # need really to check for 0?
|
|---|
| 487 | $xi = shift @prod || 0; # || 0 makes v5.005_3 happy
|
|---|
| 488 | }
|
|---|
| 489 | push @$xv, @prod;
|
|---|
| 490 | __strip_zeros($xv);
|
|---|
| 491 | $xv;
|
|---|
| 492 | }
|
|---|
| 493 |
|
|---|
| 494 | sub _div_use_mul
|
|---|
| 495 | {
|
|---|
| 496 | # ref to array, ref to array, modify first array and return remainder if
|
|---|
| 497 | # in list context
|
|---|
| 498 |
|
|---|
| 499 | # see comments in _div_use_div() for more explanations
|
|---|
| 500 |
|
|---|
| 501 | my ($c,$x,$yorg) = @_;
|
|---|
| 502 |
|
|---|
| 503 | # the general div algorithmn here is about O(N*N) and thus quite slow, so
|
|---|
| 504 | # we first check for some special cases and use shortcuts to handle them.
|
|---|
| 505 |
|
|---|
| 506 | # This works, because we store the numbers in a chunked format where each
|
|---|
| 507 | # element contains 5..7 digits (depending on system).
|
|---|
| 508 |
|
|---|
| 509 | # if both numbers have only one element:
|
|---|
| 510 | if (@$x == 1 && @$yorg == 1)
|
|---|
| 511 | {
|
|---|
| 512 | # shortcut, $yorg and $x are two small numbers
|
|---|
| 513 | if (wantarray)
|
|---|
| 514 | {
|
|---|
| 515 | my $r = [ $x->[0] % $yorg->[0] ];
|
|---|
| 516 | $x->[0] = int($x->[0] / $yorg->[0]);
|
|---|
| 517 | return ($x,$r);
|
|---|
| 518 | }
|
|---|
| 519 | else
|
|---|
| 520 | {
|
|---|
| 521 | $x->[0] = int($x->[0] / $yorg->[0]);
|
|---|
| 522 | return $x;
|
|---|
| 523 | }
|
|---|
| 524 | }
|
|---|
| 525 |
|
|---|
| 526 | # if x has more than one, but y has only one element:
|
|---|
| 527 | if (@$yorg == 1)
|
|---|
| 528 | {
|
|---|
| 529 | my $rem;
|
|---|
| 530 | $rem = _mod($c,[ @$x ],$yorg) if wantarray;
|
|---|
| 531 |
|
|---|
| 532 | # shortcut, $y is < $BASE
|
|---|
| 533 | my $j = scalar @$x; my $r = 0;
|
|---|
| 534 | my $y = $yorg->[0]; my $b;
|
|---|
| 535 | while ($j-- > 0)
|
|---|
| 536 | {
|
|---|
| 537 | $b = $r * $MBASE + $x->[$j];
|
|---|
| 538 | $x->[$j] = int($b/$y);
|
|---|
| 539 | $r = $b % $y;
|
|---|
| 540 | }
|
|---|
| 541 | pop @$x if @$x > 1 && $x->[-1] == 0; # splice up a leading zero
|
|---|
| 542 | return ($x,$rem) if wantarray;
|
|---|
| 543 | return $x;
|
|---|
| 544 | }
|
|---|
| 545 |
|
|---|
| 546 | # now x and y have more than one element
|
|---|
| 547 |
|
|---|
| 548 | # check whether y has more elements than x, if yet, the result will be 0
|
|---|
| 549 | if (@$yorg > @$x)
|
|---|
| 550 | {
|
|---|
| 551 | my $rem;
|
|---|
| 552 | $rem = [@$x] if wantarray; # make copy
|
|---|
| 553 | splice (@$x,1); # keep ref to original array
|
|---|
| 554 | $x->[0] = 0; # set to 0
|
|---|
| 555 | return ($x,$rem) if wantarray; # including remainder?
|
|---|
| 556 | return $x; # only x, which is [0] now
|
|---|
| 557 | }
|
|---|
| 558 | # check whether the numbers have the same number of elements, in that case
|
|---|
| 559 | # the result will fit into one element and can be computed efficiently
|
|---|
| 560 | if (@$yorg == @$x)
|
|---|
| 561 | {
|
|---|
| 562 | my $rem;
|
|---|
| 563 | # if $yorg has more digits than $x (it's leading element is longer than
|
|---|
| 564 | # the one from $x), the result will also be 0:
|
|---|
| 565 | if (length(int($yorg->[-1])) > length(int($x->[-1])))
|
|---|
| 566 | {
|
|---|
| 567 | $rem = [@$x] if wantarray; # make copy
|
|---|
| 568 | splice (@$x,1); # keep ref to org array
|
|---|
| 569 | $x->[0] = 0; # set to 0
|
|---|
| 570 | return ($x,$rem) if wantarray; # including remainder?
|
|---|
| 571 | return $x;
|
|---|
| 572 | }
|
|---|
| 573 | # now calculate $x / $yorg
|
|---|
| 574 | if (length(int($yorg->[-1])) == length(int($x->[-1])))
|
|---|
| 575 | {
|
|---|
| 576 | # same length, so make full compare
|
|---|
| 577 |
|
|---|
| 578 | my $a = 0; my $j = scalar @$x - 1;
|
|---|
| 579 | # manual way (abort if unequal, good for early ne)
|
|---|
| 580 | while ($j >= 0)
|
|---|
| 581 | {
|
|---|
| 582 | last if ($a = $x->[$j] - $yorg->[$j]); $j--;
|
|---|
| 583 | }
|
|---|
| 584 | # $a contains the result of the compare between X and Y
|
|---|
| 585 | # a < 0: x < y, a == 0: x == y, a > 0: x > y
|
|---|
| 586 | if ($a <= 0)
|
|---|
| 587 | {
|
|---|
| 588 | $rem = [ 0 ]; # a = 0 => x == y => rem 0
|
|---|
| 589 | $rem = [@$x] if $a != 0; # a < 0 => x < y => rem = x
|
|---|
| 590 | splice(@$x,1); # keep single element
|
|---|
| 591 | $x->[0] = 0; # if $a < 0
|
|---|
| 592 | $x->[0] = 1 if $a == 0; # $x == $y
|
|---|
| 593 | return ($x,$rem) if wantarray;
|
|---|
| 594 | return $x;
|
|---|
| 595 | }
|
|---|
| 596 | # $x >= $y, so proceed normally
|
|---|
| 597 | }
|
|---|
| 598 | }
|
|---|
| 599 |
|
|---|
| 600 | # all other cases:
|
|---|
| 601 |
|
|---|
| 602 | my $y = [ @$yorg ]; # always make copy to preserve
|
|---|
| 603 |
|
|---|
| 604 | my ($car,$bar,$prd,$dd,$xi,$yi,@q,$v2,$v1,@d,$tmp,$q,$u2,$u1,$u0);
|
|---|
| 605 |
|
|---|
| 606 | $car = $bar = $prd = 0;
|
|---|
| 607 | if (($dd = int($MBASE/($y->[-1]+1))) != 1)
|
|---|
| 608 | {
|
|---|
| 609 | for $xi (@$x)
|
|---|
| 610 | {
|
|---|
| 611 | $xi = $xi * $dd + $car;
|
|---|
| 612 | $xi -= ($car = int($xi * $RBASE)) * $MBASE; # see USE_MUL
|
|---|
| 613 | }
|
|---|
| 614 | push(@$x, $car); $car = 0;
|
|---|
| 615 | for $yi (@$y)
|
|---|
| 616 | {
|
|---|
| 617 | $yi = $yi * $dd + $car;
|
|---|
| 618 | $yi -= ($car = int($yi * $RBASE)) * $MBASE; # see USE_MUL
|
|---|
| 619 | }
|
|---|
| 620 | }
|
|---|
| 621 | else
|
|---|
| 622 | {
|
|---|
| 623 | push(@$x, 0);
|
|---|
| 624 | }
|
|---|
| 625 | @q = (); ($v2,$v1) = @$y[-2,-1];
|
|---|
| 626 | $v2 = 0 unless $v2;
|
|---|
| 627 | while ($#$x > $#$y)
|
|---|
| 628 | {
|
|---|
| 629 | ($u2,$u1,$u0) = @$x[-3..-1];
|
|---|
| 630 | $u2 = 0 unless $u2;
|
|---|
| 631 | #warn "oups v1 is 0, u0: $u0 $y->[-2] $y->[-1] l ",scalar @$y,"\n"
|
|---|
| 632 | # if $v1 == 0;
|
|---|
| 633 | $q = (($u0 == $v1) ? $MAX_VAL : int(($u0*$MBASE+$u1)/$v1));
|
|---|
| 634 | --$q while ($v2*$q > ($u0*$MBASE+$u1-$q*$v1)*$MBASE+$u2);
|
|---|
| 635 | if ($q)
|
|---|
| 636 | {
|
|---|
| 637 | ($car, $bar) = (0,0);
|
|---|
| 638 | for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
|
|---|
| 639 | {
|
|---|
| 640 | $prd = $q * $y->[$yi] + $car;
|
|---|
| 641 | $prd -= ($car = int($prd * $RBASE)) * $MBASE; # see USE_MUL
|
|---|
| 642 | $x->[$xi] += $MBASE if ($bar = (($x->[$xi] -= $prd + $bar) < 0));
|
|---|
| 643 | }
|
|---|
| 644 | if ($x->[-1] < $car + $bar)
|
|---|
| 645 | {
|
|---|
| 646 | $car = 0; --$q;
|
|---|
| 647 | for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
|
|---|
| 648 | {
|
|---|
| 649 | $x->[$xi] -= $MBASE
|
|---|
| 650 | if ($car = (($x->[$xi] += $y->[$yi] + $car) >= $MBASE));
|
|---|
| 651 | }
|
|---|
| 652 | }
|
|---|
| 653 | }
|
|---|
| 654 | pop(@$x);
|
|---|
| 655 | unshift(@q, $q);
|
|---|
| 656 | }
|
|---|
| 657 | if (wantarray)
|
|---|
| 658 | {
|
|---|
| 659 | @d = ();
|
|---|
| 660 | if ($dd != 1)
|
|---|
| 661 | {
|
|---|
| 662 | $car = 0;
|
|---|
| 663 | for $xi (reverse @$x)
|
|---|
| 664 | {
|
|---|
| 665 | $prd = $car * $MBASE + $xi;
|
|---|
| 666 | $car = $prd - ($tmp = int($prd / $dd)) * $dd; # see USE_MUL
|
|---|
| 667 | unshift(@d, $tmp);
|
|---|
| 668 | }
|
|---|
| 669 | }
|
|---|
| 670 | else
|
|---|
| 671 | {
|
|---|
| 672 | @d = @$x;
|
|---|
| 673 | }
|
|---|
| 674 | @$x = @q;
|
|---|
| 675 | my $d = \@d;
|
|---|
| 676 | __strip_zeros($x);
|
|---|
| 677 | __strip_zeros($d);
|
|---|
| 678 | return ($x,$d);
|
|---|
| 679 | }
|
|---|
| 680 | @$x = @q;
|
|---|
| 681 | __strip_zeros($x);
|
|---|
| 682 | $x;
|
|---|
| 683 | }
|
|---|
| 684 |
|
|---|
| 685 | sub _div_use_div
|
|---|
| 686 | {
|
|---|
| 687 | # ref to array, ref to array, modify first array and return remainder if
|
|---|
| 688 | # in list context
|
|---|
| 689 | my ($c,$x,$yorg) = @_;
|
|---|
| 690 |
|
|---|
| 691 | # the general div algorithmn here is about O(N*N) and thus quite slow, so
|
|---|
| 692 | # we first check for some special cases and use shortcuts to handle them.
|
|---|
| 693 |
|
|---|
| 694 | # This works, because we store the numbers in a chunked format where each
|
|---|
| 695 | # element contains 5..7 digits (depending on system).
|
|---|
| 696 |
|
|---|
| 697 | # if both numbers have only one element:
|
|---|
| 698 | if (@$x == 1 && @$yorg == 1)
|
|---|
| 699 | {
|
|---|
| 700 | # shortcut, $yorg and $x are two small numbers
|
|---|
| 701 | if (wantarray)
|
|---|
| 702 | {
|
|---|
| 703 | my $r = [ $x->[0] % $yorg->[0] ];
|
|---|
| 704 | $x->[0] = int($x->[0] / $yorg->[0]);
|
|---|
| 705 | return ($x,$r);
|
|---|
| 706 | }
|
|---|
| 707 | else
|
|---|
| 708 | {
|
|---|
| 709 | $x->[0] = int($x->[0] / $yorg->[0]);
|
|---|
| 710 | return $x;
|
|---|
| 711 | }
|
|---|
| 712 | }
|
|---|
| 713 | # if x has more than one, but y has only one element:
|
|---|
| 714 | if (@$yorg == 1)
|
|---|
| 715 | {
|
|---|
| 716 | my $rem;
|
|---|
| 717 | $rem = _mod($c,[ @$x ],$yorg) if wantarray;
|
|---|
| 718 |
|
|---|
| 719 | # shortcut, $y is < $BASE
|
|---|
| 720 | my $j = scalar @$x; my $r = 0;
|
|---|
| 721 | my $y = $yorg->[0]; my $b;
|
|---|
| 722 | while ($j-- > 0)
|
|---|
| 723 | {
|
|---|
| 724 | $b = $r * $MBASE + $x->[$j];
|
|---|
| 725 | $x->[$j] = int($b/$y);
|
|---|
| 726 | $r = $b % $y;
|
|---|
| 727 | }
|
|---|
| 728 | pop @$x if @$x > 1 && $x->[-1] == 0; # splice up a leading zero
|
|---|
| 729 | return ($x,$rem) if wantarray;
|
|---|
| 730 | return $x;
|
|---|
| 731 | }
|
|---|
| 732 | # now x and y have more than one element
|
|---|
| 733 |
|
|---|
| 734 | # check whether y has more elements than x, if yet, the result will be 0
|
|---|
| 735 | if (@$yorg > @$x)
|
|---|
| 736 | {
|
|---|
| 737 | my $rem;
|
|---|
| 738 | $rem = [@$x] if wantarray; # make copy
|
|---|
| 739 | splice (@$x,1); # keep ref to original array
|
|---|
| 740 | $x->[0] = 0; # set to 0
|
|---|
| 741 | return ($x,$rem) if wantarray; # including remainder?
|
|---|
| 742 | return $x; # only x, which is [0] now
|
|---|
| 743 | }
|
|---|
| 744 | # check whether the numbers have the same number of elements, in that case
|
|---|
| 745 | # the result will fit into one element and can be computed efficiently
|
|---|
| 746 | if (@$yorg == @$x)
|
|---|
| 747 | {
|
|---|
| 748 | my $rem;
|
|---|
| 749 | # if $yorg has more digits than $x (it's leading element is longer than
|
|---|
| 750 | # the one from $x), the result will also be 0:
|
|---|
| 751 | if (length(int($yorg->[-1])) > length(int($x->[-1])))
|
|---|
| 752 | {
|
|---|
| 753 | $rem = [@$x] if wantarray; # make copy
|
|---|
| 754 | splice (@$x,1); # keep ref to org array
|
|---|
| 755 | $x->[0] = 0; # set to 0
|
|---|
| 756 | return ($x,$rem) if wantarray; # including remainder?
|
|---|
| 757 | return $x;
|
|---|
| 758 | }
|
|---|
| 759 | # now calculate $x / $yorg
|
|---|
| 760 |
|
|---|
| 761 | if (length(int($yorg->[-1])) == length(int($x->[-1])))
|
|---|
| 762 | {
|
|---|
| 763 | # same length, so make full compare
|
|---|
| 764 |
|
|---|
| 765 | my $a = 0; my $j = scalar @$x - 1;
|
|---|
| 766 | # manual way (abort if unequal, good for early ne)
|
|---|
| 767 | while ($j >= 0)
|
|---|
| 768 | {
|
|---|
| 769 | last if ($a = $x->[$j] - $yorg->[$j]); $j--;
|
|---|
| 770 | }
|
|---|
| 771 | # $a contains the result of the compare between X and Y
|
|---|
| 772 | # a < 0: x < y, a == 0: x == y, a > 0: x > y
|
|---|
| 773 | if ($a <= 0)
|
|---|
| 774 | {
|
|---|
| 775 | $rem = [ 0 ]; # a = 0 => x == y => rem 0
|
|---|
| 776 | $rem = [@$x] if $a != 0; # a < 0 => x < y => rem = x
|
|---|
| 777 | splice(@$x,1); # keep single element
|
|---|
| 778 | $x->[0] = 0; # if $a < 0
|
|---|
| 779 | $x->[0] = 1 if $a == 0; # $x == $y
|
|---|
| 780 | return ($x,$rem) if wantarray; # including remainder?
|
|---|
| 781 | return $x;
|
|---|
| 782 | }
|
|---|
| 783 | # $x >= $y, so proceed normally
|
|---|
| 784 |
|
|---|
| 785 | }
|
|---|
| 786 | }
|
|---|
| 787 |
|
|---|
| 788 | # all other cases:
|
|---|
| 789 |
|
|---|
| 790 | my $y = [ @$yorg ]; # always make copy to preserve
|
|---|
| 791 |
|
|---|
| 792 | my ($car,$bar,$prd,$dd,$xi,$yi,@q,$v2,$v1,@d,$tmp,$q,$u2,$u1,$u0);
|
|---|
| 793 |
|
|---|
| 794 | $car = $bar = $prd = 0;
|
|---|
| 795 | if (($dd = int($MBASE/($y->[-1]+1))) != 1)
|
|---|
| 796 | {
|
|---|
| 797 | for $xi (@$x)
|
|---|
| 798 | {
|
|---|
| 799 | $xi = $xi * $dd + $car;
|
|---|
| 800 | $xi -= ($car = int($xi / $MBASE)) * $MBASE;
|
|---|
| 801 | }
|
|---|
| 802 | push(@$x, $car); $car = 0;
|
|---|
| 803 | for $yi (@$y)
|
|---|
| 804 | {
|
|---|
| 805 | $yi = $yi * $dd + $car;
|
|---|
| 806 | $yi -= ($car = int($yi / $MBASE)) * $MBASE;
|
|---|
| 807 | }
|
|---|
| 808 | }
|
|---|
| 809 | else
|
|---|
| 810 | {
|
|---|
| 811 | push(@$x, 0);
|
|---|
| 812 | }
|
|---|
| 813 |
|
|---|
| 814 | # @q will accumulate the final result, $q contains the current computed
|
|---|
| 815 | # part of the final result
|
|---|
| 816 |
|
|---|
| 817 | @q = (); ($v2,$v1) = @$y[-2,-1];
|
|---|
| 818 | $v2 = 0 unless $v2;
|
|---|
| 819 | while ($#$x > $#$y)
|
|---|
| 820 | {
|
|---|
| 821 | ($u2,$u1,$u0) = @$x[-3..-1];
|
|---|
| 822 | $u2 = 0 unless $u2;
|
|---|
| 823 | #warn "oups v1 is 0, u0: $u0 $y->[-2] $y->[-1] l ",scalar @$y,"\n"
|
|---|
| 824 | # if $v1 == 0;
|
|---|
| 825 | $q = (($u0 == $v1) ? $MAX_VAL : int(($u0*$MBASE+$u1)/$v1));
|
|---|
| 826 | --$q while ($v2*$q > ($u0*$MBASE+$u1-$q*$v1)*$MBASE+$u2);
|
|---|
| 827 | if ($q)
|
|---|
| 828 | {
|
|---|
| 829 | ($car, $bar) = (0,0);
|
|---|
| 830 | for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
|
|---|
| 831 | {
|
|---|
| 832 | $prd = $q * $y->[$yi] + $car;
|
|---|
| 833 | $prd -= ($car = int($prd / $MBASE)) * $MBASE;
|
|---|
| 834 | $x->[$xi] += $MBASE if ($bar = (($x->[$xi] -= $prd + $bar) < 0));
|
|---|
| 835 | }
|
|---|
| 836 | if ($x->[-1] < $car + $bar)
|
|---|
| 837 | {
|
|---|
| 838 | $car = 0; --$q;
|
|---|
| 839 | for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
|
|---|
| 840 | {
|
|---|
| 841 | $x->[$xi] -= $MBASE
|
|---|
| 842 | if ($car = (($x->[$xi] += $y->[$yi] + $car) >= $MBASE));
|
|---|
| 843 | }
|
|---|
| 844 | }
|
|---|
| 845 | }
|
|---|
| 846 | pop(@$x); unshift(@q, $q);
|
|---|
| 847 | }
|
|---|
| 848 | if (wantarray)
|
|---|
| 849 | {
|
|---|
| 850 | @d = ();
|
|---|
| 851 | if ($dd != 1)
|
|---|
| 852 | {
|
|---|
| 853 | $car = 0;
|
|---|
| 854 | for $xi (reverse @$x)
|
|---|
| 855 | {
|
|---|
| 856 | $prd = $car * $MBASE + $xi;
|
|---|
| 857 | $car = $prd - ($tmp = int($prd / $dd)) * $dd;
|
|---|
| 858 | unshift(@d, $tmp);
|
|---|
| 859 | }
|
|---|
| 860 | }
|
|---|
| 861 | else
|
|---|
| 862 | {
|
|---|
| 863 | @d = @$x;
|
|---|
| 864 | }
|
|---|
| 865 | @$x = @q;
|
|---|
| 866 | my $d = \@d;
|
|---|
| 867 | __strip_zeros($x);
|
|---|
| 868 | __strip_zeros($d);
|
|---|
| 869 | return ($x,$d);
|
|---|
| 870 | }
|
|---|
| 871 | @$x = @q;
|
|---|
| 872 | __strip_zeros($x);
|
|---|
| 873 | $x;
|
|---|
| 874 | }
|
|---|
| 875 |
|
|---|
| 876 | ##############################################################################
|
|---|
| 877 | # testing
|
|---|
| 878 |
|
|---|
| 879 | sub _acmp
|
|---|
| 880 | {
|
|---|
| 881 | # internal absolute post-normalized compare (ignore signs)
|
|---|
| 882 | # ref to array, ref to array, return <0, 0, >0
|
|---|
| 883 | # arrays must have at least one entry; this is not checked for
|
|---|
| 884 | my ($c,$cx,$cy) = @_;
|
|---|
| 885 |
|
|---|
| 886 | # shortcut for short numbers
|
|---|
| 887 | return (($cx->[0] <=> $cy->[0]) <=> 0)
|
|---|
| 888 | if scalar @$cx == scalar @$cy && scalar @$cx == 1;
|
|---|
| 889 |
|
|---|
| 890 | # fast comp based on number of array elements (aka pseudo-length)
|
|---|
| 891 | my $lxy = (scalar @$cx - scalar @$cy)
|
|---|
| 892 | # or length of first element if same number of elements (aka difference 0)
|
|---|
| 893 | ||
|
|---|
| 894 | # need int() here because sometimes the last element is '00018' vs '18'
|
|---|
| 895 | (length(int($cx->[-1])) - length(int($cy->[-1])));
|
|---|
| 896 | return -1 if $lxy < 0; # already differs, ret
|
|---|
| 897 | return 1 if $lxy > 0; # ditto
|
|---|
| 898 |
|
|---|
| 899 | # manual way (abort if unequal, good for early ne)
|
|---|
| 900 | my $a; my $j = scalar @$cx;
|
|---|
| 901 | while (--$j >= 0)
|
|---|
| 902 | {
|
|---|
| 903 | last if ($a = $cx->[$j] - $cy->[$j]);
|
|---|
| 904 | }
|
|---|
| 905 | $a <=> 0;
|
|---|
| 906 | }
|
|---|
| 907 |
|
|---|
| 908 | sub _len
|
|---|
| 909 | {
|
|---|
| 910 | # compute number of digits
|
|---|
| 911 |
|
|---|
| 912 | # int() because add/sub sometimes leaves strings (like '00005') instead of
|
|---|
| 913 | # '5' in this place, thus causing length() to report wrong length
|
|---|
| 914 | my $cx = $_[1];
|
|---|
| 915 |
|
|---|
| 916 | (@$cx-1)*$BASE_LEN+length(int($cx->[-1]));
|
|---|
| 917 | }
|
|---|
| 918 |
|
|---|
| 919 | sub _digit
|
|---|
| 920 | {
|
|---|
| 921 | # return the nth digit, negative values count backward
|
|---|
| 922 | # zero is rightmost, so _digit(123,0) will give 3
|
|---|
| 923 | my ($c,$x,$n) = @_;
|
|---|
| 924 |
|
|---|
| 925 | my $len = _len('',$x);
|
|---|
| 926 |
|
|---|
| 927 | $n = $len+$n if $n < 0; # -1 last, -2 second-to-last
|
|---|
| 928 | $n = abs($n); # if negative was too big
|
|---|
| 929 | $len--; $n = $len if $n > $len; # n to big?
|
|---|
| 930 |
|
|---|
| 931 | my $elem = int($n / $BASE_LEN); # which array element
|
|---|
| 932 | my $digit = $n % $BASE_LEN; # which digit in this element
|
|---|
| 933 | $elem = '0' x $BASE_LEN . @$x[$elem]; # get element padded with 0's
|
|---|
| 934 | substr($elem,-$digit-1,1);
|
|---|
| 935 | }
|
|---|
| 936 |
|
|---|
| 937 | sub _zeros
|
|---|
| 938 | {
|
|---|
| 939 | # return amount of trailing zeros in decimal
|
|---|
| 940 | # check each array elem in _m for having 0 at end as long as elem == 0
|
|---|
| 941 | # Upon finding a elem != 0, stop
|
|---|
| 942 | my $x = $_[1];
|
|---|
| 943 |
|
|---|
| 944 | return 0 if scalar @$x == 1 && $x->[0] == 0;
|
|---|
| 945 |
|
|---|
| 946 | my $zeros = 0; my $elem;
|
|---|
| 947 | foreach my $e (@$x)
|
|---|
| 948 | {
|
|---|
| 949 | if ($e != 0)
|
|---|
| 950 | {
|
|---|
| 951 | $elem = "$e"; # preserve x
|
|---|
| 952 | $elem =~ s/.*?(0*$)/$1/; # strip anything not zero
|
|---|
| 953 | $zeros *= $BASE_LEN; # elems * 5
|
|---|
| 954 | $zeros += length($elem); # count trailing zeros
|
|---|
| 955 | last; # early out
|
|---|
| 956 | }
|
|---|
| 957 | $zeros ++; # real else branch: 50% slower!
|
|---|
| 958 | }
|
|---|
| 959 | $zeros;
|
|---|
| 960 | }
|
|---|
| 961 |
|
|---|
| 962 | ##############################################################################
|
|---|
| 963 | # _is_* routines
|
|---|
| 964 |
|
|---|
| 965 | sub _is_zero
|
|---|
| 966 | {
|
|---|
| 967 | # return true if arg is zero
|
|---|
| 968 | (((scalar @{$_[1]} == 1) && ($_[1]->[0] == 0))) <=> 0;
|
|---|
| 969 | }
|
|---|
| 970 |
|
|---|
| 971 | sub _is_even
|
|---|
| 972 | {
|
|---|
| 973 | # return true if arg is even
|
|---|
| 974 | (!($_[1]->[0] & 1)) <=> 0;
|
|---|
| 975 | }
|
|---|
| 976 |
|
|---|
| 977 | sub _is_odd
|
|---|
| 978 | {
|
|---|
| 979 | # return true if arg is even
|
|---|
| 980 | (($_[1]->[0] & 1)) <=> 0;
|
|---|
| 981 | }
|
|---|
| 982 |
|
|---|
| 983 | sub _is_one
|
|---|
| 984 | {
|
|---|
| 985 | # return true if arg is one
|
|---|
| 986 | (scalar @{$_[1]} == 1) && ($_[1]->[0] == 1) <=> 0;
|
|---|
| 987 | }
|
|---|
| 988 |
|
|---|
| 989 | sub _is_two
|
|---|
| 990 | {
|
|---|
| 991 | # return true if arg is two
|
|---|
| 992 | (scalar @{$_[1]} == 1) && ($_[1]->[0] == 2) <=> 0;
|
|---|
| 993 | }
|
|---|
| 994 |
|
|---|
| 995 | sub _is_ten
|
|---|
| 996 | {
|
|---|
| 997 | # return true if arg is ten
|
|---|
| 998 | (scalar @{$_[1]} == 1) && ($_[1]->[0] == 10) <=> 0;
|
|---|
| 999 | }
|
|---|
| 1000 |
|
|---|
| 1001 | sub __strip_zeros
|
|---|
| 1002 | {
|
|---|
| 1003 | # internal normalization function that strips leading zeros from the array
|
|---|
| 1004 | # args: ref to array
|
|---|
| 1005 | my $s = shift;
|
|---|
| 1006 |
|
|---|
| 1007 | my $cnt = scalar @$s; # get count of parts
|
|---|
| 1008 | my $i = $cnt-1;
|
|---|
| 1009 | push @$s,0 if $i < 0; # div might return empty results, so fix it
|
|---|
| 1010 |
|
|---|
| 1011 | return $s if @$s == 1; # early out
|
|---|
| 1012 |
|
|---|
| 1013 | #print "strip: cnt $cnt i $i\n";
|
|---|
| 1014 | # '0', '3', '4', '0', '0',
|
|---|
| 1015 | # 0 1 2 3 4
|
|---|
| 1016 | # cnt = 5, i = 4
|
|---|
| 1017 | # i = 4
|
|---|
| 1018 | # i = 3
|
|---|
| 1019 | # => fcnt = cnt - i (5-2 => 3, cnt => 5-1 = 4, throw away from 4th pos)
|
|---|
| 1020 | # >= 1: skip first part (this can be zero)
|
|---|
| 1021 | while ($i > 0) { last if $s->[$i] != 0; $i--; }
|
|---|
| 1022 | $i++; splice @$s,$i if ($i < $cnt); # $i cant be 0
|
|---|
| 1023 | $s;
|
|---|
| 1024 | }
|
|---|
| 1025 |
|
|---|
| 1026 | ###############################################################################
|
|---|
| 1027 | # check routine to test internal state for corruptions
|
|---|
| 1028 |
|
|---|
| 1029 | sub _check
|
|---|
| 1030 | {
|
|---|
| 1031 | # used by the test suite
|
|---|
| 1032 | my $x = $_[1];
|
|---|
| 1033 |
|
|---|
| 1034 | return "$x is not a reference" if !ref($x);
|
|---|
| 1035 |
|
|---|
| 1036 | # are all parts are valid?
|
|---|
| 1037 | my $i = 0; my $j = scalar @$x; my ($e,$try);
|
|---|
| 1038 | while ($i < $j)
|
|---|
| 1039 | {
|
|---|
| 1040 | $e = $x->[$i]; $e = 'undef' unless defined $e;
|
|---|
| 1041 | $try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e)";
|
|---|
| 1042 | last if $e !~ /^[+]?[0-9]+$/;
|
|---|
| 1043 | $try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e) (stringify)";
|
|---|
| 1044 | last if "$e" !~ /^[+]?[0-9]+$/;
|
|---|
| 1045 | $try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e) (cat-stringify)";
|
|---|
| 1046 | last if '' . "$e" !~ /^[+]?[0-9]+$/;
|
|---|
| 1047 | $try = ' < 0 || >= $BASE; '."($x, $e)";
|
|---|
| 1048 | last if $e <0 || $e >= $BASE;
|
|---|
| 1049 | # this test is disabled, since new/bnorm and certain ops (like early out
|
|---|
| 1050 | # in add/sub) are allowed/expected to leave '00000' in some elements
|
|---|
| 1051 | #$try = '=~ /^00+/; '."($x, $e)";
|
|---|
| 1052 | #last if $e =~ /^00+/;
|
|---|
| 1053 | $i++;
|
|---|
| 1054 | }
|
|---|
| 1055 | return "Illegal part '$e' at pos $i (tested: $try)" if $i < $j;
|
|---|
| 1056 | 0;
|
|---|
| 1057 | }
|
|---|
| 1058 |
|
|---|
| 1059 |
|
|---|
| 1060 | ###############################################################################
|
|---|
| 1061 |
|
|---|
| 1062 | sub _mod
|
|---|
| 1063 | {
|
|---|
| 1064 | # if possible, use mod shortcut
|
|---|
| 1065 | my ($c,$x,$yo) = @_;
|
|---|
| 1066 |
|
|---|
| 1067 | # slow way since $y to big
|
|---|
| 1068 | if (scalar @$yo > 1)
|
|---|
| 1069 | {
|
|---|
| 1070 | my ($xo,$rem) = _div($c,$x,$yo);
|
|---|
| 1071 | return $rem;
|
|---|
| 1072 | }
|
|---|
| 1073 |
|
|---|
| 1074 | my $y = $yo->[0];
|
|---|
| 1075 | # both are single element arrays
|
|---|
| 1076 | if (scalar @$x == 1)
|
|---|
| 1077 | {
|
|---|
| 1078 | $x->[0] %= $y;
|
|---|
| 1079 | return $x;
|
|---|
| 1080 | }
|
|---|
| 1081 |
|
|---|
| 1082 | # @y is a single element, but @x has more than one element
|
|---|
| 1083 | my $b = $BASE % $y;
|
|---|
| 1084 | if ($b == 0)
|
|---|
| 1085 | {
|
|---|
| 1086 | # when BASE % Y == 0 then (B * BASE) % Y == 0
|
|---|
| 1087 | # (B * BASE) % $y + A % Y => A % Y
|
|---|
| 1088 | # so need to consider only last element: O(1)
|
|---|
| 1089 | $x->[0] %= $y;
|
|---|
| 1090 | }
|
|---|
| 1091 | elsif ($b == 1)
|
|---|
| 1092 | {
|
|---|
| 1093 | # else need to go through all elements: O(N), but loop is a bit simplified
|
|---|
| 1094 | my $r = 0;
|
|---|
| 1095 | foreach (@$x)
|
|---|
| 1096 | {
|
|---|
| 1097 | $r = ($r + $_) % $y; # not much faster, but heh...
|
|---|
| 1098 | #$r += $_ % $y; $r %= $y;
|
|---|
| 1099 | }
|
|---|
| 1100 | $r = 0 if $r == $y;
|
|---|
| 1101 | $x->[0] = $r;
|
|---|
| 1102 | }
|
|---|
| 1103 | else
|
|---|
| 1104 | {
|
|---|
| 1105 | # else need to go through all elements: O(N)
|
|---|
| 1106 | my $r = 0; my $bm = 1;
|
|---|
| 1107 | foreach (@$x)
|
|---|
| 1108 | {
|
|---|
| 1109 | $r = ($_ * $bm + $r) % $y;
|
|---|
| 1110 | $bm = ($bm * $b) % $y;
|
|---|
| 1111 |
|
|---|
| 1112 | #$r += ($_ % $y) * $bm;
|
|---|
| 1113 | #$bm *= $b;
|
|---|
| 1114 | #$bm %= $y;
|
|---|
| 1115 | #$r %= $y;
|
|---|
| 1116 | }
|
|---|
| 1117 | $r = 0 if $r == $y;
|
|---|
| 1118 | $x->[0] = $r;
|
|---|
| 1119 | }
|
|---|
| 1120 | splice (@$x,1); # keep one element of $x
|
|---|
| 1121 | $x;
|
|---|
| 1122 | }
|
|---|
| 1123 |
|
|---|
| 1124 | ##############################################################################
|
|---|
| 1125 | # shifts
|
|---|
| 1126 |
|
|---|
| 1127 | sub _rsft
|
|---|
| 1128 | {
|
|---|
| 1129 | my ($c,$x,$y,$n) = @_;
|
|---|
| 1130 |
|
|---|
| 1131 | if ($n != 10)
|
|---|
| 1132 | {
|
|---|
| 1133 | $n = _new($c,$n); return _div($c,$x, _pow($c,$n,$y));
|
|---|
| 1134 | }
|
|---|
| 1135 |
|
|---|
| 1136 | # shortcut (faster) for shifting by 10)
|
|---|
| 1137 | # multiples of $BASE_LEN
|
|---|
| 1138 | my $dst = 0; # destination
|
|---|
| 1139 | my $src = _num($c,$y); # as normal int
|
|---|
| 1140 | my $xlen = (@$x-1)*$BASE_LEN+length(int($x->[-1])); # len of x in digits
|
|---|
| 1141 | if ($src >= $xlen or ($src == $xlen and ! defined $x->[1]))
|
|---|
| 1142 | {
|
|---|
| 1143 | # 12345 67890 shifted right by more than 10 digits => 0
|
|---|
| 1144 | splice (@$x,1); # leave only one element
|
|---|
| 1145 | $x->[0] = 0; # set to zero
|
|---|
| 1146 | return $x;
|
|---|
| 1147 | }
|
|---|
| 1148 | my $rem = $src % $BASE_LEN; # remainder to shift
|
|---|
| 1149 | $src = int($src / $BASE_LEN); # source
|
|---|
| 1150 | if ($rem == 0)
|
|---|
| 1151 | {
|
|---|
| 1152 | splice (@$x,0,$src); # even faster, 38.4 => 39.3
|
|---|
| 1153 | }
|
|---|
| 1154 | else
|
|---|
| 1155 | {
|
|---|
| 1156 | my $len = scalar @$x - $src; # elems to go
|
|---|
| 1157 | my $vd; my $z = '0'x $BASE_LEN;
|
|---|
| 1158 | $x->[scalar @$x] = 0; # avoid || 0 test inside loop
|
|---|
| 1159 | while ($dst < $len)
|
|---|
| 1160 | {
|
|---|
| 1161 | $vd = $z.$x->[$src];
|
|---|
| 1162 | $vd = substr($vd,-$BASE_LEN,$BASE_LEN-$rem);
|
|---|
| 1163 | $src++;
|
|---|
| 1164 | $vd = substr($z.$x->[$src],-$rem,$rem) . $vd;
|
|---|
| 1165 | $vd = substr($vd,-$BASE_LEN,$BASE_LEN) if length($vd) > $BASE_LEN;
|
|---|
| 1166 | $x->[$dst] = int($vd);
|
|---|
| 1167 | $dst++;
|
|---|
| 1168 | }
|
|---|
| 1169 | splice (@$x,$dst) if $dst > 0; # kill left-over array elems
|
|---|
| 1170 | pop @$x if $x->[-1] == 0 && @$x > 1; # kill last element if 0
|
|---|
| 1171 | } # else rem == 0
|
|---|
| 1172 | $x;
|
|---|
| 1173 | }
|
|---|
| 1174 |
|
|---|
| 1175 | sub _lsft
|
|---|
| 1176 | {
|
|---|
| 1177 | my ($c,$x,$y,$n) = @_;
|
|---|
| 1178 |
|
|---|
| 1179 | if ($n != 10)
|
|---|
| 1180 | {
|
|---|
| 1181 | $n = _new($c,$n); return _mul($c,$x, _pow($c,$n,$y));
|
|---|
| 1182 | }
|
|---|
| 1183 |
|
|---|
| 1184 | # shortcut (faster) for shifting by 10) since we are in base 10eX
|
|---|
| 1185 | # multiples of $BASE_LEN:
|
|---|
| 1186 | my $src = scalar @$x; # source
|
|---|
| 1187 | my $len = _num($c,$y); # shift-len as normal int
|
|---|
| 1188 | my $rem = $len % $BASE_LEN; # remainder to shift
|
|---|
| 1189 | my $dst = $src + int($len/$BASE_LEN); # destination
|
|---|
| 1190 | my $vd; # further speedup
|
|---|
| 1191 | $x->[$src] = 0; # avoid first ||0 for speed
|
|---|
| 1192 | my $z = '0' x $BASE_LEN;
|
|---|
| 1193 | while ($src >= 0)
|
|---|
| 1194 | {
|
|---|
| 1195 | $vd = $x->[$src]; $vd = $z.$vd;
|
|---|
| 1196 | $vd = substr($vd,-$BASE_LEN+$rem,$BASE_LEN-$rem);
|
|---|
| 1197 | $vd .= $src > 0 ? substr($z.$x->[$src-1],-$BASE_LEN,$rem) : '0' x $rem;
|
|---|
| 1198 | $vd = substr($vd,-$BASE_LEN,$BASE_LEN) if length($vd) > $BASE_LEN;
|
|---|
| 1199 | $x->[$dst] = int($vd);
|
|---|
| 1200 | $dst--; $src--;
|
|---|
| 1201 | }
|
|---|
| 1202 | # set lowest parts to 0
|
|---|
| 1203 | while ($dst >= 0) { $x->[$dst--] = 0; }
|
|---|
| 1204 | # fix spurios last zero element
|
|---|
| 1205 | splice @$x,-1 if $x->[-1] == 0;
|
|---|
| 1206 | $x;
|
|---|
| 1207 | }
|
|---|
| 1208 |
|
|---|
| 1209 | sub _pow
|
|---|
| 1210 | {
|
|---|
| 1211 | # power of $x to $y
|
|---|
| 1212 | # ref to array, ref to array, return ref to array
|
|---|
| 1213 | my ($c,$cx,$cy) = @_;
|
|---|
| 1214 |
|
|---|
| 1215 | if (scalar @$cy == 1 && $cy->[0] == 0)
|
|---|
| 1216 | {
|
|---|
| 1217 | splice (@$cx,1); $cx->[0] = 1; # y == 0 => x => 1
|
|---|
| 1218 | return $cx;
|
|---|
| 1219 | }
|
|---|
| 1220 | if ((scalar @$cx == 1 && $cx->[0] == 1) || # x == 1
|
|---|
| 1221 | (scalar @$cy == 1 && $cy->[0] == 1)) # or y == 1
|
|---|
| 1222 | {
|
|---|
| 1223 | return $cx;
|
|---|
| 1224 | }
|
|---|
| 1225 | if (scalar @$cx == 1 && $cx->[0] == 0)
|
|---|
| 1226 | {
|
|---|
| 1227 | splice (@$cx,1); $cx->[0] = 0; # 0 ** y => 0 (if not y <= 0)
|
|---|
| 1228 | return $cx;
|
|---|
| 1229 | }
|
|---|
| 1230 |
|
|---|
| 1231 | my $pow2 = _one();
|
|---|
| 1232 |
|
|---|
| 1233 | my $y_bin = _as_bin($c,$cy); $y_bin =~ s/^0b//;
|
|---|
| 1234 | my $len = length($y_bin);
|
|---|
| 1235 | while (--$len > 0)
|
|---|
| 1236 | {
|
|---|
| 1237 | _mul($c,$pow2,$cx) if substr($y_bin,$len,1) eq '1'; # is odd?
|
|---|
| 1238 | _mul($c,$cx,$cx);
|
|---|
| 1239 | }
|
|---|
| 1240 |
|
|---|
| 1241 | _mul($c,$cx,$pow2);
|
|---|
| 1242 | $cx;
|
|---|
| 1243 | }
|
|---|
| 1244 |
|
|---|
| 1245 | sub _fac
|
|---|
| 1246 | {
|
|---|
| 1247 | # factorial of $x
|
|---|
| 1248 | # ref to array, return ref to array
|
|---|
| 1249 | my ($c,$cx) = @_;
|
|---|
| 1250 |
|
|---|
| 1251 | if ((@$cx == 1) && ($cx->[0] <= 2))
|
|---|
| 1252 | {
|
|---|
| 1253 | $cx->[0] ||= 1; # 0 => 1, 1 => 1, 2 => 2
|
|---|
| 1254 | return $cx;
|
|---|
| 1255 | }
|
|---|
| 1256 |
|
|---|
| 1257 | # go forward until $base is exceeded
|
|---|
| 1258 | # limit is either $x steps (steps == 100 means a result always too high) or
|
|---|
| 1259 | # $base.
|
|---|
| 1260 | my $steps = 100; $steps = $cx->[0] if @$cx == 1;
|
|---|
| 1261 | my $r = 2; my $cf = 3; my $step = 2; my $last = $r;
|
|---|
| 1262 | while ($r*$cf < $BASE && $step < $steps)
|
|---|
| 1263 | {
|
|---|
| 1264 | $last = $r; $r *= $cf++; $step++;
|
|---|
| 1265 | }
|
|---|
| 1266 | if ((@$cx == 1) && $step == $cx->[0])
|
|---|
| 1267 | {
|
|---|
| 1268 | # completely done, so keep reference to $x and return
|
|---|
| 1269 | $cx->[0] = $r;
|
|---|
| 1270 | return $cx;
|
|---|
| 1271 | }
|
|---|
| 1272 |
|
|---|
| 1273 | # now we must do the left over steps
|
|---|
| 1274 | my $n; # steps still to do
|
|---|
| 1275 | if (scalar @$cx == 1)
|
|---|
| 1276 | {
|
|---|
| 1277 | $n = $cx->[0];
|
|---|
| 1278 | }
|
|---|
| 1279 | else
|
|---|
| 1280 | {
|
|---|
| 1281 | $n = _copy($c,$cx);
|
|---|
| 1282 | }
|
|---|
| 1283 |
|
|---|
| 1284 | $cx->[0] = $last; splice (@$cx,1); # keep ref to $x
|
|---|
| 1285 | my $zero_elements = 0;
|
|---|
| 1286 |
|
|---|
| 1287 | # do left-over steps fit into a scalar?
|
|---|
| 1288 | if (ref $n eq 'ARRAY')
|
|---|
| 1289 | {
|
|---|
| 1290 | # No, so use slower inc() & cmp()
|
|---|
| 1291 | $step = [$step];
|
|---|
| 1292 | while (_acmp($step,$n) <= 0)
|
|---|
| 1293 | {
|
|---|
| 1294 | # as soon as the last element of $cx is 0, we split it up and remember
|
|---|
| 1295 | # how many zeors we got so far. The reason is that n! will accumulate
|
|---|
| 1296 | # zeros at the end rather fast.
|
|---|
| 1297 | if ($cx->[0] == 0)
|
|---|
| 1298 | {
|
|---|
| 1299 | $zero_elements ++; shift @$cx;
|
|---|
| 1300 | }
|
|---|
| 1301 | _mul($c,$cx,$step); _inc($c,$step);
|
|---|
| 1302 | }
|
|---|
| 1303 | }
|
|---|
| 1304 | else
|
|---|
| 1305 | {
|
|---|
| 1306 | # Yes, so we can speed it up slightly
|
|---|
| 1307 | while ($step <= $n)
|
|---|
| 1308 | {
|
|---|
| 1309 | # When the last element of $cx is 0, we split it up and remember
|
|---|
| 1310 | # how many we got so far. The reason is that n! will accumulate
|
|---|
| 1311 | # zeros at the end rather fast.
|
|---|
| 1312 | if ($cx->[0] == 0)
|
|---|
| 1313 | {
|
|---|
| 1314 | $zero_elements ++; shift @$cx;
|
|---|
| 1315 | }
|
|---|
| 1316 | _mul($c,$cx,[$step]); $step++;
|
|---|
| 1317 | }
|
|---|
| 1318 | }
|
|---|
| 1319 | # multiply in the zeros again
|
|---|
| 1320 | while ($zero_elements-- > 0)
|
|---|
| 1321 | {
|
|---|
| 1322 | unshift @$cx, 0;
|
|---|
| 1323 | }
|
|---|
| 1324 | $cx; # return result
|
|---|
| 1325 | }
|
|---|
| 1326 |
|
|---|
| 1327 | #############################################################################
|
|---|
| 1328 |
|
|---|
| 1329 | sub _log_int
|
|---|
| 1330 | {
|
|---|
| 1331 | # calculate integer log of $x to base $base
|
|---|
| 1332 | # ref to array, ref to array - return ref to array
|
|---|
| 1333 | my ($c,$x,$base) = @_;
|
|---|
| 1334 |
|
|---|
| 1335 | # X == 0 => NaN
|
|---|
| 1336 | return if (scalar @$x == 1 && $x->[0] == 0);
|
|---|
| 1337 | # BASE 0 or 1 => NaN
|
|---|
| 1338 | return if (scalar @$base == 1 && $base->[0] < 2);
|
|---|
| 1339 | my $cmp = _acmp($c,$x,$base); # X == BASE => 1
|
|---|
| 1340 | if ($cmp == 0)
|
|---|
| 1341 | {
|
|---|
| 1342 | splice (@$x,1); $x->[0] = 1;
|
|---|
| 1343 | return ($x,1)
|
|---|
| 1344 | }
|
|---|
| 1345 | # X < BASE
|
|---|
| 1346 | if ($cmp < 0)
|
|---|
| 1347 | {
|
|---|
| 1348 | splice (@$x,1); $x->[0] = 0;
|
|---|
| 1349 | return ($x,undef);
|
|---|
| 1350 | }
|
|---|
| 1351 |
|
|---|
| 1352 | # this trial multiplication is very fast, even for large counts (like for
|
|---|
| 1353 | # 2 ** 1024, since this still requires only 1024 very fast steps
|
|---|
| 1354 | # (multiplication of a large number by a very small number is very fast))
|
|---|
| 1355 | my $x_org = _copy($c,$x); # preserve x
|
|---|
| 1356 | splice(@$x,1); $x->[0] = 1; # keep ref to $x
|
|---|
| 1357 |
|
|---|
| 1358 | my $trial = _copy($c,$base);
|
|---|
| 1359 |
|
|---|
| 1360 | # XXX TODO this only works if $base has only one element
|
|---|
| 1361 | if (scalar @$base == 1)
|
|---|
| 1362 | {
|
|---|
| 1363 | # compute int ( length_in_base_10(X) / ( log(base) / log(10) ) )
|
|---|
| 1364 | my $len = _len($c,$x_org);
|
|---|
| 1365 | my $res = int($len / (log($base->[0]) / log(10))) || 1; # avoid $res == 0
|
|---|
| 1366 |
|
|---|
| 1367 | $x->[0] = $res;
|
|---|
| 1368 | $trial = _pow ($c, _copy($c, $base), $x);
|
|---|
| 1369 | my $a = _acmp($x,$trial,$x_org);
|
|---|
| 1370 | return ($x,1) if $a == 0;
|
|---|
| 1371 | # we now know that $res is too small
|
|---|
| 1372 | if ($res < 0)
|
|---|
| 1373 | {
|
|---|
| 1374 | _mul($c,$trial,$base); _add($c, $x, [1]);
|
|---|
| 1375 | }
|
|---|
| 1376 | else
|
|---|
| 1377 | {
|
|---|
| 1378 | # or too big
|
|---|
| 1379 | _div($c,$trial,$base); _sub($c, $x, [1]);
|
|---|
| 1380 | }
|
|---|
| 1381 | # did we now get the right result?
|
|---|
| 1382 | $a = _acmp($x,$trial,$x_org);
|
|---|
| 1383 | return ($x,1) if $a == 0; # yes, exactly
|
|---|
| 1384 | # still too big
|
|---|
| 1385 | if ($a > 0)
|
|---|
| 1386 | {
|
|---|
| 1387 | _div($c,$trial,$base); _sub($c, $x, [1]);
|
|---|
| 1388 | }
|
|---|
| 1389 | }
|
|---|
| 1390 |
|
|---|
| 1391 | # simple loop that increments $x by two in each step, possible overstepping
|
|---|
| 1392 | # the real result by one
|
|---|
| 1393 |
|
|---|
| 1394 | my $a;
|
|---|
| 1395 | my $base_mul = _mul($c, _copy($c,$base), $base);
|
|---|
| 1396 |
|
|---|
| 1397 | while (($a = _acmp($c,$trial,$x_org)) < 0)
|
|---|
| 1398 | {
|
|---|
| 1399 | _mul($c,$trial,$base_mul); _add($c, $x, [2]);
|
|---|
| 1400 | }
|
|---|
| 1401 |
|
|---|
| 1402 | my $exact = 1;
|
|---|
| 1403 | if ($a > 0)
|
|---|
| 1404 | {
|
|---|
| 1405 | # overstepped the result
|
|---|
| 1406 | _dec($c, $x);
|
|---|
| 1407 | _div($c,$trial,$base);
|
|---|
| 1408 | $a = _acmp($c,$trial,$x_org);
|
|---|
| 1409 | if ($a > 0)
|
|---|
| 1410 | {
|
|---|
| 1411 | _dec($c, $x);
|
|---|
| 1412 | }
|
|---|
| 1413 | $exact = 0 if $a != 0;
|
|---|
| 1414 | }
|
|---|
| 1415 |
|
|---|
| 1416 | ($x,$exact); # return result
|
|---|
| 1417 | }
|
|---|
| 1418 |
|
|---|
| 1419 | # for debugging:
|
|---|
| 1420 | use constant DEBUG => 0;
|
|---|
| 1421 | my $steps = 0;
|
|---|
| 1422 | sub steps { $steps };
|
|---|
| 1423 |
|
|---|
| 1424 | sub _sqrt
|
|---|
| 1425 | {
|
|---|
| 1426 | # square-root of $x in place
|
|---|
| 1427 | # Compute a guess of the result (by rule of thumb), then improve it via
|
|---|
| 1428 | # Newton's method.
|
|---|
| 1429 | my ($c,$x) = @_;
|
|---|
| 1430 |
|
|---|
| 1431 | if (scalar @$x == 1)
|
|---|
| 1432 | {
|
|---|
| 1433 | # fit's into one Perl scalar, so result can be computed directly
|
|---|
| 1434 | $x->[0] = int(sqrt($x->[0]));
|
|---|
| 1435 | return $x;
|
|---|
| 1436 | }
|
|---|
| 1437 | my $y = _copy($c,$x);
|
|---|
| 1438 | # hopefully _len/2 is < $BASE, the -1 is to always undershot the guess
|
|---|
| 1439 | # since our guess will "grow"
|
|---|
| 1440 | my $l = int((_len($c,$x)-1) / 2);
|
|---|
| 1441 |
|
|---|
| 1442 | my $lastelem = $x->[-1]; # for guess
|
|---|
| 1443 | my $elems = scalar @$x - 1;
|
|---|
| 1444 | # not enough digits, but could have more?
|
|---|
| 1445 | if ((length($lastelem) <= 3) && ($elems > 1))
|
|---|
| 1446 | {
|
|---|
| 1447 | # right-align with zero pad
|
|---|
| 1448 | my $len = length($lastelem) & 1;
|
|---|
| 1449 | print "$lastelem => " if DEBUG;
|
|---|
| 1450 | $lastelem .= substr($x->[-2] . '0' x $BASE_LEN,0,$BASE_LEN);
|
|---|
| 1451 | # former odd => make odd again, or former even to even again
|
|---|
| 1452 | $lastelem = $lastelem / 10 if (length($lastelem) & 1) != $len;
|
|---|
| 1453 | print "$lastelem\n" if DEBUG;
|
|---|
| 1454 | }
|
|---|
| 1455 |
|
|---|
| 1456 | # construct $x (instead of _lsft($c,$x,$l,10)
|
|---|
| 1457 | my $r = $l % $BASE_LEN; # 10000 00000 00000 00000 ($BASE_LEN=5)
|
|---|
| 1458 | $l = int($l / $BASE_LEN);
|
|---|
| 1459 | print "l = $l " if DEBUG;
|
|---|
| 1460 |
|
|---|
| 1461 | splice @$x,$l; # keep ref($x), but modify it
|
|---|
| 1462 |
|
|---|
| 1463 | # we make the first part of the guess not '1000...0' but int(sqrt($lastelem))
|
|---|
| 1464 | # that gives us:
|
|---|
| 1465 | # 14400 00000 => sqrt(14400) => guess first digits to be 120
|
|---|
| 1466 | # 144000 000000 => sqrt(144000) => guess 379
|
|---|
| 1467 |
|
|---|
| 1468 | print "$lastelem (elems $elems) => " if DEBUG;
|
|---|
| 1469 | $lastelem = $lastelem / 10 if ($elems & 1 == 1); # odd or even?
|
|---|
| 1470 | my $g = sqrt($lastelem); $g =~ s/\.//; # 2.345 => 2345
|
|---|
| 1471 | $r -= 1 if $elems & 1 == 0; # 70 => 7
|
|---|
| 1472 |
|
|---|
| 1473 | # padd with zeros if result is too short
|
|---|
| 1474 | $x->[$l--] = int(substr($g . '0' x $r,0,$r+1));
|
|---|
| 1475 | print "now ",$x->[-1] if DEBUG;
|
|---|
| 1476 | print " would have been ", int('1' . '0' x $r),"\n" if DEBUG;
|
|---|
| 1477 |
|
|---|
| 1478 | # If @$x > 1, we could compute the second elem of the guess, too, to create
|
|---|
| 1479 | # an even better guess. Not implemented yet. Does it improve performance?
|
|---|
| 1480 | $x->[$l--] = 0 while ($l >= 0); # all other digits of guess are zero
|
|---|
| 1481 |
|
|---|
| 1482 | print "start x= ",_str($c,$x),"\n" if DEBUG;
|
|---|
| 1483 | my $two = _two();
|
|---|
| 1484 | my $last = _zero();
|
|---|
| 1485 | my $lastlast = _zero();
|
|---|
| 1486 | $steps = 0 if DEBUG;
|
|---|
| 1487 | while (_acmp($c,$last,$x) != 0 && _acmp($c,$lastlast,$x) != 0)
|
|---|
| 1488 | {
|
|---|
| 1489 | $steps++ if DEBUG;
|
|---|
| 1490 | $lastlast = _copy($c,$last);
|
|---|
| 1491 | $last = _copy($c,$x);
|
|---|
| 1492 | _add($c,$x, _div($c,_copy($c,$y),$x));
|
|---|
| 1493 | _div($c,$x, $two );
|
|---|
| 1494 | print " x= ",_str($c,$x),"\n" if DEBUG;
|
|---|
| 1495 | }
|
|---|
| 1496 | print "\nsteps in sqrt: $steps, " if DEBUG;
|
|---|
| 1497 | _dec($c,$x) if _acmp($c,$y,_mul($c,_copy($c,$x),$x)) < 0; # overshot?
|
|---|
| 1498 | print " final ",$x->[-1],"\n" if DEBUG;
|
|---|
| 1499 | $x;
|
|---|
| 1500 | }
|
|---|
| 1501 |
|
|---|
| 1502 | sub _root
|
|---|
| 1503 | {
|
|---|
| 1504 | # take n'th root of $x in place (n >= 3)
|
|---|
| 1505 | my ($c,$x,$n) = @_;
|
|---|
| 1506 |
|
|---|
| 1507 | if (scalar @$x == 1)
|
|---|
| 1508 | {
|
|---|
| 1509 | if (scalar @$n > 1)
|
|---|
| 1510 | {
|
|---|
| 1511 | # result will always be smaller than 2 so trunc to 1 at once
|
|---|
| 1512 | $x->[0] = 1;
|
|---|
| 1513 | }
|
|---|
| 1514 | else
|
|---|
| 1515 | {
|
|---|
| 1516 | # fit's into one Perl scalar, so result can be computed directly
|
|---|
| 1517 | # cannot use int() here, because it rounds wrongly (try
|
|---|
| 1518 | # (81 ** 3) ** (1/3) to see what I mean)
|
|---|
| 1519 | #$x->[0] = int( $x->[0] ** (1 / $n->[0]) );
|
|---|
| 1520 | # round to 8 digits, then truncate result to integer
|
|---|
| 1521 | $x->[0] = int ( sprintf ("%.8f", $x->[0] ** (1 / $n->[0]) ) );
|
|---|
| 1522 | }
|
|---|
| 1523 | return $x;
|
|---|
| 1524 | }
|
|---|
| 1525 |
|
|---|
| 1526 | # we know now that X is more than one element long
|
|---|
| 1527 |
|
|---|
| 1528 | # if $n is a power of two, we can repeatedly take sqrt($X) and find the
|
|---|
| 1529 | # proper result, because sqrt(sqrt($x)) == root($x,4)
|
|---|
| 1530 | my $b = _as_bin($c,$n);
|
|---|
| 1531 | if ($b =~ /0b1(0+)$/)
|
|---|
| 1532 | {
|
|---|
| 1533 | my $count = CORE::length($1); # 0b100 => len('00') => 2
|
|---|
| 1534 | my $cnt = $count; # counter for loop
|
|---|
| 1535 | unshift (@$x, 0); # add one element, together with one
|
|---|
| 1536 | # more below in the loop this makes 2
|
|---|
| 1537 | while ($cnt-- > 0)
|
|---|
| 1538 | {
|
|---|
| 1539 | # 'inflate' $X by adding one element, basically computing
|
|---|
| 1540 | # $x * $BASE * $BASE. This gives us more $BASE_LEN digits for result
|
|---|
| 1541 | # since len(sqrt($X)) approx == len($x) / 2.
|
|---|
| 1542 | unshift (@$x, 0);
|
|---|
| 1543 | # calculate sqrt($x), $x is now one element to big, again. In the next
|
|---|
| 1544 | # round we make that two, again.
|
|---|
| 1545 | _sqrt($c,$x);
|
|---|
| 1546 | }
|
|---|
| 1547 | # $x is now one element to big, so truncate result by removing it
|
|---|
| 1548 | splice (@$x,0,1);
|
|---|
| 1549 | }
|
|---|
| 1550 | else
|
|---|
| 1551 | {
|
|---|
| 1552 | # trial computation by starting with 2,4,8,16 etc until we overstep
|
|---|
| 1553 | my $step;
|
|---|
| 1554 | my $trial = _two();
|
|---|
| 1555 |
|
|---|
| 1556 | # while still to do more than X steps
|
|---|
| 1557 | do
|
|---|
| 1558 | {
|
|---|
| 1559 | $step = _two();
|
|---|
| 1560 | while (_acmp($c, _pow($c, _copy($c, $trial), $n), $x) < 0)
|
|---|
| 1561 | {
|
|---|
| 1562 | _mul ($c, $step, [2]);
|
|---|
| 1563 | _add ($c, $trial, $step);
|
|---|
| 1564 | }
|
|---|
| 1565 |
|
|---|
| 1566 | # hit exactly?
|
|---|
| 1567 | if (_acmp($c, _pow($c, _copy($c, $trial), $n), $x) == 0)
|
|---|
| 1568 | {
|
|---|
| 1569 | @$x = @$trial; # make copy while preserving ref to $x
|
|---|
| 1570 | return $x;
|
|---|
| 1571 | }
|
|---|
| 1572 | # overstepped, so go back on step
|
|---|
| 1573 | _sub($c, $trial, $step);
|
|---|
| 1574 | } while (scalar @$step > 1 || $step->[0] > 128);
|
|---|
| 1575 |
|
|---|
| 1576 | # reset step to 2
|
|---|
| 1577 | $step = _two();
|
|---|
| 1578 | # add two, because $trial cannot be exactly the result (otherwise we would
|
|---|
| 1579 | # alrady have found it)
|
|---|
| 1580 | _add($c, $trial, $step);
|
|---|
| 1581 |
|
|---|
| 1582 | # and now add more and more (2,4,6,8,10 etc)
|
|---|
| 1583 | while (_acmp($c, _pow($c, _copy($c, $trial), $n), $x) < 0)
|
|---|
| 1584 | {
|
|---|
| 1585 | _add ($c, $trial, $step);
|
|---|
| 1586 | }
|
|---|
| 1587 |
|
|---|
| 1588 | # hit not exactly? (overstepped)
|
|---|
| 1589 | if (_acmp($c, _pow($c, _copy($c, $trial), $n), $x) > 0)
|
|---|
| 1590 | {
|
|---|
| 1591 | _dec($c,$trial);
|
|---|
| 1592 | }
|
|---|
| 1593 |
|
|---|
| 1594 | # hit not exactly? (overstepped)
|
|---|
| 1595 | # 80 too small, 81 slightly too big, 82 too big
|
|---|
| 1596 | if (_acmp($c, _pow($c, _copy($c, $trial), $n), $x) > 0)
|
|---|
| 1597 | {
|
|---|
| 1598 | _dec ($c, $trial);
|
|---|
| 1599 | }
|
|---|
| 1600 |
|
|---|
| 1601 | @$x = @$trial; # make copy while preserving ref to $x
|
|---|
| 1602 | return $x;
|
|---|
| 1603 | }
|
|---|
| 1604 | $x;
|
|---|
| 1605 | }
|
|---|
| 1606 |
|
|---|
| 1607 | ##############################################################################
|
|---|
| 1608 | # binary stuff
|
|---|
| 1609 |
|
|---|
| 1610 | sub _and
|
|---|
| 1611 | {
|
|---|
| 1612 | my ($c,$x,$y) = @_;
|
|---|
| 1613 |
|
|---|
| 1614 | # the shortcut makes equal, large numbers _really_ fast, and makes only a
|
|---|
| 1615 | # very small performance drop for small numbers (e.g. something with less
|
|---|
| 1616 | # than 32 bit) Since we optimize for large numbers, this is enabled.
|
|---|
| 1617 | return $x if _acmp($c,$x,$y) == 0; # shortcut
|
|---|
| 1618 |
|
|---|
| 1619 | my $m = _one(); my ($xr,$yr);
|
|---|
| 1620 | my $mask = $AND_MASK;
|
|---|
| 1621 |
|
|---|
| 1622 | my $x1 = $x;
|
|---|
| 1623 | my $y1 = _copy($c,$y); # make copy
|
|---|
| 1624 | $x = _zero();
|
|---|
| 1625 | my ($b,$xrr,$yrr);
|
|---|
| 1626 | use integer;
|
|---|
| 1627 | while (!_is_zero($c,$x1) && !_is_zero($c,$y1))
|
|---|
| 1628 | {
|
|---|
| 1629 | ($x1, $xr) = _div($c,$x1,$mask);
|
|---|
| 1630 | ($y1, $yr) = _div($c,$y1,$mask);
|
|---|
| 1631 |
|
|---|
| 1632 | # make ints() from $xr, $yr
|
|---|
| 1633 | # this is when the AND_BITS are greater than $BASE and is slower for
|
|---|
| 1634 | # small (<256 bits) numbers, but faster for large numbers. Disabled
|
|---|
| 1635 | # due to KISS principle
|
|---|
| 1636 |
|
|---|
| 1637 | # $b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
|
|---|
| 1638 | # $b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
|
|---|
| 1639 | # _add($c,$x, _mul($c, _new( $c, ($xrr & $yrr) ), $m) );
|
|---|
| 1640 |
|
|---|
| 1641 | # 0+ due to '&' doesn't work in strings
|
|---|
| 1642 | _add($c,$x, _mul($c, [ 0+$xr->[0] & 0+$yr->[0] ], $m) );
|
|---|
| 1643 | _mul($c,$m,$mask);
|
|---|
| 1644 | }
|
|---|
| 1645 | $x;
|
|---|
| 1646 | }
|
|---|
| 1647 |
|
|---|
| 1648 | sub _xor
|
|---|
| 1649 | {
|
|---|
| 1650 | my ($c,$x,$y) = @_;
|
|---|
| 1651 |
|
|---|
| 1652 | return _zero() if _acmp($c,$x,$y) == 0; # shortcut (see -and)
|
|---|
| 1653 |
|
|---|
| 1654 | my $m = _one(); my ($xr,$yr);
|
|---|
| 1655 | my $mask = $XOR_MASK;
|
|---|
| 1656 |
|
|---|
| 1657 | my $x1 = $x;
|
|---|
| 1658 | my $y1 = _copy($c,$y); # make copy
|
|---|
| 1659 | $x = _zero();
|
|---|
| 1660 | my ($b,$xrr,$yrr);
|
|---|
| 1661 | use integer;
|
|---|
| 1662 | while (!_is_zero($c,$x1) && !_is_zero($c,$y1))
|
|---|
| 1663 | {
|
|---|
| 1664 | ($x1, $xr) = _div($c,$x1,$mask);
|
|---|
| 1665 | ($y1, $yr) = _div($c,$y1,$mask);
|
|---|
| 1666 | # make ints() from $xr, $yr (see _and())
|
|---|
| 1667 | #$b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
|
|---|
| 1668 | #$b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
|
|---|
| 1669 | #_add($c,$x, _mul($c, _new( $c, ($xrr ^ $yrr) ), $m) );
|
|---|
| 1670 |
|
|---|
| 1671 | # 0+ due to '^' doesn't work in strings
|
|---|
| 1672 | _add($c,$x, _mul($c, [ 0+$xr->[0] ^ 0+$yr->[0] ], $m) );
|
|---|
| 1673 | _mul($c,$m,$mask);
|
|---|
| 1674 | }
|
|---|
| 1675 | # the loop stops when the shorter of the two numbers is exhausted
|
|---|
| 1676 | # the remainder of the longer one will survive bit-by-bit, so we simple
|
|---|
| 1677 | # multiply-add it in
|
|---|
| 1678 | _add($c,$x, _mul($c, $x1, $m) ) if !_is_zero($c,$x1);
|
|---|
| 1679 | _add($c,$x, _mul($c, $y1, $m) ) if !_is_zero($c,$y1);
|
|---|
| 1680 |
|
|---|
| 1681 | $x;
|
|---|
| 1682 | }
|
|---|
| 1683 |
|
|---|
| 1684 | sub _or
|
|---|
| 1685 | {
|
|---|
| 1686 | my ($c,$x,$y) = @_;
|
|---|
| 1687 |
|
|---|
| 1688 | return $x if _acmp($c,$x,$y) == 0; # shortcut (see _and)
|
|---|
| 1689 |
|
|---|
| 1690 | my $m = _one(); my ($xr,$yr);
|
|---|
| 1691 | my $mask = $OR_MASK;
|
|---|
| 1692 |
|
|---|
| 1693 | my $x1 = $x;
|
|---|
| 1694 | my $y1 = _copy($c,$y); # make copy
|
|---|
| 1695 | $x = _zero();
|
|---|
| 1696 | my ($b,$xrr,$yrr);
|
|---|
| 1697 | use integer;
|
|---|
| 1698 | while (!_is_zero($c,$x1) && !_is_zero($c,$y1))
|
|---|
| 1699 | {
|
|---|
| 1700 | ($x1, $xr) = _div($c,$x1,$mask);
|
|---|
| 1701 | ($y1, $yr) = _div($c,$y1,$mask);
|
|---|
| 1702 | # make ints() from $xr, $yr (see _and())
|
|---|
| 1703 | # $b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
|
|---|
| 1704 | # $b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
|
|---|
| 1705 | # _add($c,$x, _mul($c, _new( $c, ($xrr | $yrr) ), $m) );
|
|---|
| 1706 |
|
|---|
| 1707 | # 0+ due to '|' doesn't work in strings
|
|---|
| 1708 | _add($c,$x, _mul($c, [ 0+$xr->[0] | 0+$yr->[0] ], $m) );
|
|---|
| 1709 | _mul($c,$m,$mask);
|
|---|
| 1710 | }
|
|---|
| 1711 | # the loop stops when the shorter of the two numbers is exhausted
|
|---|
| 1712 | # the remainder of the longer one will survive bit-by-bit, so we simple
|
|---|
| 1713 | # multiply-add it in
|
|---|
| 1714 | _add($c,$x, _mul($c, $x1, $m) ) if !_is_zero($c,$x1);
|
|---|
| 1715 | _add($c,$x, _mul($c, $y1, $m) ) if !_is_zero($c,$y1);
|
|---|
| 1716 |
|
|---|
| 1717 | $x;
|
|---|
| 1718 | }
|
|---|
| 1719 |
|
|---|
| 1720 | sub _as_hex
|
|---|
| 1721 | {
|
|---|
| 1722 | # convert a decimal number to hex (ref to array, return ref to string)
|
|---|
| 1723 | my ($c,$x) = @_;
|
|---|
| 1724 |
|
|---|
| 1725 | # fit's into one element (handle also 0x0 case)
|
|---|
| 1726 | return sprintf("0x%x",$x->[0]) if @$x == 1;
|
|---|
| 1727 |
|
|---|
| 1728 | my $x1 = _copy($c,$x);
|
|---|
| 1729 |
|
|---|
| 1730 | my $es = '';
|
|---|
| 1731 | my ($xr, $h, $x10000);
|
|---|
| 1732 | if ($] >= 5.006)
|
|---|
| 1733 | {
|
|---|
| 1734 | $x10000 = [ 0x10000 ]; $h = 'h4';
|
|---|
| 1735 | }
|
|---|
| 1736 | else
|
|---|
| 1737 | {
|
|---|
| 1738 | $x10000 = [ 0x1000 ]; $h = 'h3';
|
|---|
| 1739 | }
|
|---|
| 1740 | while (@$x1 != 1 || $x1->[0] != 0) # _is_zero()
|
|---|
| 1741 | {
|
|---|
| 1742 | ($x1, $xr) = _div($c,$x1,$x10000);
|
|---|
| 1743 | $es .= unpack($h,pack('v',$xr->[0])); # XXX TODO: why pack('v',...)?
|
|---|
| 1744 | }
|
|---|
| 1745 | $es = reverse $es;
|
|---|
| 1746 | $es =~ s/^[0]+//; # strip leading zeros
|
|---|
| 1747 | '0x' . $es; # return result prepended with 0x
|
|---|
| 1748 | }
|
|---|
| 1749 |
|
|---|
| 1750 | sub _as_bin
|
|---|
| 1751 | {
|
|---|
| 1752 | # convert a decimal number to bin (ref to array, return ref to string)
|
|---|
| 1753 | my ($c,$x) = @_;
|
|---|
| 1754 |
|
|---|
| 1755 | # fit's into one element (and Perl recent enough), handle also 0b0 case
|
|---|
| 1756 | # handle zero case for older Perls
|
|---|
| 1757 | if ($] <= 5.005 && @$x == 1 && $x->[0] == 0)
|
|---|
| 1758 | {
|
|---|
| 1759 | my $t = '0b0'; return $t;
|
|---|
| 1760 | }
|
|---|
| 1761 | if (@$x == 1 && $] >= 5.006)
|
|---|
| 1762 | {
|
|---|
| 1763 | my $t = sprintf("0b%b",$x->[0]);
|
|---|
| 1764 | return $t;
|
|---|
| 1765 | }
|
|---|
| 1766 | my $x1 = _copy($c,$x);
|
|---|
| 1767 |
|
|---|
| 1768 | my $es = '';
|
|---|
| 1769 | my ($xr, $b, $x10000);
|
|---|
| 1770 | if ($] >= 5.006)
|
|---|
| 1771 | {
|
|---|
| 1772 | $x10000 = [ 0x10000 ]; $b = 'b16';
|
|---|
| 1773 | }
|
|---|
| 1774 | else
|
|---|
| 1775 | {
|
|---|
| 1776 | $x10000 = [ 0x1000 ]; $b = 'b12';
|
|---|
| 1777 | }
|
|---|
| 1778 | while (!(@$x1 == 1 && $x1->[0] == 0)) # _is_zero()
|
|---|
| 1779 | {
|
|---|
| 1780 | ($x1, $xr) = _div($c,$x1,$x10000);
|
|---|
| 1781 | $es .= unpack($b,pack('v',$xr->[0])); # XXX TODO: why pack('v',...)?
|
|---|
| 1782 | # $es .= unpack($b,$xr->[0]);
|
|---|
| 1783 | }
|
|---|
| 1784 | $es = reverse $es;
|
|---|
| 1785 | $es =~ s/^[0]+//; # strip leading zeros
|
|---|
| 1786 | '0b' . $es; # return result prepended with 0b
|
|---|
| 1787 | }
|
|---|
| 1788 |
|
|---|
| 1789 | sub _from_hex
|
|---|
| 1790 | {
|
|---|
| 1791 | # convert a hex number to decimal (ref to string, return ref to array)
|
|---|
| 1792 | my ($c,$hs) = @_;
|
|---|
| 1793 |
|
|---|
| 1794 | my $m = _new($c, 0x10000000); # 28 bit at a time (<32 bit!)
|
|---|
| 1795 | my $d = 7; # 7 digits at a time
|
|---|
| 1796 | if ($] <= 5.006)
|
|---|
| 1797 | {
|
|---|
| 1798 | # for older Perls, play safe
|
|---|
| 1799 | $m = [ 0x10000 ]; # 16 bit at a time (<32 bit!)
|
|---|
| 1800 | $d = 4; # 4 digits at a time
|
|---|
| 1801 | }
|
|---|
| 1802 |
|
|---|
| 1803 | my $mul = _one();
|
|---|
| 1804 | my $x = _zero();
|
|---|
| 1805 |
|
|---|
| 1806 | my $len = int( (length($hs)-2)/$d ); # $d digit parts, w/o the '0x'
|
|---|
| 1807 | my $val; my $i = -$d;
|
|---|
| 1808 | while ($len >= 0)
|
|---|
| 1809 | {
|
|---|
| 1810 | $val = substr($hs,$i,$d); # get hex digits
|
|---|
| 1811 | $val =~ s/^[+-]?0x// if $len == 0; # for last part only because
|
|---|
| 1812 | $val = hex($val); # hex does not like wrong chars
|
|---|
| 1813 | $i -= $d; $len --;
|
|---|
| 1814 | my $adder = [ $val ];
|
|---|
| 1815 | # if the resulting number was to big to fit into one element, create a
|
|---|
| 1816 | # two-element version (bug found by Mark Lakata - Thanx!)
|
|---|
| 1817 | if (CORE::length($val) > $BASE_LEN)
|
|---|
| 1818 | {
|
|---|
| 1819 | $adder = _new($c,$val);
|
|---|
| 1820 | }
|
|---|
| 1821 | _add ($c, $x, _mul ($c, $adder, $mul ) ) if $val != 0;
|
|---|
| 1822 | _mul ($c, $mul, $m ) if $len >= 0; # skip last mul
|
|---|
| 1823 | }
|
|---|
| 1824 | $x;
|
|---|
| 1825 | }
|
|---|
| 1826 |
|
|---|
| 1827 | sub _from_bin
|
|---|
| 1828 | {
|
|---|
| 1829 | # convert a hex number to decimal (ref to string, return ref to array)
|
|---|
| 1830 | my ($c,$bs) = @_;
|
|---|
| 1831 |
|
|---|
| 1832 | # instead of converting X (8) bit at a time, it is faster to "convert" the
|
|---|
| 1833 | # number to hex, and then call _from_hex.
|
|---|
| 1834 |
|
|---|
| 1835 | my $hs = $bs;
|
|---|
| 1836 | $hs =~ s/^[+-]?0b//; # remove sign and 0b
|
|---|
| 1837 | my $l = length($hs); # bits
|
|---|
| 1838 | $hs = '0' x (8-($l % 8)) . $hs if ($l % 8) != 0; # padd left side w/ 0
|
|---|
| 1839 | my $h = '0x' . unpack('H*', pack ('B*', $hs)); # repack as hex
|
|---|
| 1840 |
|
|---|
| 1841 | $c->_from_hex($h);
|
|---|
| 1842 | }
|
|---|
| 1843 |
|
|---|
| 1844 | ##############################################################################
|
|---|
| 1845 | # special modulus functions
|
|---|
| 1846 |
|
|---|
| 1847 | sub _modinv
|
|---|
| 1848 | {
|
|---|
| 1849 | # modular inverse
|
|---|
| 1850 | my ($c,$x,$y) = @_;
|
|---|
| 1851 |
|
|---|
| 1852 | my $u = _zero($c); my $u1 = _one($c);
|
|---|
| 1853 | my $a = _copy($c,$y); my $b = _copy($c,$x);
|
|---|
| 1854 |
|
|---|
| 1855 | # Euclid's Algorithm for bgcd(), only that we calc bgcd() ($a) and the
|
|---|
| 1856 | # result ($u) at the same time. See comments in BigInt for why this works.
|
|---|
| 1857 | my $q;
|
|---|
| 1858 | ($a, $q, $b) = ($b, _div($c,$a,$b)); # step 1
|
|---|
| 1859 | my $sign = 1;
|
|---|
| 1860 | while (!_is_zero($c,$b))
|
|---|
| 1861 | {
|
|---|
| 1862 | my $t = _add($c, # step 2:
|
|---|
| 1863 | _mul($c,_copy($c,$u1), $q) , # t = u1 * q
|
|---|
| 1864 | $u ); # + u
|
|---|
| 1865 | $u = $u1; # u = u1, u1 = t
|
|---|
| 1866 | $u1 = $t;
|
|---|
| 1867 | $sign = -$sign;
|
|---|
| 1868 | ($a, $q, $b) = ($b, _div($c,$a,$b)); # step 1
|
|---|
| 1869 | }
|
|---|
| 1870 |
|
|---|
| 1871 | # if the gcd is not 1, then return NaN
|
|---|
| 1872 | return (undef,undef) unless _is_one($c,$a);
|
|---|
| 1873 |
|
|---|
| 1874 | ($u1, $sign == 1 ? '+' : '-');
|
|---|
| 1875 | }
|
|---|
| 1876 |
|
|---|
| 1877 | sub _modpow
|
|---|
| 1878 | {
|
|---|
| 1879 | # modulus of power ($x ** $y) % $z
|
|---|
| 1880 | my ($c,$num,$exp,$mod) = @_;
|
|---|
| 1881 |
|
|---|
| 1882 | # in the trivial case,
|
|---|
| 1883 | if (_is_one($c,$mod))
|
|---|
| 1884 | {
|
|---|
| 1885 | splice @$num,0,1; $num->[0] = 0;
|
|---|
| 1886 | return $num;
|
|---|
| 1887 | }
|
|---|
| 1888 | if ((scalar @$num == 1) && (($num->[0] == 0) || ($num->[0] == 1)))
|
|---|
| 1889 | {
|
|---|
| 1890 | $num->[0] = 1;
|
|---|
| 1891 | return $num;
|
|---|
| 1892 | }
|
|---|
| 1893 |
|
|---|
| 1894 | # $num = _mod($c,$num,$mod); # this does not make it faster
|
|---|
| 1895 |
|
|---|
| 1896 | my $acc = _copy($c,$num); my $t = _one();
|
|---|
| 1897 |
|
|---|
| 1898 | my $expbin = _as_bin($c,$exp); $expbin =~ s/^0b//;
|
|---|
| 1899 | my $len = length($expbin);
|
|---|
| 1900 | while (--$len >= 0)
|
|---|
| 1901 | {
|
|---|
| 1902 | if ( substr($expbin,$len,1) eq '1') # is_odd
|
|---|
| 1903 | {
|
|---|
| 1904 | _mul($c,$t,$acc);
|
|---|
| 1905 | $t = _mod($c,$t,$mod);
|
|---|
| 1906 | }
|
|---|
| 1907 | _mul($c,$acc,$acc);
|
|---|
| 1908 | $acc = _mod($c,$acc,$mod);
|
|---|
| 1909 | }
|
|---|
| 1910 | @$num = @$t;
|
|---|
| 1911 | $num;
|
|---|
| 1912 | }
|
|---|
| 1913 |
|
|---|
| 1914 | sub _gcd
|
|---|
| 1915 | {
|
|---|
| 1916 | # greatest common divisor
|
|---|
| 1917 | my ($c,$x,$y) = @_;
|
|---|
| 1918 |
|
|---|
| 1919 | while ( (scalar @$y != 1) || ($y->[0] != 0) ) # while ($y != 0)
|
|---|
| 1920 | {
|
|---|
| 1921 | my $t = _copy($c,$y);
|
|---|
| 1922 | $y = _mod($c, $x, $y);
|
|---|
| 1923 | $x = $t;
|
|---|
| 1924 | }
|
|---|
| 1925 | $x;
|
|---|
| 1926 | }
|
|---|
| 1927 |
|
|---|
| 1928 | ##############################################################################
|
|---|
| 1929 | ##############################################################################
|
|---|
| 1930 |
|
|---|
| 1931 | 1;
|
|---|
| 1932 | __END__
|
|---|
| 1933 |
|
|---|
| 1934 | =head1 NAME
|
|---|
| 1935 |
|
|---|
| 1936 | Math::BigInt::Calc - Pure Perl module to support Math::BigInt
|
|---|
| 1937 |
|
|---|
| 1938 | =head1 SYNOPSIS
|
|---|
| 1939 |
|
|---|
| 1940 | Provides support for big integer calculations. Not intended to be used by other
|
|---|
| 1941 | modules. Other modules which sport the same functions can also be used to support
|
|---|
| 1942 | Math::BigInt, like Math::BigInt::GMP or Math::BigInt::Pari.
|
|---|
| 1943 |
|
|---|
| 1944 | =head1 DESCRIPTION
|
|---|
| 1945 |
|
|---|
| 1946 | In order to allow for multiple big integer libraries, Math::BigInt was
|
|---|
| 1947 | rewritten to use library modules for core math routines. Any module which
|
|---|
| 1948 | follows the same API as this can be used instead by using the following:
|
|---|
| 1949 |
|
|---|
| 1950 | use Math::BigInt lib => 'libname';
|
|---|
| 1951 |
|
|---|
| 1952 | 'libname' is either the long name ('Math::BigInt::Pari'), or only the short
|
|---|
| 1953 | version like 'Pari'.
|
|---|
| 1954 |
|
|---|
| 1955 | =head1 STORAGE
|
|---|
| 1956 |
|
|---|
| 1957 | =head1 METHODS
|
|---|
| 1958 |
|
|---|
| 1959 | The following functions MUST be defined in order to support the use by
|
|---|
| 1960 | Math::BigInt v1.70 or later:
|
|---|
| 1961 |
|
|---|
| 1962 | api_version() return API version, minimum 1 for v1.70
|
|---|
| 1963 | _new(string) return ref to new object from ref to decimal string
|
|---|
| 1964 | _zero() return a new object with value 0
|
|---|
| 1965 | _one() return a new object with value 1
|
|---|
| 1966 | _two() return a new object with value 2
|
|---|
| 1967 | _ten() return a new object with value 10
|
|---|
| 1968 |
|
|---|
| 1969 | _str(obj) return ref to a string representing the object
|
|---|
| 1970 | _num(obj) returns a Perl integer/floating point number
|
|---|
| 1971 | NOTE: because of Perl numeric notation defaults,
|
|---|
| 1972 | the _num'ified obj may lose accuracy due to
|
|---|
| 1973 | machine-dependend floating point size limitations
|
|---|
| 1974 |
|
|---|
| 1975 | _add(obj,obj) Simple addition of two objects
|
|---|
| 1976 | _mul(obj,obj) Multiplication of two objects
|
|---|
| 1977 | _div(obj,obj) Division of the 1st object by the 2nd
|
|---|
| 1978 | In list context, returns (result,remainder).
|
|---|
| 1979 | NOTE: this is integer math, so no
|
|---|
| 1980 | fractional part will be returned.
|
|---|
| 1981 | The second operand will be not be 0, so no need to
|
|---|
| 1982 | check for that.
|
|---|
| 1983 | _sub(obj,obj) Simple subtraction of 1 object from another
|
|---|
| 1984 | a third, optional parameter indicates that the params
|
|---|
| 1985 | are swapped. In this case, the first param needs to
|
|---|
| 1986 | be preserved, while you can destroy the second.
|
|---|
| 1987 | sub (x,y,1) => return x - y and keep x intact!
|
|---|
| 1988 | _dec(obj) decrement object by one (input is garant. to be > 0)
|
|---|
| 1989 | _inc(obj) increment object by one
|
|---|
| 1990 |
|
|---|
| 1991 |
|
|---|
| 1992 | _acmp(obj,obj) <=> operator for objects (return -1, 0 or 1)
|
|---|
| 1993 |
|
|---|
| 1994 | _len(obj) returns count of the decimal digits of the object
|
|---|
| 1995 | _digit(obj,n) returns the n'th decimal digit of object
|
|---|
| 1996 |
|
|---|
| 1997 | _is_one(obj) return true if argument is 1
|
|---|
| 1998 | _is_two(obj) return true if argument is 2
|
|---|
| 1999 | _is_ten(obj) return true if argument is 10
|
|---|
| 2000 | _is_zero(obj) return true if argument is 0
|
|---|
| 2001 | _is_even(obj) return true if argument is even (0,2,4,6..)
|
|---|
| 2002 | _is_odd(obj) return true if argument is odd (1,3,5,7..)
|
|---|
| 2003 |
|
|---|
| 2004 | _copy return a ref to a true copy of the object
|
|---|
| 2005 |
|
|---|
| 2006 | _check(obj) check whether internal representation is still intact
|
|---|
| 2007 | return 0 for ok, otherwise error message as string
|
|---|
| 2008 |
|
|---|
| 2009 | _from_hex(str) return ref to new object from ref to hexadecimal string
|
|---|
| 2010 | _from_bin(str) return ref to new object from ref to binary string
|
|---|
| 2011 |
|
|---|
| 2012 | _as_hex(str) return string containing the value as
|
|---|
| 2013 | unsigned hex string, with the '0x' prepended.
|
|---|
| 2014 | Leading zeros must be stripped.
|
|---|
| 2015 | _as_bin(str) Like as_hex, only as binary string containing only
|
|---|
| 2016 | zeros and ones. Leading zeros must be stripped and a
|
|---|
| 2017 | '0b' must be prepended.
|
|---|
| 2018 |
|
|---|
| 2019 | _rsft(obj,N,B) shift object in base B by N 'digits' right
|
|---|
| 2020 | _lsft(obj,N,B) shift object in base B by N 'digits' left
|
|---|
| 2021 |
|
|---|
| 2022 | _xor(obj1,obj2) XOR (bit-wise) object 1 with object 2
|
|---|
| 2023 | Note: XOR, AND and OR pad with zeros if size mismatches
|
|---|
| 2024 | _and(obj1,obj2) AND (bit-wise) object 1 with object 2
|
|---|
| 2025 | _or(obj1,obj2) OR (bit-wise) object 1 with object 2
|
|---|
| 2026 |
|
|---|
| 2027 | _mod(obj,obj) Return remainder of div of the 1st by the 2nd object
|
|---|
| 2028 | _sqrt(obj) return the square root of object (truncated to int)
|
|---|
| 2029 | _root(obj) return the n'th (n >= 3) root of obj (truncated to int)
|
|---|
| 2030 | _fac(obj) return factorial of object 1 (1*2*3*4..)
|
|---|
| 2031 | _pow(obj,obj) return object 1 to the power of object 2
|
|---|
| 2032 | return undef for NaN
|
|---|
| 2033 | _zeros(obj) return number of trailing decimal zeros
|
|---|
| 2034 | _modinv return inverse modulus
|
|---|
| 2035 | _modpow return modulus of power ($x ** $y) % $z
|
|---|
| 2036 | _log_int(X,N) calculate integer log() of X in base N
|
|---|
| 2037 | X >= 0, N >= 0 (return undef for NaN)
|
|---|
| 2038 | returns (RESULT, EXACT) where EXACT is:
|
|---|
| 2039 | 1 : result is exactly RESULT
|
|---|
| 2040 | 0 : result was truncated to RESULT
|
|---|
| 2041 | undef : unknown whether result is exactly RESULT
|
|---|
| 2042 | _gcd(obj,obj) return Greatest Common Divisor of two objects
|
|---|
| 2043 |
|
|---|
| 2044 | The following functions are optional, and can be defined if the underlying lib
|
|---|
| 2045 | has a fast way to do them. If undefined, Math::BigInt will use pure Perl (hence
|
|---|
| 2046 | slow) fallback routines to emulate these:
|
|---|
| 2047 |
|
|---|
| 2048 | _signed_or
|
|---|
| 2049 | _signed_and
|
|---|
| 2050 | _signed_xor
|
|---|
| 2051 |
|
|---|
| 2052 |
|
|---|
| 2053 | Input strings come in as unsigned but with prefix (i.e. as '123', '0xabc'
|
|---|
| 2054 | or '0b1101').
|
|---|
| 2055 |
|
|---|
| 2056 | So the library needs only to deal with unsigned big integers. Testing of input
|
|---|
| 2057 | parameter validity is done by the caller, so you need not worry about
|
|---|
| 2058 | underflow (f.i. in C<_sub()>, C<_dec()>) nor about division by zero or similar
|
|---|
| 2059 | cases.
|
|---|
| 2060 |
|
|---|
| 2061 | The first parameter can be modified, that includes the possibility that you
|
|---|
| 2062 | return a reference to a completely different object instead. Although keeping
|
|---|
| 2063 | the reference and just changing it's contents is prefered over creating and
|
|---|
| 2064 | returning a different reference.
|
|---|
| 2065 |
|
|---|
| 2066 | Return values are always references to objects, strings, or true/false for
|
|---|
| 2067 | comparisation routines.
|
|---|
| 2068 |
|
|---|
| 2069 | =head1 WRAP YOUR OWN
|
|---|
| 2070 |
|
|---|
| 2071 | If you want to port your own favourite c-lib for big numbers to the
|
|---|
| 2072 | Math::BigInt interface, you can take any of the already existing modules as
|
|---|
| 2073 | a rough guideline. You should really wrap up the latest BigInt and BigFloat
|
|---|
| 2074 | testsuites with your module, and replace in them any of the following:
|
|---|
| 2075 |
|
|---|
| 2076 | use Math::BigInt;
|
|---|
| 2077 |
|
|---|
| 2078 | by this:
|
|---|
| 2079 |
|
|---|
| 2080 | use Math::BigInt lib => 'yourlib';
|
|---|
| 2081 |
|
|---|
| 2082 | This way you ensure that your library really works 100% within Math::BigInt.
|
|---|
| 2083 |
|
|---|
| 2084 | =head1 LICENSE
|
|---|
| 2085 |
|
|---|
| 2086 | This program is free software; you may redistribute it and/or modify it under
|
|---|
| 2087 | the same terms as Perl itself.
|
|---|
| 2088 |
|
|---|
| 2089 | =head1 AUTHORS
|
|---|
| 2090 |
|
|---|
| 2091 | Original math code by Mark Biggar, rewritten by Tels L<http://bloodgate.com/>
|
|---|
| 2092 | in late 2000.
|
|---|
| 2093 | Seperated from BigInt and shaped API with the help of John Peacock.
|
|---|
| 2094 |
|
|---|
| 2095 | Fixed, speed-up, streamlined and enhanced by Tels 2001 - 2005.
|
|---|
| 2096 |
|
|---|
| 2097 | =head1 SEE ALSO
|
|---|
| 2098 |
|
|---|
| 2099 | L<Math::BigInt>, L<Math::BigFloat>, L<Math::BigInt::BitVect>,
|
|---|
| 2100 | L<Math::BigInt::GMP>, L<Math::BigInt::FastCalc> and L<Math::BigInt::Pari>.
|
|---|
| 2101 |
|
|---|
| 2102 | =cut
|
|---|
