home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / perl502b.zip / lib / Math / BigInt.pm < prev    next >
Text File  |  1995-01-23  |  10KB  |  350 lines
  1. package Math::BigInt;
  2.  
  3. %OVERLOAD = ( 
  4.                 # Anonymous subroutines:
  5. '+'    =>    sub {new Math::BigInt &badd},
  6. '-'    =>    sub {new Math::BigInt
  7.                $_[2]? bsub($_[1],${$_[0]}) : bsub(${$_[0]},$_[1])},
  8. '<=>'    =>    sub {new Math::BigInt
  9.                $_[2]? bcmp($_[1],${$_[0]}) : bcmp(${$_[0]},$_[1])},
  10. 'cmp'    =>    sub {new Math::BigInt
  11.                $_[2]? ($_[1] cmp ${$_[0]}) : (${$_[0]} cmp $_[1])},
  12. '*'    =>    sub {new Math::BigInt &bmul},
  13. '/'    =>    sub {new Math::BigInt 
  14.                $_[2]? scalar bdiv($_[1],${$_[0]}) :
  15.              scalar bdiv(${$_[0]},$_[1])},
  16. '%'    =>    sub {new Math::BigInt
  17.                $_[2]? bmod($_[1],${$_[0]}) : bmod(${$_[0]},$_[1])},
  18. '**'    =>    sub {new Math::BigInt
  19.                $_[2]? bpow($_[1],${$_[0]}) : bpow(${$_[0]},$_[1])},
  20. 'neg'    =>    sub {new Math::BigInt &bneg},
  21. 'abs'    =>    sub {new Math::BigInt &babs},
  22.  
  23. qw(
  24. ""    stringify
  25. 0+    numify)            # Order of arguments unsignificant
  26. );
  27.  
  28. $NaNOK=1;
  29.  
  30. sub new {
  31.   my $foo = bnorm($_[1]);
  32.   die "Not a number initialized to Math::BigInt" if !$NaNOK && $foo eq "NaN";
  33.   bless \$foo;
  34. }
  35. sub stringify { "${$_[0]}" }
  36. sub numify { 0 + "${$_[0]}" }    # Not needed, additional overhead
  37.                 # comparing to direct compilation based on
  38.                 # stringify
  39.  
  40. # arbitrary size integer math package
  41. #
  42. # by Mark Biggar
  43. #
  44. # Canonical Big integer value are strings of the form
  45. #       /^[+-]\d+$/ with leading zeros suppressed
  46. # Input values to these routines may be strings of the form
  47. #       /^\s*[+-]?[\d\s]+$/.
  48. # Examples:
  49. #   '+0'                            canonical zero value
  50. #   '   -123 123 123'               canonical value '-123123123'
  51. #   '1 23 456 7890'                 canonical value '+1234567890'
  52. # Output values always always in canonical form
  53. #
  54. # Actual math is done in an internal format consisting of an array
  55. #   whose first element is the sign (/^[+-]$/) and whose remaining 
  56. #   elements are base 100000 digits with the least significant digit first.
  57. # The string 'NaN' is used to represent the result when input arguments 
  58. #   are not numbers, as well as the result of dividing by zero
  59. #
  60. # routines provided are:
  61. #
  62. #   bneg(BINT) return BINT              negation
  63. #   babs(BINT) return BINT              absolute value
  64. #   bcmp(BINT,BINT) return CODE         compare numbers (undef,<0,=0,>0)
  65. #   badd(BINT,BINT) return BINT         addition
  66. #   bsub(BINT,BINT) return BINT         subtraction
  67. #   bmul(BINT,BINT) return BINT         multiplication
  68. #   bdiv(BINT,BINT) return (BINT,BINT)  division (quo,rem) just quo if scalar
  69. #   bmod(BINT,BINT) return BINT         modulus
  70. #   bgcd(BINT,BINT) return BINT         greatest common divisor
  71. #   bnorm(BINT) return BINT             normalization
  72. #
  73.  
  74. $zero = 0;
  75.  
  76.  
  77. # normalize string form of number.   Strip leading zeros.  Strip any
  78. #   white space and add a sign, if missing.
  79. # Strings that are not numbers result the value 'NaN'.
  80.  
  81. sub bnorm { #(num_str) return num_str
  82.     local($_) = @_;
  83.     s/\s+//g;                           # strip white space
  84.     if (s/^([+-]?)0*(\d+)$/$1$2/) {     # test if number
  85.     substr($_,$[,0) = '+' unless $1; # Add missing sign
  86.     s/^-0/+0/;
  87.     $_;
  88.     } else {
  89.     'NaN';
  90.     }
  91. }
  92.  
  93. # Convert a number from string format to internal base 100000 format.
  94. #   Assumes normalized value as input.
  95. sub internal { #(num_str) return int_num_array
  96.     local($d) = @_;
  97.     ($is,$il) = (substr($d,$[,1),length($d)-2);
  98.     substr($d,$[,1) = '';
  99.     ($is, reverse(unpack("a" . ($il%5+1) . ("a5" x ($il/5)), $d)));
  100. }
  101.  
  102. # Convert a number from internal base 100000 format to string format.
  103. #   This routine scribbles all over input array.
  104. sub external { #(int_num_array) return num_str
  105.     $es = shift;
  106.     grep($_ > 9999 || ($_ = substr('0000'.$_,-5)), @_);   # zero pad
  107.     &bnorm(join('', $es, reverse(@_)));    # reverse concat and normalize
  108. }
  109.  
  110. # Negate input value.
  111. sub bneg { #(num_str) return num_str
  112.     local($_) = &bnorm(@_);
  113.     vec($_,0,8) ^= ord('+') ^ ord('-') unless $_ eq '+0';
  114.     s/^H/N/;
  115.     $_;
  116. }
  117.  
  118. # Returns the absolute value of the input.
  119. sub babs { #(num_str) return num_str
  120.     &abs(&bnorm(@_));
  121. }
  122.  
  123. sub abs { # post-normalized abs for internal use
  124.     local($_) = @_;
  125.     s/^-/+/;
  126.     $_;
  127. }
  128.  
  129. # Compares 2 values.  Returns one of undef, <0, =0, >0. (suitable for sort)
  130. sub bcmp { #(num_str, num_str) return cond_code
  131.     local($x,$y) = (&bnorm($_[$[]),&bnorm($_[$[+1]));
  132.     if ($x eq 'NaN') {
  133.     undef;
  134.     } elsif ($y eq 'NaN') {
  135.     undef;
  136.     } else {
  137.     &cmp($x,$y);
  138.     }
  139. }
  140.  
  141. sub cmp { # post-normalized compare for internal use
  142.     local($cx, $cy) = @_;
  143.     $cx cmp $cy
  144.     &&
  145.     (
  146.     ord($cy) <=> ord($cx)
  147.     ||
  148.     ($cx cmp ',') * (length($cy) <=> length($cx) || $cy cmp $cx)
  149.     );
  150. }
  151.  
  152. sub badd { #(num_str, num_str) return num_str
  153.     local(*x, *y); ($x, $y) = (&bnorm($_[$[]),&bnorm($_[$[+1]));
  154.     if ($x eq 'NaN') {
  155.     'NaN';
  156.     } elsif ($y eq 'NaN') {
  157.     'NaN';
  158.     } else {
  159.     @x = &internal($x);             # convert to internal form
  160.     @y = &internal($y);
  161.     local($sx, $sy) = (shift @x, shift @y); # get signs
  162.     if ($sx eq $sy) {
  163.         &external($sx, &add(*x, *y)); # if same sign add
  164.     } else {
  165.         ($x, $y) = (&abs($x),&abs($y)); # make abs
  166.         if (&cmp($y,$x) > 0) {
  167.         &external($sy, &sub(*y, *x));
  168.         } else {
  169.         &external($sx, &sub(*x, *y));
  170.         }
  171.     }
  172.     }
  173. }
  174.  
  175. sub bsub { #(num_str, num_str) return num_str
  176.     &badd($_[$[],&bneg($_[$[+1]));    
  177. }
  178.  
  179. # GCD -- Euclids algorithm Knuth Vol 2 pg 296
  180. sub bgcd { #(num_str, num_str) return num_str
  181.     local($x,$y) = (&bnorm($_[$[]),&bnorm($_[$[+1]));
  182.     if ($x eq 'NaN' || $y eq 'NaN') {
  183.     'NaN';
  184.     } else {
  185.     ($x, $y) = ($y,&bmod($x,$y)) while $y ne '+0';
  186.     $x;
  187.     }
  188. }
  189.  
  190. # routine to add two base 1e5 numbers
  191. #   stolen from Knuth Vol 2 Algorithm A pg 231
  192. #   there are separate routines to add and sub as per Kunth pg 233
  193. sub add { #(int_num_array, int_num_array) return int_num_array
  194.     local(*x, *y) = @_;
  195.     $car = 0;
  196.     for $x (@x) {
  197.     last unless @y || $car;
  198.     $x -= 1e5 if $car = (($x += shift(@y) + $car) >= 1e5);
  199.     }
  200.     for $y (@y) {
  201.     last unless $car;
  202.     $y -= 1e5 if $car = (($y += $car) >= 1e5);
  203.     }
  204.     (@x, @y, $car);
  205. }
  206.  
  207. # subtract base 1e5 numbers -- stolen from Knuth Vol 2 pg 232, $x > $y
  208. sub sub { #(int_num_array, int_num_array) return int_num_array
  209.     local(*sx, *sy) = @_;
  210.     $bar = 0;
  211.     for $sx (@sx) {
  212.     last unless @y || $bar;
  213.     $sx += 1e5 if $bar = (($sx -= shift(@sy) + $bar) < 0);
  214.     }
  215.     @sx;
  216. }
  217.  
  218. # multiply two numbers -- stolen from Knuth Vol 2 pg 233
  219. sub bmul { #(num_str, num_str) return num_str
  220.     local(*x, *y); ($x, $y) = (&bnorm($_[$[]), &bnorm($_[$[+1]));
  221.     if ($x eq 'NaN') {
  222.     'NaN';
  223.     } elsif ($y eq 'NaN') {
  224.     'NaN';
  225.     } else {
  226.     @x = &internal($x);
  227.     @y = &internal($y);
  228.     &external(&mul(*x,*y));
  229.     }
  230. }
  231.  
  232. # multiply two numbers in internal representation
  233. # destroys the arguments, supposes that two arguments are different
  234. sub mul { #(*int_num_array, *int_num_array) return int_num_array
  235.     local(*x, *y) = (shift, shift);
  236.     local($signr) = (shift @x ne shift @y) ? '-' : '+';
  237.     @prod = ();
  238.     for $x (@x) {
  239.       ($car, $cty) = (0, $[);
  240.       for $y (@y) {
  241.     $prod = $x * $y + $prod[$cty] + $car;
  242.     $prod[$cty++] =
  243.       $prod - ($car = int($prod * 1e-5)) * 1e5;
  244.       }
  245.       $prod[$cty] += $car if $car;
  246.       $x = shift @prod;
  247.     }
  248.     ($signr, @x, @prod);
  249. }
  250.  
  251. # modulus
  252. sub bmod { #(num_str, num_str) return num_str
  253.     (&bdiv(@_))[$[+1];
  254. }
  255.  
  256. sub bdiv { #(dividend: num_str, divisor: num_str) return num_str
  257.     local (*x, *y); ($x, $y) = (&bnorm($_[$[]), &bnorm($_[$[+1]));
  258.     return wantarray ? ('NaN','NaN') : 'NaN'
  259.     if ($x eq 'NaN' || $y eq 'NaN' || $y eq '+0');
  260.     return wantarray ? ('+0',$x) : '+0' if (&cmp(&abs($x),&abs($y)) < 0);
  261.     @x = &internal($x); @y = &internal($y);
  262.     $srem = $y[$[];
  263.     $sr = (shift @x ne shift @y) ? '-' : '+';
  264.     $car = $bar = $prd = 0;
  265.     if (($dd = int(1e5/($y[$#y]+1))) != 1) {
  266.     for $x (@x) {
  267.         $x = $x * $dd + $car;
  268.         $x -= ($car = int($x * 1e-5)) * 1e5;
  269.     }
  270.     push(@x, $car); $car = 0;
  271.     for $y (@y) {
  272.         $y = $y * $dd + $car;
  273.         $y -= ($car = int($y * 1e-5)) * 1e5;
  274.     }
  275.     }
  276.     else {
  277.     push(@x, 0);
  278.     }
  279.     @q = (); ($v2,$v1) = @y[-2,-1];
  280.     while ($#x > $#y) {
  281.     ($u2,$u1,$u0) = @x[-3..-1];
  282.     $q = (($u0 == $v1) ? 99999 : int(($u0*1e5+$u1)/$v1));
  283.     --$q while ($v2*$q > ($u0*1e5+$u1-$q*$v1)*1e5+$u2);
  284.     if ($q) {
  285.         ($car, $bar) = (0,0);
  286.         for ($y = $[, $x = $#x-$#y+$[-1; $y <= $#y; ++$y,++$x) {
  287.         $prd = $q * $y[$y] + $car;
  288.         $prd -= ($car = int($prd * 1e-5)) * 1e5;
  289.         $x[$x] += 1e5 if ($bar = (($x[$x] -= $prd + $bar) < 0));
  290.         }
  291.         if ($x[$#x] < $car + $bar) {
  292.         $car = 0; --$q;
  293.         for ($y = $[, $x = $#x-$#y+$[-1; $y <= $#y; ++$y,++$x) {
  294.             $x[$x] -= 1e5
  295.             if ($car = (($x[$x] += $y[$y] + $car) > 1e5));
  296.         }
  297.         }   
  298.     }
  299.     pop(@x); unshift(@q, $q);
  300.     }
  301.     if (wantarray) {
  302.     @d = ();
  303.     if ($dd != 1) {
  304.         $car = 0;
  305.         for $x (reverse @x) {
  306.         $prd = $car * 1e5 + $x;
  307.         $car = $prd - ($tmp = int($prd / $dd)) * $dd;
  308.         unshift(@d, $tmp);
  309.         }
  310.     }
  311.     else {
  312.         @d = @x;
  313.     }
  314.     (&external($sr, @q), &external($srem, @d, $zero));
  315.     } else {
  316.     &external($sr, @q);
  317.     }
  318. }
  319.  
  320. # compute power of two numbers -- stolen from Knuth Vol 2 pg 233
  321. sub bpow { #(num_str, num_str) return num_str
  322.     local(*x, *y); ($x, $y) = (&bnorm($_[$[]), &bnorm($_[$[+1]));
  323.     if ($x eq 'NaN') {
  324.     'NaN';
  325.     } elsif ($y eq 'NaN') {
  326.     'NaN';
  327.     } elsif ($x eq '+1') {
  328.     '+1';
  329.     } elsif ($x eq '-1') {
  330.     &bmod($x,2) ? '-1': '+1';
  331.     } elsif ($y =~ /^-/) {
  332.     'NaN';
  333.     } elsif ($x eq '+0' && $y eq '+0') {
  334.     'NaN';
  335.     } else {
  336.     @x = &internal($x);
  337.     local(@pow2)=@x;
  338.     local(@pow)=&internal("+1");
  339.     local($y1,$res,@tmp1,@tmp2)=(1); # need tmp to send to mul
  340.     while ($y ne '+0') {
  341.       ($y,$res)=&bdiv($y,2);
  342.       if ($res ne '+0') {@tmp=@pow2; @pow=&mul(*pow,*tmp);}
  343.       if ($y ne '+0') {@tmp=@pow2;@pow2=&mul(*pow2,*tmp);}
  344.     }
  345.     &external(@pow);
  346.     }
  347. }
  348.  
  349. 1;
  350.