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