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