This is a live mirror of the Perl 5 development currently hosted at https://github.com/perl/perl5
assorted VMS test fix-ups, $Config{prefixexp} revisited
[perl5.git] / lib / Math / BigInt / Calc.pm
CommitLineData
0716bf9b
JH
1package Math::BigInt::Calc;
2
3use 5.005;
4use strict;
574bacfe 5# use warnings; # dont use warnings for older Perls
0716bf9b
JH
6
7require Exporter;
bd05a461 8use vars qw/@ISA $VERSION/;
0716bf9b
JH
9@ISA = qw(Exporter);
10
13a12e00 11$VERSION = '0.23';
0716bf9b
JH
12
13# Package to store unsigned big integers in decimal and do math with them
14
15# Internally the numbers are stored in an array with at least 1 element, no
027dc388
JH
16# leading zero parts (except the first) and in base 1eX where X is determined
17# automatically at loading time to be the maximum possible value
0716bf9b
JH
18
19# todo:
20# - fully remove funky $# stuff (maybe)
0716bf9b
JH
21
22# USE_MUL: due to problems on certain os (os390, posix-bc) "* 1e-5" is used
ee15d750
JH
23# instead of "/ 1e5" at some places, (marked with USE_MUL). Other platforms
24# BS2000, some Crays need USE_DIV instead.
bd05a461
JH
25# The BEGIN block is used to determine which of the two variants gives the
26# correct result.
0716bf9b
JH
27
28##############################################################################
29# global constants, flags and accessory
30
31# constants for easier life
32my $nan = 'NaN';
61f5c3f5 33my ($MBASE,$BASE,$RBASE,$BASE_LEN,$MAX_VAL,$BASE_LEN2,$BASE_LEN_SMALL);
394e6ffb
JH
34my ($AND_BITS,$XOR_BITS,$OR_BITS);
35my ($AND_MASK,$XOR_MASK,$OR_MASK);
61f5c3f5 36my ($LEN_CONVERT);
ee15d750
JH
37
38sub _base_len
39 {
dccbb853
JH
40 # set/get the BASE_LEN and assorted other, connected values
41 # used only be the testsuite, set is used only by the BEGIN block below
394e6ffb
JH
42 shift;
43
ee15d750
JH
44 my $b = shift;
45 if (defined $b)
46 {
61f5c3f5
T
47 # find whether we can use mul or div or none in mul()/div()
48 # (in last case reduce BASE_LEN_SMALL)
49 $BASE_LEN_SMALL = $b+1;
50 my $caught = 0;
51 while (--$BASE_LEN_SMALL > 5)
394e6ffb 52 {
61f5c3f5
T
53 $MBASE = int("1e".$BASE_LEN_SMALL);
54 $RBASE = abs('1e-'.$BASE_LEN_SMALL); # see USE_MUL
394e6ffb 55 $caught = 0;
61f5c3f5
T
56 $caught += 1 if (int($MBASE * $RBASE) != 1); # should be 1
57 $caught += 2 if (int($MBASE / $MBASE) != 1); # should be 1
394e6ffb
JH
58 last if $caught != 3;
59 }
61f5c3f5
T
60 # BASE_LEN is used for anything else than mul()/div()
61 $BASE_LEN = $BASE_LEN_SMALL;
62 $BASE_LEN = shift if (defined $_[0]); # one more arg?
ee15d750 63 $BASE = int("1e".$BASE_LEN);
61f5c3f5
T
64
65 $BASE_LEN2 = int($BASE_LEN_SMALL / 2); # for mul shortcut
66 $MBASE = int("1e".$BASE_LEN_SMALL);
67 $RBASE = abs('1e-'.$BASE_LEN_SMALL); # see USE_MUL
68 $MAX_VAL = $MBASE-1;
69 $LEN_CONVERT = 0;
70 $LEN_CONVERT = 1 if $BASE_LEN_SMALL != $BASE_LEN;
71
72 #print "BASE_LEN: $BASE_LEN MAX_VAL: $MAX_VAL BASE: $BASE RBASE: $RBASE ";
73 #print "BASE_LEN_SMALL: $BASE_LEN_SMALL MBASE: $MBASE\n";
74
394e6ffb 75 if ($caught & 1 != 0)
ee15d750
JH
76 {
77 # must USE_MUL
ee15d750
JH
78 *{_mul} = \&_mul_use_mul;
79 *{_div} = \&_div_use_mul;
80 }
394e6ffb 81 else # $caught must be 2, since it can't be 1 nor 3
ee15d750 82 {
ee15d750
JH
83 # can USE_DIV instead
84 *{_mul} = \&_mul_use_div;
85 *{_div} = \&_div_use_div;
86 }
87 }
61f5c3f5
T
88 return $BASE_LEN unless wantarray;
89 return ($BASE_LEN, $AND_BITS, $XOR_BITS, $OR_BITS, $BASE_LEN_SMALL, $MAX_VAL);
ee15d750 90 }
574bacfe
JH
91
92BEGIN
93 {
bd05a461 94 # from Daniel Pfeiffer: determine largest group of digits that is precisely
574bacfe 95 # multipliable with itself plus carry
dccbb853
JH
96 # Test now changed to expect the proper pattern, not a result off by 1 or 2
97 my ($e, $num) = 3; # lowest value we will use is 3+1-1 = 3
bd05a461
JH
98 do
99 {
100 $num = ('9' x ++$e) + 0;
394e6ffb 101 $num *= $num + 1.0;
394e6ffb
JH
102 } while ("$num" =~ /9{$e}0{$e}/); # must be a certain pattern
103 $e--; # last test failed, so retract one step
104 # the limits below brush the problems with the test above under the rug:
105 # the test should be able to find the proper $e automatically
106 $e = 5 if $^O =~ /^uts/; # UTS get's some special treatment
107 $e = 5 if $^O =~ /^unicos/; # unicos is also problematic (6 seems to work
108 # there, but we play safe)
109 $e = 8 if $e > 8; # cap, for VMS, OS/390 and other 64 bit systems
110
61f5c3f5
T
111 # determine how many digits fit into an integer and can be safely added
112 # together plus carry w/o causing an overflow
113
114 # this below detects 15 on a 64 bit system, because after that it becomes
115 # 1e16 and not 1000000 :/ I can make it detect 18, but then I get a lot of
116 # test failures. Ugh! (Tomake detect 18: uncomment lines marked with *)
117 use integer;
118 my $bi = 5; # approx. 16 bit
119 $num = int('9' x $bi);
120 # $num = 99999; # *
121 # while ( ($num+$num+1) eq '1' . '9' x $bi) # *
122 while ( int($num+$num+1) eq '1' . '9' x $bi)
123 {
124 $bi++; $num = int('9' x $bi);
125 # $bi++; $num *= 10; $num += 9; # *
126 }
127 $bi--; # back off one step
128 # by setting them equal, we ignore the findings and use the default
129 # one-size-fits-all approach from former versions
130 $bi = $e; # XXX, this should work always
131
132 __PACKAGE__->_base_len($e,$bi); # set and store
394e6ffb
JH
133
134 # find out how many bits _and, _or and _xor can take (old default = 16)
135 # I don't think anybody has yet 128 bit scalars, so let's play safe.
394e6ffb
JH
136 local $^W = 0; # don't warn about 'nonportable number'
137 $AND_BITS = 15; $XOR_BITS = 15; $OR_BITS = 15;
138
139 # find max bits, we will not go higher than numberofbits that fit into $BASE
140 # to make _and etc simpler (and faster for smaller, slower for large numbers)
141 my $max = 16;
142 while (2 ** $max < $BASE) { $max++; }
143 my ($x,$y,$z);
144 do {
145 $AND_BITS++;
146 $x = oct('0b' . '1' x $AND_BITS); $y = $x & $x;
147 $z = (2 ** $AND_BITS) - 1;
148 } while ($AND_BITS < $max && $x == $z && $y == $x);
149 $AND_BITS --; # retreat one step
150 do {
151 $XOR_BITS++;
152 $x = oct('0b' . '1' x $XOR_BITS); $y = $x ^ 0;
153 $z = (2 ** $XOR_BITS) - 1;
154 } while ($XOR_BITS < $max && $x == $z && $y == $x);
155 $XOR_BITS --; # retreat one step
156 do {
157 $OR_BITS++;
158 $x = oct('0b' . '1' x $OR_BITS); $y = $x | $x;
159 $z = (2 ** $OR_BITS) - 1;
160 } while ($OR_BITS < $max && $x == $z && $y == $x);
161 $OR_BITS --; # retreat one step
162
574bacfe
JH
163 }
164
0716bf9b 165##############################################################################
61f5c3f5
T
166# convert between the "small" and the "large" representation
167
168sub _to_large
169 {
170 # take an array in base $BASE_LEN_SMALL and convert it in-place to $BASE_LEN
171 my ($c,$x) = @_;
172
173# print "_to_large $BASE_LEN_SMALL => $BASE_LEN\n";
174
175 return $x if $LEN_CONVERT == 0 || # nothing to converconvertor
176 @$x == 1; # only one element => early out
177
178 # 12345 67890 12345 67890 contents
179 # to 3 2 1 0 index
180 # 123456 7890123 4567890 contents
181
182# # faster variant
183# my @d; my $str = '';
184# my $z = '0' x $BASE_LEN_SMALL;
185# foreach (@$x)
186# {
187# # ... . 04321 . 000321
188# $str = substr($z.$_,-$BASE_LEN_SMALL,$BASE_LEN_SMALL) . $str;
189# if (length($str) > $BASE_LEN)
190# {
191# push @d, substr($str,-$BASE_LEN,$BASE_LEN); # extract one piece
192# substr($str,-$BASE_LEN,$BASE_LEN) = ''; # remove it
193# }
194# }
195# push @d, $str if $str !~ /^0*$/; # extract last piece
196# @$x = @d;
197# $x->[-1] = int($x->[-1]); # strip leading zero
198# $x;
199
200 my $ret = "";
201 my $l = scalar @$x; # number of parts
202 $l --; $ret .= int($x->[$l]); $l--;
203 my $z = '0' x ($BASE_LEN_SMALL-1);
204 while ($l >= 0)
205 {
206 $ret .= substr($z.$x->[$l],-$BASE_LEN_SMALL);
207 $l--;
208 }
209 my $str = _new($c,\$ret); # make array
210 @$x = @$str; # clobber contents of $x
211 $x->[-1] = int($x->[-1]); # strip leading zero
212 }
213
214sub _to_small
215 {
216 # take an array in base $BASE_LEN and convert it in-place to $BASE_LEN_SMALL
217 my ($c,$x) = @_;
218
219 return $x if $LEN_CONVERT == 0; # nothing to do
220 return $x if @$x == 1 && length(int($x->[0])) <= $BASE_LEN_SMALL;
221
222 my $d = _str($c,$x);
223 my $il = length($$d)-1;
224 ## this leaves '00000' instead of int 0 and will be corrected after any op
225 # clobber contents of $x
226 @$x = reverse(unpack("a" . ($il % $BASE_LEN_SMALL+1)
227 . ("a$BASE_LEN_SMALL" x ($il / $BASE_LEN_SMALL)), $$d));
228
229 $x->[-1] = int($x->[-1]); # strip leading zero
230 }
231
232###############################################################################
0716bf9b
JH
233
234sub _new
235 {
394e6ffb 236 # (ref to string) return ref to num_array
0716bf9b
JH
237 # Convert a number from string format to internal base 100000 format.
238 # Assumes normalized value as input.
574bacfe 239 my $d = $_[1];
61f5c3f5
T
240 my $il = length($$d)-1;
241 # this leaves '00000' instead of int 0 and will be corrected after any op
242 [ reverse(unpack("a" . ($il % $BASE_LEN+1)
574bacfe 243 . ("a$BASE_LEN" x ($il / $BASE_LEN)), $$d)) ];
0716bf9b 244 }
394e6ffb
JH
245
246BEGIN
247 {
248 $AND_MASK = __PACKAGE__->_new( \( 2 ** $AND_BITS ));
249 $XOR_MASK = __PACKAGE__->_new( \( 2 ** $XOR_BITS ));
250 $OR_MASK = __PACKAGE__->_new( \( 2 ** $OR_BITS ));
251 }
0716bf9b
JH
252
253sub _zero
254 {
255 # create a zero
61f5c3f5 256 [ 0 ];
0716bf9b
JH
257 }
258
259sub _one
260 {
261 # create a one
61f5c3f5 262 [ 1 ];
0716bf9b
JH
263 }
264
027dc388
JH
265sub _two
266 {
267 # create a two (for _pow)
61f5c3f5 268 [ 2 ];
027dc388
JH
269 }
270
0716bf9b
JH
271sub _copy
272 {
61f5c3f5 273 [ @{$_[1]} ];
0716bf9b
JH
274 }
275
bd05a461
JH
276# catch and throw away
277sub import { }
278
0716bf9b
JH
279##############################################################################
280# convert back to string and number
281
282sub _str
283 {
284 # (ref to BINT) return num_str
285 # Convert number from internal base 100000 format to string format.
286 # internal format is always normalized (no leading zeros, "-0" => "+0")
574bacfe 287 my $ar = $_[1];
0716bf9b 288 my $ret = "";
61f5c3f5
T
289
290 my $l = scalar @$ar; # number of parts
291 return $nan if $l < 1; # should not happen
292
0716bf9b
JH
293 # handle first one different to strip leading zeros from it (there are no
294 # leading zero parts in internal representation)
61f5c3f5 295 $l --; $ret .= int($ar->[$l]); $l--;
0716bf9b 296 # Interestingly, the pre-padd method uses more time
574bacfe
JH
297 # the old grep variant takes longer (14 to 10 sec)
298 my $z = '0' x ($BASE_LEN-1);
0716bf9b
JH
299 while ($l >= 0)
300 {
574bacfe 301 $ret .= substr($z.$ar->[$l],-$BASE_LEN); # fastest way I could think of
0716bf9b
JH
302 $l--;
303 }
61f5c3f5 304 \$ret;
0716bf9b
JH
305 }
306
307sub _num
308 {
309 # Make a number (scalar int/float) from a BigInt object
574bacfe 310 my $x = $_[1];
0716bf9b
JH
311 return $x->[0] if scalar @$x == 1; # below $BASE
312 my $fac = 1;
313 my $num = 0;
314 foreach (@$x)
315 {
316 $num += $fac*$_; $fac *= $BASE;
317 }
61f5c3f5 318 $num;
0716bf9b
JH
319 }
320
321##############################################################################
322# actual math code
323
324sub _add
325 {
326 # (ref to int_num_array, ref to int_num_array)
574bacfe 327 # routine to add two base 1eX numbers
0716bf9b 328 # stolen from Knuth Vol 2 Algorithm A pg 231
b22b3e31 329 # there are separate routines to add and sub as per Knuth pg 233
0716bf9b
JH
330 # This routine clobbers up array x, but not y.
331
574bacfe 332 my ($c,$x,$y) = @_;
b3abae2a
JH
333
334 return $x if (@$y == 1) && $y->[0] == 0; # $x + 0 => $x
335 if ((@$x == 1) && $x->[0] == 0) # 0 + $y => $y->copy
336 {
337 # twice as slow as $x = [ @$y ], but necc. to retain $x as ref :(
338 @$x = @$y; return $x;
339 }
0716bf9b
JH
340
341 # for each in Y, add Y to X and carry. If after that, something is left in
342 # X, foreach in X add carry to X and then return X, carry
343 # Trades one "$j++" for having to shift arrays, $j could be made integer
b22b3e31 344 # but this would impose a limit to number-length of 2**32.
0716bf9b
JH
345 my $i; my $car = 0; my $j = 0;
346 for $i (@$y)
347 {
e745a66c 348 $x->[$j] -= $BASE if $car = (($x->[$j] += $i + $car) >= $BASE) ? 1 : 0;
0716bf9b
JH
349 $j++;
350 }
351 while ($car != 0)
352 {
353 $x->[$j] -= $BASE if $car = (($x->[$j] += $car) >= $BASE) ? 1 : 0; $j++;
354 }
61f5c3f5 355 $x;
e745a66c
JH
356 }
357
358sub _inc
359 {
360 # (ref to int_num_array, ref to int_num_array)
361 # routine to add 1 to a base 1eX numbers
362 # This routine clobbers up array x, but not y.
363 my ($c,$x) = @_;
364
365 for my $i (@$x)
366 {
367 return $x if (($i += 1) < $BASE); # early out
61f5c3f5 368 $i = 0; # overflow, next
e745a66c 369 }
61f5c3f5
T
370 push @$x,1 if ($x->[-1] == 0); # last overflowed, so extend
371 $x;
e745a66c
JH
372 }
373
374sub _dec
375 {
376 # (ref to int_num_array, ref to int_num_array)
377 # routine to add 1 to a base 1eX numbers
378 # This routine clobbers up array x, but not y.
379 my ($c,$x) = @_;
380
61f5c3f5 381 my $MAX = $BASE-1; # since MAX_VAL based on MBASE
e745a66c
JH
382 for my $i (@$x)
383 {
384 last if (($i -= 1) >= 0); # early out
61f5c3f5 385 $i = $MAX; # overflow, next
e745a66c
JH
386 }
387 pop @$x if $x->[-1] == 0 && @$x > 1; # last overflowed (but leave 0)
61f5c3f5 388 $x;
0716bf9b
JH
389 }
390
391sub _sub
392 {
393 # (ref to int_num_array, ref to int_num_array)
574bacfe 394 # subtract base 1eX numbers -- stolen from Knuth Vol 2 pg 232, $x > $y
b22b3e31 395 # subtract Y from X (X is always greater/equal!) by modifying x in place
574bacfe 396 my ($c,$sx,$sy,$s) = @_;
0716bf9b
JH
397
398 my $car = 0; my $i; my $j = 0;
399 if (!$s)
400 {
401 #print "case 2\n";
402 for $i (@$sx)
403 {
404 last unless defined $sy->[$j] || $car;
0716bf9b 405 $i += $BASE if $car = (($i -= ($sy->[$j] || 0) + $car) < 0); $j++;
0716bf9b
JH
406 }
407 # might leave leading zeros, so fix that
394e6ffb 408 return __strip_zeros($sx);
0716bf9b 409 }
394e6ffb
JH
410 #print "case 1 (swap)\n";
411 for $i (@$sx)
0716bf9b 412 {
394e6ffb
JH
413 last unless defined $sy->[$j] || $car;
414 $sy->[$j] += $BASE
415 if $car = (($sy->[$j] = $i-($sy->[$j]||0) - $car) < 0);
416 $j++;
0716bf9b 417 }
394e6ffb
JH
418 # might leave leading zeros, so fix that
419 __strip_zeros($sy);
0716bf9b
JH
420 }
421
ee15d750 422sub _mul_use_mul
0716bf9b
JH
423 {
424 # (BINT, BINT) return nothing
425 # multiply two numbers in internal representation
b22b3e31 426 # modifies first arg, second need not be different from first
574bacfe 427 my ($c,$xv,$yv) = @_;
dccbb853 428
b3abae2a 429 # shortcut for two very short numbers (improved by Nathan Zook)
61f5c3f5 430 # works also if xv and yv are the same reference
b3abae2a
JH
431 if ((@$xv == 1) && (@$yv == 1))
432 {
433 if (($xv->[0] *= $yv->[0]) >= $MBASE)
434 {
435 $xv->[0] = $xv->[0] - ($xv->[1] = int($xv->[0] * $RBASE)) * $MBASE;
436 };
437 return $xv;
438 }
439 # shortcut for result == 0
440 if ( ((@$xv == 1) && ($xv->[0] == 0)) ||
441 ((@$yv == 1) && ($yv->[0] == 0)) )
442 {
443 @$xv = (0);
444 return $xv;
445 }
446
0716bf9b 447 # since multiplying $x with $x fails, make copy in this case
574bacfe 448 $yv = [@$xv] if "$xv" eq "$yv"; # same references?
61f5c3f5
T
449 if ($LEN_CONVERT != 0)
450 {
451 $c->_to_small($xv); $c->_to_small($yv);
452 }
453
454 my @prod = (); my ($prod,$car,$cty,$xi,$yi);
455
0716bf9b
JH
456 for $xi (@$xv)
457 {
458 $car = 0; $cty = 0;
574bacfe
JH
459
460 # slow variant
461# for $yi (@$yv)
462# {
463# $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
464# $prod[$cty++] =
61f5c3f5 465# $prod - ($car = int($prod * RBASE)) * $MBASE; # see USE_MUL
574bacfe
JH
466# }
467# $prod[$cty] += $car if $car; # need really to check for 0?
468# $xi = shift @prod;
469
470 # faster variant
471 # looping through this if $xi == 0 is silly - so optimize it away!
472 $xi = (shift @prod || 0), next if $xi == 0;
0716bf9b
JH
473 for $yi (@$yv)
474 {
475 $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
574bacfe
JH
476## this is actually a tad slower
477## $prod = $prod[$cty]; $prod += ($car + $xi * $yi); # no ||0 here
0716bf9b 478 $prod[$cty++] =
61f5c3f5 479 $prod - ($car = int($prod * $RBASE)) * $MBASE; # see USE_MUL
0716bf9b
JH
480 }
481 $prod[$cty] += $car if $car; # need really to check for 0?
027dc388 482 $xi = shift @prod || 0; # || 0 makes v5.005_3 happy
0716bf9b 483 }
0716bf9b 484 push @$xv, @prod;
61f5c3f5
T
485 if ($LEN_CONVERT != 0)
486 {
487 $c->_to_large($yv);
488 $c->_to_large($xv);
489 }
490 else
491 {
492 __strip_zeros($xv);
493 }
494 $xv;
0716bf9b
JH
495 }
496
ee15d750
JH
497sub _mul_use_div
498 {
499 # (BINT, BINT) return nothing
500 # multiply two numbers in internal representation
501 # modifies first arg, second need not be different from first
502 my ($c,$xv,$yv) = @_;
503
b3abae2a 504 # shortcut for two very short numbers (improved by Nathan Zook)
61f5c3f5 505 # works also if xv and yv are the same reference
b3abae2a
JH
506 if ((@$xv == 1) && (@$yv == 1))
507 {
508 if (($xv->[0] *= $yv->[0]) >= $MBASE)
509 {
510 $xv->[0] =
511 $xv->[0] - ($xv->[1] = int($xv->[0] / $MBASE)) * $MBASE;
512 };
513 return $xv;
514 }
515 # shortcut for result == 0
516 if ( ((@$xv == 1) && ($xv->[0] == 0)) ||
517 ((@$yv == 1) && ($yv->[0] == 0)) )
518 {
519 @$xv = (0);
520 return $xv;
521 }
522
61f5c3f5 523
ee15d750
JH
524 # since multiplying $x with $x fails, make copy in this case
525 $yv = [@$xv] if "$xv" eq "$yv"; # same references?
61f5c3f5
T
526 if ($LEN_CONVERT != 0)
527 {
528 $c->_to_small($xv); $c->_to_small($yv);
529 }
530
531 my @prod = (); my ($prod,$car,$cty,$xi,$yi);
ee15d750
JH
532 for $xi (@$xv)
533 {
534 $car = 0; $cty = 0;
535 # looping through this if $xi == 0 is silly - so optimize it away!
536 $xi = (shift @prod || 0), next if $xi == 0;
537 for $yi (@$yv)
538 {
539 $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
540 $prod[$cty++] =
61f5c3f5 541 $prod - ($car = int($prod / $MBASE)) * $MBASE;
ee15d750
JH
542 }
543 $prod[$cty] += $car if $car; # need really to check for 0?
027dc388 544 $xi = shift @prod || 0; # || 0 makes v5.005_3 happy
ee15d750
JH
545 }
546 push @$xv, @prod;
61f5c3f5
T
547 if ($LEN_CONVERT != 0)
548 {
549 $c->_to_large($yv);
550 $c->_to_large($xv);
551 }
552 else
553 {
554 __strip_zeros($xv);
555 }
556 $xv;
ee15d750
JH
557 }
558
559sub _div_use_mul
0716bf9b 560 {
b22b3e31 561 # ref to array, ref to array, modify first array and return remainder if
0716bf9b 562 # in list context
574bacfe 563 my ($c,$x,$yorg) = @_;
0716bf9b 564
61f5c3f5
T
565 if (@$x == 1 && @$yorg == 1)
566 {
13a12e00 567 # shortcut, $yorg and $x are two small numbers
61f5c3f5
T
568 if (wantarray)
569 {
570 my $r = [ $x->[0] % $yorg->[0] ];
571 $x->[0] = int($x->[0] / $yorg->[0]);
572 return ($x,$r);
573 }
574 else
575 {
576 $x->[0] = int($x->[0] / $yorg->[0]);
577 return $x;
578 }
579 }
13a12e00
JH
580 #if (@$yorg == 1)
581 # {
582 # # shortcut, $y is < $BASE
583 #
584 # }
585
0716bf9b 586
0716bf9b 587 my $y = [ @$yorg ];
61f5c3f5
T
588 if ($LEN_CONVERT != 0)
589 {
590 $c->_to_small($x); $c->_to_small($y);
591 }
592
593 my ($car,$bar,$prd,$dd,$xi,$yi,@q,$v2,$v1,@d,$tmp,$q,$u2,$u1,$u0);
594
595 $car = $bar = $prd = 0;
596 if (($dd = int($MBASE/($y->[-1]+1))) != 1)
0716bf9b
JH
597 {
598 for $xi (@$x)
599 {
600 $xi = $xi * $dd + $car;
61f5c3f5 601 $xi -= ($car = int($xi * $RBASE)) * $MBASE; # see USE_MUL
0716bf9b
JH
602 }
603 push(@$x, $car); $car = 0;
604 for $yi (@$y)
605 {
606 $yi = $yi * $dd + $car;
61f5c3f5 607 $yi -= ($car = int($yi * $RBASE)) * $MBASE; # see USE_MUL
0716bf9b
JH
608 }
609 }
610 else
611 {
612 push(@$x, 0);
613 }
614 @q = (); ($v2,$v1) = @$y[-2,-1];
615 $v2 = 0 unless $v2;
616 while ($#$x > $#$y)
617 {
618 ($u2,$u1,$u0) = @$x[-3..-1];
619 $u2 = 0 unless $u2;
620 #warn "oups v1 is 0, u0: $u0 $y->[-2] $y->[-1] l ",scalar @$y,"\n"
621 # if $v1 == 0;
61f5c3f5
T
622 $q = (($u0 == $v1) ? $MAX_VAL : int(($u0*$MBASE+$u1)/$v1));
623 --$q while ($v2*$q > ($u0*$MBASE+$u1-$q*$v1)*$MBASE+$u2);
0716bf9b
JH
624 if ($q)
625 {
626 ($car, $bar) = (0,0);
627 for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
628 {
629 $prd = $q * $y->[$yi] + $car;
61f5c3f5
T
630 $prd -= ($car = int($prd * $RBASE)) * $MBASE; # see USE_MUL
631 $x->[$xi] += $MBASE if ($bar = (($x->[$xi] -= $prd + $bar) < 0));
0716bf9b
JH
632 }
633 if ($x->[-1] < $car + $bar)
634 {
635 $car = 0; --$q;
636 for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
637 {
61f5c3f5
T
638 $x->[$xi] -= $MBASE
639 if ($car = (($x->[$xi] += $y->[$yi] + $car) > $MBASE));
0716bf9b
JH
640 }
641 }
642 }
643 pop(@$x); unshift(@q, $q);
644 }
645 if (wantarray)
646 {
647 @d = ();
648 if ($dd != 1)
649 {
650 $car = 0;
651 for $xi (reverse @$x)
652 {
61f5c3f5 653 $prd = $car * $MBASE + $xi;
0716bf9b
JH
654 $car = $prd - ($tmp = int($prd / $dd)) * $dd; # see USE_MUL
655 unshift(@d, $tmp);
656 }
657 }
658 else
659 {
660 @d = @$x;
661 }
662 @$x = @q;
61f5c3f5
T
663 my $d = \@d;
664 if ($LEN_CONVERT != 0)
665 {
666 $c->_to_large($x); $c->_to_large($d);
667 }
668 else
669 {
670 __strip_zeros($x);
671 __strip_zeros($d);
672 }
673 return ($x,$d);
0716bf9b
JH
674 }
675 @$x = @q;
61f5c3f5
T
676 if ($LEN_CONVERT != 0)
677 {
678 $c->_to_large($x);
679 }
680 else
681 {
682 __strip_zeros($x);
683 }
684 $x;
0716bf9b
JH
685 }
686
ee15d750
JH
687sub _div_use_div
688 {
689 # ref to array, ref to array, modify first array and return remainder if
690 # in list context
ee15d750 691 my ($c,$x,$yorg) = @_;
ee15d750 692
61f5c3f5
T
693 if (@$x == 1 && @$yorg == 1)
694 {
13a12e00 695 # shortcut, $yorg and $x are two small numbers
61f5c3f5
T
696 if (wantarray)
697 {
698 my $r = [ $x->[0] % $yorg->[0] ];
699 $x->[0] = int($x->[0] / $yorg->[0]);
700 return ($x,$r);
701 }
702 else
703 {
704 $x->[0] = int($x->[0] / $yorg->[0]);
705 return $x;
706 }
707 }
13a12e00
JH
708# if (@$yorg == 1)
709# {
710# # shortcut, $y is < $BASE
711#
712# }
ee15d750 713
ee15d750 714 my $y = [ @$yorg ];
61f5c3f5
T
715 if ($LEN_CONVERT != 0)
716 {
717 $c->_to_small($x); $c->_to_small($y);
718 }
719
720 my ($car,$bar,$prd,$dd,$xi,$yi,@q,$v2,$v1,@d,$tmp,$q,$u2,$u1,$u0);
721
722 $car = $bar = $prd = 0;
723 if (($dd = int($MBASE/($y->[-1]+1))) != 1)
ee15d750
JH
724 {
725 for $xi (@$x)
726 {
727 $xi = $xi * $dd + $car;
61f5c3f5 728 $xi -= ($car = int($xi / $MBASE)) * $MBASE;
ee15d750
JH
729 }
730 push(@$x, $car); $car = 0;
731 for $yi (@$y)
732 {
733 $yi = $yi * $dd + $car;
61f5c3f5 734 $yi -= ($car = int($yi / $MBASE)) * $MBASE;
ee15d750
JH
735 }
736 }
737 else
738 {
739 push(@$x, 0);
740 }
741 @q = (); ($v2,$v1) = @$y[-2,-1];
742 $v2 = 0 unless $v2;
743 while ($#$x > $#$y)
744 {
745 ($u2,$u1,$u0) = @$x[-3..-1];
746 $u2 = 0 unless $u2;
747 #warn "oups v1 is 0, u0: $u0 $y->[-2] $y->[-1] l ",scalar @$y,"\n"
748 # if $v1 == 0;
61f5c3f5
T
749 $q = (($u0 == $v1) ? $MAX_VAL : int(($u0*$MBASE+$u1)/$v1));
750 --$q while ($v2*$q > ($u0*$MBASE+$u1-$q*$v1)*$MBASE+$u2);
ee15d750
JH
751 if ($q)
752 {
753 ($car, $bar) = (0,0);
754 for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
755 {
756 $prd = $q * $y->[$yi] + $car;
61f5c3f5
T
757 $prd -= ($car = int($prd / $MBASE)) * $MBASE;
758 $x->[$xi] += $MBASE if ($bar = (($x->[$xi] -= $prd + $bar) < 0));
ee15d750
JH
759 }
760 if ($x->[-1] < $car + $bar)
761 {
762 $car = 0; --$q;
763 for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
764 {
61f5c3f5
T
765 $x->[$xi] -= $MBASE
766 if ($car = (($x->[$xi] += $y->[$yi] + $car) > $MBASE));
ee15d750
JH
767 }
768 }
769 }
61f5c3f5 770 pop(@$x); unshift(@q, $q);
ee15d750
JH
771 }
772 if (wantarray)
773 {
774 @d = ();
775 if ($dd != 1)
776 {
777 $car = 0;
778 for $xi (reverse @$x)
779 {
61f5c3f5 780 $prd = $car * $MBASE + $xi;
ee15d750
JH
781 $car = $prd - ($tmp = int($prd / $dd)) * $dd;
782 unshift(@d, $tmp);
783 }
784 }
785 else
786 {
787 @d = @$x;
788 }
789 @$x = @q;
61f5c3f5
T
790 my $d = \@d;
791 if ($LEN_CONVERT != 0)
792 {
793 $c->_to_large($x); $c->_to_large($d);
794 }
795 else
796 {
797 __strip_zeros($x);
798 __strip_zeros($d);
799 }
800 return ($x,$d);
ee15d750
JH
801 }
802 @$x = @q;
61f5c3f5
T
803 if ($LEN_CONVERT != 0)
804 {
805 $c->_to_large($x);
806 }
807 else
808 {
809 __strip_zeros($x);
810 }
811 $x;
ee15d750
JH
812 }
813
394e6ffb
JH
814##############################################################################
815# testing
816
817sub _acmp
818 {
819 # internal absolute post-normalized compare (ignore signs)
820 # ref to array, ref to array, return <0, 0, >0
821 # arrays must have at least one entry; this is not checked for
822
823 my ($c,$cx,$cy) = @_;
824
61f5c3f5 825 # fast comp based on array elements
394e6ffb
JH
826 my $lxy = scalar @$cx - scalar @$cy;
827 return -1 if $lxy < 0; # already differs, ret
828 return 1 if $lxy > 0; # ditto
829
830 # now calculate length based on digits, not parts
831 $lxy = _len($c,$cx) - _len($c,$cy); # difference
832 return -1 if $lxy < 0;
833 return 1 if $lxy > 0;
834
835 # hm, same lengths, but same contents?
836 my $i = 0; my $a;
837 # first way takes 5.49 sec instead of 4.87, but has the early out advantage
838 # so grep is slightly faster, but more inflexible. hm. $_ instead of $k
839 # yields 5.6 instead of 5.5 sec huh?
840 # manual way (abort if unequal, good for early ne)
841 my $j = scalar @$cx - 1;
842 while ($j >= 0)
843 {
844 last if ($a = $cx->[$j] - $cy->[$j]); $j--;
845 }
846 return 1 if $a > 0;
847 return -1 if $a < 0;
61f5c3f5
T
848 0; # equal
849
394e6ffb
JH
850 # while it early aborts, it is even slower than the manual variant
851 #grep { return $a if ($a = $_ - $cy->[$i++]); } @$cx;
852 # grep way, go trough all (bad for early ne)
853 #grep { $a = $_ - $cy->[$i++]; } @$cx;
854 #return $a;
855 }
856
857sub _len
858 {
859 # compute number of digits in bigint, minus the sign
860
861 # int() because add/sub sometimes leaves strings (like '00005') instead of
862 # '5' in this place, thus causing length() to report wrong length
863 my $cx = $_[1];
864
865 return (@$cx-1)*$BASE_LEN+length(int($cx->[-1]));
866 }
867
868sub _digit
869 {
870 # return the nth digit, negative values count backward
871 # zero is rightmost, so _digit(123,0) will give 3
872 my ($c,$x,$n) = @_;
873
874 my $len = _len('',$x);
875
876 $n = $len+$n if $n < 0; # -1 last, -2 second-to-last
877 $n = abs($n); # if negative was too big
878 $len--; $n = $len if $n > $len; # n to big?
879
880 my $elem = int($n / $BASE_LEN); # which array element
881 my $digit = $n % $BASE_LEN; # which digit in this element
882 $elem = '0000'.@$x[$elem]; # get element padded with 0's
883 return substr($elem,-$digit-1,1);
884 }
885
886sub _zeros
887 {
888 # return amount of trailing zeros in decimal
889 # check each array elem in _m for having 0 at end as long as elem == 0
890 # Upon finding a elem != 0, stop
891 my $x = $_[1];
892 my $zeros = 0; my $elem;
893 foreach my $e (@$x)
894 {
895 if ($e != 0)
896 {
897 $elem = "$e"; # preserve x
898 $elem =~ s/.*?(0*$)/$1/; # strip anything not zero
899 $zeros *= $BASE_LEN; # elems * 5
61f5c3f5 900 $zeros += length($elem); # count trailing zeros
394e6ffb
JH
901 last; # early out
902 }
903 $zeros ++; # real else branch: 50% slower!
904 }
61f5c3f5 905 $zeros;
394e6ffb
JH
906 }
907
908##############################################################################
909# _is_* routines
910
911sub _is_zero
912 {
913 # return true if arg (BINT or num_str) is zero (array '+', '0')
914 my $x = $_[1];
61f5c3f5
T
915
916 (((scalar @$x == 1) && ($x->[0] == 0))) <=> 0;
394e6ffb
JH
917 }
918
919sub _is_even
920 {
921 # return true if arg (BINT or num_str) is even
922 my $x = $_[1];
61f5c3f5 923 (!($x->[0] & 1)) <=> 0;
394e6ffb
JH
924 }
925
926sub _is_odd
927 {
928 # return true if arg (BINT or num_str) is even
929 my $x = $_[1];
61f5c3f5
T
930
931 (($x->[0] & 1)) <=> 0;
394e6ffb
JH
932 }
933
934sub _is_one
935 {
936 # return true if arg (BINT or num_str) is one (array '+', '1')
937 my $x = $_[1];
61f5c3f5
T
938
939 (scalar @$x == 1) && ($x->[0] == 1) <=> 0;
394e6ffb
JH
940 }
941
942sub __strip_zeros
943 {
944 # internal normalization function that strips leading zeros from the array
945 # args: ref to array
946 my $s = shift;
947
948 my $cnt = scalar @$s; # get count of parts
949 my $i = $cnt-1;
950 push @$s,0 if $i < 0; # div might return empty results, so fix it
951
61f5c3f5
T
952 return $s if @$s == 1; # early out
953
394e6ffb
JH
954 #print "strip: cnt $cnt i $i\n";
955 # '0', '3', '4', '0', '0',
956 # 0 1 2 3 4
957 # cnt = 5, i = 4
958 # i = 4
959 # i = 3
960 # => fcnt = cnt - i (5-2 => 3, cnt => 5-1 = 4, throw away from 4th pos)
961 # >= 1: skip first part (this can be zero)
962 while ($i > 0) { last if $s->[$i] != 0; $i--; }
963 $i++; splice @$s,$i if ($i < $cnt); # $i cant be 0
964 $s;
965 }
966
967###############################################################################
968# check routine to test internal state of corruptions
969
970sub _check
971 {
972 # used by the test suite
973 my $x = $_[1];
974
975 return "$x is not a reference" if !ref($x);
976
977 # are all parts are valid?
978 my $i = 0; my $j = scalar @$x; my ($e,$try);
979 while ($i < $j)
980 {
981 $e = $x->[$i]; $e = 'undef' unless defined $e;
982 $try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e)";
983 last if $e !~ /^[+]?[0-9]+$/;
984 $try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e) (stringify)";
985 last if "$e" !~ /^[+]?[0-9]+$/;
986 $try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e) (cat-stringify)";
987 last if '' . "$e" !~ /^[+]?[0-9]+$/;
988 $try = ' < 0 || >= $BASE; '."($x, $e)";
989 last if $e <0 || $e >= $BASE;
990 # this test is disabled, since new/bnorm and certain ops (like early out
991 # in add/sub) are allowed/expected to leave '00000' in some elements
992 #$try = '=~ /^00+/; '."($x, $e)";
993 #last if $e =~ /^00+/;
994 $i++;
995 }
996 return "Illegal part '$e' at pos $i (tested: $try)" if $i < $j;
997 return 0;
998 }
999
1000
1001###############################################################################
1002###############################################################################
1003# some optional routines to make BigInt faster
1004
dccbb853
JH
1005sub _mod
1006 {
1007 # if possible, use mod shortcut
1008 my ($c,$x,$yo) = @_;
1009
1010 # slow way since $y to big
1011 if (scalar @$yo > 1)
1012 {
1013 my ($xo,$rem) = _div($c,$x,$yo);
1014 return $rem;
1015 }
1016 my $y = $yo->[0];
027dc388 1017 # both are single element arrays
dccbb853
JH
1018 if (scalar @$x == 1)
1019 {
1020 $x->[0] %= $y;
1021 return $x;
1022 }
1023
61f5c3f5 1024 # @y is single element, but @x has more than one
dccbb853
JH
1025 my $b = $BASE % $y;
1026 if ($b == 0)
1027 {
1028 # when BASE % Y == 0 then (B * BASE) % Y == 0
1029 # (B * BASE) % $y + A % Y => A % Y
1030 # so need to consider only last element: O(1)
1031 $x->[0] %= $y;
1032 }
027dc388
JH
1033 elsif ($b == 1)
1034 {
1035 # else need to go trough all elements: O(N), but loop is a bit simplified
1036 my $r = 0;
1037 foreach (@$x)
1038 {
1039 $r += $_ % $y;
1040 $r %= $y;
1041 }
1042 $r = 0 if $r == $y;
1043 $x->[0] = $r;
1044 }
dccbb853
JH
1045 else
1046 {
027dc388
JH
1047 # else need to go trough all elements: O(N)
1048 my $r = 0; my $bm = 1;
1049 foreach (@$x)
1050 {
1051 $r += ($_ % $y) * $bm;
1052 $bm *= $b;
1053 $bm %= $y;
1054 $r %= $y;
1055 }
1056 $r = 0 if $r == $y;
1057 $x->[0] = $r;
dccbb853
JH
1058 }
1059 splice (@$x,1);
61f5c3f5 1060 $x;
dccbb853
JH
1061 }
1062
0716bf9b 1063##############################################################################
574bacfe
JH
1064# shifts
1065
1066sub _rsft
1067 {
1068 my ($c,$x,$y,$n) = @_;
1069
1070 if ($n != 10)
1071 {
61f5c3f5
T
1072 $n = _new($c,\$n); return _div($c,$x, _pow($c,$n,$y));
1073 }
1074
1075 # shortcut (faster) for shifting by 10)
1076 # multiples of $BASE_LEN
1077 my $dst = 0; # destination
1078 my $src = _num($c,$y); # as normal int
1079 my $rem = $src % $BASE_LEN; # remainder to shift
1080 $src = int($src / $BASE_LEN); # source
1081 if ($rem == 0)
1082 {
1083 splice (@$x,0,$src); # even faster, 38.4 => 39.3
574bacfe
JH
1084 }
1085 else
1086 {
61f5c3f5
T
1087 my $len = scalar @$x - $src; # elems to go
1088 my $vd; my $z = '0'x $BASE_LEN;
1089 $x->[scalar @$x] = 0; # avoid || 0 test inside loop
1090 while ($dst < $len)
574bacfe 1091 {
61f5c3f5
T
1092 $vd = $z.$x->[$src];
1093 $vd = substr($vd,-$BASE_LEN,$BASE_LEN-$rem);
1094 $src++;
1095 $vd = substr($z.$x->[$src],-$rem,$rem) . $vd;
1096 $vd = substr($vd,-$BASE_LEN,$BASE_LEN) if length($vd) > $BASE_LEN;
1097 $x->[$dst] = int($vd);
1098 $dst++;
574bacfe 1099 }
61f5c3f5
T
1100 splice (@$x,$dst) if $dst > 0; # kill left-over array elems
1101 pop @$x if $x->[-1] == 0; # kill last element if 0
1102 } # else rem == 0
574bacfe
JH
1103 $x;
1104 }
1105
1106sub _lsft
1107 {
1108 my ($c,$x,$y,$n) = @_;
1109
1110 if ($n != 10)
1111 {
61f5c3f5 1112 $n = _new($c,\$n); return _mul($c,$x, _pow($c,$n,$y));
574bacfe 1113 }
61f5c3f5
T
1114
1115 # shortcut (faster) for shifting by 10) since we are in base 10eX
1116 # multiples of $BASE_LEN:
1117 my $src = scalar @$x; # source
1118 my $len = _num($c,$y); # shift-len as normal int
1119 my $rem = $len % $BASE_LEN; # remainder to shift
1120 my $dst = $src + int($len/$BASE_LEN); # destination
1121 my $vd; # further speedup
1122 $x->[$src] = 0; # avoid first ||0 for speed
1123 my $z = '0' x $BASE_LEN;
1124 while ($src >= 0)
574bacfe 1125 {
61f5c3f5
T
1126 $vd = $x->[$src]; $vd = $z.$vd;
1127 $vd = substr($vd,-$BASE_LEN+$rem,$BASE_LEN-$rem);
1128 $vd .= $src > 0 ? substr($z.$x->[$src-1],-$BASE_LEN,$rem) : '0' x $rem;
1129 $vd = substr($vd,-$BASE_LEN,$BASE_LEN) if length($vd) > $BASE_LEN;
1130 $x->[$dst] = int($vd);
1131 $dst--; $src--;
574bacfe 1132 }
61f5c3f5
T
1133 # set lowest parts to 0
1134 while ($dst >= 0) { $x->[$dst--] = 0; }
1135 # fix spurios last zero element
1136 splice @$x,-1 if $x->[-1] == 0;
574bacfe
JH
1137 $x;
1138 }
1139
027dc388
JH
1140sub _pow
1141 {
1142 # power of $x to $y
1143 # ref to array, ref to array, return ref to array
1144 my ($c,$cx,$cy) = @_;
1145
1146 my $pow2 = _one();
1147 my $two = _two();
1148 my $y1 = _copy($c,$cy);
1149 while (!_is_one($c,$y1))
1150 {
1151 _mul($c,$pow2,$cx) if _is_odd($c,$y1);
1152 _div($c,$y1,$two);
1153 _mul($c,$cx,$cx);
1154 }
1155 _mul($c,$cx,$pow2) unless _is_one($c,$pow2);
61f5c3f5 1156 $cx;
027dc388
JH
1157 }
1158
b3abae2a
JH
1159sub _fac
1160 {
1161 # factorial of $x
1162 # ref to array, return ref to array
1163 my ($c,$cx) = @_;
1164
1165 if ((@$cx == 1) && ($cx->[0] <= 2))
1166 {
1167 $cx->[0] = 1 * ($cx->[0]||1); # 0,1 => 1, 2 => 2
1168 return $cx;
1169 }
1170
1171 # go forward until $base is exceeded
1172 # limit is either $x or $base (x == 100 means as result too high)
1173 my $steps = 100; $steps = $cx->[0] if @$cx == 1;
1174 my $r = 2; my $cf = 3; my $step = 1; my $last = $r;
1175 while ($r < $BASE && $step < $steps)
1176 {
1177 $last = $r; $r *= $cf++; $step++;
1178 }
1179 if ((@$cx == 1) && ($step == $cx->[0]))
1180 {
1181 # completely done
1182 $cx = [$last];
1183 return $cx;
1184 }
1185 my $n = _copy($c,$cx);
1186 $cx = [$last];
1187
1188 #$cx = _one();
1189 while (!(@$n == 1 && $n->[0] == $step))
1190 {
1191 _mul($c,$cx,$n); _dec($c,$n);
1192 }
1193 $cx;
1194 }
1195
1196use constant DEBUG => 0;
1197
1198my $steps = 0;
1199
1200sub steps { $steps };
1201
1202sub _sqrt
0716bf9b 1203 {
394e6ffb
JH
1204 # square-root of $x
1205 # ref to array, return ref to array
1206 my ($c,$x) = @_;
0716bf9b 1207
394e6ffb
JH
1208 if (scalar @$x == 1)
1209 {
1210 # fit's into one Perl scalar
1211 $x->[0] = int(sqrt($x->[0]));
1212 return $x;
1213 }
1214 my $y = _copy($c,$x);
b3abae2a
JH
1215 # hopefully _len/2 is < $BASE, the -1 is to always undershot the guess
1216 # since our guess will "grow"
1217 my $l = int((_len($c,$x)-1) / 2);
1218
1219 my $lastelem = $x->[-1]; # for guess
1220 my $elems = scalar @$x - 1;
1221 # not enough digits, but could have more?
1222 if ((length($lastelem) <= 3) && ($elems > 1))
1223 {
1224 # right-align with zero pad
1225 my $len = length($lastelem) & 1;
1226 print "$lastelem => " if DEBUG;
1227 $lastelem .= substr($x->[-2] . '0' x $BASE_LEN,0,$BASE_LEN);
1228 # former odd => make odd again, or former even to even again
1229 $lastelem = $lastelem / 10 if (length($lastelem) & 1) != $len;
1230 print "$lastelem\n" if DEBUG;
1231 }
0716bf9b 1232
61f5c3f5
T
1233 # construct $x (instead of _lsft($c,$x,$l,10)
1234 my $r = $l % $BASE_LEN; # 10000 00000 00000 00000 ($BASE_LEN=5)
1235 $l = int($l / $BASE_LEN);
b3abae2a
JH
1236 print "l = $l " if DEBUG;
1237
1238 splice @$x,$l; # keep ref($x), but modify it
1239
1240 # we make the first part of the guess not '1000...0' but int(sqrt($lastelem))
1241 # that gives us:
1242 # 14400 00000 => sqrt(14400) => 120
1243 # 144000 000000 => sqrt(144000) => 379
1244
1245 # $x->[$l--] = int('1' . '0' x $r); # old way of guessing
1246 print "$lastelem (elems $elems) => " if DEBUG;
1247 $lastelem = $lastelem / 10 if ($elems & 1 == 1); # odd or even?
1248 my $g = sqrt($lastelem); $g =~ s/\.//; # 2.345 => 2345
1249 $r -= 1 if $elems & 1 == 0; # 70 => 7
1250
1251 # padd with zeros if result is too short
1252 $x->[$l--] = int(substr($g . '0' x $r,0,$r+1));
1253 print "now ",$x->[-1] if DEBUG;
1254 print " would have been ", int('1' . '0' x $r),"\n" if DEBUG;
1255
1256 # If @$x > 1, we could compute the second elem of the guess, too, to create
1257 # an even better guess. Not implemented yet.
1258 $x->[$l--] = 0 while ($l >= 0); # all other digits of guess are zero
61f5c3f5 1259
b3abae2a 1260 print "start x= ",${_str($c,$x)},"\n" if DEBUG;
394e6ffb
JH
1261 my $two = _two();
1262 my $last = _zero();
1263 my $lastlast = _zero();
b3abae2a 1264 $steps = 0 if DEBUG;
394e6ffb
JH
1265 while (_acmp($c,$last,$x) != 0 && _acmp($c,$lastlast,$x) != 0)
1266 {
b3abae2a 1267 $steps++ if DEBUG;
394e6ffb
JH
1268 $lastlast = _copy($c,$last);
1269 $last = _copy($c,$x);
1270 _add($c,$x, _div($c,_copy($c,$y),$x));
1271 _div($c,$x, $two );
b3abae2a 1272 print " x= ",${_str($c,$x)},"\n" if DEBUG;
394e6ffb 1273 }
b3abae2a 1274 print "\nsteps in sqrt: $steps, " if DEBUG;
394e6ffb 1275 _dec($c,$x) if _acmp($c,$y,_mul($c,_copy($c,$x),$x)) < 0; # overshot?
b3abae2a 1276 print " final ",$x->[-1],"\n" if DEBUG;
394e6ffb 1277 $x;
0716bf9b
JH
1278 }
1279
394e6ffb
JH
1280##############################################################################
1281# binary stuff
0716bf9b 1282
394e6ffb
JH
1283sub _and
1284 {
1285 my ($c,$x,$y) = @_;
0716bf9b 1286
394e6ffb
JH
1287 # the shortcut makes equal, large numbers _really_ fast, and makes only a
1288 # very small performance drop for small numbers (e.g. something with less
1289 # than 32 bit) Since we optimize for large numbers, this is enabled.
1290 return $x if _acmp($c,$x,$y) == 0; # shortcut
0716bf9b 1291
394e6ffb
JH
1292 my $m = _one(); my ($xr,$yr);
1293 my $mask = $AND_MASK;
1294
1295 my $x1 = $x;
1296 my $y1 = _copy($c,$y); # make copy
1297 $x = _zero();
1298 my ($b,$xrr,$yrr);
1299 use integer;
1300 while (!_is_zero($c,$x1) && !_is_zero($c,$y1))
1301 {
1302 ($x1, $xr) = _div($c,$x1,$mask);
1303 ($y1, $yr) = _div($c,$y1,$mask);
1304
1305 # make ints() from $xr, $yr
1306 # this is when the AND_BITS are greater tahn $BASE and is slower for
1307 # small (<256 bits) numbers, but faster for large numbers. Disabled
1308 # due to KISS principle
1309
1310# $b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
1311# $b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
1312# _add($c,$x, _mul($c, _new( $c, \($xrr & $yrr) ), $m) );
1313
61f5c3f5
T
1314 # 0+ due to '&' doesn't work in strings
1315 _add($c,$x, _mul($c, [ 0+$xr->[0] & 0+$yr->[0] ], $m) );
394e6ffb
JH
1316 _mul($c,$m,$mask);
1317 }
1318 $x;
0716bf9b
JH
1319 }
1320
394e6ffb 1321sub _xor
0716bf9b 1322 {
394e6ffb
JH
1323 my ($c,$x,$y) = @_;
1324
1325 return _zero() if _acmp($c,$x,$y) == 0; # shortcut (see -and)
1326
1327 my $m = _one(); my ($xr,$yr);
1328 my $mask = $XOR_MASK;
1329
1330 my $x1 = $x;
1331 my $y1 = _copy($c,$y); # make copy
1332 $x = _zero();
1333 my ($b,$xrr,$yrr);
1334 use integer;
1335 while (!_is_zero($c,$x1) && !_is_zero($c,$y1))
0716bf9b 1336 {
394e6ffb
JH
1337 ($x1, $xr) = _div($c,$x1,$mask);
1338 ($y1, $yr) = _div($c,$y1,$mask);
1339 # make ints() from $xr, $yr (see _and())
1340 #$b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
1341 #$b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
1342 #_add($c,$x, _mul($c, _new( $c, \($xrr ^ $yrr) ), $m) );
61f5c3f5
T
1343
1344 # 0+ due to '^' doesn't work in strings
1345 _add($c,$x, _mul($c, [ 0+$xr->[0] ^ 0+$yr->[0] ], $m) );
394e6ffb 1346 _mul($c,$m,$mask);
0716bf9b 1347 }
394e6ffb
JH
1348 # the loop stops when the shorter of the two numbers is exhausted
1349 # the remainder of the longer one will survive bit-by-bit, so we simple
1350 # multiply-add it in
1351 _add($c,$x, _mul($c, $x1, $m) ) if !_is_zero($c,$x1);
1352 _add($c,$x, _mul($c, $y1, $m) ) if !_is_zero($c,$y1);
1353
1354 $x;
0716bf9b
JH
1355 }
1356
394e6ffb 1357sub _or
0716bf9b 1358 {
394e6ffb 1359 my ($c,$x,$y) = @_;
0716bf9b 1360
394e6ffb 1361 return $x if _acmp($c,$x,$y) == 0; # shortcut (see _and)
0716bf9b 1362
394e6ffb
JH
1363 my $m = _one(); my ($xr,$yr);
1364 my $mask = $OR_MASK;
0716bf9b 1365
394e6ffb
JH
1366 my $x1 = $x;
1367 my $y1 = _copy($c,$y); # make copy
1368 $x = _zero();
1369 my ($b,$xrr,$yrr);
1370 use integer;
1371 while (!_is_zero($c,$x1) && !_is_zero($c,$y1))
1372 {
1373 ($x1, $xr) = _div($c,$x1,$mask);
1374 ($y1, $yr) = _div($c,$y1,$mask);
1375 # make ints() from $xr, $yr (see _and())
1376# $b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
1377# $b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
1378# _add($c,$x, _mul($c, _new( $c, \($xrr | $yrr) ), $m) );
1379
61f5c3f5
T
1380 # 0+ due to '|' doesn't work in strings
1381 _add($c,$x, _mul($c, [ 0+$xr->[0] | 0+$yr->[0] ], $m) );
394e6ffb
JH
1382 _mul($c,$m,$mask);
1383 }
1384 # the loop stops when the shorter of the two numbers is exhausted
1385 # the remainder of the longer one will survive bit-by-bit, so we simple
1386 # multiply-add it in
1387 _add($c,$x, _mul($c, $x1, $m) ) if !_is_zero($c,$x1);
1388 _add($c,$x, _mul($c, $y1, $m) ) if !_is_zero($c,$y1);
1389
1390 $x;
0716bf9b
JH
1391 }
1392
61f5c3f5
T
1393sub _as_hex
1394 {
1395 # convert a decimal number to hex (ref to array, return ref to string)
1396 my ($c,$x) = @_;
1397
1398 my $x1 = _copy($c,$x);
1399
1400 my $es = '';
1401 my $xr;
1402 my $x10000 = [ 0x10000 ];
1403 while (! _is_zero($c,$x1))
1404 {
1405 ($x1, $xr) = _div($c,$x1,$x10000);
1406 $es .= unpack('h4',pack('v',$xr->[0]));
1407 }
1408 $es = reverse $es;
1409 $es =~ s/^[0]+//; # strip leading zeros
1410 $es = '0x' . $es;
1411 \$es;
1412 }
1413
1414sub _as_bin
1415 {
1416 # convert a decimal number to bin (ref to array, return ref to string)
1417 my ($c,$x) = @_;
1418
1419 my $x1 = _copy($c,$x);
1420
1421 my $es = '';
1422 my $xr;
1423 my $x10000 = [ 0x10000 ];
1424 while (! _is_zero($c,$x1))
1425 {
1426 ($x1, $xr) = _div($c,$x1,$x10000);
1427 $es .= unpack('b16',pack('v',$xr->[0]));
1428 }
1429 $es = reverse $es;
1430 $es =~ s/^[0]+//; # strip leading zeros
1431 $es = '0b' . $es;
1432 \$es;
1433 }
1434
394e6ffb 1435sub _from_hex
0716bf9b 1436 {
394e6ffb
JH
1437 # convert a hex number to decimal (ref to string, return ref to array)
1438 my ($c,$hs) = @_;
0716bf9b 1439
394e6ffb
JH
1440 my $mul = _one();
1441 my $m = [ 0x10000 ]; # 16 bit at a time
1442 my $x = _zero();
0716bf9b 1443
61f5c3f5 1444 my $len = length($$hs)-2;
394e6ffb
JH
1445 $len = int($len/4); # 4-digit parts, w/o '0x'
1446 my $val; my $i = -4;
1447 while ($len >= 0)
1448 {
1449 $val = substr($$hs,$i,4);
1450 $val =~ s/^[+-]?0x// if $len == 0; # for last part only because
1451 $val = hex($val); # hex does not like wrong chars
1452 $i -= 4; $len --;
1453 _add ($c, $x, _mul ($c, [ $val ], $mul ) ) if $val != 0;
1454 _mul ($c, $mul, $m ) if $len >= 0; # skip last mul
1455 }
1456 $x;
1457 }
1458
1459sub _from_bin
0716bf9b 1460 {
394e6ffb
JH
1461 # convert a hex number to decimal (ref to string, return ref to array)
1462 my ($c,$bs) = @_;
0716bf9b 1463
13a12e00
JH
1464 # instead of converting 8 bit at a time, it is faster to convert the
1465 # number to hex, and then call _from_hex.
1466
1467 my $hs = $$bs;
1468 $hs =~ s/^[+-]?0b//; # remove sign and 0b
1469 my $l = length($hs); # bits
1470 $hs = '0' x (8-($l % 8)) . $hs if ($l % 8) != 0; # padd left side w/ 0
1471 my $h = unpack('H*', pack ('B*', $hs)); # repack as hex
1472 return $c->_from_hex(\('0x'.$h));
1473
394e6ffb
JH
1474 my $mul = _one();
1475 my $m = [ 0x100 ]; # 8 bit at a time
1476 my $x = _zero();
0716bf9b 1477
61f5c3f5 1478 my $len = length($$bs)-2;
394e6ffb
JH
1479 $len = int($len/8); # 4-digit parts, w/o '0x'
1480 my $val; my $i = -8;
1481 while ($len >= 0)
0716bf9b 1482 {
394e6ffb
JH
1483 $val = substr($$bs,$i,8);
1484 $val =~ s/^[+-]?0b// if $len == 0; # for last part only
1485
394e6ffb
JH
1486 $val = ord(pack('B8',substr('00000000'.$val,-8,8)));
1487
1488 $i -= 8; $len --;
1489 _add ($c, $x, _mul ($c, [ $val ], $mul ) ) if $val != 0;
1490 _mul ($c, $mul, $m ) if $len >= 0; # skip last mul
0716bf9b 1491 }
394e6ffb 1492 $x;
0716bf9b
JH
1493 }
1494
394e6ffb
JH
1495##############################################################################
1496##############################################################################
1497
0716bf9b
JH
14981;
1499__END__
1500
1501=head1 NAME
1502
1503Math::BigInt::Calc - Pure Perl module to support Math::BigInt
1504
1505=head1 SYNOPSIS
1506
ee15d750
JH
1507Provides support for big integer calculations. Not intended to be used by other
1508modules (except Math::BigInt::Cached). Other modules which sport the same
1509functions can also be used to support Math::Bigint, like Math::BigInt::Pari.
0716bf9b
JH
1510
1511=head1 DESCRIPTION
1512
027dc388
JH
1513In order to allow for multiple big integer libraries, Math::BigInt was
1514rewritten to use library modules for core math routines. Any module which
1515follows the same API as this can be used instead by using the following:
0716bf9b 1516
ee15d750 1517 use Math::BigInt lib => 'libname';
0716bf9b 1518
027dc388
JH
1519'libname' is either the long name ('Math::BigInt::Pari'), or only the short
1520version like 'Pari'.
1521
0716bf9b
JH
1522=head1 EXPORT
1523
027dc388
JH
1524The following functions MUST be defined in order to support the use by
1525Math::BigInt:
0716bf9b
JH
1526
1527 _new(string) return ref to new object from ref to decimal string
1528 _zero() return a new object with value 0
1529 _one() return a new object with value 1
1530
1531 _str(obj) return ref to a string representing the object
1532 _num(obj) returns a Perl integer/floating point number
1533 NOTE: because of Perl numeric notation defaults,
1534 the _num'ified obj may lose accuracy due to
1535 machine-dependend floating point size limitations
1536
1537 _add(obj,obj) Simple addition of two objects
1538 _mul(obj,obj) Multiplication of two objects
1539 _div(obj,obj) Division of the 1st object by the 2nd
b22b3e31
PN
1540 In list context, returns (result,remainder).
1541 NOTE: this is integer math, so no
1542 fractional part will be returned.
1543 _sub(obj,obj) Simple subtraction of 1 object from another
0716bf9b
JH
1544 a third, optional parameter indicates that the params
1545 are swapped. In this case, the first param needs to
1546 be preserved, while you can destroy the second.
1547 sub (x,y,1) => return x - y and keep x intact!
e745a66c
JH
1548 _dec(obj) decrement object by one (input is garant. to be > 0)
1549 _inc(obj) increment object by one
1550
0716bf9b
JH
1551
1552 _acmp(obj,obj) <=> operator for objects (return -1, 0 or 1)
1553
1554 _len(obj) returns count of the decimal digits of the object
1555 _digit(obj,n) returns the n'th decimal digit of object
1556
1557 _is_one(obj) return true if argument is +1
1558 _is_zero(obj) return true if argument is 0
1559 _is_even(obj) return true if argument is even (0,2,4,6..)
1560 _is_odd(obj) return true if argument is odd (1,3,5,7..)
1561
1562 _copy return a ref to a true copy of the object
1563
1564 _check(obj) check whether internal representation is still intact
1565 return 0 for ok, otherwise error message as string
1566
bd05a461 1567The following functions are optional, and can be defined if the underlying lib
027dc388
JH
1568has a fast way to do them. If undefined, Math::BigInt will use pure Perl (hence
1569slow) fallback routines to emulate these:
0716bf9b
JH
1570
1571 _from_hex(str) return ref to new object from ref to hexadecimal string
1572 _from_bin(str) return ref to new object from ref to binary string
1573
ee15d750
JH
1574 _as_hex(str) return ref to scalar string containing the value as
1575 unsigned hex string, with the '0x' prepended.
1576 Leading zeros must be stripped.
1577 _as_bin(str) Like as_hex, only as binary string containing only
1578 zeros and ones. Leading zeros must be stripped and a
1579 '0b' must be prepended.
1580
0716bf9b 1581 _rsft(obj,N,B) shift object in base B by N 'digits' right
dccbb853 1582 For unsupported bases B, return undef to signal failure
0716bf9b 1583 _lsft(obj,N,B) shift object in base B by N 'digits' left
dccbb853 1584 For unsupported bases B, return undef to signal failure
0716bf9b
JH
1585
1586 _xor(obj1,obj2) XOR (bit-wise) object 1 with object 2
dccbb853 1587 Note: XOR, AND and OR pad with zeros if size mismatches
0716bf9b
JH
1588 _and(obj1,obj2) AND (bit-wise) object 1 with object 2
1589 _or(obj1,obj2) OR (bit-wise) object 1 with object 2
1590
dccbb853 1591 _mod(obj,obj) Return remainder of div of the 1st by the 2nd object
394e6ffb 1592 _sqrt(obj) return the square root of object (truncate to int)
b3abae2a 1593 _fac(obj) return factorial of object 1 (1*2*3*4..)
0716bf9b
JH
1594 _pow(obj,obj) return object 1 to the power of object 2
1595 _gcd(obj,obj) return Greatest Common Divisor of two objects
1596
b22b3e31 1597 _zeros(obj) return number of trailing decimal zeros
0716bf9b 1598
b22b3e31 1599Input strings come in as unsigned but with prefix (i.e. as '123', '0xabc'
0716bf9b
JH
1600or '0b1101').
1601
b22b3e31 1602Testing of input parameter validity is done by the caller, so you need not
574bacfe
JH
1603worry about underflow (f.i. in C<_sub()>, C<_dec()>) nor about division by
1604zero or similar cases.
1605
1606The first parameter can be modified, that includes the possibility that you
1607return a reference to a completely different object instead. Although keeping
dccbb853
JH
1608the reference and just changing it's contents is prefered over creating and
1609returning a different reference.
574bacfe
JH
1610
1611Return values are always references to objects or strings. Exceptions are
1612C<_lsft()> and C<_rsft()>, which return undef if they can not shift the
027dc388
JH
1613argument. This is used to delegate shifting of bases different than the one
1614you can support back to Math::BigInt, which will use some generic code to
1615calculate the result.
574bacfe
JH
1616
1617=head1 WRAP YOUR OWN
1618
1619If you want to port your own favourite c-lib for big numbers to the
1620Math::BigInt interface, you can take any of the already existing modules as
1621a rough guideline. You should really wrap up the latest BigInt and BigFloat
bd05a461 1622testsuites with your module, and replace in them any of the following:
574bacfe
JH
1623
1624 use Math::BigInt;
1625
bd05a461 1626by this:
574bacfe
JH
1627
1628 use Math::BigInt lib => 'yourlib';
1629
1630This way you ensure that your library really works 100% within Math::BigInt.
0716bf9b
JH
1631
1632=head1 LICENSE
1633
1634This program is free software; you may redistribute it and/or modify it under
1635the same terms as Perl itself.
1636
1637=head1 AUTHORS
1638
1639Original math code by Mark Biggar, rewritten by Tels L<http://bloodgate.com/>
1640in late 2000, 2001.
1641Seperated from BigInt and shaped API with the help of John Peacock.
1642
1643=head1 SEE ALSO
1644
ee15d750
JH
1645L<Math::BigInt>, L<Math::BigFloat>, L<Math::BigInt::BitVect>,
1646L<Math::BigInt::GMP>, L<Math::BigInt::Cached> and L<Math::BigInt::Pari>.
0716bf9b
JH
1647
1648=cut