Upgrade to Math::BigInt 1.46.
[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
e745a66c 11$VERSION = '0.14';
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
16# leading zero parts (except the first) and in base 100000
17
18# todo:
19# - fully remove funky $# stuff (maybe)
0716bf9b 20
21# USE_MUL: due to problems on certain os (os390, posix-bc) "* 1e-5" is used
ee15d750 22# instead of "/ 1e5" at some places, (marked with USE_MUL). Other platforms
23# BS2000, some Crays need USE_DIV instead.
bd05a461 24# The BEGIN block is used to determine which of the two variants gives the
25# correct result.
0716bf9b 26
27##############################################################################
28# global constants, flags and accessory
29
30# constants for easier life
31my $nan = 'NaN';
ee15d750 32my ($BASE,$RBASE,$BASE_LEN,$MAX_VAL);
33
34sub _base_len
35 {
dccbb853 36 # set/get the BASE_LEN and assorted other, connected values
37 # used only be the testsuite, set is used only by the BEGIN block below
ee15d750 38 my $b = shift;
39 if (defined $b)
40 {
dccbb853 41 $b = 8 if $b > 8; # cap, for VMS, OS/390 and other 64 bit
ee15d750 42 $BASE_LEN = $b;
43 $BASE = int("1e".$BASE_LEN);
44 $RBASE = abs('1e-'.$BASE_LEN); # see USE_MUL
45 $MAX_VAL = $BASE-1;
46 # print "BASE_LEN: $BASE_LEN MAX_VAL: $MAX_VAL\n";
47 # print "int: ",int($BASE * $RBASE),"\n";
48 if (int($BASE * $RBASE) == 0) # should be 1
49 {
50 # must USE_MUL
ee15d750 51 *{_mul} = \&_mul_use_mul;
52 *{_div} = \&_div_use_mul;
53 }
54 else
55 {
ee15d750 56 # can USE_DIV instead
57 *{_mul} = \&_mul_use_div;
58 *{_div} = \&_div_use_div;
59 }
60 }
dccbb853 61 $BASE_LEN;
ee15d750 62 }
574bacfe 63
64BEGIN
65 {
bd05a461 66 # from Daniel Pfeiffer: determine largest group of digits that is precisely
574bacfe 67 # multipliable with itself plus carry
dccbb853 68 # Test now changed to expect the proper pattern, not a result off by 1 or 2
69 my ($e, $num) = 3; # lowest value we will use is 3+1-1 = 3
bd05a461 70 do
71 {
72 $num = ('9' x ++$e) + 0;
574bacfe 73 $num *= $num + 1;
dccbb853 74 # print "$num $e\n";
75 } while ("$num" =~ /9{$e}0{$e}/); # must be a certain pattern
76 # last test failed, so retract one step:
ee15d750 77 _base_len($e-1);
574bacfe 78 }
79
0716bf9b 80##############################################################################
81# create objects from various representations
82
83sub _new
84 {
85 # (string) return ref to num_array
86 # Convert a number from string format to internal base 100000 format.
87 # Assumes normalized value as input.
574bacfe 88 my $d = $_[1];
0716bf9b 89 # print "_new $d $$d\n";
90 my $il = CORE::length($$d)-1;
91 # these leaves '00000' instead of int 0 and will be corrected after any op
574bacfe 92 return [ reverse(unpack("a" . ($il % $BASE_LEN+1)
93 . ("a$BASE_LEN" x ($il / $BASE_LEN)), $$d)) ];
0716bf9b 94 }
95
96sub _zero
97 {
98 # create a zero
99 return [ 0 ];
100 }
101
102sub _one
103 {
104 # create a one
105 return [ 1 ];
106 }
107
108sub _copy
109 {
574bacfe 110 return [ @{$_[1]} ];
0716bf9b 111 }
112
bd05a461 113# catch and throw away
114sub import { }
115
0716bf9b 116##############################################################################
117# convert back to string and number
118
119sub _str
120 {
121 # (ref to BINT) return num_str
122 # Convert number from internal base 100000 format to string format.
123 # internal format is always normalized (no leading zeros, "-0" => "+0")
574bacfe 124 my $ar = $_[1];
0716bf9b 125 my $ret = "";
126 my $l = scalar @$ar; # number of parts
127 return $nan if $l < 1; # should not happen
128 # handle first one different to strip leading zeros from it (there are no
129 # leading zero parts in internal representation)
130 $l --; $ret .= $ar->[$l]; $l--;
131 # Interestingly, the pre-padd method uses more time
574bacfe 132 # the old grep variant takes longer (14 to 10 sec)
133 my $z = '0' x ($BASE_LEN-1);
0716bf9b 134 while ($l >= 0)
135 {
574bacfe 136 $ret .= substr($z.$ar->[$l],-$BASE_LEN); # fastest way I could think of
0716bf9b 137 $l--;
138 }
139 return \$ret;
140 }
141
142sub _num
143 {
144 # Make a number (scalar int/float) from a BigInt object
574bacfe 145 my $x = $_[1];
0716bf9b 146 return $x->[0] if scalar @$x == 1; # below $BASE
147 my $fac = 1;
148 my $num = 0;
149 foreach (@$x)
150 {
151 $num += $fac*$_; $fac *= $BASE;
152 }
153 return $num;
154 }
155
156##############################################################################
157# actual math code
158
159sub _add
160 {
161 # (ref to int_num_array, ref to int_num_array)
574bacfe 162 # routine to add two base 1eX numbers
0716bf9b 163 # stolen from Knuth Vol 2 Algorithm A pg 231
b22b3e31 164 # there are separate routines to add and sub as per Knuth pg 233
0716bf9b 165 # This routine clobbers up array x, but not y.
166
574bacfe 167 my ($c,$x,$y) = @_;
0716bf9b 168
169 # for each in Y, add Y to X and carry. If after that, something is left in
170 # X, foreach in X add carry to X and then return X, carry
171 # Trades one "$j++" for having to shift arrays, $j could be made integer
b22b3e31 172 # but this would impose a limit to number-length of 2**32.
0716bf9b 173 my $i; my $car = 0; my $j = 0;
174 for $i (@$y)
175 {
e745a66c 176 $x->[$j] -= $BASE if $car = (($x->[$j] += $i + $car) >= $BASE) ? 1 : 0;
0716bf9b 177 $j++;
178 }
179 while ($car != 0)
180 {
181 $x->[$j] -= $BASE if $car = (($x->[$j] += $car) >= $BASE) ? 1 : 0; $j++;
182 }
e745a66c 183 return $x;
184 }
185
186sub _inc
187 {
188 # (ref to int_num_array, ref to int_num_array)
189 # routine to add 1 to a base 1eX numbers
190 # This routine clobbers up array x, but not y.
191 my ($c,$x) = @_;
192
193 for my $i (@$x)
194 {
195 return $x if (($i += 1) < $BASE); # early out
196 $i -= $BASE;
197 }
198 if ($x->[-1] == 0) # last overflowed
199 {
200 push @$x,1; # extend
201 }
202 return $x;
203 }
204
205sub _dec
206 {
207 # (ref to int_num_array, ref to int_num_array)
208 # routine to add 1 to a base 1eX numbers
209 # This routine clobbers up array x, but not y.
210 my ($c,$x) = @_;
211
212 for my $i (@$x)
213 {
214 last if (($i -= 1) >= 0); # early out
215 $i = $MAX_VAL;
216 }
217 pop @$x if $x->[-1] == 0 && @$x > 1; # last overflowed (but leave 0)
218 return $x;
0716bf9b 219 }
220
221sub _sub
222 {
223 # (ref to int_num_array, ref to int_num_array)
574bacfe 224 # subtract base 1eX numbers -- stolen from Knuth Vol 2 pg 232, $x > $y
b22b3e31 225 # subtract Y from X (X is always greater/equal!) by modifying x in place
574bacfe 226 my ($c,$sx,$sy,$s) = @_;
0716bf9b 227
228 my $car = 0; my $i; my $j = 0;
229 if (!$s)
230 {
231 #print "case 2\n";
232 for $i (@$sx)
233 {
234 last unless defined $sy->[$j] || $car;
235 #print "x: $i y: $sy->[$j] c: $car\n";
236 $i += $BASE if $car = (($i -= ($sy->[$j] || 0) + $car) < 0); $j++;
237 #print "x: $i y: $sy->[$j-1] c: $car\n";
238 }
239 # might leave leading zeros, so fix that
240 __strip_zeros($sx);
241 return $sx;
242 }
243 else
244 {
245 #print "case 1 (swap)\n";
246 for $i (@$sx)
247 {
248 last unless defined $sy->[$j] || $car;
249 #print "$sy->[$j] $i $car => $sx->[$j]\n";
250 $sy->[$j] += $BASE
251 if $car = (($sy->[$j] = $i-($sy->[$j]||0) - $car) < 0);
252 #print "$sy->[$j] $i $car => $sy->[$j]\n";
253 $j++;
254 }
255 # might leave leading zeros, so fix that
256 __strip_zeros($sy);
257 return $sy;
258 }
259 }
260
ee15d750 261sub _mul_use_mul
0716bf9b 262 {
263 # (BINT, BINT) return nothing
264 # multiply two numbers in internal representation
b22b3e31 265 # modifies first arg, second need not be different from first
574bacfe 266 my ($c,$xv,$yv) = @_;
dccbb853 267
0716bf9b 268 my @prod = (); my ($prod,$car,$cty,$xi,$yi);
269 # since multiplying $x with $x fails, make copy in this case
574bacfe 270 $yv = [@$xv] if "$xv" eq "$yv"; # same references?
0716bf9b 271 for $xi (@$xv)
272 {
273 $car = 0; $cty = 0;
574bacfe 274
275 # slow variant
276# for $yi (@$yv)
277# {
278# $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
279# $prod[$cty++] =
280# $prod - ($car = int($prod * RBASE)) * $BASE; # see USE_MUL
281# }
282# $prod[$cty] += $car if $car; # need really to check for 0?
283# $xi = shift @prod;
284
285 # faster variant
286 # looping through this if $xi == 0 is silly - so optimize it away!
287 $xi = (shift @prod || 0), next if $xi == 0;
0716bf9b 288 for $yi (@$yv)
289 {
290 $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
574bacfe 291## this is actually a tad slower
292## $prod = $prod[$cty]; $prod += ($car + $xi * $yi); # no ||0 here
0716bf9b 293 $prod[$cty++] =
574bacfe 294 $prod - ($car = int($prod * $RBASE)) * $BASE; # see USE_MUL
0716bf9b 295 }
296 $prod[$cty] += $car if $car; # need really to check for 0?
297 $xi = shift @prod;
298 }
0716bf9b 299 push @$xv, @prod;
300 __strip_zeros($xv);
301 # normalize (handled last to save check for $y->is_zero()
302 return $xv;
303 }
304
ee15d750 305sub _mul_use_div
306 {
307 # (BINT, BINT) return nothing
308 # multiply two numbers in internal representation
309 # modifies first arg, second need not be different from first
310 my ($c,$xv,$yv) = @_;
311
312 my @prod = (); my ($prod,$car,$cty,$xi,$yi);
313 # since multiplying $x with $x fails, make copy in this case
314 $yv = [@$xv] if "$xv" eq "$yv"; # same references?
315 for $xi (@$xv)
316 {
317 $car = 0; $cty = 0;
318 # looping through this if $xi == 0 is silly - so optimize it away!
319 $xi = (shift @prod || 0), next if $xi == 0;
320 for $yi (@$yv)
321 {
322 $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
323 $prod[$cty++] =
324 $prod - ($car = int($prod / $BASE)) * $BASE;
325 }
326 $prod[$cty] += $car if $car; # need really to check for 0?
327 $xi = shift @prod;
328 }
329 push @$xv, @prod;
330 __strip_zeros($xv);
331 # normalize (handled last to save check for $y->is_zero()
332 return $xv;
333 }
334
335sub _div_use_mul
0716bf9b 336 {
b22b3e31 337 # ref to array, ref to array, modify first array and return remainder if
0716bf9b 338 # in list context
b22b3e31 339 # no longer handles sign
574bacfe 340 my ($c,$x,$yorg) = @_;
0716bf9b 341 my ($car,$bar,$prd,$dd,$xi,$yi,@q,$v2,$v1);
342
343 my (@d,$tmp,$q,$u2,$u1,$u0);
344
345 $car = $bar = $prd = 0;
346
347 my $y = [ @$yorg ];
348 if (($dd = int($BASE/($y->[-1]+1))) != 1)
349 {
350 for $xi (@$x)
351 {
352 $xi = $xi * $dd + $car;
353 $xi -= ($car = int($xi * $RBASE)) * $BASE; # see USE_MUL
354 }
355 push(@$x, $car); $car = 0;
356 for $yi (@$y)
357 {
358 $yi = $yi * $dd + $car;
359 $yi -= ($car = int($yi * $RBASE)) * $BASE; # see USE_MUL
360 }
361 }
362 else
363 {
364 push(@$x, 0);
365 }
366 @q = (); ($v2,$v1) = @$y[-2,-1];
367 $v2 = 0 unless $v2;
368 while ($#$x > $#$y)
369 {
370 ($u2,$u1,$u0) = @$x[-3..-1];
371 $u2 = 0 unless $u2;
372 #warn "oups v1 is 0, u0: $u0 $y->[-2] $y->[-1] l ",scalar @$y,"\n"
373 # if $v1 == 0;
ee15d750 374 # $q = (($u0 == $v1) ? 99999 : int(($u0*$BASE+$u1)/$v1));
375 $q = (($u0 == $v1) ? $MAX_VAL : int(($u0*$BASE+$u1)/$v1));
574bacfe 376 --$q while ($v2*$q > ($u0*$BASE+$u1-$q*$v1)*$BASE+$u2);
0716bf9b 377 if ($q)
378 {
379 ($car, $bar) = (0,0);
380 for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
381 {
382 $prd = $q * $y->[$yi] + $car;
383 $prd -= ($car = int($prd * $RBASE)) * $BASE; # see USE_MUL
574bacfe 384 $x->[$xi] += $BASE if ($bar = (($x->[$xi] -= $prd + $bar) < 0));
0716bf9b 385 }
386 if ($x->[-1] < $car + $bar)
387 {
388 $car = 0; --$q;
389 for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
390 {
574bacfe 391 $x->[$xi] -= $BASE
0716bf9b 392 if ($car = (($x->[$xi] += $y->[$yi] + $car) > $BASE));
393 }
394 }
395 }
396 pop(@$x); unshift(@q, $q);
397 }
398 if (wantarray)
399 {
400 @d = ();
401 if ($dd != 1)
402 {
403 $car = 0;
404 for $xi (reverse @$x)
405 {
406 $prd = $car * $BASE + $xi;
407 $car = $prd - ($tmp = int($prd / $dd)) * $dd; # see USE_MUL
408 unshift(@d, $tmp);
409 }
410 }
411 else
412 {
413 @d = @$x;
414 }
415 @$x = @q;
416 __strip_zeros($x);
417 __strip_zeros(\@d);
418 return ($x,\@d);
419 }
420 @$x = @q;
421 __strip_zeros($x);
422 return $x;
423 }
424
ee15d750 425sub _div_use_div
426 {
427 # ref to array, ref to array, modify first array and return remainder if
428 # in list context
429 # no longer handles sign
430 my ($c,$x,$yorg) = @_;
431 my ($car,$bar,$prd,$dd,$xi,$yi,@q,$v2,$v1);
432
433 my (@d,$tmp,$q,$u2,$u1,$u0);
434
435 $car = $bar = $prd = 0;
436
437 my $y = [ @$yorg ];
438 if (($dd = int($BASE/($y->[-1]+1))) != 1)
439 {
440 for $xi (@$x)
441 {
442 $xi = $xi * $dd + $car;
443 $xi -= ($car = int($xi / $BASE)) * $BASE;
444 }
445 push(@$x, $car); $car = 0;
446 for $yi (@$y)
447 {
448 $yi = $yi * $dd + $car;
449 $yi -= ($car = int($yi / $BASE)) * $BASE;
450 }
451 }
452 else
453 {
454 push(@$x, 0);
455 }
456 @q = (); ($v2,$v1) = @$y[-2,-1];
457 $v2 = 0 unless $v2;
458 while ($#$x > $#$y)
459 {
460 ($u2,$u1,$u0) = @$x[-3..-1];
461 $u2 = 0 unless $u2;
462 #warn "oups v1 is 0, u0: $u0 $y->[-2] $y->[-1] l ",scalar @$y,"\n"
463 # if $v1 == 0;
464 # $q = (($u0 == $v1) ? 99999 : int(($u0*$BASE+$u1)/$v1));
465 $q = (($u0 == $v1) ? $MAX_VAL : int(($u0*$BASE+$u1)/$v1));
466 --$q while ($v2*$q > ($u0*$BASE+$u1-$q*$v1)*$BASE+$u2);
467 if ($q)
468 {
469 ($car, $bar) = (0,0);
470 for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
471 {
472 $prd = $q * $y->[$yi] + $car;
473 $prd -= ($car = int($prd / $BASE)) * $BASE;
474 $x->[$xi] += $BASE if ($bar = (($x->[$xi] -= $prd + $bar) < 0));
475 }
476 if ($x->[-1] < $car + $bar)
477 {
478 $car = 0; --$q;
479 for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
480 {
481 $x->[$xi] -= $BASE
482 if ($car = (($x->[$xi] += $y->[$yi] + $car) > $BASE));
483 }
484 }
485 }
486 pop(@$x); unshift(@q, $q);
487 }
488 if (wantarray)
489 {
490 @d = ();
491 if ($dd != 1)
492 {
493 $car = 0;
494 for $xi (reverse @$x)
495 {
496 $prd = $car * $BASE + $xi;
497 $car = $prd - ($tmp = int($prd / $dd)) * $dd;
498 unshift(@d, $tmp);
499 }
500 }
501 else
502 {
503 @d = @$x;
504 }
505 @$x = @q;
506 __strip_zeros($x);
507 __strip_zeros(\@d);
508 return ($x,\@d);
509 }
510 @$x = @q;
511 __strip_zeros($x);
512 return $x;
513 }
514
dccbb853 515sub _mod
516 {
517 # if possible, use mod shortcut
518 my ($c,$x,$yo) = @_;
519
520 # slow way since $y to big
521 if (scalar @$yo > 1)
522 {
523 my ($xo,$rem) = _div($c,$x,$yo);
524 return $rem;
525 }
526 my $y = $yo->[0];
527 # both are single element
528 if (scalar @$x == 1)
529 {
530 $x->[0] %= $y;
531 return $x;
532 }
533
534 my $b = $BASE % $y;
535 if ($b == 0)
536 {
537 # when BASE % Y == 0 then (B * BASE) % Y == 0
538 # (B * BASE) % $y + A % Y => A % Y
539 # so need to consider only last element: O(1)
540 $x->[0] %= $y;
541 }
542 else
543 {
544 # else need to go trough all elemens: O(N)
545 # XXX not ready yet
546 my ($xo,$rem) = _div($c,$x,$yo);
547 return $rem;
548
549# my $i = 0; my $r = 1;
550# print "Multi: ";
551# foreach (@$x)
552# {
553# print "$_ $r $b $y\n";
554# print "\$_ % \$y = ",$_ % $y,"\n";
555# print "\$_ % \$y * \$b = ",($_ % $y) * $b,"\n";
556# $r += ($_ % $y) * $b;
557# print "$r $b $y =>";
558# $r %= $y if $r > $y;
559# print " $r\n";
560# }
561# $x->[0] = $r;
562 }
563 splice (@$x,1);
564 return $x;
565 }
566
0716bf9b 567##############################################################################
574bacfe 568# shifts
569
570sub _rsft
571 {
572 my ($c,$x,$y,$n) = @_;
573
574 if ($n != 10)
575 {
576 return; # we cant do this here, due to now _pow, so signal failure
577 }
578 else
579 {
580 # shortcut (faster) for shifting by 10)
581 # multiples of $BASE_LEN
582 my $dst = 0; # destination
583 my $src = _num($c,$y); # as normal int
dccbb853 584 my $rem = $src % $BASE_LEN; # remainder to shift
574bacfe 585 $src = int($src / $BASE_LEN); # source
586 if ($rem == 0)
587 {
588 splice (@$x,0,$src); # even faster, 38.4 => 39.3
589 }
590 else
591 {
592 my $len = scalar @$x - $src; # elems to go
593 my $vd; my $z = '0'x $BASE_LEN;
594 $x->[scalar @$x] = 0; # avoid || 0 test inside loop
595 while ($dst < $len)
596 {
597 $vd = $z.$x->[$src];
598 #print "$dst $src '$vd' ";
599 $vd = substr($vd,-$BASE_LEN,$BASE_LEN-$rem);
600 #print "'$vd' ";
601 $src++;
602 $vd = substr($z.$x->[$src],-$rem,$rem) . $vd;
603 #print "'$vd1' ";
604 #print "'$vd'\n";
605 $vd = substr($vd,-$BASE_LEN,$BASE_LEN) if length($vd) > $BASE_LEN;
606 $x->[$dst] = int($vd);
607 $dst++;
608 }
609 splice (@$x,$dst) if $dst > 0; # kill left-over array elems
610 pop @$x if $x->[-1] == 0; # kill last element if 0
611 } # else rem == 0
612 }
613 $x;
614 }
615
616sub _lsft
617 {
618 my ($c,$x,$y,$n) = @_;
619
620 if ($n != 10)
621 {
622 return; # we cant do this here, due to now _pow, so signal failure
623 }
624 else
625 {
626 # shortcut (faster) for shifting by 10) since we are in base 10eX
627 # multiples of $BASE_LEN:
628 my $src = scalar @$x; # source
629 my $len = _num($c,$y); # shift-len as normal int
dccbb853 630 my $rem = $len % $BASE_LEN; # remainder to shift
574bacfe 631 my $dst = $src + int($len/$BASE_LEN); # destination
632 my $vd; # further speedup
633 #print "src $src:",$x->[$src]||0," dst $dst:",$v->[$dst]||0," rem $rem\n";
634 $x->[$src] = 0; # avoid first ||0 for speed
635 my $z = '0' x $BASE_LEN;
636 while ($src >= 0)
637 {
638 $vd = $x->[$src]; $vd = $z.$vd;
639 #print "s $src d $dst '$vd' ";
640 $vd = substr($vd,-$BASE_LEN+$rem,$BASE_LEN-$rem);
641 #print "'$vd' ";
642 $vd .= $src > 0 ? substr($z.$x->[$src-1],-$BASE_LEN,$rem) : '0' x $rem;
643 #print "'$vd' ";
644 $vd = substr($vd,-$BASE_LEN,$BASE_LEN) if length($vd) > $BASE_LEN;
645 #print "'$vd'\n";
646 $x->[$dst] = int($vd);
647 $dst--; $src--;
648 }
649 # set lowest parts to 0
650 while ($dst >= 0) { $x->[$dst--] = 0; }
651 # fix spurios last zero element
652 splice @$x,-1 if $x->[-1] == 0;
653 #print "elems: "; my $i = 0;
654 #foreach (reverse @$v) { print "$i $_ "; $i++; } print "\n";
655 }
656 $x;
657 }
658
659##############################################################################
0716bf9b 660# testing
661
662sub _acmp
663 {
664 # internal absolute post-normalized compare (ignore signs)
665 # ref to array, ref to array, return <0, 0, >0
b22b3e31 666 # arrays must have at least one entry; this is not checked for
0716bf9b 667
574bacfe 668 my ($c,$cx, $cy) = @_;
0716bf9b 669
670 #print "$cx $cy\n";
671 my ($i,$a,$x,$y,$k);
672 # calculate length based on digits, not parts
574bacfe 673 $x = _len('',$cx); $y = _len('',$cy);
0716bf9b 674 # print "length: ",($x-$y),"\n";
574bacfe 675 my $lxy = $x - $y; # if different in length
676 return -1 if $lxy < 0;
677 return 1 if $lxy > 0;
0716bf9b 678 #print "full compare\n";
679 $i = 0; $a = 0;
680 # first way takes 5.49 sec instead of 4.87, but has the early out advantage
b22b3e31 681 # so grep is slightly faster, but more inflexible. hm. $_ instead of $k
0716bf9b 682 # yields 5.6 instead of 5.5 sec huh?
683 # manual way (abort if unequal, good for early ne)
684 my $j = scalar @$cx - 1;
685 while ($j >= 0)
686 {
687 # print "$cx->[$j] $cy->[$j] $a",$cx->[$j]-$cy->[$j],"\n";
688 last if ($a = $cx->[$j] - $cy->[$j]); $j--;
689 }
574bacfe 690 return 1 if $a > 0;
691 return -1 if $a < 0;
692 return 0; # equal
0716bf9b 693 # while it early aborts, it is even slower than the manual variant
694 #grep { return $a if ($a = $_ - $cy->[$i++]); } @$cx;
695 # grep way, go trough all (bad for early ne)
696 #grep { $a = $_ - $cy->[$i++]; } @$cx;
697 #return $a;
698 }
699
700sub _len
701 {
dccbb853 702 # compute number of digits in bigint, minus the sign
b22b3e31 703 # int() because add/sub sometimes leaves strings (like '00005') instead of
dccbb853 704 # int ('5') in this place, thus causing length() to report wrong length
574bacfe 705 my $cx = $_[1];
0716bf9b 706
574bacfe 707 return (@$cx-1)*$BASE_LEN+length(int($cx->[-1]));
0716bf9b 708 }
709
710sub _digit
711 {
712 # return the nth digit, negative values count backward
713 # zero is rightmost, so _digit(123,0) will give 3
574bacfe 714 my ($c,$x,$n) = @_;
0716bf9b 715
574bacfe 716 my $len = _len('',$x);
0716bf9b 717
718 $n = $len+$n if $n < 0; # -1 last, -2 second-to-last
719 $n = abs($n); # if negative was too big
720 $len--; $n = $len if $n > $len; # n to big?
721
574bacfe 722 my $elem = int($n / $BASE_LEN); # which array element
723 my $digit = $n % $BASE_LEN; # which digit in this element
0716bf9b 724 $elem = '0000'.@$x[$elem]; # get element padded with 0's
725 return substr($elem,-$digit-1,1);
726 }
727
728sub _zeros
729 {
730 # return amount of trailing zeros in decimal
731 # check each array elem in _m for having 0 at end as long as elem == 0
732 # Upon finding a elem != 0, stop
574bacfe 733 my $x = $_[1];
0716bf9b 734 my $zeros = 0; my $elem;
735 foreach my $e (@$x)
736 {
737 if ($e != 0)
738 {
574bacfe 739 $elem = "$e"; # preserve x
740 $elem =~ s/.*?(0*$)/$1/; # strip anything not zero
741 $zeros *= $BASE_LEN; # elems * 5
742 $zeros += CORE::length($elem); # count trailing zeros
743 last; # early out
0716bf9b 744 }
574bacfe 745 $zeros ++; # real else branch: 50% slower!
0716bf9b 746 }
747 return $zeros;
748 }
749
750##############################################################################
751# _is_* routines
752
753sub _is_zero
754 {
755 # return true if arg (BINT or num_str) is zero (array '+', '0')
574bacfe 756 my $x = $_[1];
0716bf9b 757 return (((scalar @$x == 1) && ($x->[0] == 0))) <=> 0;
758 }
759
760sub _is_even
761 {
762 # return true if arg (BINT or num_str) is even
574bacfe 763 my $x = $_[1];
0716bf9b 764 return (!($x->[0] & 1)) <=> 0;
765 }
766
767sub _is_odd
768 {
769 # return true if arg (BINT or num_str) is even
574bacfe 770 my $x = $_[1];
0716bf9b 771 return (($x->[0] & 1)) <=> 0;
772 }
773
774sub _is_one
775 {
776 # return true if arg (BINT or num_str) is one (array '+', '1')
574bacfe 777 my $x = $_[1];
0716bf9b 778 return (scalar @$x == 1) && ($x->[0] == 1) <=> 0;
779 }
780
781sub __strip_zeros
782 {
783 # internal normalization function that strips leading zeros from the array
784 # args: ref to array
0716bf9b 785 my $s = shift;
786
787 my $cnt = scalar @$s; # get count of parts
788 my $i = $cnt-1;
789 #print "strip: cnt $cnt i $i\n";
790 # '0', '3', '4', '0', '0',
791 # 0 1 2 3 4
792 # cnt = 5, i = 4
793 # i = 4
794 # i = 3
795 # => fcnt = cnt - i (5-2 => 3, cnt => 5-1 = 4, throw away from 4th pos)
796 # >= 1: skip first part (this can be zero)
797 while ($i > 0) { last if $s->[$i] != 0; $i--; }
798 $i++; splice @$s,$i if ($i < $cnt); # $i cant be 0
799 return $s;
800 }
801
802###############################################################################
803# check routine to test internal state of corruptions
804
805sub _check
806 {
bd05a461 807 # used by the test suite
574bacfe 808 my $x = $_[1];
0716bf9b 809
0716bf9b 810 return "$x is not a reference" if !ref($x);
811
812 # are all parts are valid?
813 my $i = 0; my $j = scalar @$x; my ($e,$try);
814 while ($i < $j)
815 {
816 $e = $x->[$i]; $e = 'undef' unless defined $e;
817 $try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e)";
818 last if $e !~ /^[+]?[0-9]+$/;
dccbb853 819 $try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e) (stringify)";
820 last if "$e" !~ /^[+]?[0-9]+$/;
821 $try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e) (cat-stringify)";
822 last if '' . "$e" !~ /^[+]?[0-9]+$/;
0716bf9b 823 $try = ' < 0 || >= $BASE; '."($x, $e)";
824 last if $e <0 || $e >= $BASE;
825 # this test is disabled, since new/bnorm and certain ops (like early out
826 # in add/sub) are allowed/expected to leave '00000' in some elements
827 #$try = '=~ /^00+/; '."($x, $e)";
828 #last if $e =~ /^00+/;
829 $i++;
830 }
831 return "Illegal part '$e' at pos $i (tested: $try)" if $i < $j;
832 return 0;
833 }
834
8351;
836__END__
837
838=head1 NAME
839
840Math::BigInt::Calc - Pure Perl module to support Math::BigInt
841
842=head1 SYNOPSIS
843
ee15d750 844Provides support for big integer calculations. Not intended to be used by other
845modules (except Math::BigInt::Cached). Other modules which sport the same
846functions can also be used to support Math::Bigint, like Math::BigInt::Pari.
0716bf9b 847
848=head1 DESCRIPTION
849
850In order to allow for multiple big integer libraries, Math::BigInt
851was rewritten to use library modules for core math routines. Any
852module which follows the same API as this can be used instead by
853using the following call:
854
ee15d750 855 use Math::BigInt lib => 'libname';
0716bf9b 856
857=head1 EXPORT
858
bd05a461 859The following functions MUST be defined in order to support
0716bf9b 860the use by Math::BigInt:
861
862 _new(string) return ref to new object from ref to decimal string
863 _zero() return a new object with value 0
864 _one() return a new object with value 1
865
866 _str(obj) return ref to a string representing the object
867 _num(obj) returns a Perl integer/floating point number
868 NOTE: because of Perl numeric notation defaults,
869 the _num'ified obj may lose accuracy due to
870 machine-dependend floating point size limitations
871
872 _add(obj,obj) Simple addition of two objects
873 _mul(obj,obj) Multiplication of two objects
874 _div(obj,obj) Division of the 1st object by the 2nd
b22b3e31 875 In list context, returns (result,remainder).
876 NOTE: this is integer math, so no
877 fractional part will be returned.
878 _sub(obj,obj) Simple subtraction of 1 object from another
0716bf9b 879 a third, optional parameter indicates that the params
880 are swapped. In this case, the first param needs to
881 be preserved, while you can destroy the second.
882 sub (x,y,1) => return x - y and keep x intact!
e745a66c 883 _dec(obj) decrement object by one (input is garant. to be > 0)
884 _inc(obj) increment object by one
885
0716bf9b 886
887 _acmp(obj,obj) <=> operator for objects (return -1, 0 or 1)
888
889 _len(obj) returns count of the decimal digits of the object
890 _digit(obj,n) returns the n'th decimal digit of object
891
892 _is_one(obj) return true if argument is +1
893 _is_zero(obj) return true if argument is 0
894 _is_even(obj) return true if argument is even (0,2,4,6..)
895 _is_odd(obj) return true if argument is odd (1,3,5,7..)
896
897 _copy return a ref to a true copy of the object
898
899 _check(obj) check whether internal representation is still intact
900 return 0 for ok, otherwise error message as string
901
bd05a461 902The following functions are optional, and can be defined if the underlying lib
ee15d750 903has a fast way to do them. If undefined, Math::BigInt will use a pure, but
bd05a461 904slow, Perl way as fallback to emulate these:
0716bf9b 905
906 _from_hex(str) return ref to new object from ref to hexadecimal string
907 _from_bin(str) return ref to new object from ref to binary string
908
ee15d750 909 _as_hex(str) return ref to scalar string containing the value as
910 unsigned hex string, with the '0x' prepended.
911 Leading zeros must be stripped.
912 _as_bin(str) Like as_hex, only as binary string containing only
913 zeros and ones. Leading zeros must be stripped and a
914 '0b' must be prepended.
915
0716bf9b 916 _rsft(obj,N,B) shift object in base B by N 'digits' right
dccbb853 917 For unsupported bases B, return undef to signal failure
0716bf9b 918 _lsft(obj,N,B) shift object in base B by N 'digits' left
dccbb853 919 For unsupported bases B, return undef to signal failure
0716bf9b 920
921 _xor(obj1,obj2) XOR (bit-wise) object 1 with object 2
dccbb853 922 Note: XOR, AND and OR pad with zeros if size mismatches
0716bf9b 923 _and(obj1,obj2) AND (bit-wise) object 1 with object 2
924 _or(obj1,obj2) OR (bit-wise) object 1 with object 2
925
dccbb853 926 _mod(obj,obj) Return remainder of div of the 1st by the 2nd object
0716bf9b 927 _sqrt(obj) return the square root of object
928 _pow(obj,obj) return object 1 to the power of object 2
929 _gcd(obj,obj) return Greatest Common Divisor of two objects
930
b22b3e31 931 _zeros(obj) return number of trailing decimal zeros
0716bf9b 932
b22b3e31 933Input strings come in as unsigned but with prefix (i.e. as '123', '0xabc'
0716bf9b 934or '0b1101').
935
b22b3e31 936Testing of input parameter validity is done by the caller, so you need not
574bacfe 937worry about underflow (f.i. in C<_sub()>, C<_dec()>) nor about division by
938zero or similar cases.
939
940The first parameter can be modified, that includes the possibility that you
941return a reference to a completely different object instead. Although keeping
dccbb853 942the reference and just changing it's contents is prefered over creating and
943returning a different reference.
574bacfe 944
945Return values are always references to objects or strings. Exceptions are
946C<_lsft()> and C<_rsft()>, which return undef if they can not shift the
947argument. This is used to delegate shifting of bases different than 10 back
dccbb853 948to Math::BigInt, which will use some generic code to calculate the result.
574bacfe 949
950=head1 WRAP YOUR OWN
951
952If you want to port your own favourite c-lib for big numbers to the
953Math::BigInt interface, you can take any of the already existing modules as
954a rough guideline. You should really wrap up the latest BigInt and BigFloat
bd05a461 955testsuites with your module, and replace in them any of the following:
574bacfe 956
957 use Math::BigInt;
958
bd05a461 959by this:
574bacfe 960
961 use Math::BigInt lib => 'yourlib';
962
963This way you ensure that your library really works 100% within Math::BigInt.
0716bf9b 964
965=head1 LICENSE
966
967This program is free software; you may redistribute it and/or modify it under
968the same terms as Perl itself.
969
970=head1 AUTHORS
971
972Original math code by Mark Biggar, rewritten by Tels L<http://bloodgate.com/>
973in late 2000, 2001.
974Seperated from BigInt and shaped API with the help of John Peacock.
975
976=head1 SEE ALSO
977
ee15d750 978L<Math::BigInt>, L<Math::BigFloat>, L<Math::BigInt::BitVect>,
979L<Math::BigInt::GMP>, L<Math::BigInt::Cached> and L<Math::BigInt::Pari>.
0716bf9b 980
981=cut