I don't think trying to bracket the hires time with lores
[p5sagit/p5-mst-13.2.git] / lib / bigint.pl
CommitLineData
5303340c 1package bigint;
a6d71656 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#
5303340c 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'
d1be9408 23# Output values always in canonical form
5303340c 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#
8990e307 44
1f45ae4a 45# overcome a floating point problem on certain osnames (posix-bc, os390)
46BEGIN {
47 my $x = 100000.0;
48 my $use_mult = int($x*1e-5)*1e5 == $x ? 1 : 0;
49}
50
8990e307 51$zero = 0;
52
5303340c 53\f
54# normalize string form of number. Strip leading zeros. Strip any
55# white space and add a sign, if missing.
56# Strings that are not numbers result the value 'NaN'.
8990e307 57
5303340c 58sub main'bnorm { #(num_str) return num_str
59 local($_) = @_;
60 s/\s+//g; # strip white space
61 if (s/^([+-]?)0*(\d+)$/$1$2/) { # test if number
79072805 62 substr($_,$[,0) = '+' unless $1; # Add missing sign
5303340c 63 s/^-0/+0/;
64 $_;
65 } else {
66 'NaN';
67 }
68}
69
70# Convert a number from string format to internal base 100000 format.
71# Assumes normalized value as input.
72sub internal { #(num_str) return int_num_array
73 local($d) = @_;
79072805 74 ($is,$il) = (substr($d,$[,1),length($d)-2);
75 substr($d,$[,1) = '';
5303340c 76 ($is, reverse(unpack("a" . ($il%5+1) . ("a5" x ($il/5)), $d)));
77}
78
79# Convert a number from internal base 100000 format to string format.
80# This routine scribbles all over input array.
81sub external { #(int_num_array) return num_str
82 $es = shift;
83 grep($_ > 9999 || ($_ = substr('0000'.$_,-5)), @_); # zero pad
84 &'bnorm(join('', $es, reverse(@_))); # reverse concat and normalize
85}
86
87# Negate input value.
88sub main'bneg { #(num_str) return num_str
89 local($_) = &'bnorm(@_);
90 vec($_,0,8) ^= ord('+') ^ ord('-') unless $_ eq '+0';
9d116dd7 91 s/^./N/ unless /^[-+]/; # works both in ASCII and EBCDIC
5303340c 92 $_;
93}
94
95# Returns the absolute value of the input.
96sub main'babs { #(num_str) return num_str
97 &abs(&'bnorm(@_));
98}
99
100sub abs { # post-normalized abs for internal use
101 local($_) = @_;
102 s/^-/+/;
103 $_;
104}
105\f
106# Compares 2 values. Returns one of undef, <0, =0, >0. (suitable for sort)
107sub main'bcmp { #(num_str, num_str) return cond_code
79072805 108 local($x,$y) = (&'bnorm($_[$[]),&'bnorm($_[$[+1]));
5303340c 109 if ($x eq 'NaN') {
110 undef;
111 } elsif ($y eq 'NaN') {
112 undef;
113 } else {
114 &cmp($x,$y);
115 }
116}
117
118sub cmp { # post-normalized compare for internal use
119 local($cx, $cy) = @_;
1e2e1ae8 120 return 0 if ($cx eq $cy);
121
122 local($sx, $sy) = (substr($cx, 0, 1), substr($cy, 0, 1));
123 local($ld);
124
125 if ($sx eq '+') {
126 return 1 if ($sy eq '-' || $cy eq '+0');
127 $ld = length($cx) - length($cy);
128 return $ld if ($ld);
129 return $cx cmp $cy;
130 } else { # $sx eq '-'
131 return -1 if ($sy eq '+');
132 $ld = length($cy) - length($cx);
133 return $ld if ($ld);
134 return $cy cmp $cx;
135 }
136
5303340c 137}
138
139sub main'badd { #(num_str, num_str) return num_str
79072805 140 local(*x, *y); ($x, $y) = (&'bnorm($_[$[]),&'bnorm($_[$[+1]));
5303340c 141 if ($x eq 'NaN') {
142 'NaN';
143 } elsif ($y eq 'NaN') {
144 'NaN';
145 } else {
146 @x = &internal($x); # convert to internal form
147 @y = &internal($y);
148 local($sx, $sy) = (shift @x, shift @y); # get signs
149 if ($sx eq $sy) {
150 &external($sx, &add(*x, *y)); # if same sign add
151 } else {
152 ($x, $y) = (&abs($x),&abs($y)); # make abs
153 if (&cmp($y,$x) > 0) {
154 &external($sy, &sub(*y, *x));
155 } else {
156 &external($sx, &sub(*x, *y));
157 }
158 }
159 }
160}
161
162sub main'bsub { #(num_str, num_str) return num_str
79072805 163 &'badd($_[$[],&'bneg($_[$[+1]));
5303340c 164}
165
166# GCD -- Euclids algorithm Knuth Vol 2 pg 296
167sub main'bgcd { #(num_str, num_str) return num_str
79072805 168 local($x,$y) = (&'bnorm($_[$[]),&'bnorm($_[$[+1]));
68decaef 169 if ($x eq 'NaN' || $y eq 'NaN') {
5303340c 170 'NaN';
68decaef 171 } else {
5303340c 172 ($x, $y) = ($y,&'bmod($x,$y)) while $y ne '+0';
173 $x;
174 }
175}
176\f
68decaef 177# routine to add two base 1e5 numbers
5303340c 178# stolen from Knuth Vol 2 Algorithm A pg 231
179# there are separate routines to add and sub as per Kunth pg 233
180sub add { #(int_num_array, int_num_array) return int_num_array
181 local(*x, *y) = @_;
182 $car = 0;
183 for $x (@x) {
184 last unless @y || $car;
55497cff 185 $x -= 1e5 if $car = (($x += shift(@y) + $car) >= 1e5) ? 1 : 0;
5303340c 186 }
187 for $y (@y) {
188 last unless $car;
55497cff 189 $y -= 1e5 if $car = (($y += $car) >= 1e5) ? 1 : 0;
5303340c 190 }
191 (@x, @y, $car);
192}
193
68decaef 194# subtract base 1e5 numbers -- stolen from Knuth Vol 2 pg 232, $x > $y
5303340c 195sub sub { #(int_num_array, int_num_array) return int_num_array
196 local(*sx, *sy) = @_;
197 $bar = 0;
198 for $sx (@sx) {
199 last unless @y || $bar;
e334a159 200 $sx += 1e5 if $bar = (($sx -= shift(@sy) + $bar) < 0);
5303340c 201 }
202 @sx;
203}
204
205# multiply two numbers -- stolen from Knuth Vol 2 pg 233
206sub main'bmul { #(num_str, num_str) return num_str
79072805 207 local(*x, *y); ($x, $y) = (&'bnorm($_[$[]), &'bnorm($_[$[+1]));
5303340c 208 if ($x eq 'NaN') {
209 'NaN';
210 } elsif ($y eq 'NaN') {
211 'NaN';
212 } else {
213 @x = &internal($x);
214 @y = &internal($y);
215 local($signr) = (shift @x ne shift @y) ? '-' : '+';
216 @prod = ();
217 for $x (@x) {
79072805 218 ($car, $cty) = (0, $[);
5303340c 219 for $y (@y) {
220 $prod = $x * $y + $prod[$cty] + $car;
1f45ae4a 221 if ($use_mult) {
222 $prod[$cty++] =
223 $prod - ($car = int($prod * 1e-5)) * 1e5;
224 }
225 else {
226 $prod[$cty++] =
227 $prod - ($car = int($prod / 1e5)) * 1e5;
228 }
5303340c 229 }
230 $prod[$cty] += $car if $car;
231 $x = shift @prod;
232 }
233 &external($signr, @x, @prod);
234 }
235}
236
237# modulus
238sub main'bmod { #(num_str, num_str) return num_str
79072805 239 (&'bdiv(@_))[$[+1];
5303340c 240}
241\f
242sub main'bdiv { #(dividend: num_str, divisor: num_str) return num_str
79072805 243 local (*x, *y); ($x, $y) = (&'bnorm($_[$[]), &'bnorm($_[$[+1]));
5303340c 244 return wantarray ? ('NaN','NaN') : 'NaN'
245 if ($x eq 'NaN' || $y eq 'NaN' || $y eq '+0');
246 return wantarray ? ('+0',$x) : '+0' if (&cmp(&abs($x),&abs($y)) < 0);
247 @x = &internal($x); @y = &internal($y);
79072805 248 $srem = $y[$[];
5303340c 249 $sr = (shift @x ne shift @y) ? '-' : '+';
250 $car = $bar = $prd = 0;
68decaef 251 if (($dd = int(1e5/($y[$#y]+1))) != 1) {
5303340c 252 for $x (@x) {
253 $x = $x * $dd + $car;
1f45ae4a 254 if ($use_mult) {
68decaef 255 $x -= ($car = int($x * 1e-5)) * 1e5;
1f45ae4a 256 }
257 else {
258 $x -= ($car = int($x / 1e5)) * 1e5;
259 }
5303340c 260 }
261 push(@x, $car); $car = 0;
262 for $y (@y) {
263 $y = $y * $dd + $car;
1f45ae4a 264 if ($use_mult) {
68decaef 265 $y -= ($car = int($y * 1e-5)) * 1e5;
1f45ae4a 266 }
267 else {
268 $y -= ($car = int($y / 1e5)) * 1e5;
269 }
5303340c 270 }
271 }
272 else {
273 push(@x, 0);
274 }
ed6116ce 275 @q = (); ($v2,$v1) = @y[-2,-1];
5303340c 276 while ($#x > $#y) {
ed6116ce 277 ($u2,$u1,$u0) = @x[-3..-1];
68decaef 278 $q = (($u0 == $v1) ? 99999 : int(($u0*1e5+$u1)/$v1));
279 --$q while ($v2*$q > ($u0*1e5+$u1-$q*$v1)*1e5+$u2);
5303340c 280 if ($q) {
281 ($car, $bar) = (0,0);
79072805 282 for ($y = $[, $x = $#x-$#y+$[-1; $y <= $#y; ++$y,++$x) {
5303340c 283 $prd = $q * $y[$y] + $car;
1f45ae4a 284 if ($use_mult) {
68decaef 285 $prd -= ($car = int($prd * 1e-5)) * 1e5;
1f45ae4a 286 }
287 else {
288 $prd -= ($car = int($prd / 1e5)) * 1e5;
289 }
68decaef 290 $x[$x] += 1e5 if ($bar = (($x[$x] -= $prd + $bar) < 0));
5303340c 291 }
292 if ($x[$#x] < $car + $bar) {
293 $car = 0; --$q;
79072805 294 for ($y = $[, $x = $#x-$#y+$[-1; $y <= $#y; ++$y,++$x) {
68decaef 295 $x[$x] -= 1e5
296 if ($car = (($x[$x] += $y[$y] + $car) > 1e5));
5303340c 297 }
298 }
299 }
300 pop(@x); unshift(@q, $q);
301 }
302 if (wantarray) {
303 @d = ();
304 if ($dd != 1) {
305 $car = 0;
306 for $x (reverse @x) {
68decaef 307 $prd = $car * 1e5 + $x;
5303340c 308 $car = $prd - ($tmp = int($prd / $dd)) * $dd;
309 unshift(@d, $tmp);
310 }
311 }
312 else {
313 @d = @x;
314 }
8990e307 315 (&external($sr, @q), &external($srem, @d, $zero));
5303340c 316 } else {
317 &external($sr, @q);
318 }
319}
3201;