home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / perl560.zip / lib / Math / BigInt.pm < prev    next >
Text File  |  1999-07-20  |  13KB  |  493 lines
  1. package Math::BigInt;
  2.  
  3. use overload
  4. '+'    =>    sub {new Math::BigInt &badd},
  5. '-'    =>    sub {new Math::BigInt
  6.                $_[2]? bsub($_[1],${$_[0]}) : bsub(${$_[0]},$_[1])},
  7. '<=>'    =>    sub {$_[2]? bcmp($_[1],${$_[0]}) : bcmp(${$_[0]},$_[1])},
  8. 'cmp'    =>    sub {$_[2]? ($_[1] cmp ${$_[0]}) : (${$_[0]} cmp $_[1])},
  9. '*'    =>    sub {new Math::BigInt &bmul},
  10. '/'    =>    sub {new Math::BigInt 
  11.                $_[2]? scalar bdiv($_[1],${$_[0]}) :
  12.              scalar bdiv(${$_[0]},$_[1])},
  13. '%'    =>    sub {new Math::BigInt
  14.                $_[2]? bmod($_[1],${$_[0]}) : bmod(${$_[0]},$_[1])},
  15. '**'    =>    sub {new Math::BigInt
  16.                $_[2]? bpow($_[1],${$_[0]}) : bpow(${$_[0]},$_[1])},
  17. 'neg'    =>    sub {new Math::BigInt &bneg},
  18. 'abs'    =>    sub {new Math::BigInt &babs},
  19. '<<'    =>    sub {new Math::BigInt
  20.                $_[2]? blsft($_[1],${$_[0]}) : blsft(${$_[0]},$_[1])},
  21. '>>'    =>    sub {new Math::BigInt
  22.                $_[2]? brsft($_[1],${$_[0]}) : brsft(${$_[0]},$_[1])},
  23. '&'    =>    sub {new Math::BigInt &band},
  24. '|'    =>    sub {new Math::BigInt &bior},
  25. '^'    =>    sub {new Math::BigInt &bxor},
  26. '~'    =>    sub {new Math::BigInt &bnot},
  27.  
  28. qw(
  29. ""    stringify
  30. 0+    numify)            # Order of arguments unsignificant
  31. ;
  32.  
  33. $NaNOK=1;
  34.  
  35. sub new {
  36.   my($class) = shift;
  37.   my($foo) = bnorm(shift);
  38.   die "Not a number initialized to Math::BigInt" if !$NaNOK && $foo eq "NaN";
  39.   bless \$foo, $class;
  40. }
  41. sub stringify { "${$_[0]}" }
  42. sub numify { 0 + "${$_[0]}" }    # Not needed, additional overhead
  43.                 # comparing to direct compilation based on
  44.                 # stringify
  45. sub import {
  46.   shift;
  47.   return unless @_;
  48.   die "unknown import: @_" unless @_ == 1 and $_[0] eq ':constant';
  49.   overload::constant integer => sub {Math::BigInt->new(shift)};
  50. }
  51.  
  52. $zero = 0;
  53.  
  54.  
  55. # normalize string form of number.   Strip leading zeros.  Strip any
  56. #   white space and add a sign, if missing.
  57. # Strings that are not numbers result the value 'NaN'.
  58.  
  59. sub bnorm { #(num_str) return num_str
  60.     local($_) = @_;
  61.     s/\s+//g;                           # strip white space
  62.     if (s/^([+-]?)0*(\d+)$/$1$2/) {     # test if number
  63.     substr($_,$[,0) = '+' unless $1; # Add missing sign
  64.     s/^-0/+0/;
  65.     $_;
  66.     } else {
  67.     'NaN';
  68.     }
  69. }
  70.  
  71. # Convert a number from string format to internal base 100000 format.
  72. #   Assumes normalized value as input.
  73. sub internal { #(num_str) return int_num_array
  74.     local($d) = @_;
  75.     ($is,$il) = (substr($d,$[,1),length($d)-2);
  76.     substr($d,$[,1) = '';
  77.     ($is, reverse(unpack("a" . ($il%5+1) . ("a5" x ($il/5)), $d)));
  78. }
  79.  
  80. # Convert a number from internal base 100000 format to string format.
  81. #   This routine scribbles all over input array.
  82. sub external { #(int_num_array) return num_str
  83.     $es = shift;
  84.     grep($_ > 9999 || ($_ = substr('0000'.$_,-5)), @_);   # zero pad
  85.     &bnorm(join('', $es, reverse(@_)));    # reverse concat and normalize
  86. }
  87.  
  88. # Negate input value.
  89. sub bneg { #(num_str) return num_str
  90.     local($_) = &bnorm(@_);
  91.     return $_ if $_ eq '+0' or $_ eq 'NaN';
  92.     vec($_,0,8) ^= ord('+') ^ ord('-');
  93.     $_;
  94. }
  95.  
  96. # Returns the absolute value of the input.
  97. sub babs { #(num_str) return num_str
  98.     &abs(&bnorm(@_));
  99. }
  100.  
  101. sub abs { # post-normalized abs for internal use
  102.     local($_) = @_;
  103.     s/^-/+/;
  104.     $_;
  105. }
  106.  
  107. # Compares 2 values.  Returns one of undef, <0, =0, >0. (suitable for sort)
  108. sub bcmp { #(num_str, num_str) return cond_code
  109.     local($x,$y) = (&bnorm($_[$[]),&bnorm($_[$[+1]));
  110.     if ($x eq 'NaN') {
  111.     undef;
  112.     } elsif ($y eq 'NaN') {
  113.     undef;
  114.     } else {
  115.     &cmp($x,$y) <=> 0;
  116.     }
  117. }
  118.  
  119. sub cmp { # post-normalized compare for internal use
  120.     local($cx, $cy) = @_;
  121.     
  122.     return 0 if ($cx eq $cy);
  123.  
  124.     local($sx, $sy) = (substr($cx, 0, 1), substr($cy, 0, 1));
  125.     local($ld);
  126.  
  127.     if ($sx eq '+') {
  128.       return  1 if ($sy eq '-' || $cy eq '+0');
  129.       $ld = length($cx) - length($cy);
  130.       return $ld if ($ld);
  131.       return $cx cmp $cy;
  132.     } else { # $sx eq '-'
  133.       return -1 if ($sy eq '+');
  134.       $ld = length($cy) - length($cx);
  135.       return $ld if ($ld);
  136.       return $cy cmp $cx;
  137.     }
  138. }
  139.  
  140. sub badd { #(num_str, num_str) return num_str
  141.     local(*x, *y); ($x, $y) = (&bnorm($_[$[]),&bnorm($_[$[+1]));
  142.     if ($x eq 'NaN') {
  143.     'NaN';
  144.     } elsif ($y eq 'NaN') {
  145.     'NaN';
  146.     } else {
  147.     @x = &internal($x);             # convert to internal form
  148.     @y = &internal($y);
  149.     local($sx, $sy) = (shift @x, shift @y); # get signs
  150.     if ($sx eq $sy) {
  151.         &external($sx, &add(*x, *y)); # if same sign add
  152.     } else {
  153.         ($x, $y) = (&abs($x),&abs($y)); # make abs
  154.         if (&cmp($y,$x) > 0) {
  155.         &external($sy, &sub(*y, *x));
  156.         } else {
  157.         &external($sx, &sub(*x, *y));
  158.         }
  159.     }
  160.     }
  161. }
  162.  
  163. sub bsub { #(num_str, num_str) return num_str
  164.     &badd($_[$[],&bneg($_[$[+1]));    
  165. }
  166.  
  167. # GCD -- Euclids algorithm Knuth Vol 2 pg 296
  168. sub bgcd { #(num_str, num_str) return num_str
  169.     local($x,$y) = (&bnorm($_[$[]),&bnorm($_[$[+1]));
  170.     if ($x eq 'NaN' || $y eq 'NaN') {
  171.     'NaN';
  172.     } else {
  173.     ($x, $y) = ($y,&bmod($x,$y)) while $y ne '+0';
  174.     $x;
  175.     }
  176. }
  177.  
  178. # routine to add two base 1e5 numbers
  179. #   stolen from Knuth Vol 2 Algorithm A pg 231
  180. #   there are separate routines to add and sub as per Kunth pg 233
  181. sub add { #(int_num_array, int_num_array) return int_num_array
  182.     local(*x, *y) = @_;
  183.     $car = 0;
  184.     for $x (@x) {
  185.     last unless @y || $car;
  186.     $x -= 1e5 if $car = (($x += (@y ? shift(@y) : 0) + $car) >= 1e5) ? 1 : 0;
  187.     }
  188.     for $y (@y) {
  189.     last unless $car;
  190.     $y -= 1e5 if $car = (($y += $car) >= 1e5) ? 1 : 0;
  191.     }
  192.     (@x, @y, $car);
  193. }
  194.  
  195. # subtract base 1e5 numbers -- stolen from Knuth Vol 2 pg 232, $x > $y
  196. sub sub { #(int_num_array, int_num_array) return int_num_array
  197.     local(*sx, *sy) = @_;
  198.     $bar = 0;
  199.     for $sx (@sx) {
  200.     last unless @sy || $bar;
  201.     $sx += 1e5 if $bar = (($sx -= (@sy ? shift(@sy) : 0) + $bar) < 0);
  202.     }
  203.     @sx;
  204. }
  205.  
  206. # multiply two numbers -- stolen from Knuth Vol 2 pg 233
  207. sub bmul { #(num_str, num_str) return num_str
  208.     local(*x, *y); ($x, $y) = (&bnorm($_[$[]), &bnorm($_[$[+1]));
  209.     if ($x eq 'NaN') {
  210.     'NaN';
  211.     } elsif ($y eq 'NaN') {
  212.     'NaN';
  213.     } else {
  214.     @x = &internal($x);
  215.     @y = &internal($y);
  216.     &external(&mul(*x,*y));
  217.     }
  218. }
  219.  
  220. # multiply two numbers in internal representation
  221. # destroys the arguments, supposes that two arguments are different
  222. sub mul { #(*int_num_array, *int_num_array) return int_num_array
  223.     local(*x, *y) = (shift, shift);
  224.     local($signr) = (shift @x ne shift @y) ? '-' : '+';
  225.     @prod = ();
  226.     for $x (@x) {
  227.       ($car, $cty) = (0, $[);
  228.       for $y (@y) {
  229.     $prod = $x * $y + ($prod[$cty] || 0) + $car;
  230.     $prod[$cty++] =
  231.       $prod - ($car = int($prod * 1e-5)) * 1e5;
  232.       }
  233.       $prod[$cty] += $car if $car;
  234.       $x = shift @prod;
  235.     }
  236.     ($signr, @x, @prod);
  237. }
  238.  
  239. # modulus
  240. sub bmod { #(num_str, num_str) return num_str
  241.     (&bdiv(@_))[$[+1];
  242. }
  243.  
  244. sub bdiv { #(dividend: num_str, divisor: num_str) return num_str
  245.     local (*x, *y); ($x, $y) = (&bnorm($_[$[]), &bnorm($_[$[+1]));
  246.     return wantarray ? ('NaN','NaN') : 'NaN'
  247.     if ($x eq 'NaN' || $y eq 'NaN' || $y eq '+0');
  248.     return wantarray ? ('+0',$x) : '+0' if (&cmp(&abs($x),&abs($y)) < 0);
  249.     @x = &internal($x); @y = &internal($y);
  250.     $srem = $y[$[];
  251.     $sr = (shift @x ne shift @y) ? '-' : '+';
  252.     $car = $bar = $prd = 0;
  253.     if (($dd = int(1e5/($y[$#y]+1))) != 1) {
  254.     for $x (@x) {
  255.         $x = $x * $dd + $car;
  256.         $x -= ($car = int($x * 1e-5)) * 1e5;
  257.     }
  258.     push(@x, $car); $car = 0;
  259.     for $y (@y) {
  260.         $y = $y * $dd + $car;
  261.         $y -= ($car = int($y * 1e-5)) * 1e5;
  262.     }
  263.     }
  264.     else {
  265.     push(@x, 0);
  266.     }
  267.     @q = (); ($v2,$v1) = @y[-2,-1];
  268.     $v2 = 0 unless $v2;
  269.     while ($#x > $#y) {
  270.     ($u2,$u1,$u0) = @x[-3..-1];
  271.     $u2 = 0 unless $u2;
  272.     $q = (($u0 == $v1) ? 99999 : int(($u0*1e5+$u1)/$v1));
  273.     --$q while ($v2*$q > ($u0*1e5+$u1-$q*$v1)*1e5+$u2);
  274.     if ($q) {
  275.         ($car, $bar) = (0,0);
  276.         for ($y = $[, $x = $#x-$#y+$[-1; $y <= $#y; ++$y,++$x) {
  277.         $prd = $q * $y[$y] + $car;
  278.         $prd -= ($car = int($prd * 1e-5)) * 1e5;
  279.         $x[$x] += 1e5 if ($bar = (($x[$x] -= $prd + $bar) < 0));
  280.         }
  281.         if ($x[$#x] < $car + $bar) {
  282.         $car = 0; --$q;
  283.         for ($y = $[, $x = $#x-$#y+$[-1; $y <= $#y; ++$y,++$x) {
  284.             $x[$x] -= 1e5
  285.             if ($car = (($x[$x] += $y[$y] + $car) > 1e5));
  286.         }
  287.         }   
  288.     }
  289.     pop(@x); unshift(@q, $q);
  290.     }
  291.     if (wantarray) {
  292.     @d = ();
  293.     if ($dd != 1) {
  294.         $car = 0;
  295.         for $x (reverse @x) {
  296.         $prd = $car * 1e5 + $x;
  297.         $car = $prd - ($tmp = int($prd / $dd)) * $dd;
  298.         unshift(@d, $tmp);
  299.         }
  300.     }
  301.     else {
  302.         @d = @x;
  303.     }
  304.     (&external($sr, @q), &external($srem, @d, $zero));
  305.     } else {
  306.     &external($sr, @q);
  307.     }
  308. }
  309.  
  310. # compute power of two numbers -- stolen from Knuth Vol 2 pg 233
  311. sub bpow { #(num_str, num_str) return num_str
  312.     local(*x, *y); ($x, $y) = (&bnorm($_[$[]), &bnorm($_[$[+1]));
  313.     if ($x eq 'NaN') {
  314.     'NaN';
  315.     } elsif ($y eq 'NaN') {
  316.     'NaN';
  317.     } elsif ($x eq '+1') {
  318.     '+1';
  319.     } elsif ($x eq '-1') {
  320.     &bmod($x,2) ? '-1': '+1';
  321.     } elsif ($y =~ /^-/) {
  322.     'NaN';
  323.     } elsif ($x eq '+0' && $y eq '+0') {
  324.     'NaN';
  325.     } else {
  326.     @x = &internal($x);
  327.     local(@pow2)=@x;
  328.     local(@pow)=&internal("+1");
  329.     local($y1,$res,@tmp1,@tmp2)=(1); # need tmp to send to mul
  330.     while ($y ne '+0') {
  331.       ($y,$res)=&bdiv($y,2);
  332.       if ($res ne '+0') {@tmp=@pow2; @pow=&mul(*pow,*tmp);}
  333.       if ($y ne '+0') {@tmp=@pow2;@pow2=&mul(*pow2,*tmp);}
  334.     }
  335.     &external(@pow);
  336.     }
  337. }
  338.  
  339. # compute x << y, y >= 0
  340. sub blsft { #(num_str, num_str) return num_str
  341.     &bmul($_[$[], &bpow(2, $_[$[+1]));
  342. }
  343.  
  344. # compute x >> y, y >= 0
  345. sub brsft { #(num_str, num_str) return num_str
  346.     &bdiv($_[$[], &bpow(2, $_[$[+1]));
  347. }
  348.  
  349. # compute x & y
  350. sub band { #(num_str, num_str) return num_str
  351.     local($x,$y,$r,$m,$xr,$yr) = (&bnorm($_[$[]),&bnorm($_[$[+1]),0,1);
  352.     if ($x eq 'NaN' || $y eq 'NaN') {
  353.     'NaN';
  354.     } else {
  355.     while ($x ne '+0' && $y ne '+0') {
  356.         ($x, $xr) = &bdiv($x, 0x10000);
  357.         ($y, $yr) = &bdiv($y, 0x10000);
  358.         $r = &badd(&bmul(int $xr & $yr, $m), $r);
  359.         $m = &bmul($m, 0x10000);
  360.     }
  361.     $r;
  362.     }
  363. }
  364.  
  365. # compute x | y
  366. sub bior { #(num_str, num_str) return num_str
  367.     local($x,$y,$r,$m,$xr,$yr) = (&bnorm($_[$[]),&bnorm($_[$[+1]),0,1);
  368.     if ($x eq 'NaN' || $y eq 'NaN') {
  369.     'NaN';
  370.     } else {
  371.     while ($x ne '+0' || $y ne '+0') {
  372.         ($x, $xr) = &bdiv($x, 0x10000);
  373.         ($y, $yr) = &bdiv($y, 0x10000);
  374.         $r = &badd(&bmul(int $xr | $yr, $m), $r);
  375.         $m = &bmul($m, 0x10000);
  376.     }
  377.     $r;
  378.     }
  379. }
  380.  
  381. # compute x ^ y
  382. sub bxor { #(num_str, num_str) return num_str
  383.     local($x,$y,$r,$m,$xr,$yr) = (&bnorm($_[$[]),&bnorm($_[$[+1]),0,1);
  384.     if ($x eq 'NaN' || $y eq 'NaN') {
  385.     'NaN';
  386.     } else {
  387.     while ($x ne '+0' || $y ne '+0') {
  388.         ($x, $xr) = &bdiv($x, 0x10000);
  389.         ($y, $yr) = &bdiv($y, 0x10000);
  390.         $r = &badd(&bmul(int $xr ^ $yr, $m), $r);
  391.         $m = &bmul($m, 0x10000);
  392.     }
  393.     $r;
  394.     }
  395. }
  396.  
  397. # represent ~x as twos-complement number
  398. sub bnot { #(num_str) return num_str
  399.     &bsub(-1,$_[$[]);
  400. }
  401.  
  402. 1;
  403. __END__
  404.  
  405. =head1 NAME
  406.  
  407. Math::BigInt - Arbitrary size integer math package
  408.  
  409. =head1 SYNOPSIS
  410.  
  411.   use Math::BigInt;
  412.   $i = Math::BigInt->new($string);
  413.  
  414.   $i->bneg return BINT               negation
  415.   $i->babs return BINT               absolute value
  416.   $i->bcmp(BINT) return CODE         compare numbers (undef,<0,=0,>0)
  417.   $i->badd(BINT) return BINT         addition
  418.   $i->bsub(BINT) return BINT         subtraction
  419.   $i->bmul(BINT) return BINT         multiplication
  420.   $i->bdiv(BINT) return (BINT,BINT)  division (quo,rem) just quo if scalar
  421.   $i->bmod(BINT) return BINT         modulus
  422.   $i->bgcd(BINT) return BINT         greatest common divisor
  423.   $i->bnorm return BINT              normalization
  424.   $i->blsft(BINT) return BINT        left shift
  425.   $i->brsft(BINT) return (BINT,BINT) right shift (quo,rem) just quo if scalar
  426.   $i->band(BINT) return BINT         bit-wise and
  427.   $i->bior(BINT) return BINT         bit-wise inclusive or
  428.   $i->bxor(BINT) return BINT         bit-wise exclusive or
  429.   $i->bnot return BINT               bit-wise not
  430.  
  431. =head1 DESCRIPTION
  432.  
  433. All basic math operations are overloaded if you declare your big
  434. integers as
  435.  
  436.   $i = new Math::BigInt '123 456 789 123 456 789';
  437.  
  438.  
  439. =over 2
  440.  
  441. =item Canonical notation
  442.  
  443. Big integer value are strings of the form C</^[+-]\d+$/> with leading
  444. zeros suppressed.
  445.  
  446. =item Input
  447.  
  448. Input values to these routines may be strings of the form
  449. C</^\s*[+-]?[\d\s]+$/>.
  450.  
  451. =item Output
  452.  
  453. Output values always always in canonical form
  454.  
  455. =back
  456.  
  457. Actual math is done in an internal format consisting of an array
  458. whose first element is the sign (/^[+-]$/) and whose remaining 
  459. elements are base 100000 digits with the least significant digit first.
  460. The string 'NaN' is used to represent the result when input arguments 
  461. are not numbers, as well as the result of dividing by zero.
  462.  
  463. =head1 EXAMPLES
  464.  
  465.    '+0'                            canonical zero value
  466.    '   -123 123 123'               canonical value '-123123123'
  467.    '1 23 456 7890'                 canonical value '+1234567890'
  468.  
  469.  
  470. =head1 Autocreating constants
  471.  
  472. After C<use Math::BigInt ':constant'> all the integer decimal constants
  473. in the given scope are converted to C<Math::BigInt>.  This conversion
  474. happens at compile time.
  475.  
  476. In particular
  477.  
  478.   perl -MMath::BigInt=:constant -e 'print 2**100'
  479.  
  480. print the integer value of C<2**100>.  Note that without conversion of 
  481. constants the expression 2**100 will be calculated as floating point number.
  482.  
  483. =head1 BUGS
  484.  
  485. The current version of this module is a preliminary version of the
  486. real thing that is currently (as of perl5.002) under development.
  487.  
  488. =head1 AUTHOR
  489.  
  490. Mark Biggar, overloaded interface by Ilya Zakharevich.
  491.  
  492. =cut
  493.