| 1 | package bigint;
|
|---|
| 2 | #
|
|---|
| 3 | # This library is no longer being maintained, and is included for backward
|
|---|
| 4 | # compatibility with Perl 4 programs which may require it.
|
|---|
| 5 | #
|
|---|
| 6 | # In particular, this should not be used as an example of modern Perl
|
|---|
| 7 | # programming techniques.
|
|---|
| 8 | #
|
|---|
| 9 | # Suggested alternative: Math::BigInt
|
|---|
| 10 | #
|
|---|
| 11 | # arbitrary size integer math package
|
|---|
| 12 | #
|
|---|
| 13 | # by Mark Biggar
|
|---|
| 14 | #
|
|---|
| 15 | # Canonical Big integer value are strings of the form
|
|---|
| 16 | # /^[+-]\d+$/ with leading zeros suppressed
|
|---|
| 17 | # Input values to these routines may be strings of the form
|
|---|
| 18 | # /^\s*[+-]?[\d\s]+$/.
|
|---|
| 19 | # Examples:
|
|---|
| 20 | # '+0' canonical zero value
|
|---|
| 21 | # ' -123 123 123' canonical value '-123123123'
|
|---|
| 22 | # '1 23 456 7890' canonical value '+1234567890'
|
|---|
| 23 | # Output values always in canonical form
|
|---|
| 24 | #
|
|---|
| 25 | # Actual math is done in an internal format consisting of an array
|
|---|
| 26 | # whose first element is the sign (/^[+-]$/) and whose remaining
|
|---|
| 27 | # elements are base 100000 digits with the least significant digit first.
|
|---|
| 28 | # The string 'NaN' is used to represent the result when input arguments
|
|---|
| 29 | # are not numbers, as well as the result of dividing by zero
|
|---|
| 30 | #
|
|---|
| 31 | # routines provided are:
|
|---|
| 32 | #
|
|---|
| 33 | # bneg(BINT) return BINT negation
|
|---|
| 34 | # babs(BINT) return BINT absolute value
|
|---|
| 35 | # bcmp(BINT,BINT) return CODE compare numbers (undef,<0,=0,>0)
|
|---|
| 36 | # badd(BINT,BINT) return BINT addition
|
|---|
| 37 | # bsub(BINT,BINT) return BINT subtraction
|
|---|
| 38 | # bmul(BINT,BINT) return BINT multiplication
|
|---|
| 39 | # bdiv(BINT,BINT) return (BINT,BINT) division (quo,rem) just quo if scalar
|
|---|
| 40 | # bmod(BINT,BINT) return BINT modulus
|
|---|
| 41 | # bgcd(BINT,BINT) return BINT greatest common divisor
|
|---|
| 42 | # bnorm(BINT) return BINT normalization
|
|---|
| 43 | #
|
|---|
| 44 |
|
|---|
| 45 | # overcome a floating point problem on certain osnames (posix-bc, os390)
|
|---|
| 46 | BEGIN {
|
|---|
| 47 | my $x = 100000.0;
|
|---|
| 48 | my $use_mult = int($x*1e-5)*1e5 == $x ? 1 : 0;
|
|---|
| 49 | }
|
|---|
| 50 |
|
|---|
| 51 | $zero = 0;
|
|---|
| 52 |
|
|---|
| 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 main'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 main'bneg { #(num_str) return num_str
|
|---|
| 90 | local($_) = &'bnorm(@_);
|
|---|
| 91 | vec($_,0,8) ^= ord('+') ^ ord('-') unless $_ eq '+0';
|
|---|
| 92 | s/^./N/ unless /^[-+]/; # works both in ASCII and EBCDIC
|
|---|
| 93 | $_;
|
|---|
| 94 | }
|
|---|
| 95 |
|
|---|
| 96 | # Returns the absolute value of the input.
|
|---|
| 97 | sub main'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 |
|
|---|
| 108 | # Compares 2 values. Returns one of undef, <0, =0, >0. (suitable for sort)
|
|---|
| 109 | sub main'bcmp { #(num_str, num_str) return cond_code
|
|---|
| 110 | local($x,$y) = (&'bnorm($_[$[]),&'bnorm($_[$[+1]));
|
|---|
| 111 | if ($x eq 'NaN') {
|
|---|
| 112 | undef;
|
|---|
| 113 | } elsif ($y eq 'NaN') {
|
|---|
| 114 | undef;
|
|---|
| 115 | } else {
|
|---|
| 116 | &cmp($x,$y);
|
|---|
| 117 | }
|
|---|
| 118 | }
|
|---|
| 119 |
|
|---|
| 120 | sub cmp { # post-normalized compare for internal use
|
|---|
| 121 | local($cx, $cy) = @_;
|
|---|
| 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 |
|
|---|
| 141 | sub main'badd { #(num_str, num_str) return num_str
|
|---|
| 142 | local(*x, *y); ($x, $y) = (&'bnorm($_[$[]),&'bnorm($_[$[+1]));
|
|---|
| 143 | if ($x eq 'NaN') {
|
|---|
| 144 | 'NaN';
|
|---|
| 145 | } elsif ($y eq 'NaN') {
|
|---|
| 146 | 'NaN';
|
|---|
| 147 | } else {
|
|---|
| 148 | @x = &internal($x); # convert to internal form
|
|---|
| 149 | @y = &internal($y);
|
|---|
| 150 | local($sx, $sy) = (shift @x, shift @y); # get signs
|
|---|
| 151 | if ($sx eq $sy) {
|
|---|
| 152 | &external($sx, &add(*x, *y)); # if same sign add
|
|---|
| 153 | } else {
|
|---|
| 154 | ($x, $y) = (&abs($x),&abs($y)); # make abs
|
|---|
| 155 | if (&cmp($y,$x) > 0) {
|
|---|
| 156 | &external($sy, &sub(*y, *x));
|
|---|
| 157 | } else {
|
|---|
| 158 | &external($sx, &sub(*x, *y));
|
|---|
| 159 | }
|
|---|
| 160 | }
|
|---|
| 161 | }
|
|---|
| 162 | }
|
|---|
| 163 |
|
|---|
| 164 | sub main'bsub { #(num_str, num_str) return num_str
|
|---|
| 165 | &'badd($_[$[],&'bneg($_[$[+1]));
|
|---|
| 166 | }
|
|---|
| 167 |
|
|---|
| 168 | # GCD -- Euclids algorithm Knuth Vol 2 pg 296
|
|---|
| 169 | sub main'bgcd { #(num_str, num_str) return num_str
|
|---|
| 170 | local($x,$y) = (&'bnorm($_[$[]),&'bnorm($_[$[+1]));
|
|---|
| 171 | if ($x eq 'NaN' || $y eq 'NaN') {
|
|---|
| 172 | 'NaN';
|
|---|
| 173 | } else {
|
|---|
| 174 | ($x, $y) = ($y,&'bmod($x,$y)) while $y ne '+0';
|
|---|
| 175 | $x;
|
|---|
| 176 | }
|
|---|
| 177 | }
|
|---|
| 178 | |
|---|
| 179 |
|
|---|
| 180 | # routine to add two base 1e5 numbers
|
|---|
| 181 | # stolen from Knuth Vol 2 Algorithm A pg 231
|
|---|
| 182 | # there are separate routines to add and sub as per Kunth pg 233
|
|---|
| 183 | sub add { #(int_num_array, int_num_array) return int_num_array
|
|---|
| 184 | local(*x, *y) = @_;
|
|---|
| 185 | $car = 0;
|
|---|
| 186 | for $x (@x) {
|
|---|
| 187 | last unless @y || $car;
|
|---|
| 188 | $x -= 1e5 if $car = (($x += shift(@y) + $car) >= 1e5) ? 1 : 0;
|
|---|
| 189 | }
|
|---|
| 190 | for $y (@y) {
|
|---|
| 191 | last unless $car;
|
|---|
| 192 | $y -= 1e5 if $car = (($y += $car) >= 1e5) ? 1 : 0;
|
|---|
| 193 | }
|
|---|
| 194 | (@x, @y, $car);
|
|---|
| 195 | }
|
|---|
| 196 |
|
|---|
| 197 | # subtract base 1e5 numbers -- stolen from Knuth Vol 2 pg 232, $x > $y
|
|---|
| 198 | sub sub { #(int_num_array, int_num_array) return int_num_array
|
|---|
| 199 | local(*sx, *sy) = @_;
|
|---|
| 200 | $bar = 0;
|
|---|
| 201 | for $sx (@sx) {
|
|---|
| 202 | last unless @y || $bar;
|
|---|
| 203 | $sx += 1e5 if $bar = (($sx -= shift(@sy) + $bar) < 0);
|
|---|
| 204 | }
|
|---|
| 205 | @sx;
|
|---|
| 206 | }
|
|---|
| 207 |
|
|---|
| 208 | # multiply two numbers -- stolen from Knuth Vol 2 pg 233
|
|---|
| 209 | sub main'bmul { #(num_str, num_str) return num_str
|
|---|
| 210 | local(*x, *y); ($x, $y) = (&'bnorm($_[$[]), &'bnorm($_[$[+1]));
|
|---|
| 211 | if ($x eq 'NaN') {
|
|---|
| 212 | 'NaN';
|
|---|
| 213 | } elsif ($y eq 'NaN') {
|
|---|
| 214 | 'NaN';
|
|---|
| 215 | } else {
|
|---|
| 216 | @x = &internal($x);
|
|---|
| 217 | @y = &internal($y);
|
|---|
| 218 | local($signr) = (shift @x ne shift @y) ? '-' : '+';
|
|---|
| 219 | @prod = ();
|
|---|
| 220 | for $x (@x) {
|
|---|
| 221 | ($car, $cty) = (0, $[);
|
|---|
| 222 | for $y (@y) {
|
|---|
| 223 | $prod = $x * $y + $prod[$cty] + $car;
|
|---|
| 224 | if ($use_mult) {
|
|---|
| 225 | $prod[$cty++] =
|
|---|
| 226 | $prod - ($car = int($prod * 1e-5)) * 1e5;
|
|---|
| 227 | }
|
|---|
| 228 | else {
|
|---|
| 229 | $prod[$cty++] =
|
|---|
| 230 | $prod - ($car = int($prod / 1e5)) * 1e5;
|
|---|
| 231 | }
|
|---|
| 232 | }
|
|---|
| 233 | $prod[$cty] += $car if $car;
|
|---|
| 234 | $x = shift @prod;
|
|---|
| 235 | }
|
|---|
| 236 | &external($signr, @x, @prod);
|
|---|
| 237 | }
|
|---|
| 238 | }
|
|---|
| 239 |
|
|---|
| 240 | # modulus
|
|---|
| 241 | sub main'bmod { #(num_str, num_str) return num_str
|
|---|
| 242 | (&'bdiv(@_))[$[+1];
|
|---|
| 243 | }
|
|---|
| 244 | |
|---|
| 245 |
|
|---|
| 246 | sub main'bdiv { #(dividend: num_str, divisor: num_str) return num_str
|
|---|
| 247 | local (*x, *y); ($x, $y) = (&'bnorm($_[$[]), &'bnorm($_[$[+1]));
|
|---|
| 248 | return wantarray ? ('NaN','NaN') : 'NaN'
|
|---|
| 249 | if ($x eq 'NaN' || $y eq 'NaN' || $y eq '+0');
|
|---|
| 250 | return wantarray ? ('+0',$x) : '+0' if (&cmp(&abs($x),&abs($y)) < 0);
|
|---|
| 251 | @x = &internal($x); @y = &internal($y);
|
|---|
| 252 | $srem = $y[$[];
|
|---|
| 253 | $sr = (shift @x ne shift @y) ? '-' : '+';
|
|---|
| 254 | $car = $bar = $prd = 0;
|
|---|
| 255 | if (($dd = int(1e5/($y[$#y]+1))) != 1) {
|
|---|
| 256 | for $x (@x) {
|
|---|
| 257 | $x = $x * $dd + $car;
|
|---|
| 258 | if ($use_mult) {
|
|---|
| 259 | $x -= ($car = int($x * 1e-5)) * 1e5;
|
|---|
| 260 | }
|
|---|
| 261 | else {
|
|---|
| 262 | $x -= ($car = int($x / 1e5)) * 1e5;
|
|---|
| 263 | }
|
|---|
| 264 | }
|
|---|
| 265 | push(@x, $car); $car = 0;
|
|---|
| 266 | for $y (@y) {
|
|---|
| 267 | $y = $y * $dd + $car;
|
|---|
| 268 | if ($use_mult) {
|
|---|
| 269 | $y -= ($car = int($y * 1e-5)) * 1e5;
|
|---|
| 270 | }
|
|---|
| 271 | else {
|
|---|
| 272 | $y -= ($car = int($y / 1e5)) * 1e5;
|
|---|
| 273 | }
|
|---|
| 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 | if ($use_mult) {
|
|---|
| 289 | $prd -= ($car = int($prd * 1e-5)) * 1e5;
|
|---|
| 290 | }
|
|---|
| 291 | else {
|
|---|
| 292 | $prd -= ($car = int($prd / 1e5)) * 1e5;
|
|---|
| 293 | }
|
|---|
| 294 | $x[$x] += 1e5 if ($bar = (($x[$x] -= $prd + $bar) < 0));
|
|---|
| 295 | }
|
|---|
| 296 | if ($x[$#x] < $car + $bar) {
|
|---|
| 297 | $car = 0; --$q;
|
|---|
| 298 | for ($y = $[, $x = $#x-$#y+$[-1; $y <= $#y; ++$y,++$x) {
|
|---|
| 299 | $x[$x] -= 1e5
|
|---|
| 300 | if ($car = (($x[$x] += $y[$y] + $car) > 1e5));
|
|---|
| 301 | }
|
|---|
| 302 | }
|
|---|
| 303 | }
|
|---|
| 304 | pop(@x); unshift(@q, $q);
|
|---|
| 305 | }
|
|---|
| 306 | if (wantarray) {
|
|---|
| 307 | @d = ();
|
|---|
| 308 | if ($dd != 1) {
|
|---|
| 309 | $car = 0;
|
|---|
| 310 | for $x (reverse @x) {
|
|---|
| 311 | $prd = $car * 1e5 + $x;
|
|---|
| 312 | $car = $prd - ($tmp = int($prd / $dd)) * $dd;
|
|---|
| 313 | unshift(@d, $tmp);
|
|---|
| 314 | }
|
|---|
| 315 | }
|
|---|
| 316 | else {
|
|---|
| 317 | @d = @x;
|
|---|
| 318 | }
|
|---|
| 319 | (&external($sr, @q), &external($srem, @d, $zero));
|
|---|
| 320 | } else {
|
|---|
| 321 | &external($sr, @q);
|
|---|
| 322 | }
|
|---|
| 323 | }
|
|---|
| 324 | 1;
|
|---|
