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