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     $cx cmp $cy
 110     &&
 111     (
 112         ord($cy) <=> ord($cx)
 113         ||
 114         ($cx cmp ',') * (length($cy) <=> length($cx) || $cy cmp $cx)
 115     );
 116 }
 117
 118 sub badd { #(num_str, num_str) return num_str
 119     local(*x, *y); ($x, $y) = (&bnorm($_[$[]),&bnorm($_[$[+1]));
 120     if ($x eq 'NaN') {
 121         'NaN';
 122     } elsif ($y eq 'NaN') {
 123         'NaN';
 124     } else {
 125         @x = &internal($x);             # convert to internal form
 126         @y = &internal($y);
 127         local($sx, $sy) = (shift @x, shift @y); # get signs
 128         if ($sx eq $sy) {
 129             &external($sx, &add(*x, *y)); # if same sign add
 130         } else {
 131             ($x, $y) = (&abs($x),&abs($y)); # make abs
 132             if (&cmp($y,$x) > 0) {
 133                 &external($sy, &sub(*y, *x));
 134             } else {
 135                 &external($sx, &sub(*x, *y));
 136             }
 137         }
 138     }
 139 }
 140
 141 sub bsub { #(num_str, num_str) return num_str
 142     &badd($_[$[],&bneg($_[$[+1]));
 143 }
 144
 145 # GCD -- Euclids algorithm Knuth Vol 2 pg 296
 146 sub bgcd { #(num_str, num_str) return num_str
 147     local($x,$y) = (&bnorm($_[$[]),&bnorm($_[$[+1]));
 148     if ($x eq 'NaN' || $y eq 'NaN') {
 149         'NaN';
 150     } else {
 151         ($x, $y) = ($y,&bmod($x,$y)) while $y ne '+0';
 152         $x;
 153     }
 154 }
 155
 156 # routine to add two base 1e5 numbers
 157 #   stolen from Knuth Vol 2 Algorithm A pg 231
 158 #   there are separate routines to add and sub as per Kunth pg 233
 159 sub add { #(int_num_array, int_num_array) return int_num_array
 160     local(*x, *y) = @_;
 161     $car = 0;
 162     for $x (@x) {
 163         last unless @y || $car;
 164         $x -= 1e5 if $car = (($x += shift(@y) + $car) >= 1e5);
 165     }
 166     for $y (@y) {
 167         last unless $car;
 168         $y -= 1e5 if $car = (($y += $car) >= 1e5);
 169     }
 170     (@x, @y, $car);
 171 }
 172
 173 # subtract base 1e5 numbers -- stolen from Knuth Vol 2 pg 232, $x > $y
 174 sub sub { #(int_num_array, int_num_array) return int_num_array
 175     local(*sx, *sy) = @_;
 176     $bar = 0;
 177     for $sx (@sx) {
 178         last unless @y || $bar;
 179         $sx += 1e5 if $bar = (($sx -= shift(@sy) + $bar) < 0);
 180     }
 181     @sx;
 182 }
 183
 184 # multiply two numbers -- stolen from Knuth Vol 2 pg 233
 185 sub bmul { #(num_str, num_str) return num_str
 186     local(*x, *y); ($x, $y) = (&bnorm($_[$[]), &bnorm($_[$[+1]));
 187     if ($x eq 'NaN') {
 188         'NaN';
 189     } elsif ($y eq 'NaN') {
 190         'NaN';
 191     } else {
 192         @x = &internal($x);
 193         @y = &internal($y);
 194         &external(&mul(*x,*y));
 195     }
 196 }
 197
 198 # multiply two numbers in internal representation
 199 # destroys the arguments, supposes that two arguments are different
 200 sub mul { #(*int_num_array, *int_num_array) return int_num_array
 201     local(*x, *y) = (shift, shift);
 202     local($signr) = (shift @x ne shift @y) ? '-' : '+';
 203     @prod = ();
 204     for $x (@x) {
 205       ($car, $cty) = (0, $[);
 206       for $y (@y) {
 207         $prod = $x * $y + $prod[$cty] + $car;
 208         $prod[$cty++] =
 209           $prod - ($car = int($prod * 1e-5)) * 1e5;
 210       }
 211       $prod[$cty] += $car if $car;
 212       $x = shift @prod;
 213     }
 214     ($signr, @x, @prod);
 215 }
 216
 217 # modulus
 218 sub bmod { #(num_str, num_str) return num_str
 219     (&bdiv(@_))[$[+1];
 220 }
 221
 222 sub bdiv { #(dividend: num_str, divisor: num_str) return num_str
 223     local (*x, *y); ($x, $y) = (&bnorm($_[$[]), &bnorm($_[$[+1]));
 224     return wantarray ? ('NaN','NaN') : 'NaN'
 225         if ($x eq 'NaN' || $y eq 'NaN' || $y eq '+0');
 226     return wantarray ? ('+0',$x) : '+0' if (&cmp(&abs($x),&abs($y)) < 0);
 227     @x = &internal($x); @y = &internal($y);
 228     $srem = $y[$[];
 229     $sr = (shift @x ne shift @y) ? '-' : '+';
 230     $car = $bar = $prd = 0;
 231     if (($dd = int(1e5/($y[$#y]+1))) != 1) {
 232         for $x (@x) {
 233             $x = $x * $dd + $car;
 234             $x -= ($car = int($x * 1e-5)) * 1e5;
 235         }
 236         push(@x, $car); $car = 0;
 237         for $y (@y) {
 238             $y = $y * $dd + $car;
 239             $y -= ($car = int($y * 1e-5)) * 1e5;
 240         }
 241     }
 242     else {
 243         push(@x, 0);
 244     }
 245     @q = (); ($v2,$v1) = @y[-2,-1];
 246     while ($#x > $#y) {
 247         ($u2,$u1,$u0) = @x[-3..-1];
 248         $q = (($u0 == $v1) ? 99999 : int(($u0*1e5+$u1)/$v1));
 249         --$q while ($v2*$q > ($u0*1e5+$u1-$q*$v1)*1e5+$u2);
 250         if ($q) {
 251             ($car, $bar) = (0,0);
 252             for ($y = $[, $x = $#x-$#y+$[-1; $y <= $#y; ++$y,++$x) {
 253                 $prd = $q * $y[$y] + $car;
 254                 $prd -= ($car = int($prd * 1e-5)) * 1e5;
 255                 $x[$x] += 1e5 if ($bar = (($x[$x] -= $prd + $bar) < 0));
 256             }
 257             if ($x[$#x] < $car + $bar) {
 258                 $car = 0; --$q;
 259                 for ($y = $[, $x = $#x-$#y+$[-1; $y <= $#y; ++$y,++$x) {
 260                     $x[$x] -= 1e5
 261                         if ($car = (($x[$x] += $y[$y] + $car) > 1e5));
 262                 }
 263             }
 264         }
 265         pop(@x); unshift(@q, $q);
 266     }
 267     if (wantarray) {
 268         @d = ();
 269         if ($dd != 1) {
 270             $car = 0;
 271             for $x (reverse @x) {
 272                 $prd = $car * 1e5 + $x;
 273                 $car = $prd - ($tmp = int($prd / $dd)) * $dd;
 274                 unshift(@d, $tmp);
 275             }
 276         }
 277         else {
 278             @d = @x;
 279         }
 280         (&external($sr, @q), &external($srem, @d, $zero));
 281     } else {
 282         &external($sr, @q);
 283     }
 284 }
 285
 286 # compute power of two numbers -- stolen from Knuth Vol 2 pg 233
 287 sub bpow { #(num_str, num_str) return num_str
 288     local(*x, *y); ($x, $y) = (&bnorm($_[$[]), &bnorm($_[$[+1]));
 289     if ($x eq 'NaN') {
 290         'NaN';
 291     } elsif ($y eq 'NaN') {
 292         'NaN';
 293     } elsif ($x eq '+1') {
 294         '+1';
 295     } elsif ($x eq '-1') {
 296         &bmod($x,2) ? '-1': '+1';
 297     } elsif ($y =~ /^-/) {
 298         'NaN';
 299     } elsif ($x eq '+0' && $y eq '+0') {
 300         'NaN';
 301     } else {
 302         @x = &internal($x);
 303         local(@pow2)=@x;
 304         local(@pow)=&internal("+1");
 305         local($y1,$res,@tmp1,@tmp2)=(1); # need tmp to send to mul
 306         while ($y ne '+0') {
 307           ($y,$res)=&bdiv($y,2);
 308           if ($res ne '+0') {@tmp=@pow2; @pow=&mul(*pow,*tmp);}
 309           if ($y ne '+0') {@tmp=@pow2;@pow2=&mul(*pow2,*tmp);}
 310         }
 311         &external(@pow);
 312     }
 313 }
 314
 315 1;
 316 __END__
 317
 318 =head1 NAME
 319
 320 Math::BigInt - Arbitrary size integer math package
 321
 322 =head1 SYNOPSIS
 323
 324   use Math::BigInt;
 325   $i = Math::BigInt->new($string);
 326
 327   $i->bneg return BINT               negation
 328   $i->babs return BINT               absolute value
 329   $i->bcmp(BINT) return CODE         compare numbers (undef,<0,=0,>0)
 330   $i->badd(BINT) return BINT         addition
 331   $i->bsub(BINT) return BINT         subtraction
 332   $i->bmul(BINT) return BINT         multiplication
 333   $i->bdiv(BINT) return (BINT,BINT)  division (quo,rem) just quo if scalar
 334   $i->bmod(BINT) return BINT         modulus
 335   $i->bgcd(BINT) return BINT         greatest common divisor
 336   $i->bnorm return BINT              normalization
 337
 338 =head1 DESCRIPTION
 339
 340 All basic math operations are overloaded if you declare your big
 341 integers as
 342
 343   $i = new Math::BigInt '123 456 789 123 456 789';
 344
 345
 346 =over 2
 347
 348 =item Canonical notation
 349
 350 Big integer value are strings of the form C</^[+-]\d+$/> with leading
 351 zeros suppressed.
 352
 353 =item Input
 354
 355 Input values to these routines may be strings of the form
 356 C</^\s*[+-]?[\d\s]+$/>.
 357
 358 =item Output
 359
 360 Output values always always in canonical form
 361
 362 =back
 363
 364 Actual math is done in an internal format consisting of an array
 365 whose first element is the sign (/^[+-]$/) and whose remaining
 366 elements are base 100000 digits with the least significant digit first.
 367 The string 'NaN' is used to represent the result when input arguments
 368 are not numbers, as well as the result of dividing by zero.
 369
 370 =head1 EXAMPLES
 371
 372    '+0'                            canonical zero value
 373    '   -123 123 123'               canonical value '-123123123'
 374    '1 23 456 7890'                 canonical value '+1234567890'
 375
 376
 377 =head1 BUGS
 378
 379 The current version of this module is a preliminary version of the
 380 real thing that is currently (as of perl5.002) under development.
 381
 382 =head1 AUTHOR
 383
 384 Mark Biggar, overloaded interface by Ilya Zakharevich.
 385
 386 =cut