Upgrade to Math::BigInt 1.51.
[p5sagit/p5-mst-13.2.git] / lib / Math / BigInt / Calc.pm
CommitLineData
0716bf9b 1package Math::BigInt::Calc;
2
3use 5.005;
4use strict;
574bacfe 5# use warnings; # dont use warnings for older Perls
0716bf9b 6
7require Exporter;
bd05a461 8use vars qw/@ISA $VERSION/;
0716bf9b 9@ISA = qw(Exporter);
10
b3abae2a 11$VERSION = '0.22';
0716bf9b 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 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 18
19# todo:
20# - fully remove funky $# stuff (maybe)
0716bf9b 21
22# USE_MUL: due to problems on certain os (os390, posix-bc) "* 1e-5" is used
ee15d750 23# instead of "/ 1e5" at some places, (marked with USE_MUL). Other platforms
24# BS2000, some Crays need USE_DIV instead.
bd05a461 25# The BEGIN block is used to determine which of the two variants gives the
26# correct result.
0716bf9b 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 34my ($AND_BITS,$XOR_BITS,$OR_BITS);
35my ($AND_MASK,$XOR_MASK,$OR_MASK);
61f5c3f5 36my ($LEN_CONVERT);
ee15d750 37
38sub _base_len
39 {
dccbb853 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 42 shift;
43
ee15d750 44 my $b = shift;
45 if (defined $b)
46 {
61f5c3f5 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 53 $MBASE = int("1e".$BASE_LEN_SMALL);
54 $RBASE = abs('1e-'.$BASE_LEN_SMALL); # see USE_MUL
394e6ffb 55 $caught = 0;
61f5c3f5 56 $caught += 1 if (int($MBASE * $RBASE) != 1); # should be 1
57 $caught += 2 if (int($MBASE / $MBASE) != 1); # should be 1
394e6ffb 58 last if $caught != 3;
59 }
61f5c3f5 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 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 76 {
77 # must USE_MUL
ee15d750 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 83 # can USE_DIV instead
84 *{_mul} = \&_mul_use_div;
85 *{_div} = \&_div_use_div;
86 }
87 }
61f5c3f5 88 return $BASE_LEN unless wantarray;
89 return ($BASE_LEN, $AND_BITS, $XOR_BITS, $OR_BITS, $BASE_LEN_SMALL, $MAX_VAL);
ee15d750 90 }
574bacfe 91
92BEGIN
93 {
bd05a461 94 # from Daniel Pfeiffer: determine largest group of digits that is precisely
574bacfe 95 # multipliable with itself plus carry
dccbb853 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 98 do
99 {
100 $num = ('9' x ++$e) + 0;
394e6ffb 101 $num *= $num + 1.0;
394e6ffb 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 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 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 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 163 }
164
0716bf9b 165##############################################################################
61f5c3f5 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 233
234sub _new
235 {
394e6ffb 236 # (ref to string) return ref to num_array
0716bf9b 237 # Convert a number from string format to internal base 100000 format.
238 # Assumes normalized value as input.
574bacfe 239 my $d = $_[1];
61f5c3f5 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 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 252
253sub _zero
254 {
255 # create a zero
61f5c3f5 256 [ 0 ];
0716bf9b 257 }
258
259sub _one
260 {
261 # create a one
61f5c3f5 262 [ 1 ];
0716bf9b 263 }
264
027dc388 265sub _two
266 {
267 # create a two (for _pow)
61f5c3f5 268 [ 2 ];
027dc388 269 }
270
0716bf9b 271sub _copy
272 {
61f5c3f5 273 [ @{$_[1]} ];
0716bf9b 274 }
275
bd05a461 276# catch and throw away
277sub import { }
278
0716bf9b 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 289
290 my $l = scalar @$ar; # number of parts
291 return $nan if $l < 1; # should not happen
292
0716bf9b 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 297 # the old grep variant takes longer (14 to 10 sec)
298 my $z = '0' x ($BASE_LEN-1);
0716bf9b 299 while ($l >= 0)
300 {
574bacfe 301 $ret .= substr($z.$ar->[$l],-$BASE_LEN); # fastest way I could think of
0716bf9b 302 $l--;
303 }
61f5c3f5 304 \$ret;
0716bf9b 305 }
306
307sub _num
308 {
309 # Make a number (scalar int/float) from a BigInt object
574bacfe 310 my $x = $_[1];
0716bf9b 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 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 330 # This routine clobbers up array x, but not y.
331
574bacfe 332 my ($c,$x,$y) = @_;
b3abae2a 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 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 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 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 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 370 push @$x,1 if ($x->[-1] == 0); # last overflowed, so extend
371 $x;
e745a66c 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 382 for my $i (@$x)
383 {
384 last if (($i -= 1) >= 0); # early out
61f5c3f5 385 $i = $MAX; # overflow, next
e745a66c 386 }
387 pop @$x if $x->[-1] == 0 && @$x > 1; # last overflowed (but leave 0)
61f5c3f5 388 $x;
0716bf9b 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 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 406 }
407 # might leave leading zeros, so fix that
394e6ffb 408 return __strip_zeros($sx);
0716bf9b 409 }
394e6ffb 410 #print "case 1 (swap)\n";
411 for $i (@$sx)
0716bf9b 412 {
394e6ffb 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 418 # might leave leading zeros, so fix that
419 __strip_zeros($sy);
0716bf9b 420 }
421
ee15d750 422sub _mul_use_mul
0716bf9b 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 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 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 456 for $xi (@$xv)
457 {
458 $car = 0; $cty = 0;
574bacfe 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 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 473 for $yi (@$yv)
474 {
475 $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
574bacfe 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 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 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 495 }
496
ee15d750 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 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 524 # since multiplying $x with $x fails, make copy in this case
525 $yv = [@$xv] if "$xv" eq "$yv"; # same references?
61f5c3f5 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 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 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 545 }
546 push @$xv, @prod;
61f5c3f5 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 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 565 if (@$x == 1 && @$yorg == 1)
566 {
567 # shortcut, $y is smaller than $x
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 }
0716bf9b 580
0716bf9b 581 my $y = [ @$yorg ];
61f5c3f5 582 if ($LEN_CONVERT != 0)
583 {
584 $c->_to_small($x); $c->_to_small($y);
585 }
586
587 my ($car,$bar,$prd,$dd,$xi,$yi,@q,$v2,$v1,@d,$tmp,$q,$u2,$u1,$u0);
588
589 $car = $bar = $prd = 0;
590 if (($dd = int($MBASE/($y->[-1]+1))) != 1)
0716bf9b 591 {
592 for $xi (@$x)
593 {
594 $xi = $xi * $dd + $car;
61f5c3f5 595 $xi -= ($car = int($xi * $RBASE)) * $MBASE; # see USE_MUL
0716bf9b 596 }
597 push(@$x, $car); $car = 0;
598 for $yi (@$y)
599 {
600 $yi = $yi * $dd + $car;
61f5c3f5 601 $yi -= ($car = int($yi * $RBASE)) * $MBASE; # see USE_MUL
0716bf9b 602 }
603 }
604 else
605 {
606 push(@$x, 0);
607 }
608 @q = (); ($v2,$v1) = @$y[-2,-1];
609 $v2 = 0 unless $v2;
610 while ($#$x > $#$y)
611 {
612 ($u2,$u1,$u0) = @$x[-3..-1];
613 $u2 = 0 unless $u2;
614 #warn "oups v1 is 0, u0: $u0 $y->[-2] $y->[-1] l ",scalar @$y,"\n"
615 # if $v1 == 0;
61f5c3f5 616 $q = (($u0 == $v1) ? $MAX_VAL : int(($u0*$MBASE+$u1)/$v1));
617 --$q while ($v2*$q > ($u0*$MBASE+$u1-$q*$v1)*$MBASE+$u2);
0716bf9b 618 if ($q)
619 {
620 ($car, $bar) = (0,0);
621 for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
622 {
623 $prd = $q * $y->[$yi] + $car;
61f5c3f5 624 $prd -= ($car = int($prd * $RBASE)) * $MBASE; # see USE_MUL
625 $x->[$xi] += $MBASE if ($bar = (($x->[$xi] -= $prd + $bar) < 0));
0716bf9b 626 }
627 if ($x->[-1] < $car + $bar)
628 {
629 $car = 0; --$q;
630 for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
631 {
61f5c3f5 632 $x->[$xi] -= $MBASE
633 if ($car = (($x->[$xi] += $y->[$yi] + $car) > $MBASE));
0716bf9b 634 }
635 }
636 }
637 pop(@$x); unshift(@q, $q);
638 }
639 if (wantarray)
640 {
641 @d = ();
642 if ($dd != 1)
643 {
644 $car = 0;
645 for $xi (reverse @$x)
646 {
61f5c3f5 647 $prd = $car * $MBASE + $xi;
0716bf9b 648 $car = $prd - ($tmp = int($prd / $dd)) * $dd; # see USE_MUL
649 unshift(@d, $tmp);
650 }
651 }
652 else
653 {
654 @d = @$x;
655 }
656 @$x = @q;
61f5c3f5 657 my $d = \@d;
658 if ($LEN_CONVERT != 0)
659 {
660 $c->_to_large($x); $c->_to_large($d);
661 }
662 else
663 {
664 __strip_zeros($x);
665 __strip_zeros($d);
666 }
667 return ($x,$d);
0716bf9b 668 }
669 @$x = @q;
61f5c3f5 670 if ($LEN_CONVERT != 0)
671 {
672 $c->_to_large($x);
673 }
674 else
675 {
676 __strip_zeros($x);
677 }
678 $x;
0716bf9b 679 }
680
ee15d750 681sub _div_use_div
682 {
683 # ref to array, ref to array, modify first array and return remainder if
684 # in list context
ee15d750 685 my ($c,$x,$yorg) = @_;
ee15d750 686
61f5c3f5 687 if (@$x == 1 && @$yorg == 1)
688 {
689 # shortcut, $y is smaller than $x
690 if (wantarray)
691 {
692 my $r = [ $x->[0] % $yorg->[0] ];
693 $x->[0] = int($x->[0] / $yorg->[0]);
694 return ($x,$r);
695 }
696 else
697 {
698 $x->[0] = int($x->[0] / $yorg->[0]);
699 return $x;
700 }
701 }
ee15d750 702
ee15d750 703 my $y = [ @$yorg ];
61f5c3f5 704 if ($LEN_CONVERT != 0)
705 {
706 $c->_to_small($x); $c->_to_small($y);
707 }
708
709 my ($car,$bar,$prd,$dd,$xi,$yi,@q,$v2,$v1,@d,$tmp,$q,$u2,$u1,$u0);
710
711 $car = $bar = $prd = 0;
712 if (($dd = int($MBASE/($y->[-1]+1))) != 1)
ee15d750 713 {
714 for $xi (@$x)
715 {
716 $xi = $xi * $dd + $car;
61f5c3f5 717 $xi -= ($car = int($xi / $MBASE)) * $MBASE;
ee15d750 718 }
719 push(@$x, $car); $car = 0;
720 for $yi (@$y)
721 {
722 $yi = $yi * $dd + $car;
61f5c3f5 723 $yi -= ($car = int($yi / $MBASE)) * $MBASE;
ee15d750 724 }
725 }
726 else
727 {
728 push(@$x, 0);
729 }
730 @q = (); ($v2,$v1) = @$y[-2,-1];
731 $v2 = 0 unless $v2;
732 while ($#$x > $#$y)
733 {
734 ($u2,$u1,$u0) = @$x[-3..-1];
735 $u2 = 0 unless $u2;
736 #warn "oups v1 is 0, u0: $u0 $y->[-2] $y->[-1] l ",scalar @$y,"\n"
737 # if $v1 == 0;
61f5c3f5 738 $q = (($u0 == $v1) ? $MAX_VAL : int(($u0*$MBASE+$u1)/$v1));
739 --$q while ($v2*$q > ($u0*$MBASE+$u1-$q*$v1)*$MBASE+$u2);
ee15d750 740 if ($q)
741 {
742 ($car, $bar) = (0,0);
743 for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
744 {
745 $prd = $q * $y->[$yi] + $car;
61f5c3f5 746 $prd -= ($car = int($prd / $MBASE)) * $MBASE;
747 $x->[$xi] += $MBASE if ($bar = (($x->[$xi] -= $prd + $bar) < 0));
ee15d750 748 }
749 if ($x->[-1] < $car + $bar)
750 {
751 $car = 0; --$q;
752 for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
753 {
61f5c3f5 754 $x->[$xi] -= $MBASE
755 if ($car = (($x->[$xi] += $y->[$yi] + $car) > $MBASE));
ee15d750 756 }
757 }
758 }
61f5c3f5 759 pop(@$x); unshift(@q, $q);
ee15d750 760 }
761 if (wantarray)
762 {
763 @d = ();
764 if ($dd != 1)
765 {
766 $car = 0;
767 for $xi (reverse @$x)
768 {
61f5c3f5 769 $prd = $car * $MBASE + $xi;
ee15d750 770 $car = $prd - ($tmp = int($prd / $dd)) * $dd;
771 unshift(@d, $tmp);
772 }
773 }
774 else
775 {
776 @d = @$x;
777 }
778 @$x = @q;
61f5c3f5 779 my $d = \@d;
780 if ($LEN_CONVERT != 0)
781 {
782 $c->_to_large($x); $c->_to_large($d);
783 }
784 else
785 {
786 __strip_zeros($x);
787 __strip_zeros($d);
788 }
789 return ($x,$d);
ee15d750 790 }
791 @$x = @q;
61f5c3f5 792 if ($LEN_CONVERT != 0)
793 {
794 $c->_to_large($x);
795 }
796 else
797 {
798 __strip_zeros($x);
799 }
800 $x;
ee15d750 801 }
802
394e6ffb 803##############################################################################
804# testing
805
806sub _acmp
807 {
808 # internal absolute post-normalized compare (ignore signs)
809 # ref to array, ref to array, return <0, 0, >0
810 # arrays must have at least one entry; this is not checked for
811
812 my ($c,$cx,$cy) = @_;
813
61f5c3f5 814 # fast comp based on array elements
394e6ffb 815 my $lxy = scalar @$cx - scalar @$cy;
816 return -1 if $lxy < 0; # already differs, ret
817 return 1 if $lxy > 0; # ditto
818
819 # now calculate length based on digits, not parts
820 $lxy = _len($c,$cx) - _len($c,$cy); # difference
821 return -1 if $lxy < 0;
822 return 1 if $lxy > 0;
823
824 # hm, same lengths, but same contents?
825 my $i = 0; my $a;
826 # first way takes 5.49 sec instead of 4.87, but has the early out advantage
827 # so grep is slightly faster, but more inflexible. hm. $_ instead of $k
828 # yields 5.6 instead of 5.5 sec huh?
829 # manual way (abort if unequal, good for early ne)
830 my $j = scalar @$cx - 1;
831 while ($j >= 0)
832 {
833 last if ($a = $cx->[$j] - $cy->[$j]); $j--;
834 }
835 return 1 if $a > 0;
836 return -1 if $a < 0;
61f5c3f5 837 0; # equal
838
394e6ffb 839 # while it early aborts, it is even slower than the manual variant
840 #grep { return $a if ($a = $_ - $cy->[$i++]); } @$cx;
841 # grep way, go trough all (bad for early ne)
842 #grep { $a = $_ - $cy->[$i++]; } @$cx;
843 #return $a;
844 }
845
846sub _len
847 {
848 # compute number of digits in bigint, minus the sign
849
850 # int() because add/sub sometimes leaves strings (like '00005') instead of
851 # '5' in this place, thus causing length() to report wrong length
852 my $cx = $_[1];
853
854 return (@$cx-1)*$BASE_LEN+length(int($cx->[-1]));
855 }
856
857sub _digit
858 {
859 # return the nth digit, negative values count backward
860 # zero is rightmost, so _digit(123,0) will give 3
861 my ($c,$x,$n) = @_;
862
863 my $len = _len('',$x);
864
865 $n = $len+$n if $n < 0; # -1 last, -2 second-to-last
866 $n = abs($n); # if negative was too big
867 $len--; $n = $len if $n > $len; # n to big?
868
869 my $elem = int($n / $BASE_LEN); # which array element
870 my $digit = $n % $BASE_LEN; # which digit in this element
871 $elem = '0000'.@$x[$elem]; # get element padded with 0's
872 return substr($elem,-$digit-1,1);
873 }
874
875sub _zeros
876 {
877 # return amount of trailing zeros in decimal
878 # check each array elem in _m for having 0 at end as long as elem == 0
879 # Upon finding a elem != 0, stop
880 my $x = $_[1];
881 my $zeros = 0; my $elem;
882 foreach my $e (@$x)
883 {
884 if ($e != 0)
885 {
886 $elem = "$e"; # preserve x
887 $elem =~ s/.*?(0*$)/$1/; # strip anything not zero
888 $zeros *= $BASE_LEN; # elems * 5
61f5c3f5 889 $zeros += length($elem); # count trailing zeros
394e6ffb 890 last; # early out
891 }
892 $zeros ++; # real else branch: 50% slower!
893 }
61f5c3f5 894 $zeros;
394e6ffb 895 }
896
897##############################################################################
898# _is_* routines
899
900sub _is_zero
901 {
902 # return true if arg (BINT or num_str) is zero (array '+', '0')
903 my $x = $_[1];
61f5c3f5 904
905 (((scalar @$x == 1) && ($x->[0] == 0))) <=> 0;
394e6ffb 906 }
907
908sub _is_even
909 {
910 # return true if arg (BINT or num_str) is even
911 my $x = $_[1];
61f5c3f5 912 (!($x->[0] & 1)) <=> 0;
394e6ffb 913 }
914
915sub _is_odd
916 {
917 # return true if arg (BINT or num_str) is even
918 my $x = $_[1];
61f5c3f5 919
920 (($x->[0] & 1)) <=> 0;
394e6ffb 921 }
922
923sub _is_one
924 {
925 # return true if arg (BINT or num_str) is one (array '+', '1')
926 my $x = $_[1];
61f5c3f5 927
928 (scalar @$x == 1) && ($x->[0] == 1) <=> 0;
394e6ffb 929 }
930
931sub __strip_zeros
932 {
933 # internal normalization function that strips leading zeros from the array
934 # args: ref to array
935 my $s = shift;
936
937 my $cnt = scalar @$s; # get count of parts
938 my $i = $cnt-1;
939 push @$s,0 if $i < 0; # div might return empty results, so fix it
940
61f5c3f5 941 return $s if @$s == 1; # early out
942
394e6ffb 943 #print "strip: cnt $cnt i $i\n";
944 # '0', '3', '4', '0', '0',
945 # 0 1 2 3 4
946 # cnt = 5, i = 4
947 # i = 4
948 # i = 3
949 # => fcnt = cnt - i (5-2 => 3, cnt => 5-1 = 4, throw away from 4th pos)
950 # >= 1: skip first part (this can be zero)
951 while ($i > 0) { last if $s->[$i] != 0; $i--; }
952 $i++; splice @$s,$i if ($i < $cnt); # $i cant be 0
953 $s;
954 }
955
956###############################################################################
957# check routine to test internal state of corruptions
958
959sub _check
960 {
961 # used by the test suite
962 my $x = $_[1];
963
964 return "$x is not a reference" if !ref($x);
965
966 # are all parts are valid?
967 my $i = 0; my $j = scalar @$x; my ($e,$try);
968 while ($i < $j)
969 {
970 $e = $x->[$i]; $e = 'undef' unless defined $e;
971 $try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e)";
972 last if $e !~ /^[+]?[0-9]+$/;
973 $try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e) (stringify)";
974 last if "$e" !~ /^[+]?[0-9]+$/;
975 $try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e) (cat-stringify)";
976 last if '' . "$e" !~ /^[+]?[0-9]+$/;
977 $try = ' < 0 || >= $BASE; '."($x, $e)";
978 last if $e <0 || $e >= $BASE;
979 # this test is disabled, since new/bnorm and certain ops (like early out
980 # in add/sub) are allowed/expected to leave '00000' in some elements
981 #$try = '=~ /^00+/; '."($x, $e)";
982 #last if $e =~ /^00+/;
983 $i++;
984 }
985 return "Illegal part '$e' at pos $i (tested: $try)" if $i < $j;
986 return 0;
987 }
988
989
990###############################################################################
991###############################################################################
992# some optional routines to make BigInt faster
993
dccbb853 994sub _mod
995 {
996 # if possible, use mod shortcut
997 my ($c,$x,$yo) = @_;
998
999 # slow way since $y to big
1000 if (scalar @$yo > 1)
1001 {
1002 my ($xo,$rem) = _div($c,$x,$yo);
1003 return $rem;
1004 }
1005 my $y = $yo->[0];
027dc388 1006 # both are single element arrays
dccbb853 1007 if (scalar @$x == 1)
1008 {
1009 $x->[0] %= $y;
1010 return $x;
1011 }
1012
61f5c3f5 1013 # @y is single element, but @x has more than one
dccbb853 1014 my $b = $BASE % $y;
1015 if ($b == 0)
1016 {
1017 # when BASE % Y == 0 then (B * BASE) % Y == 0
1018 # (B * BASE) % $y + A % Y => A % Y
1019 # so need to consider only last element: O(1)
1020 $x->[0] %= $y;
1021 }
027dc388 1022 elsif ($b == 1)
1023 {
1024 # else need to go trough all elements: O(N), but loop is a bit simplified
1025 my $r = 0;
1026 foreach (@$x)
1027 {
1028 $r += $_ % $y;
1029 $r %= $y;
1030 }
1031 $r = 0 if $r == $y;
1032 $x->[0] = $r;
1033 }
dccbb853 1034 else
1035 {
027dc388 1036 # else need to go trough all elements: O(N)
1037 my $r = 0; my $bm = 1;
1038 foreach (@$x)
1039 {
1040 $r += ($_ % $y) * $bm;
1041 $bm *= $b;
1042 $bm %= $y;
1043 $r %= $y;
1044 }
1045 $r = 0 if $r == $y;
1046 $x->[0] = $r;
dccbb853 1047 }
1048 splice (@$x,1);
61f5c3f5 1049 $x;
dccbb853 1050 }
1051
0716bf9b 1052##############################################################################
574bacfe 1053# shifts
1054
1055sub _rsft
1056 {
1057 my ($c,$x,$y,$n) = @_;
1058
1059 if ($n != 10)
1060 {
61f5c3f5 1061 $n = _new($c,\$n); return _div($c,$x, _pow($c,$n,$y));
1062 }
1063
1064 # shortcut (faster) for shifting by 10)
1065 # multiples of $BASE_LEN
1066 my $dst = 0; # destination
1067 my $src = _num($c,$y); # as normal int
1068 my $rem = $src % $BASE_LEN; # remainder to shift
1069 $src = int($src / $BASE_LEN); # source
1070 if ($rem == 0)
1071 {
1072 splice (@$x,0,$src); # even faster, 38.4 => 39.3
574bacfe 1073 }
1074 else
1075 {
61f5c3f5 1076 my $len = scalar @$x - $src; # elems to go
1077 my $vd; my $z = '0'x $BASE_LEN;
1078 $x->[scalar @$x] = 0; # avoid || 0 test inside loop
1079 while ($dst < $len)
574bacfe 1080 {
61f5c3f5 1081 $vd = $z.$x->[$src];
1082 $vd = substr($vd,-$BASE_LEN,$BASE_LEN-$rem);
1083 $src++;
1084 $vd = substr($z.$x->[$src],-$rem,$rem) . $vd;
1085 $vd = substr($vd,-$BASE_LEN,$BASE_LEN) if length($vd) > $BASE_LEN;
1086 $x->[$dst] = int($vd);
1087 $dst++;
574bacfe 1088 }
61f5c3f5 1089 splice (@$x,$dst) if $dst > 0; # kill left-over array elems
1090 pop @$x if $x->[-1] == 0; # kill last element if 0
1091 } # else rem == 0
574bacfe 1092 $x;
1093 }
1094
1095sub _lsft
1096 {
1097 my ($c,$x,$y,$n) = @_;
1098
1099 if ($n != 10)
1100 {
61f5c3f5 1101 $n = _new($c,\$n); return _mul($c,$x, _pow($c,$n,$y));
574bacfe 1102 }
61f5c3f5 1103
1104 # shortcut (faster) for shifting by 10) since we are in base 10eX
1105 # multiples of $BASE_LEN:
1106 my $src = scalar @$x; # source
1107 my $len = _num($c,$y); # shift-len as normal int
1108 my $rem = $len % $BASE_LEN; # remainder to shift
1109 my $dst = $src + int($len/$BASE_LEN); # destination
1110 my $vd; # further speedup
1111 $x->[$src] = 0; # avoid first ||0 for speed
1112 my $z = '0' x $BASE_LEN;
1113 while ($src >= 0)
574bacfe 1114 {
61f5c3f5 1115 $vd = $x->[$src]; $vd = $z.$vd;
1116 $vd = substr($vd,-$BASE_LEN+$rem,$BASE_LEN-$rem);
1117 $vd .= $src > 0 ? substr($z.$x->[$src-1],-$BASE_LEN,$rem) : '0' x $rem;
1118 $vd = substr($vd,-$BASE_LEN,$BASE_LEN) if length($vd) > $BASE_LEN;
1119 $x->[$dst] = int($vd);
1120 $dst--; $src--;
574bacfe 1121 }
61f5c3f5 1122 # set lowest parts to 0
1123 while ($dst >= 0) { $x->[$dst--] = 0; }
1124 # fix spurios last zero element
1125 splice @$x,-1 if $x->[-1] == 0;
574bacfe 1126 $x;
1127 }
1128
027dc388 1129sub _pow
1130 {
1131 # power of $x to $y
1132 # ref to array, ref to array, return ref to array
1133 my ($c,$cx,$cy) = @_;
1134
1135 my $pow2 = _one();
1136 my $two = _two();
1137 my $y1 = _copy($c,$cy);
1138 while (!_is_one($c,$y1))
1139 {
1140 _mul($c,$pow2,$cx) if _is_odd($c,$y1);
1141 _div($c,$y1,$two);
1142 _mul($c,$cx,$cx);
1143 }
1144 _mul($c,$cx,$pow2) unless _is_one($c,$pow2);
61f5c3f5 1145 $cx;
027dc388 1146 }
1147
b3abae2a 1148sub _fac
1149 {
1150 # factorial of $x
1151 # ref to array, return ref to array
1152 my ($c,$cx) = @_;
1153
1154 if ((@$cx == 1) && ($cx->[0] <= 2))
1155 {
1156 $cx->[0] = 1 * ($cx->[0]||1); # 0,1 => 1, 2 => 2
1157 return $cx;
1158 }
1159
1160 # go forward until $base is exceeded
1161 # limit is either $x or $base (x == 100 means as result too high)
1162 my $steps = 100; $steps = $cx->[0] if @$cx == 1;
1163 my $r = 2; my $cf = 3; my $step = 1; my $last = $r;
1164 while ($r < $BASE && $step < $steps)
1165 {
1166 $last = $r; $r *= $cf++; $step++;
1167 }
1168 if ((@$cx == 1) && ($step == $cx->[0]))
1169 {
1170 # completely done
1171 $cx = [$last];
1172 return $cx;
1173 }
1174 my $n = _copy($c,$cx);
1175 $cx = [$last];
1176
1177 #$cx = _one();
1178 while (!(@$n == 1 && $n->[0] == $step))
1179 {
1180 _mul($c,$cx,$n); _dec($c,$n);
1181 }
1182 $cx;
1183 }
1184
1185use constant DEBUG => 0;
1186
1187my $steps = 0;
1188
1189sub steps { $steps };
1190
1191sub _sqrt
0716bf9b 1192 {
394e6ffb 1193 # square-root of $x
1194 # ref to array, return ref to array
1195 my ($c,$x) = @_;
0716bf9b 1196
394e6ffb 1197 if (scalar @$x == 1)
1198 {
1199 # fit's into one Perl scalar
1200 $x->[0] = int(sqrt($x->[0]));
1201 return $x;
1202 }
1203 my $y = _copy($c,$x);
b3abae2a 1204 # hopefully _len/2 is < $BASE, the -1 is to always undershot the guess
1205 # since our guess will "grow"
1206 my $l = int((_len($c,$x)-1) / 2);
1207
1208 my $lastelem = $x->[-1]; # for guess
1209 my $elems = scalar @$x - 1;
1210 # not enough digits, but could have more?
1211 if ((length($lastelem) <= 3) && ($elems > 1))
1212 {
1213 # right-align with zero pad
1214 my $len = length($lastelem) & 1;
1215 print "$lastelem => " if DEBUG;
1216 $lastelem .= substr($x->[-2] . '0' x $BASE_LEN,0,$BASE_LEN);
1217 # former odd => make odd again, or former even to even again
1218 $lastelem = $lastelem / 10 if (length($lastelem) & 1) != $len;
1219 print "$lastelem\n" if DEBUG;
1220 }
0716bf9b 1221
61f5c3f5 1222 # construct $x (instead of _lsft($c,$x,$l,10)
1223 my $r = $l % $BASE_LEN; # 10000 00000 00000 00000 ($BASE_LEN=5)
1224 $l = int($l / $BASE_LEN);
b3abae2a 1225 print "l = $l " if DEBUG;
1226
1227 splice @$x,$l; # keep ref($x), but modify it
1228
1229 # we make the first part of the guess not '1000...0' but int(sqrt($lastelem))
1230 # that gives us:
1231 # 14400 00000 => sqrt(14400) => 120
1232 # 144000 000000 => sqrt(144000) => 379
1233
1234 # $x->[$l--] = int('1' . '0' x $r); # old way of guessing
1235 print "$lastelem (elems $elems) => " if DEBUG;
1236 $lastelem = $lastelem / 10 if ($elems & 1 == 1); # odd or even?
1237 my $g = sqrt($lastelem); $g =~ s/\.//; # 2.345 => 2345
1238 $r -= 1 if $elems & 1 == 0; # 70 => 7
1239
1240 # padd with zeros if result is too short
1241 $x->[$l--] = int(substr($g . '0' x $r,0,$r+1));
1242 print "now ",$x->[-1] if DEBUG;
1243 print " would have been ", int('1' . '0' x $r),"\n" if DEBUG;
1244
1245 # If @$x > 1, we could compute the second elem of the guess, too, to create
1246 # an even better guess. Not implemented yet.
1247 $x->[$l--] = 0 while ($l >= 0); # all other digits of guess are zero
61f5c3f5 1248
b3abae2a 1249 print "start x= ",${_str($c,$x)},"\n" if DEBUG;
394e6ffb 1250 my $two = _two();
1251 my $last = _zero();
1252 my $lastlast = _zero();
b3abae2a 1253 $steps = 0 if DEBUG;
394e6ffb 1254 while (_acmp($c,$last,$x) != 0 && _acmp($c,$lastlast,$x) != 0)
1255 {
b3abae2a 1256 $steps++ if DEBUG;
394e6ffb 1257 $lastlast = _copy($c,$last);
1258 $last = _copy($c,$x);
1259 _add($c,$x, _div($c,_copy($c,$y),$x));
1260 _div($c,$x, $two );
b3abae2a 1261 print " x= ",${_str($c,$x)},"\n" if DEBUG;
394e6ffb 1262 }
b3abae2a 1263 print "\nsteps in sqrt: $steps, " if DEBUG;
394e6ffb 1264 _dec($c,$x) if _acmp($c,$y,_mul($c,_copy($c,$x),$x)) < 0; # overshot?
b3abae2a 1265 print " final ",$x->[-1],"\n" if DEBUG;
394e6ffb 1266 $x;
0716bf9b 1267 }
1268
394e6ffb 1269##############################################################################
1270# binary stuff
0716bf9b 1271
394e6ffb 1272sub _and
1273 {
1274 my ($c,$x,$y) = @_;
0716bf9b 1275
394e6ffb 1276 # the shortcut makes equal, large numbers _really_ fast, and makes only a
1277 # very small performance drop for small numbers (e.g. something with less
1278 # than 32 bit) Since we optimize for large numbers, this is enabled.
1279 return $x if _acmp($c,$x,$y) == 0; # shortcut
0716bf9b 1280
394e6ffb 1281 my $m = _one(); my ($xr,$yr);
1282 my $mask = $AND_MASK;
1283
1284 my $x1 = $x;
1285 my $y1 = _copy($c,$y); # make copy
1286 $x = _zero();
1287 my ($b,$xrr,$yrr);
1288 use integer;
1289 while (!_is_zero($c,$x1) && !_is_zero($c,$y1))
1290 {
1291 ($x1, $xr) = _div($c,$x1,$mask);
1292 ($y1, $yr) = _div($c,$y1,$mask);
1293
1294 # make ints() from $xr, $yr
1295 # this is when the AND_BITS are greater tahn $BASE and is slower for
1296 # small (<256 bits) numbers, but faster for large numbers. Disabled
1297 # due to KISS principle
1298
1299# $b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
1300# $b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
1301# _add($c,$x, _mul($c, _new( $c, \($xrr & $yrr) ), $m) );
1302
61f5c3f5 1303 # 0+ due to '&' doesn't work in strings
1304 _add($c,$x, _mul($c, [ 0+$xr->[0] & 0+$yr->[0] ], $m) );
394e6ffb 1305 _mul($c,$m,$mask);
1306 }
1307 $x;
0716bf9b 1308 }
1309
394e6ffb 1310sub _xor
0716bf9b 1311 {
394e6ffb 1312 my ($c,$x,$y) = @_;
1313
1314 return _zero() if _acmp($c,$x,$y) == 0; # shortcut (see -and)
1315
1316 my $m = _one(); my ($xr,$yr);
1317 my $mask = $XOR_MASK;
1318
1319 my $x1 = $x;
1320 my $y1 = _copy($c,$y); # make copy
1321 $x = _zero();
1322 my ($b,$xrr,$yrr);
1323 use integer;
1324 while (!_is_zero($c,$x1) && !_is_zero($c,$y1))
0716bf9b 1325 {
394e6ffb 1326 ($x1, $xr) = _div($c,$x1,$mask);
1327 ($y1, $yr) = _div($c,$y1,$mask);
1328 # make ints() from $xr, $yr (see _and())
1329 #$b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
1330 #$b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
1331 #_add($c,$x, _mul($c, _new( $c, \($xrr ^ $yrr) ), $m) );
61f5c3f5 1332
1333 # 0+ due to '^' doesn't work in strings
1334 _add($c,$x, _mul($c, [ 0+$xr->[0] ^ 0+$yr->[0] ], $m) );
394e6ffb 1335 _mul($c,$m,$mask);
0716bf9b 1336 }
394e6ffb 1337 # the loop stops when the shorter of the two numbers is exhausted
1338 # the remainder of the longer one will survive bit-by-bit, so we simple
1339 # multiply-add it in
1340 _add($c,$x, _mul($c, $x1, $m) ) if !_is_zero($c,$x1);
1341 _add($c,$x, _mul($c, $y1, $m) ) if !_is_zero($c,$y1);
1342
1343 $x;
0716bf9b 1344 }
1345
394e6ffb 1346sub _or
0716bf9b 1347 {
394e6ffb 1348 my ($c,$x,$y) = @_;
0716bf9b 1349
394e6ffb 1350 return $x if _acmp($c,$x,$y) == 0; # shortcut (see _and)
0716bf9b 1351
394e6ffb 1352 my $m = _one(); my ($xr,$yr);
1353 my $mask = $OR_MASK;
0716bf9b 1354
394e6ffb 1355 my $x1 = $x;
1356 my $y1 = _copy($c,$y); # make copy
1357 $x = _zero();
1358 my ($b,$xrr,$yrr);
1359 use integer;
1360 while (!_is_zero($c,$x1) && !_is_zero($c,$y1))
1361 {
1362 ($x1, $xr) = _div($c,$x1,$mask);
1363 ($y1, $yr) = _div($c,$y1,$mask);
1364 # make ints() from $xr, $yr (see _and())
1365# $b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
1366# $b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
1367# _add($c,$x, _mul($c, _new( $c, \($xrr | $yrr) ), $m) );
1368
61f5c3f5 1369 # 0+ due to '|' doesn't work in strings
1370 _add($c,$x, _mul($c, [ 0+$xr->[0] | 0+$yr->[0] ], $m) );
394e6ffb 1371 _mul($c,$m,$mask);
1372 }
1373 # the loop stops when the shorter of the two numbers is exhausted
1374 # the remainder of the longer one will survive bit-by-bit, so we simple
1375 # multiply-add it in
1376 _add($c,$x, _mul($c, $x1, $m) ) if !_is_zero($c,$x1);
1377 _add($c,$x, _mul($c, $y1, $m) ) if !_is_zero($c,$y1);
1378
1379 $x;
0716bf9b 1380 }
1381
61f5c3f5 1382sub _as_hex
1383 {
1384 # convert a decimal number to hex (ref to array, return ref to string)
1385 my ($c,$x) = @_;
1386
1387 my $x1 = _copy($c,$x);
1388
1389 my $es = '';
1390 my $xr;
1391 my $x10000 = [ 0x10000 ];
1392 while (! _is_zero($c,$x1))
1393 {
1394 ($x1, $xr) = _div($c,$x1,$x10000);
1395 $es .= unpack('h4',pack('v',$xr->[0]));
1396 }
1397 $es = reverse $es;
1398 $es =~ s/^[0]+//; # strip leading zeros
1399 $es = '0x' . $es;
1400 \$es;
1401 }
1402
1403sub _as_bin
1404 {
1405 # convert a decimal number to bin (ref to array, return ref to string)
1406 my ($c,$x) = @_;
1407
1408 my $x1 = _copy($c,$x);
1409
1410 my $es = '';
1411 my $xr;
1412 my $x10000 = [ 0x10000 ];
1413 while (! _is_zero($c,$x1))
1414 {
1415 ($x1, $xr) = _div($c,$x1,$x10000);
1416 $es .= unpack('b16',pack('v',$xr->[0]));
1417 }
1418 $es = reverse $es;
1419 $es =~ s/^[0]+//; # strip leading zeros
1420 $es = '0b' . $es;
1421 \$es;
1422 }
1423
394e6ffb 1424sub _from_hex
0716bf9b 1425 {
394e6ffb 1426 # convert a hex number to decimal (ref to string, return ref to array)
1427 my ($c,$hs) = @_;
0716bf9b 1428
394e6ffb 1429 my $mul = _one();
1430 my $m = [ 0x10000 ]; # 16 bit at a time
1431 my $x = _zero();
0716bf9b 1432
61f5c3f5 1433 my $len = length($$hs)-2;
394e6ffb 1434 $len = int($len/4); # 4-digit parts, w/o '0x'
1435 my $val; my $i = -4;
1436 while ($len >= 0)
1437 {
1438 $val = substr($$hs,$i,4);
1439 $val =~ s/^[+-]?0x// if $len == 0; # for last part only because
1440 $val = hex($val); # hex does not like wrong chars
1441 $i -= 4; $len --;
1442 _add ($c, $x, _mul ($c, [ $val ], $mul ) ) if $val != 0;
1443 _mul ($c, $mul, $m ) if $len >= 0; # skip last mul
1444 }
1445 $x;
1446 }
1447
1448sub _from_bin
0716bf9b 1449 {
394e6ffb 1450 # convert a hex number to decimal (ref to string, return ref to array)
1451 my ($c,$bs) = @_;
0716bf9b 1452
394e6ffb 1453 my $mul = _one();
1454 my $m = [ 0x100 ]; # 8 bit at a time
1455 my $x = _zero();
0716bf9b 1456
61f5c3f5 1457 my $len = length($$bs)-2;
394e6ffb 1458 $len = int($len/8); # 4-digit parts, w/o '0x'
1459 my $val; my $i = -8;
1460 while ($len >= 0)
0716bf9b 1461 {
394e6ffb 1462 $val = substr($$bs,$i,8);
1463 $val =~ s/^[+-]?0b// if $len == 0; # for last part only
1464
394e6ffb 1465 $val = ord(pack('B8',substr('00000000'.$val,-8,8)));
1466
1467 $i -= 8; $len --;
1468 _add ($c, $x, _mul ($c, [ $val ], $mul ) ) if $val != 0;
1469 _mul ($c, $mul, $m ) if $len >= 0; # skip last mul
0716bf9b 1470 }
394e6ffb 1471 $x;
0716bf9b 1472 }
1473
394e6ffb 1474##############################################################################
1475##############################################################################
1476
0716bf9b 14771;
1478__END__
1479
1480=head1 NAME
1481
1482Math::BigInt::Calc - Pure Perl module to support Math::BigInt
1483
1484=head1 SYNOPSIS
1485
ee15d750 1486Provides support for big integer calculations. Not intended to be used by other
1487modules (except Math::BigInt::Cached). Other modules which sport the same
1488functions can also be used to support Math::Bigint, like Math::BigInt::Pari.
0716bf9b 1489
1490=head1 DESCRIPTION
1491
027dc388 1492In order to allow for multiple big integer libraries, Math::BigInt was
1493rewritten to use library modules for core math routines. Any module which
1494follows the same API as this can be used instead by using the following:
0716bf9b 1495
ee15d750 1496 use Math::BigInt lib => 'libname';
0716bf9b 1497
027dc388 1498'libname' is either the long name ('Math::BigInt::Pari'), or only the short
1499version like 'Pari'.
1500
0716bf9b 1501=head1 EXPORT
1502
027dc388 1503The following functions MUST be defined in order to support the use by
1504Math::BigInt:
0716bf9b 1505
1506 _new(string) return ref to new object from ref to decimal string
1507 _zero() return a new object with value 0
1508 _one() return a new object with value 1
1509
1510 _str(obj) return ref to a string representing the object
1511 _num(obj) returns a Perl integer/floating point number
1512 NOTE: because of Perl numeric notation defaults,
1513 the _num'ified obj may lose accuracy due to
1514 machine-dependend floating point size limitations
1515
1516 _add(obj,obj) Simple addition of two objects
1517 _mul(obj,obj) Multiplication of two objects
1518 _div(obj,obj) Division of the 1st object by the 2nd
b22b3e31 1519 In list context, returns (result,remainder).
1520 NOTE: this is integer math, so no
1521 fractional part will be returned.
1522 _sub(obj,obj) Simple subtraction of 1 object from another
0716bf9b 1523 a third, optional parameter indicates that the params
1524 are swapped. In this case, the first param needs to
1525 be preserved, while you can destroy the second.
1526 sub (x,y,1) => return x - y and keep x intact!
e745a66c 1527 _dec(obj) decrement object by one (input is garant. to be > 0)
1528 _inc(obj) increment object by one
1529
0716bf9b 1530
1531 _acmp(obj,obj) <=> operator for objects (return -1, 0 or 1)
1532
1533 _len(obj) returns count of the decimal digits of the object
1534 _digit(obj,n) returns the n'th decimal digit of object
1535
1536 _is_one(obj) return true if argument is +1
1537 _is_zero(obj) return true if argument is 0
1538 _is_even(obj) return true if argument is even (0,2,4,6..)
1539 _is_odd(obj) return true if argument is odd (1,3,5,7..)
1540
1541 _copy return a ref to a true copy of the object
1542
1543 _check(obj) check whether internal representation is still intact
1544 return 0 for ok, otherwise error message as string
1545
bd05a461 1546The following functions are optional, and can be defined if the underlying lib
027dc388 1547has a fast way to do them. If undefined, Math::BigInt will use pure Perl (hence
1548slow) fallback routines to emulate these:
0716bf9b 1549
1550 _from_hex(str) return ref to new object from ref to hexadecimal string
1551 _from_bin(str) return ref to new object from ref to binary string
1552
ee15d750 1553 _as_hex(str) return ref to scalar string containing the value as
1554 unsigned hex string, with the '0x' prepended.
1555 Leading zeros must be stripped.
1556 _as_bin(str) Like as_hex, only as binary string containing only
1557 zeros and ones. Leading zeros must be stripped and a
1558 '0b' must be prepended.
1559
0716bf9b 1560 _rsft(obj,N,B) shift object in base B by N 'digits' right
dccbb853 1561 For unsupported bases B, return undef to signal failure
0716bf9b 1562 _lsft(obj,N,B) shift object in base B by N 'digits' left
dccbb853 1563 For unsupported bases B, return undef to signal failure
0716bf9b 1564
1565 _xor(obj1,obj2) XOR (bit-wise) object 1 with object 2
dccbb853 1566 Note: XOR, AND and OR pad with zeros if size mismatches
0716bf9b 1567 _and(obj1,obj2) AND (bit-wise) object 1 with object 2
1568 _or(obj1,obj2) OR (bit-wise) object 1 with object 2
1569
dccbb853 1570 _mod(obj,obj) Return remainder of div of the 1st by the 2nd object
394e6ffb 1571 _sqrt(obj) return the square root of object (truncate to int)
b3abae2a 1572 _fac(obj) return factorial of object 1 (1*2*3*4..)
0716bf9b 1573 _pow(obj,obj) return object 1 to the power of object 2
1574 _gcd(obj,obj) return Greatest Common Divisor of two objects
1575
b22b3e31 1576 _zeros(obj) return number of trailing decimal zeros
0716bf9b 1577
b22b3e31 1578Input strings come in as unsigned but with prefix (i.e. as '123', '0xabc'
0716bf9b 1579or '0b1101').
1580
b22b3e31 1581Testing of input parameter validity is done by the caller, so you need not
574bacfe 1582worry about underflow (f.i. in C<_sub()>, C<_dec()>) nor about division by
1583zero or similar cases.
1584
1585The first parameter can be modified, that includes the possibility that you
1586return a reference to a completely different object instead. Although keeping
dccbb853 1587the reference and just changing it's contents is prefered over creating and
1588returning a different reference.
574bacfe 1589
1590Return values are always references to objects or strings. Exceptions are
1591C<_lsft()> and C<_rsft()>, which return undef if they can not shift the
027dc388 1592argument. This is used to delegate shifting of bases different than the one
1593you can support back to Math::BigInt, which will use some generic code to
1594calculate the result.
574bacfe 1595
1596=head1 WRAP YOUR OWN
1597
1598If you want to port your own favourite c-lib for big numbers to the
1599Math::BigInt interface, you can take any of the already existing modules as
1600a rough guideline. You should really wrap up the latest BigInt and BigFloat
bd05a461 1601testsuites with your module, and replace in them any of the following:
574bacfe 1602
1603 use Math::BigInt;
1604
bd05a461 1605by this:
574bacfe 1606
1607 use Math::BigInt lib => 'yourlib';
1608
1609This way you ensure that your library really works 100% within Math::BigInt.
0716bf9b 1610
1611=head1 LICENSE
1612
1613This program is free software; you may redistribute it and/or modify it under
1614the same terms as Perl itself.
1615
1616=head1 AUTHORS
1617
1618Original math code by Mark Biggar, rewritten by Tels L<http://bloodgate.com/>
1619in late 2000, 2001.
1620Seperated from BigInt and shaped API with the help of John Peacock.
1621
1622=head1 SEE ALSO
1623
ee15d750 1624L<Math::BigInt>, L<Math::BigFloat>, L<Math::BigInt::BitVect>,
1625L<Math::BigInt::GMP>, L<Math::BigInt::Cached> and L<Math::BigInt::Pari>.
0716bf9b 1626
1627=cut