1package Math::BigInt::Lib; 2 3use 5.006001; 4use strict; 5use warnings; 6 7our $VERSION = '1.999811'; 8 9use Carp; 10 11use overload 12 13 # overload key: with_assign 14 15 '+' => sub { 16 my $class = ref $_[0]; 17 my $x = $class -> _copy($_[0]); 18 my $y = ref($_[1]) ? $_[1] : $class -> _new($_[1]); 19 return $class -> _add($x, $y); 20 }, 21 22 '-' => sub { 23 my $class = ref $_[0]; 24 my ($x, $y); 25 if ($_[2]) { # if swapped 26 $y = $_[0]; 27 $x = ref($_[1]) ? $_[1] : $class -> _new($_[1]); 28 } else { 29 $x = $class -> _copy($_[0]); 30 $y = ref($_[1]) ? $_[1] : $class -> _new($_[1]); 31 } 32 return $class -> _sub($x, $y); 33 }, 34 35 '*' => sub { 36 my $class = ref $_[0]; 37 my $x = $class -> _copy($_[0]); 38 my $y = ref($_[1]) ? $_[1] : $class -> _new($_[1]); 39 return $class -> _mul($x, $y); 40 }, 41 42 '/' => sub { 43 my $class = ref $_[0]; 44 my ($x, $y); 45 if ($_[2]) { # if swapped 46 $y = $_[0]; 47 $x = ref($_[1]) ? $_[1] : $class -> _new($_[1]); 48 } else { 49 $x = $class -> _copy($_[0]); 50 $y = ref($_[1]) ? $_[1] : $class -> _new($_[1]); 51 } 52 return $class -> _div($x, $y); 53 }, 54 55 '%' => sub { 56 my $class = ref $_[0]; 57 my ($x, $y); 58 if ($_[2]) { # if swapped 59 $y = $_[0]; 60 $x = ref($_[1]) ? $_[1] : $class -> _new($_[1]); 61 } else { 62 $x = $class -> _copy($_[0]); 63 $y = ref($_[1]) ? $_[1] : $class -> _new($_[1]); 64 } 65 return $class -> _mod($x, $y); 66 }, 67 68 '**' => sub { 69 my $class = ref $_[0]; 70 my ($x, $y); 71 if ($_[2]) { # if swapped 72 $y = $_[0]; 73 $x = ref($_[1]) ? $_[1] : $class -> _new($_[1]); 74 } else { 75 $x = $class -> _copy($_[0]); 76 $y = ref($_[1]) ? $_[1] : $class -> _new($_[1]); 77 } 78 return $class -> _pow($x, $y); 79 }, 80 81 '<<' => sub { 82 my $class = ref $_[0]; 83 my ($x, $y); 84 if ($_[2]) { # if swapped 85 $y = $class -> _num($_[0]); 86 $x = ref($_[1]) ? $_[1] : $class -> _new($_[1]); 87 } else { 88 $x = $_[0]; 89 $y = ref($_[1]) ? $class -> _num($_[1]) : $_[1]; 90 } 91 return $class -> _blsft($x, $y); 92 }, 93 94 '>>' => sub { 95 my $class = ref $_[0]; 96 my ($x, $y); 97 if ($_[2]) { # if swapped 98 $y = $_[0]; 99 $x = ref($_[1]) ? $_[1] : $class -> _new($_[1]); 100 } else { 101 $x = $class -> _copy($_[0]); 102 $y = ref($_[1]) ? $_[1] : $class -> _new($_[1]); 103 } 104 return $class -> _brsft($x, $y); 105 }, 106 107 # overload key: num_comparison 108 109 '<' => sub { 110 my $class = ref $_[0]; 111 my ($x, $y); 112 if ($_[2]) { # if swapped 113 $y = $_[0]; 114 $x = ref($_[1]) ? $_[1] : $class -> _new($_[1]); 115 } else { 116 $x = $class -> _copy($_[0]); 117 $y = ref($_[1]) ? $_[1] : $class -> _new($_[1]); 118 } 119 return $class -> _acmp($x, $y) < 0; 120 }, 121 122 '<=' => sub { 123 my $class = ref $_[0]; 124 my ($x, $y); 125 if ($_[2]) { # if swapped 126 $y = $_[0]; 127 $x = ref($_[1]) ? $_[1] : $class -> _new($_[1]); 128 } else { 129 $x = $class -> _copy($_[0]); 130 $y = ref($_[1]) ? $_[1] : $class -> _new($_[1]); 131 } 132 return $class -> _acmp($x, $y) <= 0; 133 }, 134 135 '>' => sub { 136 my $class = ref $_[0]; 137 my ($x, $y); 138 if ($_[2]) { # if swapped 139 $y = $_[0]; 140 $x = ref($_[1]) ? $_[1] : $class -> _new($_[1]); 141 } else { 142 $x = $class -> _copy($_[0]); 143 $y = ref($_[1]) ? $_[1] : $class -> _new($_[1]); 144 } 145 return $class -> _acmp($x, $y) > 0; 146 }, 147 148 '>=' => sub { 149 my $class = ref $_[0]; 150 my ($x, $y); 151 if ($_[2]) { # if swapped 152 $y = $_[0]; 153 $x = ref($_[1]) ? $_[1] : $class -> _new($_[1]); 154 } else { 155 $x = $class -> _copy($_[0]); 156 $y = ref($_[1]) ? $_[1] : $class -> _new($_[1]); 157 } 158 return $class -> _acmp($x, $y) >= 0; 159 }, 160 161 '==' => sub { 162 my $class = ref $_[0]; 163 my $x = $class -> _copy($_[0]); 164 my $y = ref($_[1]) ? $_[1] : $class -> _new($_[1]); 165 return $class -> _acmp($x, $y) == 0; 166 }, 167 168 '!=' => sub { 169 my $class = ref $_[0]; 170 my $x = $class -> _copy($_[0]); 171 my $y = ref($_[1]) ? $_[1] : $class -> _new($_[1]); 172 return $class -> _acmp($x, $y) != 0; 173 }, 174 175 # overload key: 3way_comparison 176 177 '<=>' => sub { 178 my $class = ref $_[0]; 179 my ($x, $y); 180 if ($_[2]) { # if swapped 181 $y = $_[0]; 182 $x = ref($_[1]) ? $_[1] : $class -> _new($_[1]); 183 } else { 184 $x = $class -> _copy($_[0]); 185 $y = ref($_[1]) ? $_[1] : $class -> _new($_[1]); 186 } 187 return $class -> _acmp($x, $y); 188 }, 189 190 # overload key: binary 191 192 '&' => sub { 193 my $class = ref $_[0]; 194 my ($x, $y); 195 if ($_[2]) { # if swapped 196 $y = $_[0]; 197 $x = ref($_[1]) ? $_[1] : $class -> _new($_[1]); 198 } else { 199 $x = $class -> _copy($_[0]); 200 $y = ref($_[1]) ? $_[1] : $class -> _new($_[1]); 201 } 202 return $class -> _and($x, $y); 203 }, 204 205 '|' => sub { 206 my $class = ref $_[0]; 207 my ($x, $y); 208 if ($_[2]) { # if swapped 209 $y = $_[0]; 210 $x = ref($_[1]) ? $_[1] : $class -> _new($_[1]); 211 } else { 212 $x = $class -> _copy($_[0]); 213 $y = ref($_[1]) ? $_[1] : $class -> _new($_[1]); 214 } 215 return $class -> _or($x, $y); 216 }, 217 218 '^' => sub { 219 my $class = ref $_[0]; 220 my ($x, $y); 221 if ($_[2]) { # if swapped 222 $y = $_[0]; 223 $x = ref($_[1]) ? $_[1] : $class -> _new($_[1]); 224 } else { 225 $x = $class -> _copy($_[0]); 226 $y = ref($_[1]) ? $_[1] : $class -> _new($_[1]); 227 } 228 return $class -> _xor($x, $y); 229 }, 230 231 # overload key: func 232 233 'abs' => sub { $_[0] }, 234 235 'sqrt' => sub { 236 my $class = ref $_[0]; 237 return $class -> _sqrt($class -> _copy($_[0])); 238 }, 239 240 'int' => sub { $_[0] }, 241 242 # overload key: conversion 243 244 'bool' => sub { ref($_[0]) -> _is_zero($_[0]) ? '' : 1; }, 245 246 '""' => sub { ref($_[0]) -> _str($_[0]); }, 247 248 '0+' => sub { ref($_[0]) -> _num($_[0]); }, 249 250 '=' => sub { ref($_[0]) -> _copy($_[0]); }, 251 252 ; 253 254# Do we need api_version() at all, now that we have a virtual parent class that 255# will provide any missing methods? Fixme! 256 257sub api_version () { 258 croak "@{[(caller 0)[3]]} method not implemented"; 259} 260 261sub _new { 262 croak "@{[(caller 0)[3]]} method not implemented"; 263} 264 265sub _zero { 266 my $class = shift; 267 return $class -> _new("0"); 268} 269 270sub _one { 271 my $class = shift; 272 return $class -> _new("1"); 273} 274 275sub _two { 276 my $class = shift; 277 return $class -> _new("2"); 278 279} 280sub _ten { 281 my $class = shift; 282 return $class -> _new("10"); 283} 284 285sub _1ex { 286 my ($class, $exp) = @_; 287 $exp = $class -> _num($exp) if ref($exp); 288 return $class -> _new("1" . ("0" x $exp)); 289} 290 291sub _copy { 292 my ($class, $x) = @_; 293 return $class -> _new($class -> _str($x)); 294} 295 296# catch and throw away 297sub import { } 298 299############################################################################## 300# convert back to string and number 301 302sub _str { 303 # Convert number from internal base 1eN format to string format. Internal 304 # format is always normalized, i.e., no leading zeros. 305 croak "@{[(caller 0)[3]]} method not implemented"; 306} 307 308sub _num { 309 my ($class, $x) = @_; 310 0 + $class -> _str($x); 311} 312 313############################################################################## 314# actual math code 315 316sub _add { 317 croak "@{[(caller 0)[3]]} method not implemented"; 318} 319 320sub _sub { 321 croak "@{[(caller 0)[3]]} method not implemented"; 322} 323 324sub _mul { 325 my ($class, $x, $y) = @_; 326 my $sum = $class -> _zero(); 327 my $i = $class -> _zero(); 328 while ($class -> _acmp($i, $y) < 0) { 329 $sum = $class -> _add($sum, $x); 330 $i = $class -> _inc($i); 331 } 332 return $sum; 333} 334 335sub _div { 336 my ($class, $x, $y) = @_; 337 338 croak "@{[(caller 0)[3]]} requires non-zero divisor" 339 if $class -> _is_zero($y); 340 341 my $r = $class -> _copy($x); 342 my $q = $class -> _zero(); 343 while ($class -> _acmp($r, $y) >= 0) { 344 $q = $class -> _inc($q); 345 $r = $class -> _sub($r, $y); 346 } 347 348 return $q, $r if wantarray; 349 return $q; 350} 351 352sub _inc { 353 my ($class, $x) = @_; 354 $class -> _add($x, $class -> _one()); 355} 356 357sub _dec { 358 my ($class, $x) = @_; 359 $class -> _sub($x, $class -> _one()); 360} 361 362############################################################################## 363# testing 364 365sub _acmp { 366 # Compare two (absolute) values. Return -1, 0, or 1. 367 my ($class, $x, $y) = @_; 368 my $xstr = $class -> _str($x); 369 my $ystr = $class -> _str($y); 370 371 length($xstr) <=> length($ystr) || $xstr cmp $ystr; 372} 373 374sub _len { 375 my ($class, $x) = @_; 376 CORE::length($class -> _str($x)); 377} 378 379sub _alen { 380 my ($class, $x) = @_; 381 $class -> _len($x); 382} 383 384sub _digit { 385 my ($class, $x, $n) = @_; 386 substr($class ->_str($x), -($n+1), 1); 387} 388 389sub _zeros { 390 my ($class, $x) = @_; 391 my $str = $class -> _str($x); 392 $str =~ /[^0](0*)\z/ ? CORE::length($1) : 0; 393} 394 395############################################################################## 396# _is_* routines 397 398sub _is_zero { 399 # return true if arg is zero 400 my ($class, $x) = @_; 401 $class -> _str($x) == 0; 402} 403 404sub _is_even { 405 # return true if arg is even 406 my ($class, $x) = @_; 407 substr($class -> _str($x), -1, 1) % 2 == 0; 408} 409 410sub _is_odd { 411 # return true if arg is odd 412 my ($class, $x) = @_; 413 substr($class -> _str($x), -1, 1) % 2 != 0; 414} 415 416sub _is_one { 417 # return true if arg is one 418 my ($class, $x) = @_; 419 $class -> _str($x) == 1; 420} 421 422sub _is_two { 423 # return true if arg is two 424 my ($class, $x) = @_; 425 $class -> _str($x) == 2; 426} 427 428sub _is_ten { 429 # return true if arg is ten 430 my ($class, $x) = @_; 431 $class -> _str($x) == 10; 432} 433 434############################################################################### 435# check routine to test internal state for corruptions 436 437sub _check { 438 # used by the test suite 439 my ($class, $x) = @_; 440 return "Input is undefined" unless defined $x; 441 return "$x is not a reference" unless ref($x); 442 return 0; 443} 444 445############################################################################### 446 447sub _mod { 448 # modulus 449 my ($class, $x, $y) = @_; 450 451 croak "@{[(caller 0)[3]]} requires non-zero second operand" 452 if $class -> _is_zero($y); 453 454 if ($class -> can('_div')) { 455 $x = $class -> _copy($x); 456 my ($q, $r) = $class -> _div($x, $y); 457 return $r; 458 } else { 459 my $r = $class -> _copy($x); 460 while ($class -> _acmp($r, $y) >= 0) { 461 $r = $class -> _sub($r, $y); 462 } 463 return $r; 464 } 465} 466 467############################################################################## 468# shifts 469 470sub _rsft { 471 my ($class, $x, $n, $b) = @_; 472 $b = $class -> _new($b) unless ref $b; 473 return scalar $class -> _div($x, $class -> _pow($class -> _copy($b), $n)); 474} 475 476sub _lsft { 477 my ($class, $x, $n, $b) = @_; 478 $b = $class -> _new($b) unless ref $b; 479 return $class -> _mul($x, $class -> _pow($class -> _copy($b), $n)); 480} 481 482sub _pow { 483 # power of $x to $y 484 my ($class, $x, $y) = @_; 485 486 if ($class -> _is_zero($y)) { 487 return $class -> _one(); # y == 0 => x => 1 488 } 489 490 if (($class -> _is_one($x)) || # x == 1 491 ($class -> _is_one($y))) # or y == 1 492 { 493 return $x; 494 } 495 496 if ($class -> _is_zero($x)) { 497 return $class -> _zero(); # 0 ** y => 0 (if not y <= 0) 498 } 499 500 my $pow2 = $class -> _one(); 501 502 my $y_bin = $class -> _as_bin($y); 503 $y_bin =~ s/^0b//; 504 my $len = length($y_bin); 505 506 while (--$len > 0) { 507 $pow2 = $class -> _mul($pow2, $x) if substr($y_bin, $len, 1) eq '1'; 508 $x = $class -> _mul($x, $x); 509 } 510 511 $x = $class -> _mul($x, $pow2); 512 return $x; 513} 514 515sub _nok { 516 # Return binomial coefficient (n over k). 517 my ($class, $n, $k) = @_; 518 519 # If k > n/2, or, equivalently, 2*k > n, compute nok(n, k) as 520 # nok(n, n-k), to minimize the number if iterations in the loop. 521 522 { 523 my $twok = $class -> _mul($class -> _two(), $class -> _copy($k)); 524 if ($class -> _acmp($twok, $n) > 0) { 525 $k = $class -> _sub($class -> _copy($n), $k); 526 } 527 } 528 529 # Example: 530 # 531 # / 7 \ 7! 1*2*3*4 * 5*6*7 5 * 6 * 7 532 # | | = --------- = --------------- = --------- = ((5 * 6) / 2 * 7) / 3 533 # \ 3 / (7-3)! 3! 1*2*3*4 * 1*2*3 1 * 2 * 3 534 # 535 # Equivalently, _nok(11, 5) is computed as 536 # 537 # (((((((7 * 8) / 2) * 9) / 3) * 10) / 4) * 11) / 5 538 539 if ($class -> _is_zero($k)) { 540 return $class -> _one(); 541 } 542 543 # Make a copy of the original n, in case the subclass modifies n in-place. 544 545 my $n_orig = $class -> _copy($n); 546 547 # n = 5, f = 6, d = 2 (cf. example above) 548 549 $n = $class -> _sub($n, $k); 550 $n = $class -> _inc($n); 551 552 my $f = $class -> _copy($n); 553 $f = $class -> _inc($f); 554 555 my $d = $class -> _two(); 556 557 # while f <= n (the original n, that is) ... 558 559 while ($class -> _acmp($f, $n_orig) <= 0) { 560 $n = $class -> _mul($n, $f); 561 $n = $class -> _div($n, $d); 562 $f = $class -> _inc($f); 563 $d = $class -> _inc($d); 564 } 565 566 return $n; 567} 568 569sub _fac { 570 # factorial 571 my ($class, $x) = @_; 572 573 my $two = $class -> _two(); 574 575 if ($class -> _acmp($x, $two) < 0) { 576 return $class -> _one(); 577 } 578 579 my $i = $class -> _copy($x); 580 while ($class -> _acmp($i, $two) > 0) { 581 $i = $class -> _dec($i); 582 $x = $class -> _mul($x, $i); 583 } 584 585 return $x; 586} 587 588sub _dfac { 589 # double factorial 590 my ($class, $x) = @_; 591 592 my $two = $class -> _two(); 593 594 if ($class -> _acmp($x, $two) < 0) { 595 return $class -> _one(); 596 } 597 598 my $i = $class -> _copy($x); 599 while ($class -> _acmp($i, $two) > 0) { 600 $i = $class -> _sub($i, $two); 601 $x = $class -> _mul($x, $i); 602 } 603 604 return $x; 605} 606 607sub _log_int { 608 # calculate integer log of $x to base $base 609 # calculate integer log of $x to base $base 610 # ref to array, ref to array - return ref to array 611 my ($class, $x, $base) = @_; 612 613 # X == 0 => NaN 614 return if $class -> _is_zero($x); 615 616 $base = $class -> _new(2) unless defined($base); 617 $base = $class -> _new($base) unless ref($base); 618 619 # BASE 0 or 1 => NaN 620 return if $class -> _is_zero($base) || $class -> _is_one($base); 621 622 # X == 1 => 0 (is exact) 623 if ($class -> _is_one($x)) { 624 return $class -> _zero(), 1; 625 } 626 627 my $cmp = $class -> _acmp($x, $base); 628 629 # X == BASE => 1 (is exact) 630 if ($cmp == 0) { 631 return $class -> _one(), 1; 632 } 633 634 # 1 < X < BASE => 0 (is truncated) 635 if ($cmp < 0) { 636 return $class -> _zero(), 0; 637 } 638 639 my $y; 640 641 # log(x) / log(b) = log(xm * 10^xe) / log(bm * 10^be) 642 # = (log(xm) + xe*(log(10))) / (log(bm) + be*log(10)) 643 644 { 645 my $x_str = $class -> _str($x); 646 my $b_str = $class -> _str($base); 647 my $xm = "." . $x_str; 648 my $bm = "." . $b_str; 649 my $xe = length($x_str); 650 my $be = length($b_str); 651 my $log10 = log(10); 652 my $guess = int((log($xm) + $xe * $log10) / (log($bm) + $be * $log10)); 653 $y = $class -> _new($guess); 654 } 655 656 my $trial = $class -> _pow($class -> _copy($base), $y); 657 my $acmp = $class -> _acmp($trial, $x); 658 659 # Did we get the exact result? 660 661 return $y, 1 if $acmp == 0; 662 663 # Too small? 664 665 while ($acmp < 0) { 666 $trial = $class -> _mul($trial, $base); 667 $y = $class -> _inc($y); 668 $acmp = $class -> _acmp($trial, $x); 669 } 670 671 # Too big? 672 673 while ($acmp > 0) { 674 $trial = $class -> _div($trial, $base); 675 $y = $class -> _dec($y); 676 $acmp = $class -> _acmp($trial, $x); 677 } 678 679 return $y, 1 if $acmp == 0; # result is exact 680 return $y, 0; # result is too small 681} 682 683sub _sqrt { 684 # square-root of $y in place 685 my ($class, $y) = @_; 686 687 return $y if $class -> _is_zero($y); 688 689 my $y_str = $class -> _str($y); 690 my $y_len = length($y_str); 691 692 # Compute the guess $x. 693 694 my $xm; 695 my $xe; 696 if ($y_len % 2 == 0) { 697 $xm = sqrt("." . $y_str); 698 $xe = $y_len / 2; 699 $xm = sprintf "%.0f", int($xm * 1e15); 700 $xe -= 15; 701 } else { 702 $xm = sqrt(".0" . $y_str); 703 $xe = ($y_len + 1) / 2; 704 $xm = sprintf "%.0f", int($xm * 1e16); 705 $xe -= 16; 706 } 707 708 my $x; 709 if ($xe < 0) { 710 $x = substr $xm, 0, length($xm) + $xe; 711 } else { 712 $x = $xm . ("0" x $xe); 713 } 714 715 $x = $class -> _new($x); 716 717 # Newton's method for computing square root of y 718 # 719 # x(i+1) = x(i) - f(x(i)) / f'(x(i)) 720 # = x(i) - (x(i)^2 - y) / (2 * x(i)) # use if x(i)^2 > y 721 # = y(i) + (y - x(i)^2) / (2 * x(i)) # use if x(i)^2 < y 722 723 # Determine if x, our guess, is too small, correct, or too large. 724 725 my $xsq = $class -> _mul($class -> _copy($x), $x); # x(i)^2 726 my $acmp = $class -> _acmp($xsq, $y); # x(i)^2 <=> y 727 728 # Only assign a value to this variable if we will be using it. 729 730 my $two; 731 $two = $class -> _two() if $acmp != 0; 732 733 # If x is too small, do one iteration of Newton's method. Since the 734 # function f(x) = x^2 - y is concave and monotonically increasing, the next 735 # guess for x will either be correct or too large. 736 737 if ($acmp < 0) { 738 739 # x(i+1) = x(i) + (y - x(i)^2) / (2 * x(i)) 740 741 my $numer = $class -> _sub($class -> _copy($y), $xsq); # y - x(i)^2 742 my $denom = $class -> _mul($class -> _copy($two), $x); # 2 * x(i) 743 my $delta = $class -> _div($numer, $denom); 744 745 unless ($class -> _is_zero($delta)) { 746 $x = $class -> _add($x, $delta); 747 $xsq = $class -> _mul($class -> _copy($x), $x); # x(i)^2 748 $acmp = $class -> _acmp($xsq, $y); # x(i)^2 <=> y 749 } 750 } 751 752 # If our guess for x is too large, apply Newton's method repeatedly until 753 # we either have got the correct value, or the delta is zero. 754 755 while ($acmp > 0) { 756 757 # x(i+1) = x(i) - (x(i)^2 - y) / (2 * x(i)) 758 759 my $numer = $class -> _sub($xsq, $y); # x(i)^2 - y 760 my $denom = $class -> _mul($class -> _copy($two), $x); # 2 * x(i) 761 my $delta = $class -> _div($numer, $denom); 762 last if $class -> _is_zero($delta); 763 764 $x = $class -> _sub($x, $delta); 765 $xsq = $class -> _mul($class -> _copy($x), $x); # x(i)^2 766 $acmp = $class -> _acmp($xsq, $y); # x(i)^2 <=> y 767 } 768 769 # When the delta is zero, our value for x might still be too large. We 770 # require that the outout is either exact or too small (i.e., rounded down 771 # to the nearest integer), so do a final check. 772 773 while ($acmp > 0) { 774 $x = $class -> _dec($x); 775 $xsq = $class -> _mul($class -> _copy($x), $x); # x(i)^2 776 $acmp = $class -> _acmp($xsq, $y); # x(i)^2 <=> y 777 } 778 779 return $x; 780} 781 782sub _root { 783 my ($class, $y, $n) = @_; 784 785 return $y if $class -> _is_zero($y) || $class -> _is_one($y) || 786 $class -> _is_one($n); 787 788 # If y <= n, the result is always (truncated to) 1. 789 790 return $class -> _one() if $class -> _acmp($y, $n) <= 0; 791 792 # Compute the initial guess x of y^(1/n). When n is large, Newton's method 793 # converges slowly if the "guess" (initial value) is poor, so we need a 794 # good guess. It the guess is too small, the next guess will be too large, 795 # and from then on all guesses are too large. 796 797 my $DEBUG = 0; 798 799 # Split y into mantissa and exponent in base 10, so that 800 # 801 # y = xm * 10^xe, where 0 < xm < 1 and xe is an integer 802 803 my $y_str = $class -> _str($y); 804 my $ym = "." . $y_str; 805 my $ye = length($y_str); 806 807 # From this compute the approximate base 10 logarithm of y 808 # 809 # log_10(y) = log_10(ym) + log_10(ye^10) 810 # = log(ym)/log(10) + ye 811 812 my $log10y = log($ym) / log(10) + $ye; 813 814 # And from this compute the approximate base 10 logarithm of x, where 815 # x = y^(1/n) 816 # 817 # log_10(x) = log_10(y)/n 818 819 my $log10x = $log10y / $class -> _num($n); 820 821 # From this compute xm and xe, the mantissa and exponent (in base 10) of x, 822 # where 1 < xm <= 10 and xe is an integer. 823 824 my $xe = int $log10x; 825 my $xm = 10 ** ($log10x - $xe); 826 827 # Scale the mantissa and exponent to increase the integer part of ym, which 828 # gives us better accuracy. 829 830 if ($DEBUG) { 831 print "\n"; 832 print "y_str = $y_str\n"; 833 print "ym = $ym\n"; 834 print "ye = $ye\n"; 835 print "log10y = $log10y\n"; 836 print "log10x = $log10x\n"; 837 print "xm = $xm\n"; 838 print "xe = $xe\n"; 839 } 840 841 my $d = $xe < 15 ? $xe : 15; 842 $xm *= 10 ** $d; 843 $xe -= $d; 844 845 if ($DEBUG) { 846 print "\n"; 847 print "xm = $xm\n"; 848 print "xe = $xe\n"; 849 } 850 851 # If the mantissa is not an integer, round up to nearest integer, and then 852 # convert the number to a string. It is important to always round up due to 853 # how Newton's method behaves in this case. If the initial guess is too 854 # small, the next guess will be too large, after which every succeeding 855 # guess converges the correct value from above. Now, if the initial guess 856 # is too small and n is large, the next guess will be much too large and 857 # require a large number of iterations to get close to the solution. 858 # Because of this, we are likely to find the solution faster if we make 859 # sure the initial guess is not too small. 860 861 my $xm_int = int($xm); 862 my $x_str = sprintf '%.0f', $xm > $xm_int ? $xm_int + 1 : $xm_int; 863 $x_str .= "0" x $xe; 864 865 my $x = $class -> _new($x_str); 866 867 if ($DEBUG) { 868 print "xm = $xm\n"; 869 print "xe = $xe\n"; 870 print "\n"; 871 print "x_str = $x_str (initial guess)\n"; 872 print "\n"; 873 } 874 875 # Use Newton's method for computing n'th root of y. 876 # 877 # x(i+1) = x(i) - f(x(i)) / f'(x(i)) 878 # = x(i) - (x(i)^n - y) / (n * x(i)^(n-1)) # use if x(i)^n > y 879 # = x(i) + (y - x(i)^n) / (n * x(i)^(n-1)) # use if x(i)^n < y 880 881 # Determine if x, our guess, is too small, correct, or too large. Rather 882 # than computing x(i)^n and x(i)^(n-1) directly, compute x(i)^(n-1) and 883 # then the same value multiplied by x. 884 885 my $nm1 = $class -> _dec($class -> _copy($n)); # n-1 886 my $xpownm1 = $class -> _pow($class -> _copy($x), $nm1); # x(i)^(n-1) 887 my $xpown = $class -> _mul($class -> _copy($xpownm1), $x); # x(i)^n 888 my $acmp = $class -> _acmp($xpown, $y); # x(i)^n <=> y 889 890 if ($DEBUG) { 891 print "\n"; 892 print "x = ", $class -> _str($x), "\n"; 893 print "x^n = ", $class -> _str($xpown), "\n"; 894 print "y = ", $class -> _str($y), "\n"; 895 print "acmp = $acmp\n"; 896 } 897 898 # If x is too small, do one iteration of Newton's method. Since the 899 # function f(x) = x^n - y is concave and monotonically increasing, the next 900 # guess for x will either be correct or too large. 901 902 if ($acmp < 0) { 903 904 # x(i+1) = x(i) + (y - x(i)^n) / (n * x(i)^(n-1)) 905 906 my $numer = $class -> _sub($class -> _copy($y), $xpown); # y - x(i)^n 907 my $denom = $class -> _mul($class -> _copy($n), $xpownm1); # n * x(i)^(n-1) 908 my $delta = $class -> _div($numer, $denom); 909 910 if ($DEBUG) { 911 print "\n"; 912 print "numer = ", $class -> _str($numer), "\n"; 913 print "denom = ", $class -> _str($denom), "\n"; 914 print "delta = ", $class -> _str($delta), "\n"; 915 } 916 917 unless ($class -> _is_zero($delta)) { 918 $x = $class -> _add($x, $delta); 919 $xpownm1 = $class -> _pow($class -> _copy($x), $nm1); # x(i)^(n-1) 920 $xpown = $class -> _mul($class -> _copy($xpownm1), $x); # x(i)^n 921 $acmp = $class -> _acmp($xpown, $y); # x(i)^n <=> y 922 923 if ($DEBUG) { 924 print "\n"; 925 print "x = ", $class -> _str($x), "\n"; 926 print "x^n = ", $class -> _str($xpown), "\n"; 927 print "y = ", $class -> _str($y), "\n"; 928 print "acmp = $acmp\n"; 929 } 930 } 931 } 932 933 # If our guess for x is too large, apply Newton's method repeatedly until 934 # we either have got the correct value, or the delta is zero. 935 936 while ($acmp > 0) { 937 938 # x(i+1) = x(i) - (x(i)^n - y) / (n * x(i)^(n-1)) 939 940 my $numer = $class -> _sub($class -> _copy($xpown), $y); # x(i)^n - y 941 my $denom = $class -> _mul($class -> _copy($n), $xpownm1); # n * x(i)^(n-1) 942 943 if ($DEBUG) { 944 print "numer = ", $class -> _str($numer), "\n"; 945 print "denom = ", $class -> _str($denom), "\n"; 946 } 947 948 my $delta = $class -> _div($numer, $denom); 949 950 if ($DEBUG) { 951 print "delta = ", $class -> _str($delta), "\n"; 952 } 953 954 last if $class -> _is_zero($delta); 955 956 $x = $class -> _sub($x, $delta); 957 $xpownm1 = $class -> _pow($class -> _copy($x), $nm1); # x(i)^(n-1) 958 $xpown = $class -> _mul($class -> _copy($xpownm1), $x); # x(i)^n 959 $acmp = $class -> _acmp($xpown, $y); # x(i)^n <=> y 960 961 if ($DEBUG) { 962 print "\n"; 963 print "x = ", $class -> _str($x), "\n"; 964 print "x^n = ", $class -> _str($xpown), "\n"; 965 print "y = ", $class -> _str($y), "\n"; 966 print "acmp = $acmp\n"; 967 } 968 } 969 970 # When the delta is zero, our value for x might still be too large. We 971 # require that the outout is either exact or too small (i.e., rounded down 972 # to the nearest integer), so do a final check. 973 974 while ($acmp > 0) { 975 $x = $class -> _dec($x); 976 $xpown = $class -> _pow($class -> _copy($x), $n); # x(i)^n 977 $acmp = $class -> _acmp($xpown, $y); # x(i)^n <=> y 978 } 979 980 return $x; 981} 982 983############################################################################## 984# binary stuff 985 986sub _and { 987 my ($class, $x, $y) = @_; 988 989 return $x if $class -> _acmp($x, $y) == 0; 990 991 my $m = $class -> _one(); 992 my $mask = $class -> _new("32768"); 993 994 my ($xr, $yr); # remainders after division 995 996 my $xc = $class -> _copy($x); 997 my $yc = $class -> _copy($y); 998 my $z = $class -> _zero(); 999 1000 until ($class -> _is_zero($xc) || $class -> _is_zero($yc)) { 1001 ($xc, $xr) = $class -> _div($xc, $mask); 1002 ($yc, $yr) = $class -> _div($yc, $mask); 1003 my $bits = $class -> _new($class -> _num($xr) & $class -> _num($yr)); 1004 $z = $class -> _add($z, $class -> _mul($bits, $m)); 1005 $m = $class -> _mul($m, $mask); 1006 } 1007 1008 return $z; 1009} 1010 1011sub _xor { 1012 my ($class, $x, $y) = @_; 1013 1014 return $class -> _zero() if $class -> _acmp($x, $y) == 0; 1015 1016 my $m = $class -> _one(); 1017 my $mask = $class -> _new("32768"); 1018 1019 my ($xr, $yr); # remainders after division 1020 1021 my $xc = $class -> _copy($x); 1022 my $yc = $class -> _copy($y); 1023 my $z = $class -> _zero(); 1024 1025 until ($class -> _is_zero($xc) || $class -> _is_zero($yc)) { 1026 ($xc, $xr) = $class -> _div($xc, $mask); 1027 ($yc, $yr) = $class -> _div($yc, $mask); 1028 my $bits = $class -> _new($class -> _num($xr) ^ $class -> _num($yr)); 1029 $z = $class -> _add($z, $class -> _mul($bits, $m)); 1030 $m = $class -> _mul($m, $mask); 1031 } 1032 1033 # The loop above stops when the smallest of the two numbers is exhausted. 1034 # The remainder of the longer one will survive bit-by-bit, so we simple 1035 # multiply-add it in. 1036 1037 $z = $class -> _add($z, $class -> _mul($xc, $m)) 1038 unless $class -> _is_zero($xc); 1039 $z = $class -> _add($z, $class -> _mul($yc, $m)) 1040 unless $class -> _is_zero($yc); 1041 1042 return $z; 1043} 1044 1045sub _or { 1046 my ($class, $x, $y) = @_; 1047 1048 return $x if $class -> _acmp($x, $y) == 0; # shortcut (see _and) 1049 1050 my $m = $class -> _one(); 1051 my $mask = $class -> _new("32768"); 1052 1053 my ($xr, $yr); # remainders after division 1054 1055 my $xc = $class -> _copy($x); 1056 my $yc = $class -> _copy($y); 1057 my $z = $class -> _zero(); 1058 1059 until ($class -> _is_zero($xc) || $class -> _is_zero($yc)) { 1060 ($xc, $xr) = $class -> _div($xc, $mask); 1061 ($yc, $yr) = $class -> _div($yc, $mask); 1062 my $bits = $class -> _new($class -> _num($xr) | $class -> _num($yr)); 1063 $z = $class -> _add($z, $class -> _mul($bits, $m)); 1064 $m = $class -> _mul($m, $mask); 1065 } 1066 1067 # The loop above stops when the smallest of the two numbers is exhausted. 1068 # The remainder of the longer one will survive bit-by-bit, so we simple 1069 # multiply-add it in. 1070 1071 $z = $class -> _add($z, $class -> _mul($xc, $m)) 1072 unless $class -> _is_zero($xc); 1073 $z = $class -> _add($z, $class -> _mul($yc, $m)) 1074 unless $class -> _is_zero($yc); 1075 1076 return $z; 1077} 1078 1079sub _to_bin { 1080 # convert the number to a string of binary digits without prefix 1081 my ($class, $x) = @_; 1082 my $str = ''; 1083 my $tmp = $class -> _copy($x); 1084 my $chunk = $class -> _new("16777216"); # 2^24 = 24 binary digits 1085 my $rem; 1086 until ($class -> _acmp($tmp, $chunk) < 0) { 1087 ($tmp, $rem) = $class -> _div($tmp, $chunk); 1088 $str = sprintf("%024b", $class -> _num($rem)) . $str; 1089 } 1090 unless ($class -> _is_zero($tmp)) { 1091 $str = sprintf("%b", $class -> _num($tmp)) . $str; 1092 } 1093 return length($str) ? $str : '0'; 1094} 1095 1096sub _to_oct { 1097 # convert the number to a string of octal digits without prefix 1098 my ($class, $x) = @_; 1099 my $str = ''; 1100 my $tmp = $class -> _copy($x); 1101 my $chunk = $class -> _new("16777216"); # 2^24 = 8 octal digits 1102 my $rem; 1103 until ($class -> _acmp($tmp, $chunk) < 0) { 1104 ($tmp, $rem) = $class -> _div($tmp, $chunk); 1105 $str = sprintf("%08o", $class -> _num($rem)) . $str; 1106 } 1107 unless ($class -> _is_zero($tmp)) { 1108 $str = sprintf("%o", $class -> _num($tmp)) . $str; 1109 } 1110 return length($str) ? $str : '0'; 1111} 1112 1113sub _to_hex { 1114 # convert the number to a string of hexadecimal digits without prefix 1115 my ($class, $x) = @_; 1116 my $str = ''; 1117 my $tmp = $class -> _copy($x); 1118 my $chunk = $class -> _new("16777216"); # 2^24 = 6 hexadecimal digits 1119 my $rem; 1120 until ($class -> _acmp($tmp, $chunk) < 0) { 1121 ($tmp, $rem) = $class -> _div($tmp, $chunk); 1122 $str = sprintf("%06x", $class -> _num($rem)) . $str; 1123 } 1124 unless ($class -> _is_zero($tmp)) { 1125 $str = sprintf("%x", $class -> _num($tmp)) . $str; 1126 } 1127 return length($str) ? $str : '0'; 1128} 1129 1130sub _as_bin { 1131 # convert the number to a string of binary digits with prefix 1132 my ($class, $x) = @_; 1133 return '0b' . $class -> _to_bin($x); 1134} 1135 1136sub _as_oct { 1137 # convert the number to a string of octal digits with prefix 1138 my ($class, $x) = @_; 1139 return '0' . $class -> _to_oct($x); # yes, 0 becomes "00" 1140} 1141 1142sub _as_hex { 1143 # convert the number to a string of hexadecimal digits with prefix 1144 my ($class, $x) = @_; 1145 return '0x' . $class -> _to_hex($x); 1146} 1147 1148sub _to_bytes { 1149 # convert the number to a string of bytes 1150 my ($class, $x) = @_; 1151 my $str = ''; 1152 my $tmp = $class -> _copy($x); 1153 my $chunk = $class -> _new("65536"); 1154 my $rem; 1155 until ($class -> _is_zero($tmp)) { 1156 ($tmp, $rem) = $class -> _div($tmp, $chunk); 1157 $str = pack('n', $class -> _num($rem)) . $str; 1158 } 1159 $str =~ s/^\0+//; 1160 return length($str) ? $str : "\x00"; 1161} 1162 1163*_as_bytes = \&_to_bytes; 1164 1165sub _from_hex { 1166 # Convert a string of hexadecimal digits to a number. 1167 1168 my ($class, $hex) = @_; 1169 $hex =~ s/^0[xX]//; 1170 1171 # Find the largest number of hexadecimal digits that we can safely use with 1172 # 32 bit integers. There are 4 bits pr hexadecimal digit, and we use only 1173 # 31 bits to play safe. This gives us int(31 / 4) = 7. 1174 1175 my $len = length $hex; 1176 my $rem = 1 + ($len - 1) % 7; 1177 1178 # Do the first chunk. 1179 1180 my $ret = $class -> _new(int hex substr $hex, 0, $rem); 1181 return $ret if $rem == $len; 1182 1183 # Do the remaining chunks, if any. 1184 1185 my $shift = $class -> _new(1 << (4 * 7)); 1186 for (my $offset = $rem ; $offset < $len ; $offset += 7) { 1187 my $part = int hex substr $hex, $offset, 7; 1188 $ret = $class -> _mul($ret, $shift); 1189 $ret = $class -> _add($ret, $class -> _new($part)); 1190 } 1191 1192 return $ret; 1193} 1194 1195sub _from_oct { 1196 # Convert a string of octal digits to a number. 1197 1198 my ($class, $oct) = @_; 1199 1200 # Find the largest number of octal digits that we can safely use with 32 1201 # bit integers. There are 3 bits pr octal digit, and we use only 31 bits to 1202 # play safe. This gives us int(31 / 3) = 10. 1203 1204 my $len = length $oct; 1205 my $rem = 1 + ($len - 1) % 10; 1206 1207 # Do the first chunk. 1208 1209 my $ret = $class -> _new(int oct substr $oct, 0, $rem); 1210 return $ret if $rem == $len; 1211 1212 # Do the remaining chunks, if any. 1213 1214 my $shift = $class -> _new(1 << (3 * 10)); 1215 for (my $offset = $rem ; $offset < $len ; $offset += 10) { 1216 my $part = int oct substr $oct, $offset, 10; 1217 $ret = $class -> _mul($ret, $shift); 1218 $ret = $class -> _add($ret, $class -> _new($part)); 1219 } 1220 1221 return $ret; 1222} 1223 1224sub _from_bin { 1225 # Convert a string of binary digits to a number. 1226 1227 my ($class, $bin) = @_; 1228 $bin =~ s/^0[bB]//; 1229 1230 # The largest number of binary digits that we can safely use with 32 bit 1231 # integers is 31. We use only 31 bits to play safe. 1232 1233 my $len = length $bin; 1234 my $rem = 1 + ($len - 1) % 31; 1235 1236 # Do the first chunk. 1237 1238 my $ret = $class -> _new(int oct '0b' . substr $bin, 0, $rem); 1239 return $ret if $rem == $len; 1240 1241 # Do the remaining chunks, if any. 1242 1243 my $shift = $class -> _new(1 << 31); 1244 for (my $offset = $rem ; $offset < $len ; $offset += 31) { 1245 my $part = int oct '0b' . substr $bin, $offset, 31; 1246 $ret = $class -> _mul($ret, $shift); 1247 $ret = $class -> _add($ret, $class -> _new($part)); 1248 } 1249 1250 return $ret; 1251} 1252 1253sub _from_bytes { 1254 # convert string of bytes to a number 1255 my ($class, $str) = @_; 1256 my $x = $class -> _zero(); 1257 my $base = $class -> _new("256"); 1258 my $n = length($str); 1259 for (my $i = 0 ; $i < $n ; ++$i) { 1260 $x = $class -> _mul($x, $base); 1261 my $byteval = $class -> _new(unpack 'C', substr($str, $i, 1)); 1262 $x = $class -> _add($x, $byteval); 1263 } 1264 return $x; 1265} 1266 1267############################################################################## 1268# special modulus functions 1269 1270sub _modinv { 1271 # modular multiplicative inverse 1272 my ($class, $x, $y) = @_; 1273 1274 # modulo zero 1275 if ($class -> _is_zero($y)) { 1276 return (undef, undef); 1277 } 1278 1279 # modulo one 1280 if ($class -> _is_one($y)) { 1281 return ($class -> _zero(), '+'); 1282 } 1283 1284 my $u = $class -> _zero(); 1285 my $v = $class -> _one(); 1286 my $a = $class -> _copy($y); 1287 my $b = $class -> _copy($x); 1288 1289 # Euclid's Algorithm for bgcd(). 1290 1291 my $q; 1292 my $sign = 1; 1293 { 1294 ($a, $q, $b) = ($b, $class -> _div($a, $b)); 1295 last if $class -> _is_zero($b); 1296 1297 my $vq = $class -> _mul($class -> _copy($v), $q); 1298 my $t = $class -> _add($vq, $u); 1299 $u = $v; 1300 $v = $t; 1301 $sign = -$sign; 1302 redo; 1303 } 1304 1305 # if the gcd is not 1, there exists no modular multiplicative inverse 1306 return (undef, undef) unless $class -> _is_one($a); 1307 1308 ($v, $sign == 1 ? '+' : '-'); 1309} 1310 1311sub _modpow { 1312 # modulus of power ($x ** $y) % $z 1313 my ($class, $num, $exp, $mod) = @_; 1314 1315 # a^b (mod 1) = 0 for all a and b 1316 if ($class -> _is_one($mod)) { 1317 return $class -> _zero(); 1318 } 1319 1320 # 0^a (mod m) = 0 if m != 0, a != 0 1321 # 0^0 (mod m) = 1 if m != 0 1322 if ($class -> _is_zero($num)) { 1323 return $class -> _is_zero($exp) ? $class -> _one() 1324 : $class -> _zero(); 1325 } 1326 1327 # $num = $class -> _mod($num, $mod); # this does not make it faster 1328 1329 my $acc = $class -> _copy($num); 1330 my $t = $class -> _one(); 1331 1332 my $expbin = $class -> _as_bin($exp); 1333 $expbin =~ s/^0b//; 1334 my $len = length($expbin); 1335 1336 while (--$len >= 0) { 1337 if (substr($expbin, $len, 1) eq '1') { 1338 $t = $class -> _mul($t, $acc); 1339 $t = $class -> _mod($t, $mod); 1340 } 1341 $acc = $class -> _mul($acc, $acc); 1342 $acc = $class -> _mod($acc, $mod); 1343 } 1344 return $t; 1345} 1346 1347sub _gcd { 1348 # Greatest common divisor. 1349 1350 my ($class, $x, $y) = @_; 1351 1352 # gcd(0, 0) = 0 1353 # gcd(0, a) = a, if a != 0 1354 1355 if ($class -> _acmp($x, $y) == 0) { 1356 return $class -> _copy($x); 1357 } 1358 1359 if ($class -> _is_zero($x)) { 1360 if ($class -> _is_zero($y)) { 1361 return $class -> _zero(); 1362 } else { 1363 return $class -> _copy($y); 1364 } 1365 } else { 1366 if ($class -> _is_zero($y)) { 1367 return $class -> _copy($x); 1368 } else { 1369 1370 # Until $y is zero ... 1371 1372 $x = $class -> _copy($x); 1373 until ($class -> _is_zero($y)) { 1374 1375 # Compute remainder. 1376 1377 $x = $class -> _mod($x, $y); 1378 1379 # Swap $x and $y. 1380 1381 my $tmp = $x; 1382 $x = $class -> _copy($y); 1383 $y = $tmp; 1384 } 1385 1386 return $x; 1387 } 1388 } 1389} 1390 1391sub _lcm { 1392 # Least common multiple. 1393 1394 my ($class, $x, $y) = @_; 1395 1396 # lcm(0, x) = 0 for all x 1397 1398 return $class -> _zero() 1399 if ($class -> _is_zero($x) || 1400 $class -> _is_zero($y)); 1401 1402 my $gcd = $class -> _gcd($class -> _copy($x), $y); 1403 $x = $class -> _div($x, $gcd); 1404 $x = $class -> _mul($x, $y); 1405 return $x; 1406} 1407 1408sub _lucas { 1409 my ($class, $n) = @_; 1410 1411 $n = $class -> _num($n) if ref $n; 1412 1413 # In list context, use lucas(n) = lucas(n-1) + lucas(n-2) 1414 1415 if (wantarray) { 1416 my @y; 1417 1418 push @y, $class -> _two(); 1419 return @y if $n == 0; 1420 1421 push @y, $class -> _one(); 1422 return @y if $n == 1; 1423 1424 for (my $i = 2 ; $i <= $n ; ++ $i) { 1425 $y[$i] = $class -> _add($class -> _copy($y[$i - 1]), $y[$i - 2]); 1426 } 1427 1428 return @y; 1429 } 1430 1431 require Scalar::Util; 1432 1433 # In scalar context use that lucas(n) = fib(n-1) + fib(n+1). 1434 # 1435 # Remember that _fib() behaves differently in scalar context and list 1436 # context, so we must add scalar() to get the desired behaviour. 1437 1438 return $class -> _two() if $n == 0; 1439 1440 return $class -> _add(scalar $class -> _fib($n - 1), 1441 scalar $class -> _fib($n + 1)); 1442} 1443 1444sub _fib { 1445 my ($class, $n) = @_; 1446 1447 $n = $class -> _num($n) if ref $n; 1448 1449 # In list context, use fib(n) = fib(n-1) + fib(n-2) 1450 1451 if (wantarray) { 1452 my @y; 1453 1454 push @y, $class -> _zero(); 1455 return @y if $n == 0; 1456 1457 push @y, $class -> _one(); 1458 return @y if $n == 1; 1459 1460 for (my $i = 2 ; $i <= $n ; ++ $i) { 1461 $y[$i] = $class -> _add($class -> _copy($y[$i - 1]), $y[$i - 2]); 1462 } 1463 1464 return @y; 1465 } 1466 1467 # In scalar context use a fast algorithm that is much faster than the 1468 # recursive algorith used in list context. 1469 1470 my $cache = {}; 1471 my $two = $class -> _two(); 1472 my $fib; 1473 1474 $fib = sub { 1475 my $n = shift; 1476 return $class -> _zero() if $n <= 0; 1477 return $class -> _one() if $n <= 2; 1478 return $cache -> {$n} if exists $cache -> {$n}; 1479 1480 my $k = int($n / 2); 1481 my $a = $fib -> ($k + 1); 1482 my $b = $fib -> ($k); 1483 my $y; 1484 1485 if ($n % 2 == 1) { 1486 # a*a + b*b 1487 $y = $class -> _add($class -> _mul($class -> _copy($a), $a), 1488 $class -> _mul($class -> _copy($b), $b)); 1489 } else { 1490 # (2*a - b)*b 1491 $y = $class -> _mul($class -> _sub($class -> _mul( 1492 $class -> _copy($two), $a), $b), $b); 1493 } 1494 1495 $cache -> {$n} = $y; 1496 return $y; 1497 }; 1498 1499 return $fib -> ($n); 1500} 1501 1502############################################################################## 1503############################################################################## 1504 15051; 1506 1507__END__ 1508 1509=pod 1510 1511=head1 NAME 1512 1513Math::BigInt::Lib - virtual parent class for Math::BigInt libraries 1514 1515=head1 SYNOPSIS 1516 1517 # In the backend library for Math::BigInt et al. 1518 1519 package Math::BigInt::MyBackend; 1520 1521 use Math::BigInt::lib; 1522 our @ISA = qw< Math::BigInt::lib >; 1523 1524 sub _new { ... } 1525 sub _str { ... } 1526 sub _add { ... } 1527 str _sub { ... } 1528 ... 1529 1530 # In your main program. 1531 1532 use Math::BigInt lib => 'MyBackend'; 1533 1534=head1 DESCRIPTION 1535 1536This module provides support for big integer calculations. It is not intended 1537to be used directly, but rather as a parent class for backend libraries used by 1538Math::BigInt, Math::BigFloat, Math::BigRat, and related modules. 1539 1540Other backend libraries include Math::BigInt::Calc, Math::BigInt::FastCalc, 1541Math::BigInt::GMP, and Math::BigInt::Pari. 1542 1543In order to allow for multiple big integer libraries, Math::BigInt was 1544rewritten to use a plug-in library for core math routines. Any module which 1545conforms to the API can be used by Math::BigInt by using this in your program: 1546 1547 use Math::BigInt lib => 'libname'; 1548 1549'libname' is either the long name, like 'Math::BigInt::Pari', or only the short 1550version, like 'Pari'. 1551 1552=head2 General Notes 1553 1554A library only needs to deal with unsigned big integers. Testing of input 1555parameter validity is done by the caller, so there is no need to worry about 1556underflow (e.g., in C<_sub()> and C<_dec()>) or about division by zero (e.g., 1557in C<_div()> and C<_mod()>)) or similar cases. 1558 1559Some libraries use methods that don't modify their argument, and some libraries 1560don't even use objects, but rather unblessed references. Because of this, 1561liberary methods are always called as class methods, not instance methods: 1562 1563 $x = Class -> method($x, $y); # like this 1564 $x = $x -> method($y); # not like this ... 1565 $x -> method($y); # ... or like this 1566 1567And with boolean methods 1568 1569 $bool = Class -> method($x, $y); # like this 1570 $bool = $x -> method($y); # not like this 1571 1572Return values are always objects, strings, Perl scalars, or true/false for 1573comparison routines. 1574 1575=head3 API version 1576 1577=over 4 1578 1579=item CLASS-E<gt>api_version() 1580 1581Return API version as a Perl scalar, 1 for Math::BigInt v1.70, 2 for 1582Math::BigInt v1.83. 1583 1584This method is no longer used. Methods that are not implemented by a subclass 1585will be inherited from this class. 1586 1587=back 1588 1589=head3 Constructors 1590 1591The following methods are mandatory: _new(), _str(), _add(), and _sub(). 1592However, computations will be very slow without _mul() and _div(). 1593 1594=over 4 1595 1596=item CLASS-E<gt>_new(STR) 1597 1598Convert a string representing an unsigned decimal number to an object 1599representing the same number. The input is normalized, i.e., it matches 1600C<^(0|[1-9]\d*)$>. 1601 1602=item CLASS-E<gt>_zero() 1603 1604Return an object representing the number zero. 1605 1606=item CLASS-E<gt>_one() 1607 1608Return an object representing the number one. 1609 1610=item CLASS-E<gt>_two() 1611 1612Return an object representing the number two. 1613 1614=item CLASS-E<gt>_ten() 1615 1616Return an object representing the number ten. 1617 1618=item CLASS-E<gt>_from_bin(STR) 1619 1620Return an object given a string representing a binary number. The input has a 1621'0b' prefix and matches the regular expression C<^0[bB](0|1[01]*)$>. 1622 1623=item CLASS-E<gt>_from_oct(STR) 1624 1625Return an object given a string representing an octal number. The input has a 1626'0' prefix and matches the regular expression C<^0[1-7]*$>. 1627 1628=item CLASS-E<gt>_from_hex(STR) 1629 1630Return an object given a string representing a hexadecimal number. The input 1631has a '0x' prefix and matches the regular expression 1632C<^0x(0|[1-9a-fA-F][\da-fA-F]*)$>. 1633 1634=item CLASS-E<gt>_from_bytes(STR) 1635 1636Returns an object given a byte string representing the number. The byte string 1637is in big endian byte order, so the two-byte input string "\x01\x00" should 1638give an output value representing the number 256. 1639 1640=back 1641 1642=head3 Mathematical functions 1643 1644=over 4 1645 1646=item CLASS-E<gt>_add(OBJ1, OBJ2) 1647 1648Returns the result of adding OBJ2 to OBJ1. 1649 1650=item CLASS-E<gt>_mul(OBJ1, OBJ2) 1651 1652Returns the result of multiplying OBJ2 and OBJ1. 1653 1654=item CLASS-E<gt>_div(OBJ1, OBJ2) 1655 1656In scalar context, returns the quotient after dividing OBJ1 by OBJ2 and 1657truncating the result to an integer. In list context, return the quotient and 1658the remainder. 1659 1660=item CLASS-E<gt>_sub(OBJ1, OBJ2, FLAG) 1661 1662=item CLASS-E<gt>_sub(OBJ1, OBJ2) 1663 1664Returns the result of subtracting OBJ2 by OBJ1. If C<flag> is false or omitted, 1665OBJ1 might be modified. If C<flag> is true, OBJ2 might be modified. 1666 1667=item CLASS-E<gt>_dec(OBJ) 1668 1669Returns the result after decrementing OBJ by one. 1670 1671=item CLASS-E<gt>_inc(OBJ) 1672 1673Returns the result after incrementing OBJ by one. 1674 1675=item CLASS-E<gt>_mod(OBJ1, OBJ2) 1676 1677Returns OBJ1 modulo OBJ2, i.e., the remainder after dividing OBJ1 by OBJ2. 1678 1679=item CLASS-E<gt>_sqrt(OBJ) 1680 1681Returns the square root of OBJ, truncated to an integer. 1682 1683=item CLASS-E<gt>_root(OBJ, N) 1684 1685Returns the Nth root of OBJ, truncated to an integer. 1686 1687=item CLASS-E<gt>_fac(OBJ) 1688 1689Returns the factorial of OBJ, i.e., the product of all positive integers up to 1690and including OBJ. 1691 1692=item CLASS-E<gt>_dfac(OBJ) 1693 1694Returns the double factorial of OBJ. If OBJ is an even integer, returns the 1695product of all positive, even integers up to and including OBJ, i.e., 16962*4*6*...*OBJ. If OBJ is an odd integer, returns the product of all positive, 1697odd integers, i.e., 1*3*5*...*OBJ. 1698 1699=item CLASS-E<gt>_pow(OBJ1, OBJ2) 1700 1701Returns OBJ1 raised to the power of OBJ2. By convention, 0**0 = 1. 1702 1703=item CLASS-E<gt>_modinv(OBJ1, OBJ2) 1704 1705Returns the modular multiplicative inverse, i.e., return OBJ3 so that 1706 1707 (OBJ3 * OBJ1) % OBJ2 = 1 % OBJ2 1708 1709The result is returned as two arguments. If the modular multiplicative inverse 1710does not exist, both arguments are undefined. Otherwise, the arguments are a 1711number (object) and its sign ("+" or "-"). 1712 1713The output value, with its sign, must either be a positive value in the range 17141,2,...,OBJ2-1 or the same value subtracted OBJ2. For instance, if the input 1715arguments are objects representing the numbers 7 and 5, the method must either 1716return an object representing the number 3 and a "+" sign, since (3*7) % 5 = 1 1717% 5, or an object representing the number 2 and a "-" sign, since (-2*7) % 5 = 1 1718% 5. 1719 1720=item CLASS-E<gt>_modpow(OBJ1, OBJ2, OBJ3) 1721 1722Returns the modular exponentiation, i.e., (OBJ1 ** OBJ2) % OBJ3. 1723 1724=item CLASS-E<gt>_rsft(OBJ, N, B) 1725 1726Returns the result after shifting OBJ N digits to thee right in base B. This is 1727equivalent to performing integer division by B**N and discarding the remainder, 1728except that it might be much faster. 1729 1730For instance, if the object $obj represents the hexadecimal number 0xabcde, 1731then C<_rsft($obj, 2, 16)> returns an object representing the number 0xabc. The 1732"remainer", 0xde, is discarded and not returned. 1733 1734=item CLASS-E<gt>_lsft(OBJ, N, B) 1735 1736Returns the result after shifting OBJ N digits to the left in base B. This is 1737equivalent to multiplying by B**N, except that it might be much faster. 1738 1739=item CLASS-E<gt>_log_int(OBJ, B) 1740 1741Returns the logarithm of OBJ to base BASE truncted to an integer. This method 1742has two output arguments, the OBJECT and a STATUS. The STATUS is Perl scalar; 1743it is 1 if OBJ is the exact result, 0 if the result was truncted to give OBJ, 1744and undef if it is unknown whether OBJ is the exact result. 1745 1746=item CLASS-E<gt>_gcd(OBJ1, OBJ2) 1747 1748Returns the greatest common divisor of OBJ1 and OBJ2. 1749 1750=item CLASS-E<gt>_lcm(OBJ1, OBJ2) 1751 1752Return the least common multiple of OBJ1 and OBJ2. 1753 1754=item CLASS-E<gt>_fib(OBJ) 1755 1756In scalar context, returns the nth Fibonacci number: _fib(0) returns 0, _fib(1) 1757returns 1, _fib(2) returns 1, _fib(3) returns 2 etc. In list context, returns 1758the Fibonacci numbers from F(0) to F(n): 0, 1, 1, 2, 3, 5, 8, 13, 21, 34, ... 1759 1760=item CLASS-E<gt>_lucas(OBJ) 1761 1762In scalar context, returns the nth Lucas number: _lucas(0) returns 2, _lucas(1) 1763returns 1, _lucas(2) returns 3, etc. In list context, returns the Lucas numbers 1764from L(0) to L(n): 2, 1, 3, 4, 7, 11, 18, 29,47, 76, ... 1765 1766=back 1767 1768=head3 Bitwise operators 1769 1770=over 4 1771 1772=item CLASS-E<gt>_and(OBJ1, OBJ2) 1773 1774Returns bitwise and. 1775 1776=item CLASS-E<gt>_or(OBJ1, OBJ2) 1777 1778Return bitwise or. 1779 1780=item CLASS-E<gt>_xor(OBJ1, OBJ2) 1781 1782Return bitwise exclusive or. 1783 1784=back 1785 1786=head3 Boolean operators 1787 1788=over 4 1789 1790=item CLASS-E<gt>_is_zero(OBJ) 1791 1792Returns a true value if OBJ is zero, and false value otherwise. 1793 1794=item CLASS-E<gt>_is_one(OBJ) 1795 1796Returns a true value if OBJ is one, and false value otherwise. 1797 1798=item CLASS-E<gt>_is_two(OBJ) 1799 1800Returns a true value if OBJ is two, and false value otherwise. 1801 1802=item CLASS-E<gt>_is_ten(OBJ) 1803 1804Returns a true value if OBJ is ten, and false value otherwise. 1805 1806=item CLASS-E<gt>_is_even(OBJ) 1807 1808Return a true value if OBJ is an even integer, and a false value otherwise. 1809 1810=item CLASS-E<gt>_is_odd(OBJ) 1811 1812Return a true value if OBJ is an even integer, and a false value otherwise. 1813 1814=item CLASS-E<gt>_acmp(OBJ1, OBJ2) 1815 1816Compare OBJ1 and OBJ2 and return -1, 0, or 1, if OBJ1 is numerically less than, 1817equal to, or larger than OBJ2, respectively. 1818 1819=back 1820 1821=head3 String conversion 1822 1823=over 4 1824 1825=item CLASS-E<gt>_str(OBJ) 1826 1827Returns a string representing OBJ in decimal notation. The returned string 1828should have no leading zeros, i.e., it should match C<^(0|[1-9]\d*)$>. 1829 1830=item CLASS-E<gt>_to_bin(OBJ) 1831 1832Returns the binary string representation of OBJ. 1833 1834=item CLASS-E<gt>_to_oct(OBJ) 1835 1836Returns the octal string representation of the number. 1837 1838=item CLASS-E<gt>_to_hex(OBJ) 1839 1840Returns the hexadecimal string representation of the number. 1841 1842=item CLASS-E<gt>_to_bytes(OBJ) 1843 1844Returns a byte string representation of OBJ. The byte string is in big endian 1845byte order, so if OBJ represents the number 256, the output should be the 1846two-byte string "\x01\x00". 1847 1848=item CLASS-E<gt>_as_bin(OBJ) 1849 1850Like C<_to_bin()> but with a '0b' prefix. 1851 1852=item CLASS-E<gt>_as_oct(OBJ) 1853 1854Like C<_to_oct()> but with a '0' prefix. 1855 1856=item CLASS-E<gt>_as_hex(OBJ) 1857 1858Like C<_to_hex()> but with a '0x' prefix. 1859 1860=item CLASS-E<gt>_as_bytes(OBJ) 1861 1862This is an alias to C<_to_bytes()>. 1863 1864=back 1865 1866=head3 Numeric conversion 1867 1868=over 4 1869 1870=item CLASS-E<gt>_num(OBJ) 1871 1872Returns a Perl scalar number representing the number OBJ as close as 1873possible. Since Perl scalars have limited precision, the returned value might 1874not be exactly the same as OBJ. 1875 1876=back 1877 1878=head3 Miscellaneous 1879 1880=over 4 1881 1882=item CLASS-E<gt>_copy(OBJ) 1883 1884Returns a true copy OBJ. 1885 1886=item CLASS-E<gt>_len(OBJ) 1887 1888Returns the number of the decimal digits in OBJ. The output is a Perl scalar. 1889 1890=item CLASS-E<gt>_zeros(OBJ) 1891 1892Returns the number of trailing decimal zeros. The output is a Perl scalar. The 1893number zero has no trailing decimal zeros. 1894 1895=item CLASS-E<gt>_digit(OBJ, N) 1896 1897Returns the Nth digit in OBJ as a Perl scalar. N is a Perl scalar, where zero 1898refers to the rightmost (least significant) digit, and negative values count 1899from the left (most significant digit). If $obj represents the number 123, then 1900 1901 CLASS->_digit($obj, 0) # returns 3 1902 CLASS->_digit($obj, 1) # returns 2 1903 CLASS->_digit($obj, 2) # returns 1 1904 CLASS->_digit($obj, -1) # returns 1 1905 1906=item CLASS-E<gt>_check(OBJ) 1907 1908Returns true if the object is invalid and false otherwise. Preferably, the true 1909value is a string describing the problem with the object. This is a check 1910routine to test the internal state of the object for corruption. 1911 1912=item CLASS-E<gt>_set(OBJ) 1913 1914xxx 1915 1916=back 1917 1918=head2 API version 2 1919 1920The following methods are required for an API version of 2 or greater. 1921 1922=head3 Constructors 1923 1924=over 4 1925 1926=item CLASS-E<gt>_1ex(N) 1927 1928Return an object representing the number 10**N where N E<gt>= 0 is a Perl 1929scalar. 1930 1931=back 1932 1933=head3 Mathematical functions 1934 1935=over 4 1936 1937=item CLASS-E<gt>_nok(OBJ1, OBJ2) 1938 1939Return the binomial coefficient OBJ1 over OBJ1. 1940 1941=back 1942 1943=head3 Miscellaneous 1944 1945=over 4 1946 1947=item CLASS-E<gt>_alen(OBJ) 1948 1949Return the approximate number of decimal digits of the object. The output is a 1950Perl scalar. 1951 1952=back 1953 1954=head2 API optional methods 1955 1956The following methods are optional, and can be defined if the underlying lib 1957has a fast way to do them. If undefined, Math::BigInt will use pure Perl (hence 1958slow) fallback routines to emulate these: 1959 1960=head3 Signed bitwise operators. 1961 1962=over 4 1963 1964=item CLASS-E<gt>_signed_or(OBJ1, OBJ2, SIGN1, SIGN2) 1965 1966Return the signed bitwise or. 1967 1968=item CLASS-E<gt>_signed_and(OBJ1, OBJ2, SIGN1, SIGN2) 1969 1970Return the signed bitwise and. 1971 1972=item CLASS-E<gt>_signed_xor(OBJ1, OBJ2, SIGN1, SIGN2) 1973 1974Return the signed bitwise exclusive or. 1975 1976=back 1977 1978=head1 WRAP YOUR OWN 1979 1980If you want to port your own favourite C library for big numbers to the 1981Math::BigInt interface, you can take any of the already existing modules as a 1982rough guideline. You should really wrap up the latest Math::BigInt and 1983Math::BigFloat testsuites with your module, and replace in them any of the 1984following: 1985 1986 use Math::BigInt; 1987 1988by this: 1989 1990 use Math::BigInt lib => 'yourlib'; 1991 1992This way you ensure that your library really works 100% within Math::BigInt. 1993 1994=head1 BUGS 1995 1996Please report any bugs or feature requests to 1997C<bug-math-bigint at rt.cpan.org>, or through the web interface at 1998L<https://rt.cpan.org/Ticket/Create.html?Queue=Math-BigInt> 1999(requires login). 2000We will be notified, and then you'll automatically be notified of progress on 2001your bug as I make changes. 2002 2003=head1 SUPPORT 2004 2005You can find documentation for this module with the perldoc command. 2006 2007 perldoc Math::BigInt::Calc 2008 2009You can also look for information at: 2010 2011=over 4 2012 2013=item * RT: CPAN's request tracker 2014 2015L<https://rt.cpan.org/Public/Dist/Display.html?Name=Math-BigInt> 2016 2017=item * AnnoCPAN: Annotated CPAN documentation 2018 2019L<http://annocpan.org/dist/Math-BigInt> 2020 2021=item * CPAN Ratings 2022 2023L<http://cpanratings.perl.org/dist/Math-BigInt> 2024 2025=item * Search CPAN 2026 2027L<http://search.cpan.org/dist/Math-BigInt/> 2028 2029=item * CPAN Testers Matrix 2030 2031L<http://matrix.cpantesters.org/?dist=Math-BigInt> 2032 2033=item * The Bignum mailing list 2034 2035=over 4 2036 2037=item * Post to mailing list 2038 2039C<bignum at lists.scsys.co.uk> 2040 2041=item * View mailing list 2042 2043L<http://lists.scsys.co.uk/pipermail/bignum/> 2044 2045=item * Subscribe/Unsubscribe 2046 2047L<http://lists.scsys.co.uk/cgi-bin/mailman/listinfo/bignum> 2048 2049=back 2050 2051=back 2052 2053=head1 LICENSE 2054 2055This program is free software; you may redistribute it and/or modify it under 2056the same terms as Perl itself. 2057 2058=head1 AUTHOR 2059 2060Peter John Acklam, E<lt>pjacklam@online.noE<gt> 2061 2062Code and documentation based on the Math::BigInt::Calc module by Tels 2063E<lt>nospam-abuse@bloodgate.comE<gt> 2064 2065=head1 SEE ALSO 2066 2067L<Math::BigInt>, L<Math::BigInt::Calc>, L<Math::BigInt::GMP>, 2068L<Math::BigInt::FastCalc> and L<Math::BigInt::Pari>. 2069 2070=cut 2071