projects / p5sagit/p5-mst-13.2.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | history | HEAD
Math-BigInt v1.49 released
[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
61f5c3f5 11$VERSION = '0.20';
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) = @_;
0716bf9b 333
334 # for each in Y, add Y to X and carry. If after that, something is left in
335 # X, foreach in X add carry to X and then return X, carry
336 # Trades one "$j++" for having to shift arrays, $j could be made integer
b22b3e31 337 # but this would impose a limit to number-length of 2**32.
0716bf9b 338 my $i; my $car = 0; my $j = 0;
339 for $i (@$y)
340 {
e745a66c 341 $x->[$j] -= $BASE if $car = (($x->[$j] += $i + $car) >= $BASE) ? 1 : 0;
0716bf9b 342 $j++;
343 }
344 while ($car != 0)
345 {
346 $x->[$j] -= $BASE if $car = (($x->[$j] += $car) >= $BASE) ? 1 : 0; $j++;
347 }
61f5c3f5 348 $x;
e745a66c 349 }
350
351sub _inc
352 {
353 # (ref to int_num_array, ref to int_num_array)
354 # routine to add 1 to a base 1eX numbers
355 # This routine clobbers up array x, but not y.
356 my ($c,$x) = @_;
357
358 for my $i (@$x)
359 {
360 return $x if (($i += 1) < $BASE); # early out
61f5c3f5 361 $i = 0; # overflow, next
e745a66c 362 }
61f5c3f5 363 push @$x,1 if ($x->[-1] == 0); # last overflowed, so extend
364 $x;
e745a66c 365 }
366
367sub _dec
368 {
369 # (ref to int_num_array, ref to int_num_array)
370 # routine to add 1 to a base 1eX numbers
371 # This routine clobbers up array x, but not y.
372 my ($c,$x) = @_;
373
61f5c3f5 374 my $MAX = $BASE-1; # since MAX_VAL based on MBASE
e745a66c 375 for my $i (@$x)
376 {
377 last if (($i -= 1) >= 0); # early out
61f5c3f5 378 $i = $MAX; # overflow, next
e745a66c 379 }
380 pop @$x if $x->[-1] == 0 && @$x > 1; # last overflowed (but leave 0)
61f5c3f5 381 $x;
0716bf9b 382 }
383
384sub _sub
385 {
386 # (ref to int_num_array, ref to int_num_array)
574bacfe 387 # subtract base 1eX numbers -- stolen from Knuth Vol 2 pg 232, $x > $y
b22b3e31 388 # subtract Y from X (X is always greater/equal!) by modifying x in place
574bacfe 389 my ($c,$sx,$sy,$s) = @_;
0716bf9b 390
391 my $car = 0; my $i; my $j = 0;
392 if (!$s)
393 {
394 #print "case 2\n";
395 for $i (@$sx)
396 {
397 last unless defined $sy->[$j] || $car;
0716bf9b 398 $i += $BASE if $car = (($i -= ($sy->[$j] || 0) + $car) < 0); $j++;
0716bf9b 399 }
400 # might leave leading zeros, so fix that
394e6ffb 401 return __strip_zeros($sx);
0716bf9b 402 }
394e6ffb 403 #print "case 1 (swap)\n";
404 for $i (@$sx)
0716bf9b 405 {
394e6ffb 406 last unless defined $sy->[$j] || $car;
407 $sy->[$j] += $BASE
408 if $car = (($sy->[$j] = $i-($sy->[$j]||0) - $car) < 0);
409 $j++;
0716bf9b 410 }
394e6ffb 411 # might leave leading zeros, so fix that
412 __strip_zeros($sy);
0716bf9b 413 }
414
ee15d750 415sub _mul_use_mul
0716bf9b 416 {
417 # (BINT, BINT) return nothing
418 # multiply two numbers in internal representation
b22b3e31 419 # modifies first arg, second need not be different from first
574bacfe 420 my ($c,$xv,$yv) = @_;
dccbb853 421
394e6ffb 422 # shortcut for two very short numbers
423 # +0 since part maybe string '00001' from new()
61f5c3f5 424 # works also if xv and yv are the same reference
394e6ffb 425 if ((@$xv == 1) && (@$yv == 1)
426 && (length($xv->[0]+0) <= $BASE_LEN2)
427 && (length($yv->[0]+0) <= $BASE_LEN2))
428 {
429 $xv->[0] *= $yv->[0];
430 return $xv;
431 }
432
0716bf9b 433 # since multiplying $x with $x fails, make copy in this case
574bacfe 434 $yv = [@$xv] if "$xv" eq "$yv"; # same references?
61f5c3f5 435 if ($LEN_CONVERT != 0)
436 {
437 $c->_to_small($xv); $c->_to_small($yv);
438 }
439
440 my @prod = (); my ($prod,$car,$cty,$xi,$yi);
441
0716bf9b 442 for $xi (@$xv)
443 {
444 $car = 0; $cty = 0;
574bacfe 445
446 # slow variant
447# for $yi (@$yv)
448# {
449# $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
450# $prod[$cty++] =
61f5c3f5 451# $prod - ($car = int($prod * RBASE)) * $MBASE; # see USE_MUL
574bacfe 452# }
453# $prod[$cty] += $car if $car; # need really to check for 0?
454# $xi = shift @prod;
455
456 # faster variant
457 # looping through this if $xi == 0 is silly - so optimize it away!
458 $xi = (shift @prod || 0), next if $xi == 0;
0716bf9b 459 for $yi (@$yv)
460 {
461 $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
574bacfe 462## this is actually a tad slower
463## $prod = $prod[$cty]; $prod += ($car + $xi * $yi); # no ||0 here
0716bf9b 464 $prod[$cty++] =
61f5c3f5 465 $prod - ($car = int($prod * $RBASE)) * $MBASE; # see USE_MUL
0716bf9b 466 }
467 $prod[$cty] += $car if $car; # need really to check for 0?
027dc388 468 $xi = shift @prod || 0; # || 0 makes v5.005_3 happy
0716bf9b 469 }
0716bf9b 470 push @$xv, @prod;
61f5c3f5 471 if ($LEN_CONVERT != 0)
472 {
473 $c->_to_large($yv);
474 $c->_to_large($xv);
475 }
476 else
477 {
478 __strip_zeros($xv);
479 }
480 $xv;
0716bf9b 481 }
482
ee15d750 483sub _mul_use_div
484 {
485 # (BINT, BINT) return nothing
486 # multiply two numbers in internal representation
487 # modifies first arg, second need not be different from first
488 my ($c,$xv,$yv) = @_;
489
394e6ffb 490 # shortcut for two very short numbers
491 # +0 since part maybe string '00001' from new()
61f5c3f5 492 # works also if xv and yv are the same reference
394e6ffb 493 if ((@$xv == 1) && (@$yv == 1)
494 && (length($xv->[0]+0) <= $BASE_LEN2)
495 && (length($yv->[0]+0) <= $BASE_LEN2))
496 {
497 $xv->[0] *= $yv->[0];
498 return $xv;
499 }
61f5c3f5 500
ee15d750 501 # since multiplying $x with $x fails, make copy in this case
502 $yv = [@$xv] if "$xv" eq "$yv"; # same references?
61f5c3f5 503 if ($LEN_CONVERT != 0)
504 {
505 $c->_to_small($xv); $c->_to_small($yv);
506 }
507
508 my @prod = (); my ($prod,$car,$cty,$xi,$yi);
ee15d750 509 for $xi (@$xv)
510 {
511 $car = 0; $cty = 0;
512 # looping through this if $xi == 0 is silly - so optimize it away!
513 $xi = (shift @prod || 0), next if $xi == 0;
514 for $yi (@$yv)
515 {
516 $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
517 $prod[$cty++] =
61f5c3f5 518 $prod - ($car = int($prod / $MBASE)) * $MBASE;
ee15d750 519 }
520 $prod[$cty] += $car if $car; # need really to check for 0?
027dc388 521 $xi = shift @prod || 0; # || 0 makes v5.005_3 happy
ee15d750 522 }
523 push @$xv, @prod;
61f5c3f5 524 if ($LEN_CONVERT != 0)
525 {
526 $c->_to_large($yv);
527 $c->_to_large($xv);
528 }
529 else
530 {
531 __strip_zeros($xv);
532 }
533 $xv;
ee15d750 534 }
535
536sub _div_use_mul
0716bf9b 537 {
b22b3e31 538 # ref to array, ref to array, modify first array and return remainder if
0716bf9b 539 # in list context
574bacfe 540 my ($c,$x,$yorg) = @_;
0716bf9b 541
61f5c3f5 542 if (@$x == 1 && @$yorg == 1)
543 {
544 # shortcut, $y is smaller than $x
545 if (wantarray)
546 {
547 my $r = [ $x->[0] % $yorg->[0] ];
548 $x->[0] = int($x->[0] / $yorg->[0]);
549 return ($x,$r);
550 }
551 else
552 {
553 $x->[0] = int($x->[0] / $yorg->[0]);
554 return $x;
555 }
556 }
0716bf9b 557
0716bf9b 558 my $y = [ @$yorg ];
61f5c3f5 559 if ($LEN_CONVERT != 0)
560 {
561 $c->_to_small($x); $c->_to_small($y);
562 }
563
564 my ($car,$bar,$prd,$dd,$xi,$yi,@q,$v2,$v1,@d,$tmp,$q,$u2,$u1,$u0);
565
566 $car = $bar = $prd = 0;
567 if (($dd = int($MBASE/($y->[-1]+1))) != 1)
0716bf9b 568 {
569 for $xi (@$x)
570 {
571 $xi = $xi * $dd + $car;
61f5c3f5 572 $xi -= ($car = int($xi * $RBASE)) * $MBASE; # see USE_MUL
0716bf9b 573 }
574 push(@$x, $car); $car = 0;
575 for $yi (@$y)
576 {
577 $yi = $yi * $dd + $car;
61f5c3f5 578 $yi -= ($car = int($yi * $RBASE)) * $MBASE; # see USE_MUL
0716bf9b 579 }
580 }
581 else
582 {
583 push(@$x, 0);
584 }
585 @q = (); ($v2,$v1) = @$y[-2,-1];
586 $v2 = 0 unless $v2;
587 while ($#$x > $#$y)
588 {
589 ($u2,$u1,$u0) = @$x[-3..-1];
590 $u2 = 0 unless $u2;
591 #warn "oups v1 is 0, u0: $u0 $y->[-2] $y->[-1] l ",scalar @$y,"\n"
592 # if $v1 == 0;
61f5c3f5 593 $q = (($u0 == $v1) ? $MAX_VAL : int(($u0*$MBASE+$u1)/$v1));
594 --$q while ($v2*$q > ($u0*$MBASE+$u1-$q*$v1)*$MBASE+$u2);
0716bf9b 595 if ($q)
596 {
597 ($car, $bar) = (0,0);
598 for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
599 {
600 $prd = $q * $y->[$yi] + $car;
61f5c3f5 601 $prd -= ($car = int($prd * $RBASE)) * $MBASE; # see USE_MUL
602 $x->[$xi] += $MBASE if ($bar = (($x->[$xi] -= $prd + $bar) < 0));
0716bf9b 603 }
604 if ($x->[-1] < $car + $bar)
605 {
606 $car = 0; --$q;
607 for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
608 {
61f5c3f5 609 $x->[$xi] -= $MBASE
610 if ($car = (($x->[$xi] += $y->[$yi] + $car) > $MBASE));
0716bf9b 611 }
612 }
613 }
614 pop(@$x); unshift(@q, $q);
615 }
616 if (wantarray)
617 {
618 @d = ();
619 if ($dd != 1)
620 {
621 $car = 0;
622 for $xi (reverse @$x)
623 {
61f5c3f5 624 $prd = $car * $MBASE + $xi;
0716bf9b 625 $car = $prd - ($tmp = int($prd / $dd)) * $dd; # see USE_MUL
626 unshift(@d, $tmp);
627 }
628 }
629 else
630 {
631 @d = @$x;
632 }
633 @$x = @q;
61f5c3f5 634 my $d = \@d;
635 if ($LEN_CONVERT != 0)
636 {
637 $c->_to_large($x); $c->_to_large($d);
638 }
639 else
640 {
641 __strip_zeros($x);
642 __strip_zeros($d);
643 }
644 return ($x,$d);
0716bf9b 645 }
646 @$x = @q;
61f5c3f5 647 if ($LEN_CONVERT != 0)
648 {
649 $c->_to_large($x);
650 }
651 else
652 {
653 __strip_zeros($x);
654 }
655 $x;
0716bf9b 656 }
657
ee15d750 658sub _div_use_div
659 {
660 # ref to array, ref to array, modify first array and return remainder if
661 # in list context
ee15d750 662 my ($c,$x,$yorg) = @_;
ee15d750 663
61f5c3f5 664 if (@$x == 1 && @$yorg == 1)
665 {
666 # shortcut, $y is smaller than $x
667 if (wantarray)
668 {
669 my $r = [ $x->[0] % $yorg->[0] ];
670 $x->[0] = int($x->[0] / $yorg->[0]);
671 return ($x,$r);
672 }
673 else
674 {
675 $x->[0] = int($x->[0] / $yorg->[0]);
676 return $x;
677 }
678 }
ee15d750 679
ee15d750 680 my $y = [ @$yorg ];
61f5c3f5 681 if ($LEN_CONVERT != 0)
682 {
683 $c->_to_small($x); $c->_to_small($y);
684 }
685
686 my ($car,$bar,$prd,$dd,$xi,$yi,@q,$v2,$v1,@d,$tmp,$q,$u2,$u1,$u0);
687
688 $car = $bar = $prd = 0;
689 if (($dd = int($MBASE/($y->[-1]+1))) != 1)
ee15d750 690 {
691 for $xi (@$x)
692 {
693 $xi = $xi * $dd + $car;
61f5c3f5 694 $xi -= ($car = int($xi / $MBASE)) * $MBASE;
ee15d750 695 }
696 push(@$x, $car); $car = 0;
697 for $yi (@$y)
698 {
699 $yi = $yi * $dd + $car;
61f5c3f5 700 $yi -= ($car = int($yi / $MBASE)) * $MBASE;
ee15d750 701 }
702 }
703 else
704 {
705 push(@$x, 0);
706 }
707 @q = (); ($v2,$v1) = @$y[-2,-1];
708 $v2 = 0 unless $v2;
709 while ($#$x > $#$y)
710 {
711 ($u2,$u1,$u0) = @$x[-3..-1];
712 $u2 = 0 unless $u2;
713 #warn "oups v1 is 0, u0: $u0 $y->[-2] $y->[-1] l ",scalar @$y,"\n"
714 # if $v1 == 0;
61f5c3f5 715 $q = (($u0 == $v1) ? $MAX_VAL : int(($u0*$MBASE+$u1)/$v1));
716 --$q while ($v2*$q > ($u0*$MBASE+$u1-$q*$v1)*$MBASE+$u2);
ee15d750 717 if ($q)
718 {
719 ($car, $bar) = (0,0);
720 for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
721 {
722 $prd = $q * $y->[$yi] + $car;
61f5c3f5 723 $prd -= ($car = int($prd / $MBASE)) * $MBASE;
724 $x->[$xi] += $MBASE if ($bar = (($x->[$xi] -= $prd + $bar) < 0));
ee15d750 725 }
726 if ($x->[-1] < $car + $bar)
727 {
728 $car = 0; --$q;
729 for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
730 {
61f5c3f5 731 $x->[$xi] -= $MBASE
732 if ($car = (($x->[$xi] += $y->[$yi] + $car) > $MBASE));
ee15d750 733 }
734 }
735 }
61f5c3f5 736 pop(@$x); unshift(@q, $q);
ee15d750 737 }
738 if (wantarray)
739 {
740 @d = ();
741 if ($dd != 1)
742 {
743 $car = 0;
744 for $xi (reverse @$x)
745 {
61f5c3f5 746 $prd = $car * $MBASE + $xi;
ee15d750 747 $car = $prd - ($tmp = int($prd / $dd)) * $dd;
748 unshift(@d, $tmp);
749 }
750 }
751 else
752 {
753 @d = @$x;
754 }
755 @$x = @q;
61f5c3f5 756 my $d = \@d;
757 if ($LEN_CONVERT != 0)
758 {
759 $c->_to_large($x); $c->_to_large($d);
760 }
761 else
762 {
763 __strip_zeros($x);
764 __strip_zeros($d);
765 }
766 return ($x,$d);
ee15d750 767 }
768 @$x = @q;
61f5c3f5 769 if ($LEN_CONVERT != 0)
770 {
771 $c->_to_large($x);
772 }
773 else
774 {
775 __strip_zeros($x);
776 }
777 $x;
ee15d750 778 }
779
394e6ffb 780##############################################################################
781# testing
782
783sub _acmp
784 {
785 # internal absolute post-normalized compare (ignore signs)
786 # ref to array, ref to array, return <0, 0, >0
787 # arrays must have at least one entry; this is not checked for
788
789 my ($c,$cx,$cy) = @_;
790
61f5c3f5 791 # fast comp based on array elements
394e6ffb 792 my $lxy = scalar @$cx - scalar @$cy;
793 return -1 if $lxy < 0; # already differs, ret
794 return 1 if $lxy > 0; # ditto
795
796 # now calculate length based on digits, not parts
797 $lxy = _len($c,$cx) - _len($c,$cy); # difference
798 return -1 if $lxy < 0;
799 return 1 if $lxy > 0;
800
801 # hm, same lengths, but same contents?
802 my $i = 0; my $a;
803 # first way takes 5.49 sec instead of 4.87, but has the early out advantage
804 # so grep is slightly faster, but more inflexible. hm. $_ instead of $k
805 # yields 5.6 instead of 5.5 sec huh?
806 # manual way (abort if unequal, good for early ne)
807 my $j = scalar @$cx - 1;
808 while ($j >= 0)
809 {
810 last if ($a = $cx->[$j] - $cy->[$j]); $j--;
811 }
812 return 1 if $a > 0;
813 return -1 if $a < 0;
61f5c3f5 814 0; # equal
815
394e6ffb 816 # while it early aborts, it is even slower than the manual variant
817 #grep { return $a if ($a = $_ - $cy->[$i++]); } @$cx;
818 # grep way, go trough all (bad for early ne)
819 #grep { $a = $_ - $cy->[$i++]; } @$cx;
820 #return $a;
821 }
822
823sub _len
824 {
825 # compute number of digits in bigint, minus the sign
826
827 # int() because add/sub sometimes leaves strings (like '00005') instead of
828 # '5' in this place, thus causing length() to report wrong length
829 my $cx = $_[1];
830
831 return (@$cx-1)*$BASE_LEN+length(int($cx->[-1]));
832 }
833
834sub _digit
835 {
836 # return the nth digit, negative values count backward
837 # zero is rightmost, so _digit(123,0) will give 3
838 my ($c,$x,$n) = @_;
839
840 my $len = _len('',$x);
841
842 $n = $len+$n if $n < 0; # -1 last, -2 second-to-last
843 $n = abs($n); # if negative was too big
844 $len--; $n = $len if $n > $len; # n to big?
845
846 my $elem = int($n / $BASE_LEN); # which array element
847 my $digit = $n % $BASE_LEN; # which digit in this element
848 $elem = '0000'.@$x[$elem]; # get element padded with 0's
849 return substr($elem,-$digit-1,1);
850 }
851
852sub _zeros
853 {
854 # return amount of trailing zeros in decimal
855 # check each array elem in _m for having 0 at end as long as elem == 0
856 # Upon finding a elem != 0, stop
857 my $x = $_[1];
858 my $zeros = 0; my $elem;
859 foreach my $e (@$x)
860 {
861 if ($e != 0)
862 {
863 $elem = "$e"; # preserve x
864 $elem =~ s/.*?(0*$)/$1/; # strip anything not zero
865 $zeros *= $BASE_LEN; # elems * 5
61f5c3f5 866 $zeros += length($elem); # count trailing zeros
394e6ffb 867 last; # early out
868 }
869 $zeros ++; # real else branch: 50% slower!
870 }
61f5c3f5 871 $zeros;
394e6ffb 872 }
873
874##############################################################################
875# _is_* routines
876
877sub _is_zero
878 {
879 # return true if arg (BINT or num_str) is zero (array '+', '0')
880 my $x = $_[1];
61f5c3f5 881
882 (((scalar @$x == 1) && ($x->[0] == 0))) <=> 0;
394e6ffb 883 }
884
885sub _is_even
886 {
887 # return true if arg (BINT or num_str) is even
888 my $x = $_[1];
61f5c3f5 889 (!($x->[0] & 1)) <=> 0;
394e6ffb 890 }
891
892sub _is_odd
893 {
894 # return true if arg (BINT or num_str) is even
895 my $x = $_[1];
61f5c3f5 896
897 (($x->[0] & 1)) <=> 0;
394e6ffb 898 }
899
900sub _is_one
901 {
902 # return true if arg (BINT or num_str) is one (array '+', '1')
903 my $x = $_[1];
61f5c3f5 904
905 (scalar @$x == 1) && ($x->[0] == 1) <=> 0;
394e6ffb 906 }
907
908sub __strip_zeros
909 {
910 # internal normalization function that strips leading zeros from the array
911 # args: ref to array
912 my $s = shift;
913
914 my $cnt = scalar @$s; # get count of parts
915 my $i = $cnt-1;
916 push @$s,0 if $i < 0; # div might return empty results, so fix it
917
61f5c3f5 918 return $s if @$s == 1; # early out
919
394e6ffb 920 #print "strip: cnt $cnt i $i\n";
921 # '0', '3', '4', '0', '0',
922 # 0 1 2 3 4
923 # cnt = 5, i = 4
924 # i = 4
925 # i = 3
926 # => fcnt = cnt - i (5-2 => 3, cnt => 5-1 = 4, throw away from 4th pos)
927 # >= 1: skip first part (this can be zero)
928 while ($i > 0) { last if $s->[$i] != 0; $i--; }
929 $i++; splice @$s,$i if ($i < $cnt); # $i cant be 0
930 $s;
931 }
932
933###############################################################################
934# check routine to test internal state of corruptions
935
936sub _check
937 {
938 # used by the test suite
939 my $x = $_[1];
940
941 return "$x is not a reference" if !ref($x);
942
943 # are all parts are valid?
944 my $i = 0; my $j = scalar @$x; my ($e,$try);
945 while ($i < $j)
946 {
947 $e = $x->[$i]; $e = 'undef' unless defined $e;
948 $try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e)";
949 last if $e !~ /^[+]?[0-9]+$/;
950 $try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e) (stringify)";
951 last if "$e" !~ /^[+]?[0-9]+$/;
952 $try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e) (cat-stringify)";
953 last if '' . "$e" !~ /^[+]?[0-9]+$/;
954 $try = ' < 0 || >= $BASE; '."($x, $e)";
955 last if $e <0 || $e >= $BASE;
956 # this test is disabled, since new/bnorm and certain ops (like early out
957 # in add/sub) are allowed/expected to leave '00000' in some elements
958 #$try = '=~ /^00+/; '."($x, $e)";
959 #last if $e =~ /^00+/;
960 $i++;
961 }
962 return "Illegal part '$e' at pos $i (tested: $try)" if $i < $j;
963 return 0;
964 }
965
966
967###############################################################################
968###############################################################################
969# some optional routines to make BigInt faster
970
dccbb853 971sub _mod
972 {
973 # if possible, use mod shortcut
974 my ($c,$x,$yo) = @_;
975
976 # slow way since $y to big
977 if (scalar @$yo > 1)
978 {
979 my ($xo,$rem) = _div($c,$x,$yo);
980 return $rem;
981 }
982 my $y = $yo->[0];
027dc388 983 # both are single element arrays
dccbb853 984 if (scalar @$x == 1)
985 {
986 $x->[0] %= $y;
987 return $x;
988 }
989
61f5c3f5 990 # @y is single element, but @x has more than one
dccbb853 991 my $b = $BASE % $y;
992 if ($b == 0)
993 {
994 # when BASE % Y == 0 then (B * BASE) % Y == 0
995 # (B * BASE) % $y + A % Y => A % Y
996 # so need to consider only last element: O(1)
997 $x->[0] %= $y;
998 }
027dc388 999 elsif ($b == 1)
1000 {
1001 # else need to go trough all elements: O(N), but loop is a bit simplified
1002 my $r = 0;
1003 foreach (@$x)
1004 {
1005 $r += $_ % $y;
1006 $r %= $y;
1007 }
1008 $r = 0 if $r == $y;
1009 $x->[0] = $r;
1010 }
dccbb853 1011 else
1012 {
027dc388 1013 # else need to go trough all elements: O(N)
1014 my $r = 0; my $bm = 1;
1015 foreach (@$x)
1016 {
1017 $r += ($_ % $y) * $bm;
1018 $bm *= $b;
1019 $bm %= $y;
1020 $r %= $y;
1021 }
1022 $r = 0 if $r == $y;
1023 $x->[0] = $r;
dccbb853 1024 }
1025 splice (@$x,1);
61f5c3f5 1026 $x;
dccbb853 1027 }
1028
0716bf9b 1029##############################################################################
574bacfe 1030# shifts
1031
1032sub _rsft
1033 {
1034 my ($c,$x,$y,$n) = @_;
1035
1036 if ($n != 10)
1037 {
61f5c3f5 1038 $n = _new($c,\$n); return _div($c,$x, _pow($c,$n,$y));
1039 }
1040
1041 # shortcut (faster) for shifting by 10)
1042 # multiples of $BASE_LEN
1043 my $dst = 0; # destination
1044 my $src = _num($c,$y); # as normal int
1045 my $rem = $src % $BASE_LEN; # remainder to shift
1046 $src = int($src / $BASE_LEN); # source
1047 if ($rem == 0)
1048 {
1049 splice (@$x,0,$src); # even faster, 38.4 => 39.3
574bacfe 1050 }
1051 else
1052 {
61f5c3f5 1053 my $len = scalar @$x - $src; # elems to go
1054 my $vd; my $z = '0'x $BASE_LEN;
1055 $x->[scalar @$x] = 0; # avoid || 0 test inside loop
1056 while ($dst < $len)
574bacfe 1057 {
61f5c3f5 1058 $vd = $z.$x->[$src];
1059 $vd = substr($vd,-$BASE_LEN,$BASE_LEN-$rem);
1060 $src++;
1061 $vd = substr($z.$x->[$src],-$rem,$rem) . $vd;
1062 $vd = substr($vd,-$BASE_LEN,$BASE_LEN) if length($vd) > $BASE_LEN;
1063 $x->[$dst] = int($vd);
1064 $dst++;
574bacfe 1065 }
61f5c3f5 1066 splice (@$x,$dst) if $dst > 0; # kill left-over array elems
1067 pop @$x if $x->[-1] == 0; # kill last element if 0
1068 } # else rem == 0
574bacfe 1069 $x;
1070 }
1071
1072sub _lsft
1073 {
1074 my ($c,$x,$y,$n) = @_;
1075
1076 if ($n != 10)
1077 {
61f5c3f5 1078 $n = _new($c,\$n); return _mul($c,$x, _pow($c,$n,$y));
574bacfe 1079 }
61f5c3f5 1080
1081 # shortcut (faster) for shifting by 10) since we are in base 10eX
1082 # multiples of $BASE_LEN:
1083 my $src = scalar @$x; # source
1084 my $len = _num($c,$y); # shift-len as normal int
1085 my $rem = $len % $BASE_LEN; # remainder to shift
1086 my $dst = $src + int($len/$BASE_LEN); # destination
1087 my $vd; # further speedup
1088 $x->[$src] = 0; # avoid first ||0 for speed
1089 my $z = '0' x $BASE_LEN;
1090 while ($src >= 0)
574bacfe 1091 {
61f5c3f5 1092 $vd = $x->[$src]; $vd = $z.$vd;
1093 $vd = substr($vd,-$BASE_LEN+$rem,$BASE_LEN-$rem);
1094 $vd .= $src > 0 ? substr($z.$x->[$src-1],-$BASE_LEN,$rem) : '0' x $rem;
1095 $vd = substr($vd,-$BASE_LEN,$BASE_LEN) if length($vd) > $BASE_LEN;
1096 $x->[$dst] = int($vd);
1097 $dst--; $src--;
574bacfe 1098 }
61f5c3f5 1099 # set lowest parts to 0
1100 while ($dst >= 0) { $x->[$dst--] = 0; }
1101 # fix spurios last zero element
1102 splice @$x,-1 if $x->[-1] == 0;
574bacfe 1103 $x;
1104 }
1105
027dc388 1106sub _pow
1107 {
1108 # power of $x to $y
1109 # ref to array, ref to array, return ref to array
1110 my ($c,$cx,$cy) = @_;
1111
1112 my $pow2 = _one();
1113 my $two = _two();
1114 my $y1 = _copy($c,$cy);
1115 while (!_is_one($c,$y1))
1116 {
1117 _mul($c,$pow2,$cx) if _is_odd($c,$y1);
1118 _div($c,$y1,$two);
1119 _mul($c,$cx,$cx);
1120 }
1121 _mul($c,$cx,$pow2) unless _is_one($c,$pow2);
61f5c3f5 1122 $cx;
027dc388 1123 }
1124
61f5c3f5 1125sub _sqrt1
0716bf9b 1126 {
394e6ffb 1127 # square-root of $x
1128 # ref to array, return ref to array
1129 my ($c,$x) = @_;
0716bf9b 1130
394e6ffb 1131 if (scalar @$x == 1)
1132 {
1133 # fit's into one Perl scalar
1134 $x->[0] = int(sqrt($x->[0]));
1135 return $x;
1136 }
1137 my $y = _copy($c,$x);
61f5c3f5 1138 my $l = _len($c,$x) / 2; # hopefully _len/2 is < $BASE
1139 # my $l2 = [ _len($c,$x) / 2 ]; # old way: hopefully _len/2 is < $BASE
0716bf9b 1140
394e6ffb 1141 splice @$x,0; $x->[0] = 1; # keep ref($x), but modify it
0716bf9b 1142
61f5c3f5 1143 # old way
1144 # _lsft($c,$x,$l2,10);
0716bf9b 1145
61f5c3f5 1146 # construct $x (instead of _lsft($c,$x,$l,10)
1147 my $r = $l % $BASE_LEN; # 10000 00000 00000 00000 ($BASE_LEN=5)
1148 $l = int($l / $BASE_LEN);
1149 $x->[$l--] = int('1' . '0' x $r);
1150 $x->[$l--] = 0 while ($l >= 0);
1151
394e6ffb 1152 my $two = _two();
1153 my $last = _zero();
1154 my $lastlast = _zero();
1155 while (_acmp($c,$last,$x) != 0 && _acmp($c,$lastlast,$x) != 0)
1156 {
1157 $lastlast = _copy($c,$last);
1158 $last = _copy($c,$x);
1159 _add($c,$x, _div($c,_copy($c,$y),$x));
1160 _div($c,$x, $two );
1161 }
1162 _dec($c,$x) if _acmp($c,$y,_mul($c,_copy($c,$x),$x)) < 0; # overshot?
1163 $x;
0716bf9b 1164 }
1165
394e6ffb 1166##############################################################################
1167# binary stuff
0716bf9b 1168
394e6ffb 1169sub _and
1170 {
1171 my ($c,$x,$y) = @_;
0716bf9b 1172
394e6ffb 1173 # the shortcut makes equal, large numbers _really_ fast, and makes only a
1174 # very small performance drop for small numbers (e.g. something with less
1175 # than 32 bit) Since we optimize for large numbers, this is enabled.
1176 return $x if _acmp($c,$x,$y) == 0; # shortcut
0716bf9b 1177
394e6ffb 1178 my $m = _one(); my ($xr,$yr);
1179 my $mask = $AND_MASK;
1180
1181 my $x1 = $x;
1182 my $y1 = _copy($c,$y); # make copy
1183 $x = _zero();
1184 my ($b,$xrr,$yrr);
1185 use integer;
1186 while (!_is_zero($c,$x1) && !_is_zero($c,$y1))
1187 {
1188 ($x1, $xr) = _div($c,$x1,$mask);
1189 ($y1, $yr) = _div($c,$y1,$mask);
1190
1191 # make ints() from $xr, $yr
1192 # this is when the AND_BITS are greater tahn $BASE and is slower for
1193 # small (<256 bits) numbers, but faster for large numbers. Disabled
1194 # due to KISS principle
1195
1196# $b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
1197# $b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
1198# _add($c,$x, _mul($c, _new( $c, \($xrr & $yrr) ), $m) );
1199
61f5c3f5 1200 # 0+ due to '&' doesn't work in strings
1201 _add($c,$x, _mul($c, [ 0+$xr->[0] & 0+$yr->[0] ], $m) );
394e6ffb 1202 _mul($c,$m,$mask);
1203 }
1204 $x;
0716bf9b 1205 }
1206
394e6ffb 1207sub _xor
0716bf9b 1208 {
394e6ffb 1209 my ($c,$x,$y) = @_;
1210
1211 return _zero() if _acmp($c,$x,$y) == 0; # shortcut (see -and)
1212
1213 my $m = _one(); my ($xr,$yr);
1214 my $mask = $XOR_MASK;
1215
1216 my $x1 = $x;
1217 my $y1 = _copy($c,$y); # make copy
1218 $x = _zero();
1219 my ($b,$xrr,$yrr);
1220 use integer;
1221 while (!_is_zero($c,$x1) && !_is_zero($c,$y1))
0716bf9b 1222 {
394e6ffb 1223 ($x1, $xr) = _div($c,$x1,$mask);
1224 ($y1, $yr) = _div($c,$y1,$mask);
1225 # make ints() from $xr, $yr (see _and())
1226 #$b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
1227 #$b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
1228 #_add($c,$x, _mul($c, _new( $c, \($xrr ^ $yrr) ), $m) );
61f5c3f5 1229
1230 # 0+ due to '^' doesn't work in strings
1231 _add($c,$x, _mul($c, [ 0+$xr->[0] ^ 0+$yr->[0] ], $m) );
394e6ffb 1232 _mul($c,$m,$mask);
0716bf9b 1233 }
394e6ffb 1234 # the loop stops when the shorter of the two numbers is exhausted
1235 # the remainder of the longer one will survive bit-by-bit, so we simple
1236 # multiply-add it in
1237 _add($c,$x, _mul($c, $x1, $m) ) if !_is_zero($c,$x1);
1238 _add($c,$x, _mul($c, $y1, $m) ) if !_is_zero($c,$y1);
1239
1240 $x;
0716bf9b 1241 }
1242
394e6ffb 1243sub _or
0716bf9b 1244 {
394e6ffb 1245 my ($c,$x,$y) = @_;
0716bf9b 1246
394e6ffb 1247 return $x if _acmp($c,$x,$y) == 0; # shortcut (see _and)
0716bf9b 1248
394e6ffb 1249 my $m = _one(); my ($xr,$yr);
1250 my $mask = $OR_MASK;
0716bf9b 1251
394e6ffb 1252 my $x1 = $x;
1253 my $y1 = _copy($c,$y); # make copy
1254 $x = _zero();
1255 my ($b,$xrr,$yrr);
1256 use integer;
1257 while (!_is_zero($c,$x1) && !_is_zero($c,$y1))
1258 {
1259 ($x1, $xr) = _div($c,$x1,$mask);
1260 ($y1, $yr) = _div($c,$y1,$mask);
1261 # make ints() from $xr, $yr (see _and())
1262# $b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
1263# $b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
1264# _add($c,$x, _mul($c, _new( $c, \($xrr | $yrr) ), $m) );
1265
61f5c3f5 1266 # 0+ due to '|' doesn't work in strings
1267 _add($c,$x, _mul($c, [ 0+$xr->[0] | 0+$yr->[0] ], $m) );
394e6ffb 1268 _mul($c,$m,$mask);
1269 }
1270 # the loop stops when the shorter of the two numbers is exhausted
1271 # the remainder of the longer one will survive bit-by-bit, so we simple
1272 # multiply-add it in
1273 _add($c,$x, _mul($c, $x1, $m) ) if !_is_zero($c,$x1);
1274 _add($c,$x, _mul($c, $y1, $m) ) if !_is_zero($c,$y1);
1275
1276 $x;
0716bf9b 1277 }
1278
61f5c3f5 1279sub _as_hex
1280 {
1281 # convert a decimal number to hex (ref to array, return ref to string)
1282 my ($c,$x) = @_;
1283
1284 my $x1 = _copy($c,$x);
1285
1286 my $es = '';
1287 my $xr;
1288 my $x10000 = [ 0x10000 ];
1289 while (! _is_zero($c,$x1))
1290 {
1291 ($x1, $xr) = _div($c,$x1,$x10000);
1292 $es .= unpack('h4',pack('v',$xr->[0]));
1293 }
1294 $es = reverse $es;
1295 $es =~ s/^[0]+//; # strip leading zeros
1296 $es = '0x' . $es;
1297 \$es;
1298 }
1299
1300sub _as_bin
1301 {
1302 # convert a decimal number to bin (ref to array, return ref to string)
1303 my ($c,$x) = @_;
1304
1305 my $x1 = _copy($c,$x);
1306
1307 my $es = '';
1308 my $xr;
1309 my $x10000 = [ 0x10000 ];
1310 while (! _is_zero($c,$x1))
1311 {
1312 ($x1, $xr) = _div($c,$x1,$x10000);
1313 $es .= unpack('b16',pack('v',$xr->[0]));
1314 }
1315 $es = reverse $es;
1316 $es =~ s/^[0]+//; # strip leading zeros
1317 $es = '0b' . $es;
1318 \$es;
1319 }
1320
394e6ffb 1321sub _from_hex
0716bf9b 1322 {
394e6ffb 1323 # convert a hex number to decimal (ref to string, return ref to array)
1324 my ($c,$hs) = @_;
0716bf9b 1325
394e6ffb 1326 my $mul = _one();
1327 my $m = [ 0x10000 ]; # 16 bit at a time
1328 my $x = _zero();
0716bf9b 1329
61f5c3f5 1330 my $len = length($$hs)-2;
394e6ffb 1331 $len = int($len/4); # 4-digit parts, w/o '0x'
1332 my $val; my $i = -4;
1333 while ($len >= 0)
1334 {
1335 $val = substr($$hs,$i,4);
1336 $val =~ s/^[+-]?0x// if $len == 0; # for last part only because
1337 $val = hex($val); # hex does not like wrong chars
1338 $i -= 4; $len --;
1339 _add ($c, $x, _mul ($c, [ $val ], $mul ) ) if $val != 0;
1340 _mul ($c, $mul, $m ) if $len >= 0; # skip last mul
1341 }
1342 $x;
1343 }
1344
1345sub _from_bin
0716bf9b 1346 {
394e6ffb 1347 # convert a hex number to decimal (ref to string, return ref to array)
1348 my ($c,$bs) = @_;
0716bf9b 1349
394e6ffb 1350 my $mul = _one();
1351 my $m = [ 0x100 ]; # 8 bit at a time
1352 my $x = _zero();
0716bf9b 1353
61f5c3f5 1354 my $len = length($$bs)-2;
394e6ffb 1355 $len = int($len/8); # 4-digit parts, w/o '0x'
1356 my $val; my $i = -8;
1357 while ($len >= 0)
0716bf9b 1358 {
394e6ffb 1359 $val = substr($$bs,$i,8);
1360 $val =~ s/^[+-]?0b// if $len == 0; # for last part only
1361
394e6ffb 1362 $val = ord(pack('B8',substr('00000000'.$val,-8,8)));
1363
1364 $i -= 8; $len --;
1365 _add ($c, $x, _mul ($c, [ $val ], $mul ) ) if $val != 0;
1366 _mul ($c, $mul, $m ) if $len >= 0; # skip last mul
0716bf9b 1367 }
394e6ffb 1368 $x;
0716bf9b 1369 }
1370
394e6ffb 1371##############################################################################
1372##############################################################################
1373
0716bf9b 13741;
1375__END__
1376
1377=head1 NAME
1378
1379Math::BigInt::Calc - Pure Perl module to support Math::BigInt
1380
1381=head1 SYNOPSIS
1382
ee15d750 1383Provides support for big integer calculations. Not intended to be used by other
1384modules (except Math::BigInt::Cached). Other modules which sport the same
1385functions can also be used to support Math::Bigint, like Math::BigInt::Pari.
0716bf9b 1386
1387=head1 DESCRIPTION
1388
027dc388 1389In order to allow for multiple big integer libraries, Math::BigInt was
1390rewritten to use library modules for core math routines. Any module which
1391follows the same API as this can be used instead by using the following:
0716bf9b 1392
ee15d750 1393 use Math::BigInt lib => 'libname';
0716bf9b 1394
027dc388 1395'libname' is either the long name ('Math::BigInt::Pari'), or only the short
1396version like 'Pari'.
1397
0716bf9b 1398=head1 EXPORT
1399
027dc388 1400The following functions MUST be defined in order to support the use by
1401Math::BigInt:
0716bf9b 1402
1403 _new(string) return ref to new object from ref to decimal string
1404 _zero() return a new object with value 0
1405 _one() return a new object with value 1
1406
1407 _str(obj) return ref to a string representing the object
1408 _num(obj) returns a Perl integer/floating point number
1409 NOTE: because of Perl numeric notation defaults,
1410 the _num'ified obj may lose accuracy due to
1411 machine-dependend floating point size limitations
1412
1413 _add(obj,obj) Simple addition of two objects
1414 _mul(obj,obj) Multiplication of two objects
1415 _div(obj,obj) Division of the 1st object by the 2nd
b22b3e31 1416 In list context, returns (result,remainder).
1417 NOTE: this is integer math, so no
1418 fractional part will be returned.
1419 _sub(obj,obj) Simple subtraction of 1 object from another
0716bf9b 1420 a third, optional parameter indicates that the params
1421 are swapped. In this case, the first param needs to
1422 be preserved, while you can destroy the second.
1423 sub (x,y,1) => return x - y and keep x intact!
e745a66c 1424 _dec(obj) decrement object by one (input is garant. to be > 0)
1425 _inc(obj) increment object by one
1426
0716bf9b 1427
1428 _acmp(obj,obj) <=> operator for objects (return -1, 0 or 1)
1429
1430 _len(obj) returns count of the decimal digits of the object
1431 _digit(obj,n) returns the n'th decimal digit of object
1432
1433 _is_one(obj) return true if argument is +1
1434 _is_zero(obj) return true if argument is 0
1435 _is_even(obj) return true if argument is even (0,2,4,6..)
1436 _is_odd(obj) return true if argument is odd (1,3,5,7..)
1437
1438 _copy return a ref to a true copy of the object
1439
1440 _check(obj) check whether internal representation is still intact
1441 return 0 for ok, otherwise error message as string
1442
bd05a461 1443The following functions are optional, and can be defined if the underlying lib
027dc388 1444has a fast way to do them. If undefined, Math::BigInt will use pure Perl (hence
1445slow) fallback routines to emulate these:
0716bf9b 1446
1447 _from_hex(str) return ref to new object from ref to hexadecimal string
1448 _from_bin(str) return ref to new object from ref to binary string
1449
ee15d750 1450 _as_hex(str) return ref to scalar string containing the value as
1451 unsigned hex string, with the '0x' prepended.
1452 Leading zeros must be stripped.
1453 _as_bin(str) Like as_hex, only as binary string containing only
1454 zeros and ones. Leading zeros must be stripped and a
1455 '0b' must be prepended.
1456
0716bf9b 1457 _rsft(obj,N,B) shift object in base B by N 'digits' right
dccbb853 1458 For unsupported bases B, return undef to signal failure
0716bf9b 1459 _lsft(obj,N,B) shift object in base B by N 'digits' left
dccbb853 1460 For unsupported bases B, return undef to signal failure
0716bf9b 1461
1462 _xor(obj1,obj2) XOR (bit-wise) object 1 with object 2
dccbb853 1463 Note: XOR, AND and OR pad with zeros if size mismatches
0716bf9b 1464 _and(obj1,obj2) AND (bit-wise) object 1 with object 2
1465 _or(obj1,obj2) OR (bit-wise) object 1 with object 2
1466
dccbb853 1467 _mod(obj,obj) Return remainder of div of the 1st by the 2nd object
394e6ffb 1468 _sqrt(obj) return the square root of object (truncate to int)
0716bf9b 1469 _pow(obj,obj) return object 1 to the power of object 2
1470 _gcd(obj,obj) return Greatest Common Divisor of two objects
1471
b22b3e31 1472 _zeros(obj) return number of trailing decimal zeros
0716bf9b 1473
b22b3e31 1474Input strings come in as unsigned but with prefix (i.e. as '123', '0xabc'
0716bf9b 1475or '0b1101').
1476
b22b3e31 1477Testing of input parameter validity is done by the caller, so you need not
574bacfe 1478worry about underflow (f.i. in C<_sub()>, C<_dec()>) nor about division by
1479zero or similar cases.
1480
1481The first parameter can be modified, that includes the possibility that you
1482return a reference to a completely different object instead. Although keeping
dccbb853 1483the reference and just changing it's contents is prefered over creating and
1484returning a different reference.
574bacfe 1485
1486Return values are always references to objects or strings. Exceptions are
1487C<_lsft()> and C<_rsft()>, which return undef if they can not shift the
027dc388 1488argument. This is used to delegate shifting of bases different than the one
1489you can support back to Math::BigInt, which will use some generic code to
1490calculate the result.
574bacfe 1491
1492=head1 WRAP YOUR OWN
1493
1494If you want to port your own favourite c-lib for big numbers to the
1495Math::BigInt interface, you can take any of the already existing modules as
1496a rough guideline. You should really wrap up the latest BigInt and BigFloat
bd05a461 1497testsuites with your module, and replace in them any of the following:
574bacfe 1498
1499 use Math::BigInt;
1500
bd05a461 1501by this:
574bacfe 1502
1503 use Math::BigInt lib => 'yourlib';
1504
1505This way you ensure that your library really works 100% within Math::BigInt.
0716bf9b 1506
1507=head1 LICENSE
1508
1509This program is free software; you may redistribute it and/or modify it under
1510the same terms as Perl itself.
1511
1512=head1 AUTHORS
1513
1514Original math code by Mark Biggar, rewritten by Tels L<http://bloodgate.com/>
1515in late 2000, 2001.
1516Seperated from BigInt and shaped API with the help of John Peacock.
1517
1518=head1 SEE ALSO
1519
ee15d750 1520L<Math::BigInt>, L<Math::BigFloat>, L<Math::BigInt::BitVect>,
1521L<Math::BigInt::GMP>, L<Math::BigInt::Cached> and L<Math::BigInt::Pari>.
0716bf9b 1522
1523=cut
Unnamed repository; edit this file to name it for gitweb.