lib/Math/BigInt.pm

   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 += (@y ? shift(@y) : 0) + $car) >= 1e5) ? 1 : 0;
 175     }
 176     for $y (@y) {
 177         last unless $car;
 178         $y -= 1e5 if $car = (($y += $car) >= 1e5) ? 1 : 0;
 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 @sy || $bar;
 189         $sx += 1e5 if $bar = (($sx -= (@sy ? shift(@sy) : 0) + $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] || 0) + $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