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