This is a live mirror of the Perl 5 development currently hosted at https://github.com/perl/perl5
Upgrade Math::BigInt from 1.999701 to 1.999704
[perl5.git] / cpan / Math-BigInt / lib / Math / BigInt / Calc.pm
CommitLineData
0716bf9b
JH
1package Math::BigInt::Calc;
2
0d71d61a 3use 5.006002;
0716bf9b 4use strict;
ada8209b 5# use warnings; # do not use warnings for older Perls
0716bf9b 6
06ce15ad 7our $VERSION = '1.999704';
0716bf9b
JH
8
9# Package to store unsigned big integers in decimal and do math with them
10
11# Internally the numbers are stored in an array with at least 1 element, no
027dc388
JH
12# leading zero parts (except the first) and in base 1eX where X is determined
13# automatically at loading time to be the maximum possible value
0716bf9b
JH
14
15# todo:
a87115f0 16# - fully remove funky $# stuff in div() (maybe - that code scares me...)
0716bf9b
JH
17
18# USE_MUL: due to problems on certain os (os390, posix-bc) "* 1e-5" is used
ee15d750
JH
19# instead of "/ 1e5" at some places, (marked with USE_MUL). Other platforms
20# BS2000, some Crays need USE_DIV instead.
bd05a461
JH
21# The BEGIN block is used to determine which of the two variants gives the
22# correct result.
0716bf9b 23
990fb837 24# Beware of things like:
2ebb273f 25# $i = $i * $y + $car; $car = int($i / $BASE); $i = $i % $BASE;
ada8209b 26# This works on x86, but fails on ARM (SA1100, iPAQ) due to who knows what
990fb837 27# reasons. So, use this instead (slower, but correct):
2ebb273f 28# $i = $i * $y + $car; $car = int($i / $BASE); $i -= $BASE * $car;
990fb837 29
0716bf9b
JH
30##############################################################################
31# global constants, flags and accessory
9b924220 32
2ebb273f 33# announce that we are compatible with MBI v1.83 and up
50109ad0 34sub api_version () { 2; }
0716bf9b
JH
35
36# constants for easier life
2ebb273f 37my ($BASE,$BASE_LEN,$RBASE,$MAX_VAL);
394e6ffb
JH
38my ($AND_BITS,$XOR_BITS,$OR_BITS);
39my ($AND_MASK,$XOR_MASK,$OR_MASK);
ee15d750
JH
40
41sub _base_len
42 {
50109ad0
RGS
43 # Set/get the BASE_LEN and assorted other, connected values.
44 # Used only by the testsuite, the set variant is used only by the BEGIN
45 # block below:
394e6ffb
JH
46 shift;
47
8d8376c6 48 my ($b, $int) = @_;
ee15d750
JH
49 if (defined $b)
50 {
2ebb273f
T
51 # avoid redefinitions
52 undef &_mul;
53 undef &_div;
54
08a3f4a9 55 if ($] >= 5.008 && $int && $b > 7)
2ebb273f
T
56 {
57 $BASE_LEN = $b;
58 *_mul = \&_mul_use_div_64;
59 *_div = \&_div_use_div_64;
60 $BASE = int("1e".$BASE_LEN);
61 $MAX_VAL = $BASE-1;
62 return $BASE_LEN unless wantarray;
f800d0e1 63 return ($BASE_LEN, $BASE, $AND_BITS, $XOR_BITS, $OR_BITS, $BASE_LEN, $MAX_VAL,);
2ebb273f
T
64 }
65
66 # find whether we can use mul or div in mul()/div()
67 $BASE_LEN = $b+1;
61f5c3f5 68 my $caught = 0;
2ebb273f 69 while (--$BASE_LEN > 5)
394e6ffb 70 {
2ebb273f
T
71 $BASE = int("1e".$BASE_LEN);
72 $RBASE = abs('1e-'.$BASE_LEN); # see USE_MUL
394e6ffb 73 $caught = 0;
2ebb273f
T
74 $caught += 1 if (int($BASE * $RBASE) != 1); # should be 1
75 $caught += 2 if (int($BASE / $BASE) != 1); # should be 1
394e6ffb
JH
76 last if $caught != 3;
77 }
ee15d750 78 $BASE = int("1e".$BASE_LEN);
2ebb273f
T
79 $RBASE = abs('1e-'.$BASE_LEN); # see USE_MUL
80 $MAX_VAL = $BASE-1;
b68b7ab1 81
50109ad0
RGS
82 # ($caught & 1) != 0 => cannot use MUL
83 # ($caught & 2) != 0 => cannot use DIV
990fb837 84 if ($caught == 2) # 2
ee15d750 85 {
990fb837 86 # must USE_MUL since we cannot use DIV
2ebb273f
T
87 *_mul = \&_mul_use_mul;
88 *_div = \&_div_use_mul;
ee15d750 89 }
990fb837 90 else # 0 or 1
ee15d750 91 {
ee15d750 92 # can USE_DIV instead
2ebb273f
T
93 *_mul = \&_mul_use_div;
94 *_div = \&_div_use_div;
ee15d750
JH
95 }
96 }
61f5c3f5 97 return $BASE_LEN unless wantarray;
f800d0e1 98 return ($BASE_LEN, $BASE, $AND_BITS, $XOR_BITS, $OR_BITS, $BASE_LEN, $MAX_VAL);
ee15d750 99 }
574bacfe 100
03874afe
T
101sub _new
102 {
103 # (ref to string) return ref to num_array
104 # Convert a number from string format (without sign) to internal base
105 # 1ex format. Assumes normalized value as input.
106 my $il = length($_[1])-1;
107
108 # < BASE_LEN due len-1 above
109 return [ int($_[1]) ] if $il < $BASE_LEN; # shortcut for short numbers
110
111 # this leaves '00000' instead of int 0 and will be corrected after any op
112 [ reverse(unpack("a" . ($il % $BASE_LEN+1)
113 . ("a$BASE_LEN" x ($il / $BASE_LEN)), $_[1])) ];
114 }
115
574bacfe
JH
116BEGIN
117 {
bd05a461 118 # from Daniel Pfeiffer: determine largest group of digits that is precisely
574bacfe 119 # multipliable with itself plus carry
dccbb853
JH
120 # Test now changed to expect the proper pattern, not a result off by 1 or 2
121 my ($e, $num) = 3; # lowest value we will use is 3+1-1 = 3
bd05a461
JH
122 do
123 {
124 $num = ('9' x ++$e) + 0;
394e6ffb 125 $num *= $num + 1.0;
394e6ffb
JH
126 } while ("$num" =~ /9{$e}0{$e}/); # must be a certain pattern
127 $e--; # last test failed, so retract one step
128 # the limits below brush the problems with the test above under the rug:
129 # the test should be able to find the proper $e automatically
130 $e = 5 if $^O =~ /^uts/; # UTS get's some special treatment
131 $e = 5 if $^O =~ /^unicos/; # unicos is also problematic (6 seems to work
132 # there, but we play safe)
2ebb273f 133
8d8376c6
T
134 my $int = 0;
135 if ($e > 7)
136 {
137 use integer;
138 my $e1 = 7;
139 $num = 7;
140 do
141 {
142 $num = ('9' x ++$e1) + 0;
143 $num *= $num + 1;
144 } while ("$num" =~ /9{$e1}0{$e1}/); # must be a certain pattern
145 $e1--; # last test failed, so retract one step
146 if ($e1 > 7)
147 {
148 $int = 1; $e = $e1;
149 }
150 }
151
152 __PACKAGE__->_base_len($e,$int); # set and store
394e6ffb 153
b68b7ab1 154 use integer;
394e6ffb
JH
155 # find out how many bits _and, _or and _xor can take (old default = 16)
156 # I don't think anybody has yet 128 bit scalars, so let's play safe.
394e6ffb 157 local $^W = 0; # don't warn about 'nonportable number'
c38b2de2 158 $AND_BITS = 15; $XOR_BITS = 15; $OR_BITS = 15;
394e6ffb
JH
159
160 # find max bits, we will not go higher than numberofbits that fit into $BASE
161 # to make _and etc simpler (and faster for smaller, slower for large numbers)
162 my $max = 16;
163 while (2 ** $max < $BASE) { $max++; }
1ddff52a
JH
164 {
165 no integer;
166 $max = 16 if $] < 5.006; # older Perls might not take >16 too well
167 }
394e6ffb
JH
168 my ($x,$y,$z);
169 do {
170 $AND_BITS++;
0617f807 171 $x = CORE::oct('0b' . '1' x $AND_BITS); $y = $x & $x;
394e6ffb
JH
172 $z = (2 ** $AND_BITS) - 1;
173 } while ($AND_BITS < $max && $x == $z && $y == $x);
174 $AND_BITS --; # retreat one step
175 do {
176 $XOR_BITS++;
0617f807 177 $x = CORE::oct('0b' . '1' x $XOR_BITS); $y = $x ^ 0;
394e6ffb
JH
178 $z = (2 ** $XOR_BITS) - 1;
179 } while ($XOR_BITS < $max && $x == $z && $y == $x);
180 $XOR_BITS --; # retreat one step
181 do {
182 $OR_BITS++;
0617f807 183 $x = CORE::oct('0b' . '1' x $OR_BITS); $y = $x | $x;
394e6ffb
JH
184 $z = (2 ** $OR_BITS) - 1;
185 } while ($OR_BITS < $max && $x == $z && $y == $x);
186 $OR_BITS --; # retreat one step
187
9b924220
RGS
188 $AND_MASK = __PACKAGE__->_new( ( 2 ** $AND_BITS ));
189 $XOR_MASK = __PACKAGE__->_new( ( 2 ** $XOR_BITS ));
190 $OR_MASK = __PACKAGE__->_new( ( 2 ** $OR_BITS ));
50109ad0 191
9681bfa6 192 # We can compute the approximate length no faster than the real length:
50109ad0 193 *_alen = \&_len;
394e6ffb 194 }
0716bf9b 195
03874afe
T
196###############################################################################
197
0716bf9b
JH
198sub _zero
199 {
200 # create a zero
61f5c3f5 201 [ 0 ];
0716bf9b
JH
202 }
203
204sub _one
205 {
206 # create a one
61f5c3f5 207 [ 1 ];
0716bf9b
JH
208 }
209
027dc388
JH
210sub _two
211 {
1ddff52a 212 # create a two (used internally for shifting)
61f5c3f5 213 [ 2 ];
027dc388
JH
214 }
215
9b924220
RGS
216sub _ten
217 {
218 # create a 10 (used internally for shifting)
219 [ 10 ];
220 }
221
50109ad0
RGS
222sub _1ex
223 {
224 # create a 1Ex
225 my $rem = $_[1] % $BASE_LEN; # remainder
226 my $parts = $_[1] / $BASE_LEN; # parts
227
228 # 000000, 000000, 100
229 [ (0) x $parts, '1' . ('0' x $rem) ];
230 }
231
0716bf9b
JH
232sub _copy
233 {
091c87b1 234 # make a true copy
61f5c3f5 235 [ @{$_[1]} ];
0716bf9b
JH
236 }
237
bd05a461
JH
238# catch and throw away
239sub import { }
240
0716bf9b
JH
241##############################################################################
242# convert back to string and number
243
244sub _str
245 {
246 # (ref to BINT) return num_str
247 # Convert number from internal base 100000 format to string format.
248 # internal format is always normalized (no leading zeros, "-0" => "+0")
574bacfe 249 my $ar = $_[1];
61f5c3f5 250
b68b7ab1
T
251 my $l = scalar @$ar; # number of parts
252 if ($l < 1) # should not happen
253 {
254 require Carp;
255 Carp::croak("$_[1] has no elements");
256 }
61f5c3f5 257
b68b7ab1 258 my $ret = "";
0716bf9b
JH
259 # handle first one different to strip leading zeros from it (there are no
260 # leading zero parts in internal representation)
61f5c3f5 261 $l --; $ret .= int($ar->[$l]); $l--;
0716bf9b 262 # Interestingly, the pre-padd method uses more time
091c87b1 263 # the old grep variant takes longer (14 vs. 10 sec)
574bacfe 264 my $z = '0' x ($BASE_LEN-1);
0716bf9b
JH
265 while ($l >= 0)
266 {
574bacfe 267 $ret .= substr($z.$ar->[$l],-$BASE_LEN); # fastest way I could think of
0716bf9b
JH
268 $l--;
269 }
9b924220 270 $ret;
0716bf9b
JH
271 }
272
273sub _num
274 {
e5dcd524
PJA
275 # Make a Perl scalar number (int/float) from a BigInt object.
276 my $x = $_[1];
9b924220 277
e5dcd524
PJA
278 return 0 + $x->[0] if scalar @$x == 1; # below $BASE
279
280 # Start with the most significant element and work towards the least
281 # significant element. Avoid multiplying "inf" (which happens if the number
282 # overflows) with "0" (if there are zero elements in $x) since this gives
283 # "nan" which propagates to the output.
284
285 my $num = 0;
286 for (my $i = $#$x ; $i >= 0 ; --$i) {
287 $num *= $BASE;
288 $num += $x -> [$i];
0716bf9b 289 }
e5dcd524 290 return $num;
0716bf9b
JH
291 }
292
293##############################################################################
294# actual math code
295
296sub _add
297 {
298 # (ref to int_num_array, ref to int_num_array)
574bacfe 299 # routine to add two base 1eX numbers
0716bf9b 300 # stolen from Knuth Vol 2 Algorithm A pg 231
b22b3e31 301 # there are separate routines to add and sub as per Knuth pg 233
db2e1fb3 302 # This routine modifies array x, but not y.
0716bf9b 303
574bacfe 304 my ($c,$x,$y) = @_;
b3abae2a
JH
305
306 return $x if (@$y == 1) && $y->[0] == 0; # $x + 0 => $x
307 if ((@$x == 1) && $x->[0] == 0) # 0 + $y => $y->copy
308 {
7b29e1e6 309 # twice as slow as $x = [ @$y ], but nec. to retain $x as ref :(
b3abae2a
JH
310 @$x = @$y; return $x;
311 }
0716bf9b
JH
312
313 # for each in Y, add Y to X and carry. If after that, something is left in
314 # X, foreach in X add carry to X and then return X, carry
091c87b1 315 # Trades one "$j++" for having to shift arrays
0716bf9b
JH
316 my $i; my $car = 0; my $j = 0;
317 for $i (@$y)
318 {
e745a66c 319 $x->[$j] -= $BASE if $car = (($x->[$j] += $i + $car) >= $BASE) ? 1 : 0;
0716bf9b
JH
320 $j++;
321 }
322 while ($car != 0)
323 {
324 $x->[$j] -= $BASE if $car = (($x->[$j] += $car) >= $BASE) ? 1 : 0; $j++;
325 }
61f5c3f5 326 $x;
e745a66c
JH
327 }
328
329sub _inc
330 {
331 # (ref to int_num_array, ref to int_num_array)
091c87b1 332 # Add 1 to $x, modify $x in place
e745a66c
JH
333 my ($c,$x) = @_;
334
335 for my $i (@$x)
336 {
337 return $x if (($i += 1) < $BASE); # early out
61f5c3f5 338 $i = 0; # overflow, next
e745a66c 339 }
ae161977 340 push @$x,1 if (($x->[-1] || 0) == 0); # last overflowed, so extend
61f5c3f5 341 $x;
e745a66c
JH
342 }
343
344sub _dec
345 {
346 # (ref to int_num_array, ref to int_num_array)
091c87b1 347 # Sub 1 from $x, modify $x in place
e745a66c
JH
348 my ($c,$x) = @_;
349
2ebb273f 350 my $MAX = $BASE-1; # since MAX_VAL based on BASE
e745a66c
JH
351 for my $i (@$x)
352 {
353 last if (($i -= 1) >= 0); # early out
091c87b1 354 $i = $MAX; # underflow, next
e745a66c 355 }
091c87b1 356 pop @$x if $x->[-1] == 0 && @$x > 1; # last underflowed (but leave 0)
61f5c3f5 357 $x;
0716bf9b
JH
358 }
359
360sub _sub
361 {
9393ace2 362 # (ref to int_num_array, ref to int_num_array, swap)
574bacfe 363 # subtract base 1eX numbers -- stolen from Knuth Vol 2 pg 232, $x > $y
56b9c951 364 # subtract Y from X by modifying x in place
574bacfe 365 my ($c,$sx,$sy,$s) = @_;
0716bf9b
JH
366
367 my $car = 0; my $i; my $j = 0;
368 if (!$s)
369 {
0716bf9b
JH
370 for $i (@$sx)
371 {
372 last unless defined $sy->[$j] || $car;
0716bf9b 373 $i += $BASE if $car = (($i -= ($sy->[$j] || 0) + $car) < 0); $j++;
0716bf9b
JH
374 }
375 # might leave leading zeros, so fix that
394e6ffb 376 return __strip_zeros($sx);
0716bf9b 377 }
394e6ffb 378 for $i (@$sx)
0716bf9b 379 {
07d34614 380 # we can't do an early out if $x is < than $y, since we
56b9c951
JH
381 # need to copy the high chunks from $y. Found by Bob Mathews.
382 #last unless defined $sy->[$j] || $car;
394e6ffb
JH
383 $sy->[$j] += $BASE
384 if $car = (($sy->[$j] = $i-($sy->[$j]||0) - $car) < 0);
385 $j++;
0716bf9b 386 }
394e6ffb
JH
387 # might leave leading zeros, so fix that
388 __strip_zeros($sy);
0716bf9b
JH
389 }
390
ee15d750 391sub _mul_use_mul
0716bf9b 392 {
9393ace2 393 # (ref to int_num_array, ref to int_num_array)
0716bf9b 394 # multiply two numbers in internal representation
b22b3e31 395 # modifies first arg, second need not be different from first
574bacfe 396 my ($c,$xv,$yv) = @_;
dccbb853 397
990fb837 398 if (@$yv == 1)
b3abae2a 399 {
990fb837
RGS
400 # shortcut for two very short numbers (improved by Nathan Zook)
401 # works also if xv and yv are the same reference, and handles also $x == 0
402 if (@$xv == 1)
403 {
2ebb273f 404 if (($xv->[0] *= $yv->[0]) >= $BASE)
990fb837 405 {
2ebb273f 406 $xv->[0] = $xv->[0] - ($xv->[1] = int($xv->[0] * $RBASE)) * $BASE;
990fb837
RGS
407 };
408 return $xv;
409 }
410 # $x * 0 => 0
411 if ($yv->[0] == 0)
412 {
413 @$xv = (0);
414 return $xv;
415 }
416 # multiply a large number a by a single element one, so speed up
417 my $y = $yv->[0]; my $car = 0;
418 foreach my $i (@$xv)
419 {
2ebb273f 420 $i = $i * $y + $car; $car = int($i * $RBASE); $i -= $car * $BASE;
990fb837
RGS
421 }
422 push @$xv, $car if $car != 0;
b3abae2a
JH
423 return $xv;
424 }
990fb837
RGS
425 # shortcut for result $x == 0 => result = 0
426 return $xv if ( ((@$xv == 1) && ($xv->[0] == 0)) );
b3abae2a 427
0716bf9b 428 # since multiplying $x with $x fails, make copy in this case
d614cd8b 429 $yv = [@$xv] if $xv == $yv; # same references?
9393ace2 430
61f5c3f5
T
431 my @prod = (); my ($prod,$car,$cty,$xi,$yi);
432
0716bf9b
JH
433 for $xi (@$xv)
434 {
435 $car = 0; $cty = 0;
574bacfe
JH
436
437 # slow variant
438# for $yi (@$yv)
439# {
440# $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
441# $prod[$cty++] =
2ebb273f 442# $prod - ($car = int($prod * RBASE)) * $BASE; # see USE_MUL
574bacfe
JH
443# }
444# $prod[$cty] += $car if $car; # need really to check for 0?
445# $xi = shift @prod;
446
447 # faster variant
448 # looping through this if $xi == 0 is silly - so optimize it away!
449 $xi = (shift @prod || 0), next if $xi == 0;
0716bf9b
JH
450 for $yi (@$yv)
451 {
452 $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
574bacfe
JH
453## this is actually a tad slower
454## $prod = $prod[$cty]; $prod += ($car + $xi * $yi); # no ||0 here
0716bf9b 455 $prod[$cty++] =
2ebb273f 456 $prod - ($car = int($prod * $RBASE)) * $BASE; # see USE_MUL
0716bf9b
JH
457 }
458 $prod[$cty] += $car if $car; # need really to check for 0?
027dc388 459 $xi = shift @prod || 0; # || 0 makes v5.005_3 happy
0716bf9b 460 }
0716bf9b 461 push @$xv, @prod;
50109ad0
RGS
462 # can't have leading zeros
463# __strip_zeros($xv);
61f5c3f5 464 $xv;
0716bf9b
JH
465 }
466
2ebb273f
T
467sub _mul_use_div_64
468 {
469 # (ref to int_num_array, ref to int_num_array)
470 # multiply two numbers in internal representation
471 # modifies first arg, second need not be different from first
472 # works for 64 bit integer with "use integer"
473 my ($c,$xv,$yv) = @_;
474
475 use integer;
476 if (@$yv == 1)
477 {
478 # shortcut for two small numbers, also handles $x == 0
479 if (@$xv == 1)
480 {
481 # shortcut for two very short numbers (improved by Nathan Zook)
482 # works also if xv and yv are the same reference, and handles also $x == 0
483 if (($xv->[0] *= $yv->[0]) >= $BASE)
484 {
485 $xv->[0] =
486 $xv->[0] - ($xv->[1] = $xv->[0] / $BASE) * $BASE;
487 };
488 return $xv;
489 }
490 # $x * 0 => 0
491 if ($yv->[0] == 0)
492 {
493 @$xv = (0);
494 return $xv;
495 }
496 # multiply a large number a by a single element one, so speed up
497 my $y = $yv->[0]; my $car = 0;
498 foreach my $i (@$xv)
499 {
500 #$i = $i * $y + $car; $car = $i / $BASE; $i -= $car * $BASE;
501 $i = $i * $y + $car; $i -= ($car = $i / $BASE) * $BASE;
502 }
503 push @$xv, $car if $car != 0;
504 return $xv;
505 }
506 # shortcut for result $x == 0 => result = 0
507 return $xv if ( ((@$xv == 1) && ($xv->[0] == 0)) );
508
509 # since multiplying $x with $x fails, make copy in this case
510 $yv = [@$xv] if $xv == $yv; # same references?
511
512 my @prod = (); my ($prod,$car,$cty,$xi,$yi);
513 for $xi (@$xv)
514 {
515 $car = 0; $cty = 0;
516 # looping through this if $xi == 0 is silly - so optimize it away!
517 $xi = (shift @prod || 0), next if $xi == 0;
518 for $yi (@$yv)
519 {
520 $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
521 $prod[$cty++] = $prod - ($car = $prod / $BASE) * $BASE;
522 }
523 $prod[$cty] += $car if $car; # need really to check for 0?
524 $xi = shift @prod || 0; # || 0 makes v5.005_3 happy
525 }
526 push @$xv, @prod;
527 $xv;
528 }
529
ee15d750
JH
530sub _mul_use_div
531 {
9393ace2 532 # (ref to int_num_array, ref to int_num_array)
ee15d750
JH
533 # multiply two numbers in internal representation
534 # modifies first arg, second need not be different from first
535 my ($c,$xv,$yv) = @_;
50109ad0 536
990fb837 537 if (@$yv == 1)
b3abae2a 538 {
990fb837
RGS
539 # shortcut for two small numbers, also handles $x == 0
540 if (@$xv == 1)
541 {
542 # shortcut for two very short numbers (improved by Nathan Zook)
543 # works also if xv and yv are the same reference, and handles also $x == 0
2ebb273f 544 if (($xv->[0] *= $yv->[0]) >= $BASE)
990fb837
RGS
545 {
546 $xv->[0] =
2ebb273f 547 $xv->[0] - ($xv->[1] = int($xv->[0] / $BASE)) * $BASE;
990fb837
RGS
548 };
549 return $xv;
550 }
551 # $x * 0 => 0
552 if ($yv->[0] == 0)
553 {
554 @$xv = (0);
555 return $xv;
556 }
557 # multiply a large number a by a single element one, so speed up
558 my $y = $yv->[0]; my $car = 0;
559 foreach my $i (@$xv)
560 {
2ebb273f
T
561 $i = $i * $y + $car; $car = int($i / $BASE); $i -= $car * $BASE;
562 # This (together with use integer;) does not work on 32-bit Perls
563 #$i = $i * $y + $car; $i -= ($car = $i / $BASE) * $BASE;
990fb837
RGS
564 }
565 push @$xv, $car if $car != 0;
b3abae2a
JH
566 return $xv;
567 }
990fb837
RGS
568 # shortcut for result $x == 0 => result = 0
569 return $xv if ( ((@$xv == 1) && ($xv->[0] == 0)) );
b3abae2a 570
ee15d750 571 # since multiplying $x with $x fails, make copy in this case
d614cd8b 572 $yv = [@$xv] if $xv == $yv; # same references?
9393ace2 573
61f5c3f5 574 my @prod = (); my ($prod,$car,$cty,$xi,$yi);
ee15d750
JH
575 for $xi (@$xv)
576 {
577 $car = 0; $cty = 0;
578 # looping through this if $xi == 0 is silly - so optimize it away!
579 $xi = (shift @prod || 0), next if $xi == 0;
580 for $yi (@$yv)
581 {
582 $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
2ebb273f 583 $prod[$cty++] = $prod - ($car = int($prod / $BASE)) * $BASE;
ee15d750
JH
584 }
585 $prod[$cty] += $car if $car; # need really to check for 0?
027dc388 586 $xi = shift @prod || 0; # || 0 makes v5.005_3 happy
ee15d750
JH
587 }
588 push @$xv, @prod;
50109ad0
RGS
589 # can't have leading zeros
590# __strip_zeros($xv);
61f5c3f5 591 $xv;
ee15d750
JH
592 }
593
594sub _div_use_mul
0716bf9b 595 {
b22b3e31 596 # ref to array, ref to array, modify first array and return remainder if
0716bf9b 597 # in list context
aef458a0
JH
598
599 # see comments in _div_use_div() for more explanations
600
574bacfe 601 my ($c,$x,$yorg) = @_;
aef458a0 602
9681bfa6 603 # the general div algorithm here is about O(N*N) and thus quite slow, so
aef458a0 604 # we first check for some special cases and use shortcuts to handle them.
0716bf9b 605
aef458a0
JH
606 # This works, because we store the numbers in a chunked format where each
607 # element contains 5..7 digits (depending on system).
608
609 # if both numbers have only one element:
61f5c3f5
T
610 if (@$x == 1 && @$yorg == 1)
611 {
13a12e00 612 # shortcut, $yorg and $x are two small numbers
61f5c3f5
T
613 if (wantarray)
614 {
615 my $r = [ $x->[0] % $yorg->[0] ];
616 $x->[0] = int($x->[0] / $yorg->[0]);
617 return ($x,$r);
618 }
619 else
620 {
621 $x->[0] = int($x->[0] / $yorg->[0]);
622 return $x;
623 }
624 }
aef458a0
JH
625
626 # if x has more than one, but y has only one element:
28df3e88
JH
627 if (@$yorg == 1)
628 {
629 my $rem;
630 $rem = _mod($c,[ @$x ],$yorg) if wantarray;
13a12e00 631
28df3e88
JH
632 # shortcut, $y is < $BASE
633 my $j = scalar @$x; my $r = 0;
634 my $y = $yorg->[0]; my $b;
635 while ($j-- > 0)
636 {
2ebb273f 637 $b = $r * $BASE + $x->[$j];
28df3e88
JH
638 $x->[$j] = int($b/$y);
639 $r = $b % $y;
640 }
641 pop @$x if @$x > 1 && $x->[-1] == 0; # splice up a leading zero
642 return ($x,$rem) if wantarray;
643 return $x;
644 }
0716bf9b 645
aef458a0
JH
646 # now x and y have more than one element
647
648 # check whether y has more elements than x, if yet, the result will be 0
649 if (@$yorg > @$x)
650 {
651 my $rem;
652 $rem = [@$x] if wantarray; # make copy
653 splice (@$x,1); # keep ref to original array
654 $x->[0] = 0; # set to 0
655 return ($x,$rem) if wantarray; # including remainder?
656 return $x; # only x, which is [0] now
657 }
658 # check whether the numbers have the same number of elements, in that case
659 # the result will fit into one element and can be computed efficiently
660 if (@$yorg == @$x)
661 {
662 my $rem;
663 # if $yorg has more digits than $x (it's leading element is longer than
664 # the one from $x), the result will also be 0:
665 if (length(int($yorg->[-1])) > length(int($x->[-1])))
666 {
667 $rem = [@$x] if wantarray; # make copy
668 splice (@$x,1); # keep ref to org array
669 $x->[0] = 0; # set to 0
670 return ($x,$rem) if wantarray; # including remainder?
671 return $x;
672 }
673 # now calculate $x / $yorg
674 if (length(int($yorg->[-1])) == length(int($x->[-1])))
675 {
b68b7ab1
T
676 # same length, so make full compare
677
aef458a0
JH
678 my $a = 0; my $j = scalar @$x - 1;
679 # manual way (abort if unequal, good for early ne)
680 while ($j >= 0)
681 {
682 last if ($a = $x->[$j] - $yorg->[$j]); $j--;
683 }
684 # $a contains the result of the compare between X and Y
b68b7ab1 685 # a < 0: x < y, a == 0: x == y, a > 0: x > y
aef458a0
JH
686 if ($a <= 0)
687 {
b68b7ab1
T
688 $rem = [ 0 ]; # a = 0 => x == y => rem 0
689 $rem = [@$x] if $a != 0; # a < 0 => x < y => rem = x
690 splice(@$x,1); # keep single element
691 $x->[0] = 0; # if $a < 0
692 $x->[0] = 1 if $a == 0; # $x == $y
7596a890
RGS
693 return ($x,$rem) if wantarray;
694 return $x;
aef458a0 695 }
b68b7ab1 696 # $x >= $y, so proceed normally
aef458a0
JH
697 }
698 }
699
700 # all other cases:
701
d614cd8b 702 my $y = [ @$yorg ]; # always make copy to preserve
61f5c3f5
T
703
704 my ($car,$bar,$prd,$dd,$xi,$yi,@q,$v2,$v1,@d,$tmp,$q,$u2,$u1,$u0);
705
706 $car = $bar = $prd = 0;
2ebb273f 707 if (($dd = int($BASE/($y->[-1]+1))) != 1)
0716bf9b
JH
708 {
709 for $xi (@$x)
710 {
711 $xi = $xi * $dd + $car;
2ebb273f 712 $xi -= ($car = int($xi * $RBASE)) * $BASE; # see USE_MUL
0716bf9b
JH
713 }
714 push(@$x, $car); $car = 0;
715 for $yi (@$y)
716 {
717 $yi = $yi * $dd + $car;
2ebb273f 718 $yi -= ($car = int($yi * $RBASE)) * $BASE; # see USE_MUL
0716bf9b
JH
719 }
720 }
721 else
722 {
723 push(@$x, 0);
724 }
725 @q = (); ($v2,$v1) = @$y[-2,-1];
726 $v2 = 0 unless $v2;
727 while ($#$x > $#$y)
728 {
729 ($u2,$u1,$u0) = @$x[-3..-1];
730 $u2 = 0 unless $u2;
731 #warn "oups v1 is 0, u0: $u0 $y->[-2] $y->[-1] l ",scalar @$y,"\n"
732 # if $v1 == 0;
2ebb273f
T
733 $q = (($u0 == $v1) ? $MAX_VAL : int(($u0*$BASE+$u1)/$v1));
734 --$q while ($v2*$q > ($u0*$BASE+$u1-$q*$v1)*$BASE+$u2);
0716bf9b
JH
735 if ($q)
736 {
737 ($car, $bar) = (0,0);
738 for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
739 {
740 $prd = $q * $y->[$yi] + $car;
2ebb273f
T
741 $prd -= ($car = int($prd * $RBASE)) * $BASE; # see USE_MUL
742 $x->[$xi] += $BASE if ($bar = (($x->[$xi] -= $prd + $bar) < 0));
0716bf9b
JH
743 }
744 if ($x->[-1] < $car + $bar)
745 {
746 $car = 0; --$q;
747 for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
748 {
2ebb273f
T
749 $x->[$xi] -= $BASE
750 if ($car = (($x->[$xi] += $y->[$yi] + $car) >= $BASE));
0716bf9b
JH
751 }
752 }
753 }
aef458a0
JH
754 pop(@$x);
755 unshift(@q, $q);
0716bf9b
JH
756 }
757 if (wantarray)
758 {
759 @d = ();
760 if ($dd != 1)
761 {
762 $car = 0;
763 for $xi (reverse @$x)
764 {
2ebb273f 765 $prd = $car * $BASE + $xi;
0716bf9b
JH
766 $car = $prd - ($tmp = int($prd / $dd)) * $dd; # see USE_MUL
767 unshift(@d, $tmp);
768 }
769 }
770 else
771 {
772 @d = @$x;
773 }
774 @$x = @q;
61f5c3f5 775 my $d = \@d;
990fb837
RGS
776 __strip_zeros($x);
777 __strip_zeros($d);
61f5c3f5 778 return ($x,$d);
0716bf9b
JH
779 }
780 @$x = @q;
990fb837 781 __strip_zeros($x);
61f5c3f5 782 $x;
0716bf9b
JH
783 }
784
2ebb273f
T
785sub _div_use_div_64
786 {
787 # ref to array, ref to array, modify first array and return remainder if
788 # in list context
789 # This version works on 64 bit integers
790 my ($c,$x,$yorg) = @_;
791
792 use integer;
9681bfa6 793 # the general div algorithm here is about O(N*N) and thus quite slow, so
2ebb273f
T
794 # we first check for some special cases and use shortcuts to handle them.
795
796 # This works, because we store the numbers in a chunked format where each
797 # element contains 5..7 digits (depending on system).
798
799 # if both numbers have only one element:
800 if (@$x == 1 && @$yorg == 1)
801 {
802 # shortcut, $yorg and $x are two small numbers
803 if (wantarray)
804 {
805 my $r = [ $x->[0] % $yorg->[0] ];
806 $x->[0] = int($x->[0] / $yorg->[0]);
807 return ($x,$r);
808 }
809 else
810 {
811 $x->[0] = int($x->[0] / $yorg->[0]);
812 return $x;
813 }
814 }
815 # if x has more than one, but y has only one element:
816 if (@$yorg == 1)
817 {
818 my $rem;
819 $rem = _mod($c,[ @$x ],$yorg) if wantarray;
820
821 # shortcut, $y is < $BASE
822 my $j = scalar @$x; my $r = 0;
823 my $y = $yorg->[0]; my $b;
824 while ($j-- > 0)
825 {
826 $b = $r * $BASE + $x->[$j];
827 $x->[$j] = int($b/$y);
828 $r = $b % $y;
829 }
830 pop @$x if @$x > 1 && $x->[-1] == 0; # splice up a leading zero
831 return ($x,$rem) if wantarray;
832 return $x;
833 }
834 # now x and y have more than one element
835
836 # check whether y has more elements than x, if yet, the result will be 0
837 if (@$yorg > @$x)
838 {
839 my $rem;
840 $rem = [@$x] if wantarray; # make copy
841 splice (@$x,1); # keep ref to original array
842 $x->[0] = 0; # set to 0
843 return ($x,$rem) if wantarray; # including remainder?
844 return $x; # only x, which is [0] now
845 }
846 # check whether the numbers have the same number of elements, in that case
847 # the result will fit into one element and can be computed efficiently
848 if (@$yorg == @$x)
849 {
850 my $rem;
851 # if $yorg has more digits than $x (it's leading element is longer than
852 # the one from $x), the result will also be 0:
853 if (length(int($yorg->[-1])) > length(int($x->[-1])))
854 {
855 $rem = [@$x] if wantarray; # make copy
856 splice (@$x,1); # keep ref to org array
857 $x->[0] = 0; # set to 0
858 return ($x,$rem) if wantarray; # including remainder?
859 return $x;
860 }
861 # now calculate $x / $yorg
862
863 if (length(int($yorg->[-1])) == length(int($x->[-1])))
864 {
865 # same length, so make full compare
866
867 my $a = 0; my $j = scalar @$x - 1;
868 # manual way (abort if unequal, good for early ne)
869 while ($j >= 0)
870 {
871 last if ($a = $x->[$j] - $yorg->[$j]); $j--;
872 }
873 # $a contains the result of the compare between X and Y
874 # a < 0: x < y, a == 0: x == y, a > 0: x > y
875 if ($a <= 0)
876 {
877 $rem = [ 0 ]; # a = 0 => x == y => rem 0
878 $rem = [@$x] if $a != 0; # a < 0 => x < y => rem = x
879 splice(@$x,1); # keep single element
880 $x->[0] = 0; # if $a < 0
881 $x->[0] = 1 if $a == 0; # $x == $y
882 return ($x,$rem) if wantarray; # including remainder?
883 return $x;
884 }
885 # $x >= $y, so proceed normally
886
887 }
888 }
889
890 # all other cases:
891
892 my $y = [ @$yorg ]; # always make copy to preserve
893
894 my ($car,$bar,$prd,$dd,$xi,$yi,@q,$v2,$v1,@d,$tmp,$q,$u2,$u1,$u0);
895
896 $car = $bar = $prd = 0;
897 if (($dd = int($BASE/($y->[-1]+1))) != 1)
898 {
899 for $xi (@$x)
900 {
901 $xi = $xi * $dd + $car;
902 $xi -= ($car = int($xi / $BASE)) * $BASE;
903 }
904 push(@$x, $car); $car = 0;
905 for $yi (@$y)
906 {
907 $yi = $yi * $dd + $car;
908 $yi -= ($car = int($yi / $BASE)) * $BASE;
909 }
910 }
911 else
912 {
913 push(@$x, 0);
914 }
915
916 # @q will accumulate the final result, $q contains the current computed
917 # part of the final result
918
919 @q = (); ($v2,$v1) = @$y[-2,-1];
920 $v2 = 0 unless $v2;
921 while ($#$x > $#$y)
922 {
923 ($u2,$u1,$u0) = @$x[-3..-1];
924 $u2 = 0 unless $u2;
925 #warn "oups v1 is 0, u0: $u0 $y->[-2] $y->[-1] l ",scalar @$y,"\n"
926 # if $v1 == 0;
927 $q = (($u0 == $v1) ? $MAX_VAL : int(($u0*$BASE+$u1)/$v1));
928 --$q while ($v2*$q > ($u0*$BASE+$u1-$q*$v1)*$BASE+$u2);
929 if ($q)
930 {
931 ($car, $bar) = (0,0);
932 for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
933 {
934 $prd = $q * $y->[$yi] + $car;
935 $prd -= ($car = int($prd / $BASE)) * $BASE;
936 $x->[$xi] += $BASE if ($bar = (($x->[$xi] -= $prd + $bar) < 0));
937 }
938 if ($x->[-1] < $car + $bar)
939 {
940 $car = 0; --$q;
941 for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
942 {
943 $x->[$xi] -= $BASE
944 if ($car = (($x->[$xi] += $y->[$yi] + $car) >= $BASE));
945 }
946 }
947 }
948 pop(@$x); unshift(@q, $q);
949 }
950 if (wantarray)
951 {
952 @d = ();
953 if ($dd != 1)
954 {
955 $car = 0;
956 for $xi (reverse @$x)
957 {
958 $prd = $car * $BASE + $xi;
959 $car = $prd - ($tmp = int($prd / $dd)) * $dd;
960 unshift(@d, $tmp);
961 }
962 }
963 else
964 {
965 @d = @$x;
966 }
967 @$x = @q;
968 my $d = \@d;
969 __strip_zeros($x);
970 __strip_zeros($d);
971 return ($x,$d);
972 }
973 @$x = @q;
974 __strip_zeros($x);
975 $x;
976 }
977
ee15d750
JH
978sub _div_use_div
979 {
980 # ref to array, ref to array, modify first array and return remainder if
981 # in list context
ee15d750 982 my ($c,$x,$yorg) = @_;
ee15d750 983
9681bfa6 984 # the general div algorithm here is about O(N*N) and thus quite slow, so
990fb837
RGS
985 # we first check for some special cases and use shortcuts to handle them.
986
987 # This works, because we store the numbers in a chunked format where each
988 # element contains 5..7 digits (depending on system).
989
990 # if both numbers have only one element:
61f5c3f5
T
991 if (@$x == 1 && @$yorg == 1)
992 {
13a12e00 993 # shortcut, $yorg and $x are two small numbers
61f5c3f5
T
994 if (wantarray)
995 {
996 my $r = [ $x->[0] % $yorg->[0] ];
997 $x->[0] = int($x->[0] / $yorg->[0]);
998 return ($x,$r);
999 }
1000 else
1001 {
1002 $x->[0] = int($x->[0] / $yorg->[0]);
1003 return $x;
1004 }
1005 }
990fb837 1006 # if x has more than one, but y has only one element:
28df3e88
JH
1007 if (@$yorg == 1)
1008 {
1009 my $rem;
1010 $rem = _mod($c,[ @$x ],$yorg) if wantarray;
1011
1012 # shortcut, $y is < $BASE
1013 my $j = scalar @$x; my $r = 0;
1014 my $y = $yorg->[0]; my $b;
1015 while ($j-- > 0)
1016 {
2ebb273f 1017 $b = $r * $BASE + $x->[$j];
28df3e88
JH
1018 $x->[$j] = int($b/$y);
1019 $r = $b % $y;
1020 }
1021 pop @$x if @$x > 1 && $x->[-1] == 0; # splice up a leading zero
1022 return ($x,$rem) if wantarray;
1023 return $x;
1024 }
990fb837 1025 # now x and y have more than one element
ee15d750 1026
990fb837
RGS
1027 # check whether y has more elements than x, if yet, the result will be 0
1028 if (@$yorg > @$x)
61f5c3f5 1029 {
990fb837
RGS
1030 my $rem;
1031 $rem = [@$x] if wantarray; # make copy
1032 splice (@$x,1); # keep ref to original array
1033 $x->[0] = 0; # set to 0
1034 return ($x,$rem) if wantarray; # including remainder?
aef458a0 1035 return $x; # only x, which is [0] now
61f5c3f5 1036 }
990fb837
RGS
1037 # check whether the numbers have the same number of elements, in that case
1038 # the result will fit into one element and can be computed efficiently
1039 if (@$yorg == @$x)
1040 {
1041 my $rem;
1042 # if $yorg has more digits than $x (it's leading element is longer than
1043 # the one from $x), the result will also be 0:
1044 if (length(int($yorg->[-1])) > length(int($x->[-1])))
1045 {
1046 $rem = [@$x] if wantarray; # make copy
1047 splice (@$x,1); # keep ref to org array
1048 $x->[0] = 0; # set to 0
1049 return ($x,$rem) if wantarray; # including remainder?
1050 return $x;
1051 }
1052 # now calculate $x / $yorg
091c87b1 1053
990fb837
RGS
1054 if (length(int($yorg->[-1])) == length(int($x->[-1])))
1055 {
b68b7ab1
T
1056 # same length, so make full compare
1057
990fb837
RGS
1058 my $a = 0; my $j = scalar @$x - 1;
1059 # manual way (abort if unequal, good for early ne)
1060 while ($j >= 0)
1061 {
1062 last if ($a = $x->[$j] - $yorg->[$j]); $j--;
1063 }
aef458a0 1064 # $a contains the result of the compare between X and Y
b68b7ab1 1065 # a < 0: x < y, a == 0: x == y, a > 0: x > y
990fb837
RGS
1066 if ($a <= 0)
1067 {
b68b7ab1
T
1068 $rem = [ 0 ]; # a = 0 => x == y => rem 0
1069 $rem = [@$x] if $a != 0; # a < 0 => x < y => rem = x
aef458a0 1070 splice(@$x,1); # keep single element
990fb837 1071 $x->[0] = 0; # if $a < 0
b68b7ab1 1072 $x->[0] = 1 if $a == 0; # $x == $y
7596a890
RGS
1073 return ($x,$rem) if wantarray; # including remainder?
1074 return $x;
990fb837 1075 }
aef458a0 1076 # $x >= $y, so proceed normally
b68b7ab1 1077
990fb837 1078 }
990fb837
RGS
1079 }
1080
1081 # all other cases:
1082
1083 my $y = [ @$yorg ]; # always make copy to preserve
61f5c3f5
T
1084
1085 my ($car,$bar,$prd,$dd,$xi,$yi,@q,$v2,$v1,@d,$tmp,$q,$u2,$u1,$u0);
1086
1087 $car = $bar = $prd = 0;
2ebb273f 1088 if (($dd = int($BASE/($y->[-1]+1))) != 1)
ee15d750
JH
1089 {
1090 for $xi (@$x)
1091 {
1092 $xi = $xi * $dd + $car;
2ebb273f 1093 $xi -= ($car = int($xi / $BASE)) * $BASE;
ee15d750
JH
1094 }
1095 push(@$x, $car); $car = 0;
1096 for $yi (@$y)
1097 {
1098 $yi = $yi * $dd + $car;
2ebb273f 1099 $yi -= ($car = int($yi / $BASE)) * $BASE;
ee15d750
JH
1100 }
1101 }
1102 else
1103 {
1104 push(@$x, 0);
1105 }
aef458a0
JH
1106
1107 # @q will accumulate the final result, $q contains the current computed
1108 # part of the final result
1109
ee15d750
JH
1110 @q = (); ($v2,$v1) = @$y[-2,-1];
1111 $v2 = 0 unless $v2;
1112 while ($#$x > $#$y)
1113 {
1114 ($u2,$u1,$u0) = @$x[-3..-1];
1115 $u2 = 0 unless $u2;
1116 #warn "oups v1 is 0, u0: $u0 $y->[-2] $y->[-1] l ",scalar @$y,"\n"
1117 # if $v1 == 0;
2ebb273f
T
1118 $q = (($u0 == $v1) ? $MAX_VAL : int(($u0*$BASE+$u1)/$v1));
1119 --$q while ($v2*$q > ($u0*$BASE+$u1-$q*$v1)*$BASE+$u2);
ee15d750
JH
1120 if ($q)
1121 {
1122 ($car, $bar) = (0,0);
1123 for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
1124 {
1125 $prd = $q * $y->[$yi] + $car;
2ebb273f
T
1126 $prd -= ($car = int($prd / $BASE)) * $BASE;
1127 $x->[$xi] += $BASE if ($bar = (($x->[$xi] -= $prd + $bar) < 0));
ee15d750
JH
1128 }
1129 if ($x->[-1] < $car + $bar)
1130 {
1131 $car = 0; --$q;
1132 for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
1133 {
2ebb273f
T
1134 $x->[$xi] -= $BASE
1135 if ($car = (($x->[$xi] += $y->[$yi] + $car) >= $BASE));
ee15d750
JH
1136 }
1137 }
1138 }
61f5c3f5 1139 pop(@$x); unshift(@q, $q);
ee15d750
JH
1140 }
1141 if (wantarray)
1142 {
1143 @d = ();
1144 if ($dd != 1)
1145 {
1146 $car = 0;
1147 for $xi (reverse @$x)
1148 {
2ebb273f 1149 $prd = $car * $BASE + $xi;
ee15d750
JH
1150 $car = $prd - ($tmp = int($prd / $dd)) * $dd;
1151 unshift(@d, $tmp);
1152 }
1153 }
1154 else
1155 {
1156 @d = @$x;
1157 }
1158 @$x = @q;
61f5c3f5 1159 my $d = \@d;
990fb837
RGS
1160 __strip_zeros($x);
1161 __strip_zeros($d);
61f5c3f5 1162 return ($x,$d);
ee15d750
JH
1163 }
1164 @$x = @q;
990fb837 1165 __strip_zeros($x);
61f5c3f5 1166 $x;
ee15d750
JH
1167 }
1168
394e6ffb
JH
1169##############################################################################
1170# testing
1171
1172sub _acmp
1173 {
1174 # internal absolute post-normalized compare (ignore signs)
1175 # ref to array, ref to array, return <0, 0, >0
1176 # arrays must have at least one entry; this is not checked for
394e6ffb 1177 my ($c,$cx,$cy) = @_;
091c87b1
T
1178
1179 # shortcut for short numbers
1180 return (($cx->[0] <=> $cy->[0]) <=> 0)
1181 if scalar @$cx == scalar @$cy && scalar @$cx == 1;
394e6ffb 1182
f9a08e12 1183 # fast comp based on number of array elements (aka pseudo-length)
091c87b1
T
1184 my $lxy = (scalar @$cx - scalar @$cy)
1185 # or length of first element if same number of elements (aka difference 0)
1186 ||
1187 # need int() here because sometimes the last element is '00018' vs '18'
1188 (length(int($cx->[-1])) - length(int($cy->[-1])));
394e6ffb
JH
1189 return -1 if $lxy < 0; # already differs, ret
1190 return 1 if $lxy > 0; # ditto
56d9de68 1191
394e6ffb 1192 # manual way (abort if unequal, good for early ne)
091c87b1
T
1193 my $a; my $j = scalar @$cx;
1194 while (--$j >= 0)
9393ace2 1195 {
091c87b1 1196 last if ($a = $cx->[$j] - $cy->[$j]);
9393ace2 1197 }
091c87b1 1198 $a <=> 0;
394e6ffb
JH
1199 }
1200
1201sub _len
1202 {
50109ad0 1203 # compute number of digits in base 10
394e6ffb
JH
1204
1205 # int() because add/sub sometimes leaves strings (like '00005') instead of
1206 # '5' in this place, thus causing length() to report wrong length
1207 my $cx = $_[1];
1208
56d9de68 1209 (@$cx-1)*$BASE_LEN+length(int($cx->[-1]));
394e6ffb
JH
1210 }
1211
1212sub _digit
1213 {
84f8f1ae
PJA
1214 # Return the nth digit. Zero is rightmost, so _digit(123,0) gives 3.
1215 # Negative values count from the left, so _digit(123, -1) gives 1.
394e6ffb
JH
1216 my ($c,$x,$n) = @_;
1217
1218 my $len = _len('',$x);
1219
84f8f1ae
PJA
1220 $n += $len if $n < 0; # -1 last, -2 second-to-last
1221 return "0" if $n < 0 || $n >= $len; # return 0 for digits out of range
1222
1223 my $elem = int($n / $BASE_LEN); # which array element
1224 my $digit = $n % $BASE_LEN; # which digit in this element
1225 substr("$x->[$elem]", -$digit-1, 1);
394e6ffb
JH
1226 }
1227
1228sub _zeros
1229 {
1230 # return amount of trailing zeros in decimal
1231 # check each array elem in _m for having 0 at end as long as elem == 0
1232 # Upon finding a elem != 0, stop
1233 my $x = $_[1];
9b924220
RGS
1234
1235 return 0 if scalar @$x == 1 && $x->[0] == 0;
1236
394e6ffb
JH
1237 my $zeros = 0; my $elem;
1238 foreach my $e (@$x)
1239 {
1240 if ($e != 0)
1241 {
1242 $elem = "$e"; # preserve x
1243 $elem =~ s/.*?(0*$)/$1/; # strip anything not zero
1244 $zeros *= $BASE_LEN; # elems * 5
61f5c3f5 1245 $zeros += length($elem); # count trailing zeros
394e6ffb
JH
1246 last; # early out
1247 }
1248 $zeros ++; # real else branch: 50% slower!
1249 }
61f5c3f5 1250 $zeros;
394e6ffb
JH
1251 }
1252
1253##############################################################################
1254# _is_* routines
1255
1256sub _is_zero
1257 {
9b924220
RGS
1258 # return true if arg is zero
1259 (((scalar @{$_[1]} == 1) && ($_[1]->[0] == 0))) <=> 0;
394e6ffb
JH
1260 }
1261
1262sub _is_even
1263 {
9b924220
RGS
1264 # return true if arg is even
1265 (!($_[1]->[0] & 1)) <=> 0;
394e6ffb
JH
1266 }
1267
1268sub _is_odd
1269 {
d5735945
PF
1270 # return true if arg is odd
1271 (($_[1]->[0] & 1)) <=> 0;
394e6ffb
JH
1272 }
1273
1274sub _is_one
1275 {
9b924220
RGS
1276 # return true if arg is one
1277 (scalar @{$_[1]} == 1) && ($_[1]->[0] == 1) <=> 0;
1278 }
1279
1280sub _is_two
1281 {
1282 # return true if arg is two
1283 (scalar @{$_[1]} == 1) && ($_[1]->[0] == 2) <=> 0;
1284 }
61f5c3f5 1285
9b924220
RGS
1286sub _is_ten
1287 {
1288 # return true if arg is ten
1289 (scalar @{$_[1]} == 1) && ($_[1]->[0] == 10) <=> 0;
394e6ffb
JH
1290 }
1291
1292sub __strip_zeros
1293 {
1294 # internal normalization function that strips leading zeros from the array
1295 # args: ref to array
1296 my $s = shift;
1297
1298 my $cnt = scalar @$s; # get count of parts
1299 my $i = $cnt-1;
1300 push @$s,0 if $i < 0; # div might return empty results, so fix it
1301
61f5c3f5
T
1302 return $s if @$s == 1; # early out
1303
394e6ffb
JH
1304 #print "strip: cnt $cnt i $i\n";
1305 # '0', '3', '4', '0', '0',
1306 # 0 1 2 3 4
1307 # cnt = 5, i = 4
1308 # i = 4
1309 # i = 3
1310 # => fcnt = cnt - i (5-2 => 3, cnt => 5-1 = 4, throw away from 4th pos)
1311 # >= 1: skip first part (this can be zero)
1312 while ($i > 0) { last if $s->[$i] != 0; $i--; }
1313 $i++; splice @$s,$i if ($i < $cnt); # $i cant be 0
1314 $s;
1315 }
1316
1317###############################################################################
3a427a11 1318# check routine to test internal state for corruptions
394e6ffb
JH
1319
1320sub _check
1321 {
1322 # used by the test suite
1323 my $x = $_[1];
1324
1325 return "$x is not a reference" if !ref($x);
1326
1327 # are all parts are valid?
1328 my $i = 0; my $j = scalar @$x; my ($e,$try);
1329 while ($i < $j)
1330 {
1331 $e = $x->[$i]; $e = 'undef' unless defined $e;
1332 $try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e)";
1333 last if $e !~ /^[+]?[0-9]+$/;
1334 $try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e) (stringify)";
1335 last if "$e" !~ /^[+]?[0-9]+$/;
1336 $try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e) (cat-stringify)";
1337 last if '' . "$e" !~ /^[+]?[0-9]+$/;
1338 $try = ' < 0 || >= $BASE; '."($x, $e)";
1339 last if $e <0 || $e >= $BASE;
1340 # this test is disabled, since new/bnorm and certain ops (like early out
1341 # in add/sub) are allowed/expected to leave '00000' in some elements
1342 #$try = '=~ /^00+/; '."($x, $e)";
1343 #last if $e =~ /^00+/;
1344 $i++;
1345 }
1346 return "Illegal part '$e' at pos $i (tested: $try)" if $i < $j;
3a427a11 1347 0;
394e6ffb
JH
1348 }
1349
1350
1351###############################################################################
394e6ffb 1352
dccbb853
JH
1353sub _mod
1354 {
1355 # if possible, use mod shortcut
1356 my ($c,$x,$yo) = @_;
1357
de1ac46b 1358 # slow way since $y too big
dccbb853
JH
1359 if (scalar @$yo > 1)
1360 {
1361 my ($xo,$rem) = _div($c,$x,$yo);
de1ac46b
PJA
1362 @$x = @$rem;
1363 return $x;
dccbb853 1364 }
aef458a0 1365
dccbb853 1366 my $y = $yo->[0];
de1ac46b
PJA
1367
1368 # if both are single element arrays
dccbb853
JH
1369 if (scalar @$x == 1)
1370 {
1371 $x->[0] %= $y;
1372 return $x;
1373 }
1374
de1ac46b 1375 # if @$x has more than one element, but @$y is a single element
dccbb853
JH
1376 my $b = $BASE % $y;
1377 if ($b == 0)
1378 {
1379 # when BASE % Y == 0 then (B * BASE) % Y == 0
1380 # (B * BASE) % $y + A % Y => A % Y
1381 # so need to consider only last element: O(1)
1382 $x->[0] %= $y;
1383 }
027dc388
JH
1384 elsif ($b == 1)
1385 {
de1ac46b
PJA
1386 # else need to go through all elements in @$x: O(N), but loop is a bit
1387 # simplified
027dc388
JH
1388 my $r = 0;
1389 foreach (@$x)
1390 {
28df3e88
JH
1391 $r = ($r + $_) % $y; # not much faster, but heh...
1392 #$r += $_ % $y; $r %= $y;
027dc388
JH
1393 }
1394 $r = 0 if $r == $y;
1395 $x->[0] = $r;
1396 }
dccbb853
JH
1397 else
1398 {
de1ac46b
PJA
1399 # else need to go through all elements in @$x: O(N)
1400 my $r = 0;
1401 my $bm = 1;
027dc388
JH
1402 foreach (@$x)
1403 {
28df3e88
JH
1404 $r = ($_ * $bm + $r) % $y;
1405 $bm = ($bm * $b) % $y;
1406
1407 #$r += ($_ % $y) * $bm;
1408 #$bm *= $b;
1409 #$bm %= $y;
1410 #$r %= $y;
027dc388
JH
1411 }
1412 $r = 0 if $r == $y;
1413 $x->[0] = $r;
dccbb853 1414 }
de1ac46b
PJA
1415 @$x = $x->[0]; # keep one element of @$x
1416 return $x;
dccbb853
JH
1417 }
1418
0716bf9b 1419##############################################################################
574bacfe
JH
1420# shifts
1421
1422sub _rsft
1423 {
1424 my ($c,$x,$y,$n) = @_;
1425
1426 if ($n != 10)
1427 {
9b924220 1428 $n = _new($c,$n); return _div($c,$x, _pow($c,$n,$y));
61f5c3f5
T
1429 }
1430
1431 # shortcut (faster) for shifting by 10)
1432 # multiples of $BASE_LEN
1433 my $dst = 0; # destination
1434 my $src = _num($c,$y); # as normal int
56d9de68 1435 my $xlen = (@$x-1)*$BASE_LEN+length(int($x->[-1])); # len of x in digits
90d1b129 1436 if ($src >= $xlen or ($src == $xlen and ! defined $x->[1]))
56d9de68
T
1437 {
1438 # 12345 67890 shifted right by more than 10 digits => 0
1439 splice (@$x,1); # leave only one element
1440 $x->[0] = 0; # set to zero
1441 return $x;
1442 }
61f5c3f5
T
1443 my $rem = $src % $BASE_LEN; # remainder to shift
1444 $src = int($src / $BASE_LEN); # source
1445 if ($rem == 0)
1446 {
1447 splice (@$x,0,$src); # even faster, 38.4 => 39.3
574bacfe
JH
1448 }
1449 else
1450 {
61f5c3f5
T
1451 my $len = scalar @$x - $src; # elems to go
1452 my $vd; my $z = '0'x $BASE_LEN;
1453 $x->[scalar @$x] = 0; # avoid || 0 test inside loop
1454 while ($dst < $len)
574bacfe 1455 {
61f5c3f5
T
1456 $vd = $z.$x->[$src];
1457 $vd = substr($vd,-$BASE_LEN,$BASE_LEN-$rem);
1458 $src++;
1459 $vd = substr($z.$x->[$src],-$rem,$rem) . $vd;
1460 $vd = substr($vd,-$BASE_LEN,$BASE_LEN) if length($vd) > $BASE_LEN;
1461 $x->[$dst] = int($vd);
1462 $dst++;
574bacfe 1463 }
61f5c3f5 1464 splice (@$x,$dst) if $dst > 0; # kill left-over array elems
56b9c951 1465 pop @$x if $x->[-1] == 0 && @$x > 1; # kill last element if 0
61f5c3f5 1466 } # else rem == 0
574bacfe
JH
1467 $x;
1468 }
1469
1470sub _lsft
1471 {
1472 my ($c,$x,$y,$n) = @_;
1473
1474 if ($n != 10)
1475 {
9b924220 1476 $n = _new($c,$n); return _mul($c,$x, _pow($c,$n,$y));
574bacfe 1477 }
61f5c3f5
T
1478
1479 # shortcut (faster) for shifting by 10) since we are in base 10eX
1480 # multiples of $BASE_LEN:
1481 my $src = scalar @$x; # source
1482 my $len = _num($c,$y); # shift-len as normal int
1483 my $rem = $len % $BASE_LEN; # remainder to shift
1484 my $dst = $src + int($len/$BASE_LEN); # destination
1485 my $vd; # further speedup
1486 $x->[$src] = 0; # avoid first ||0 for speed
1487 my $z = '0' x $BASE_LEN;
1488 while ($src >= 0)
574bacfe 1489 {
61f5c3f5
T
1490 $vd = $x->[$src]; $vd = $z.$vd;
1491 $vd = substr($vd,-$BASE_LEN+$rem,$BASE_LEN-$rem);
1492 $vd .= $src > 0 ? substr($z.$x->[$src-1],-$BASE_LEN,$rem) : '0' x $rem;
1493 $vd = substr($vd,-$BASE_LEN,$BASE_LEN) if length($vd) > $BASE_LEN;
1494 $x->[$dst] = int($vd);
1495 $dst--; $src--;
574bacfe 1496 }
61f5c3f5
T
1497 # set lowest parts to 0
1498 while ($dst >= 0) { $x->[$dst--] = 0; }
c4a6f826 1499 # fix spurious last zero element
61f5c3f5 1500 splice @$x,-1 if $x->[-1] == 0;
574bacfe
JH
1501 $x;
1502 }
1503
027dc388
JH
1504sub _pow
1505 {
1506 # power of $x to $y
1507 # ref to array, ref to array, return ref to array
1508 my ($c,$cx,$cy) = @_;
1509
b282a552
T
1510 if (scalar @$cy == 1 && $cy->[0] == 0)
1511 {
1512 splice (@$cx,1); $cx->[0] = 1; # y == 0 => x => 1
1513 return $cx;
1514 }
1515 if ((scalar @$cx == 1 && $cx->[0] == 1) || # x == 1
1516 (scalar @$cy == 1 && $cy->[0] == 1)) # or y == 1
1517 {
1518 return $cx;
1519 }
1520 if (scalar @$cx == 1 && $cx->[0] == 0)
1521 {
1522 splice (@$cx,1); $cx->[0] = 0; # 0 ** y => 0 (if not y <= 0)
1523 return $cx;
1524 }
1525
027dc388 1526 my $pow2 = _one();
1ddff52a 1527
9b924220 1528 my $y_bin = _as_bin($c,$cy); $y_bin =~ s/^0b//;
1ddff52a
JH
1529 my $len = length($y_bin);
1530 while (--$len > 0)
027dc388 1531 {
1ddff52a 1532 _mul($c,$pow2,$cx) if substr($y_bin,$len,1) eq '1'; # is odd?
027dc388
JH
1533 _mul($c,$cx,$cx);
1534 }
1ddff52a
JH
1535
1536 _mul($c,$cx,$pow2);
61f5c3f5 1537 $cx;
027dc388
JH
1538 }
1539
de1ac46b
PJA
1540sub _nok {
1541 # Return binomial coefficient (n over k).
1542 # Given refs to arrays, return ref to array.
1543 # First input argument is modified.
50109ad0 1544
de1ac46b
PJA
1545 my ($c, $n, $k) = @_;
1546
1547 # If k > n/2, or, equivalently, 2*k > n, compute nok(n, k) as
1548 # nok(n, n-k), to minimize the number if iterations in the loop.
50109ad0 1549
50109ad0 1550 {
de1ac46b
PJA
1551 my $twok = _mul($c, _two($c), _copy($c, $k)); # 2 * k
1552 if (_acmp($c, $twok, $n) > 0) { # if 2*k > n
1553 $k = _sub($c, _copy($c, $n), $k); # k = n - k
1554 }
50109ad0 1555 }
de1ac46b
PJA
1556
1557 # Example:
1558 #
1559 # / 7 \ 7! 1*2*3*4 * 5*6*7 5 * 6 * 7 6 7
1560 # | | = --------- = --------------- = --------- = 5 * - * -
1561 # \ 3 / (7-3)! 3! 1*2*3*4 * 1*2*3 1 * 2 * 3 2 3
1562
1563 if (_is_zero($c, $k)) {
1564 @$n = 1;
50109ad0 1565 }
de1ac46b
PJA
1566
1567 else {
1568
ada8209b 1569 # Make a copy of the original n, since we'll be modifying n in-place.
de1ac46b
PJA
1570
1571 my $n_orig = _copy($c, $n);
1572
1573 # n = 5, f = 6, d = 2 (cf. example above)
1574
1575 _sub($c, $n, $k);
1576 _inc($c, $n);
1577
1578 my $f = _copy($c, $n);
1579 _inc($c, $f);
1580
1581 my $d = _two($c);
1582
1583 # while f <= n (the original n, that is) ...
1584
1585 while (_acmp($c, $f, $n_orig) <= 0) {
1586
1587 # n = (n * f / d) == 5 * 6 / 2 (cf. example above)
1588
1589 _mul($c, $n, $f);
1590 _div($c, $n, $d);
1591
1592 # f = 7, d = 3 (cf. example above)
1593
1594 _inc($c, $f);
1595 _inc($c, $d);
1596 }
1597
1598 }
1599
1600 return $n;
1601}
50109ad0
RGS
1602
1603my @factorials = (
1604 1,
1605 1,
1606 2,
1607 2*3,
1608 2*3*4,
1609 2*3*4*5,
1610 2*3*4*5*6,
1611 2*3*4*5*6*7,
1612);
1613
b3abae2a
JH
1614sub _fac
1615 {
1616 # factorial of $x
1617 # ref to array, return ref to array
1618 my ($c,$cx) = @_;
1619
50109ad0 1620 if ((@$cx == 1) && ($cx->[0] <= 7))
b3abae2a 1621 {
50109ad0 1622 $cx->[0] = $factorials[$cx->[0]]; # 0 => 1, 1 => 1, 2 => 2 etc.
b3abae2a
JH
1623 return $cx;
1624 }
1625
50109ad0
RGS
1626 if ((@$cx == 1) && # we do this only if $x >= 12 and $x <= 7000
1627 ($cx->[0] >= 12 && $cx->[0] < 7000))
1628 {
1629
1630 # Calculate (k-j) * (k-j+1) ... k .. (k+j-1) * (k + j)
1631 # See http://blogten.blogspot.com/2007/01/calculating-n.html
1632 # The above series can be expressed as factors:
1633 # k * k - (j - i) * 2
1634 # We cache k*k, and calculate (j * j) as the sum of the first j odd integers
1635
1636 # This will not work when N exceeds the storage of a Perl scalar, however,
1637 # in this case the algorithm would be way to slow to terminate, anyway.
1638
1639 # As soon as the last element of $cx is 0, we split it up and remember
1640 # how many zeors we got so far. The reason is that n! will accumulate
1641 # zeros at the end rather fast.
1642 my $zero_elements = 0;
1643
1644 # If n is even, set n = n -1
1645 my $k = _num($c,$cx); my $even = 1;
1646 if (($k & 1) == 0)
1647 {
1648 $even = $k; $k --;
1649 }
1650 # set k to the center point
1651 $k = ($k + 1) / 2;
1652# print "k $k even: $even\n";
1653 # now calculate k * k
1654 my $k2 = $k * $k;
1655 my $odd = 1; my $sum = 1;
1656 my $i = $k - 1;
1657 # keep reference to x
1658 my $new_x = _new($c, $k * $even);
1659 @$cx = @$new_x;
1660 if ($cx->[0] == 0)
1661 {
1662 $zero_elements ++; shift @$cx;
1663 }
1664# print STDERR "x = ", _str($c,$cx),"\n";
1665 my $BASE2 = int(sqrt($BASE))-1;
1666 my $j = 1;
1667 while ($j <= $i)
1668 {
1669 my $m = ($k2 - $sum); $odd += 2; $sum += $odd; $j++;
1670 while ($j <= $i && ($m < $BASE2) && (($k2 - $sum) < $BASE2))
1671 {
1672 $m *= ($k2 - $sum);
1673 $odd += 2; $sum += $odd; $j++;
1674# print STDERR "\n k2 $k2 m $m sum $sum odd $odd\n"; sleep(1);
1675 }
1676 if ($m < $BASE)
1677 {
1678 _mul($c,$cx,[$m]);
1679 }
1680 else
1681 {
1682 _mul($c,$cx,$c->_new($m));
1683 }
1684 if ($cx->[0] == 0)
1685 {
1686 $zero_elements ++; shift @$cx;
1687 }
1688# print STDERR "Calculate $k2 - $sum = $m (x = ", _str($c,$cx),")\n";
1689 }
1690 # multiply in the zeros again
1691 unshift @$cx, (0) x $zero_elements;
1692 return $cx;
1693 }
1694
b3abae2a 1695 # go forward until $base is exceeded
091c87b1
T
1696 # limit is either $x steps (steps == 100 means a result always too high) or
1697 # $base.
b3abae2a 1698 my $steps = 100; $steps = $cx->[0] if @$cx == 1;
091c87b1
T
1699 my $r = 2; my $cf = 3; my $step = 2; my $last = $r;
1700 while ($r*$cf < $BASE && $step < $steps)
b3abae2a
JH
1701 {
1702 $last = $r; $r *= $cf++; $step++;
1703 }
091c87b1 1704 if ((@$cx == 1) && $step == $cx->[0])
b3abae2a 1705 {
091c87b1
T
1706 # completely done, so keep reference to $x and return
1707 $cx->[0] = $r;
b3abae2a
JH
1708 return $cx;
1709 }
091c87b1 1710
990fb837 1711 # now we must do the left over steps
091c87b1
T
1712 my $n; # steps still to do
1713 if (scalar @$cx == 1)
1714 {
1715 $n = $cx->[0];
1716 }
1717 else
1718 {
1719 $n = _copy($c,$cx);
1720 }
b3abae2a 1721
50109ad0
RGS
1722 # Set $cx to the last result below $BASE (but keep ref to $x)
1723 $cx->[0] = $last; splice (@$cx,1);
1724 # As soon as the last element of $cx is 0, we split it up and remember
1725 # how many zeors we got so far. The reason is that n! will accumulate
1726 # zeros at the end rather fast.
990fb837 1727 my $zero_elements = 0;
091c87b1
T
1728
1729 # do left-over steps fit into a scalar?
1730 if (ref $n eq 'ARRAY')
b3abae2a 1731 {
091c87b1 1732 # No, so use slower inc() & cmp()
50109ad0
RGS
1733 # ($n is at least $BASE here)
1734 my $base_2 = int(sqrt($BASE)) - 1;
1735 #print STDERR "base_2: $base_2\n";
1736 while ($step < $base_2)
1737 {
1738 if ($cx->[0] == 0)
1739 {
1740 $zero_elements ++; shift @$cx;
1741 }
1742 my $b = $step * ($step + 1); $step += 2;
1743 _mul($c,$cx,[$b]);
1744 }
091c87b1 1745 $step = [$step];
50109ad0 1746 while (_acmp($c,$step,$n) <= 0)
990fb837
RGS
1747 {
1748 if ($cx->[0] == 0)
1749 {
1750 $zero_elements ++; shift @$cx;
1751 }
091c87b1 1752 _mul($c,$cx,$step); _inc($c,$step);
990fb837 1753 }
990fb837 1754 }
091c87b1 1755 else
990fb837 1756 {
091c87b1 1757 # Yes, so we can speed it up slightly
50109ad0
RGS
1758
1759# print "# left over steps $n\n";
1760
1761 my $base_4 = int(sqrt(sqrt($BASE))) - 2;
1762 #print STDERR "base_4: $base_4\n";
1763 my $n4 = $n - 4;
1764 while ($step < $n4 && $step < $base_4)
990fb837 1765 {
091c87b1
T
1766 if ($cx->[0] == 0)
1767 {
1768 $zero_elements ++; shift @$cx;
1769 }
50109ad0
RGS
1770 my $b = $step * ($step + 1); $step += 2; $b *= $step * ($step + 1); $step += 2;
1771 _mul($c,$cx,[$b]);
1772 }
1773 my $base_2 = int(sqrt($BASE)) - 1;
1774 my $n2 = $n - 2;
1775 #print STDERR "base_2: $base_2\n";
1776 while ($step < $n2 && $step < $base_2)
1777 {
1778 if ($cx->[0] == 0)
1779 {
1780 $zero_elements ++; shift @$cx;
1781 }
1782 my $b = $step * ($step + 1); $step += 2;
1783 _mul($c,$cx,[$b]);
1784 }
1785 # do what's left over
1786 while ($step <= $n)
1787 {
091c87b1 1788 _mul($c,$cx,[$step]); $step++;
50109ad0
RGS
1789 if ($cx->[0] == 0)
1790 {
1791 $zero_elements ++; shift @$cx;
1792 }
990fb837 1793 }
990fb837
RGS
1794 }
1795 # multiply in the zeros again
50109ad0 1796 unshift @$cx, (0) x $zero_elements;
091c87b1
T
1797 $cx; # return result
1798 }
1799
9b924220
RGS
1800#############################################################################
1801
091c87b1
T
1802sub _log_int
1803 {
1804 # calculate integer log of $x to base $base
1805 # ref to array, ref to array - return ref to array
1806 my ($c,$x,$base) = @_;
1807
1808 # X == 0 => NaN
1809 return if (scalar @$x == 1 && $x->[0] == 0);
1810 # BASE 0 or 1 => NaN
1811 return if (scalar @$base == 1 && $base->[0] < 2);
b282a552 1812 my $cmp = _acmp($c,$x,$base); # X == BASE => 1
091c87b1
T
1813 if ($cmp == 0)
1814 {
1815 splice (@$x,1); $x->[0] = 1;
8df1e0a2 1816 return ($x,1)
091c87b1
T
1817 }
1818 # X < BASE
1819 if ($cmp < 0)
1820 {
1821 splice (@$x,1); $x->[0] = 0;
8df1e0a2 1822 return ($x,undef);
091c87b1
T
1823 }
1824
091c87b1 1825 my $x_org = _copy($c,$x); # preserve x
8df1e0a2
T
1826 splice(@$x,1); $x->[0] = 1; # keep ref to $x
1827
86b76201
T
1828 # Compute a guess for the result based on:
1829 # $guess = int ( length_in_base_10(X) / ( log(base) / log(10) ) )
2ebb273f
T
1830 my $len = _len($c,$x_org);
1831 my $log = log($base->[-1]) / log(10);
b282a552 1832
2ebb273f
T
1833 # for each additional element in $base, we add $BASE_LEN to the result,
1834 # based on the observation that log($BASE,10) is BASE_LEN and
1835 # log(x*y) == log(x) + log(y):
1836 $log += ((scalar @$base)-1) * $BASE_LEN;
b282a552 1837
2ebb273f
T
1838 # calculate now a guess based on the values obtained above:
1839 my $res = int($len / $log);
1840
1841 $x->[0] = $res;
1842 my $trial = _pow ($c, _copy($c, $base), $x);
1843 my $a = _acmp($c,$trial,$x_org);
1844
1845# print STDERR "# trial ", _str($c,$x)," was: $a (0 = exact, -1 too small, +1 too big)\n";
1846
1847 # found an exact result?
1848 return ($x,1) if $a == 0;
1849
1850 if ($a > 0)
1851 {
1852 # or too big
1853 _div($c,$trial,$base); _dec($c, $x);
1854 while (($a = _acmp($c,$trial,$x_org)) > 0)
b282a552 1855 {
2ebb273f
T
1856# print STDERR "# big _log_int at ", _str($c,$x), "\n";
1857 _div($c,$trial,$base); _dec($c, $x);
b282a552 1858 }
2ebb273f
T
1859 # result is now exact (a == 0), or too small (a < 0)
1860 return ($x, $a == 0 ? 1 : 0);
1861 }
1862
1863 # else: result was to small
1864 _mul($c,$trial,$base);
1865
1866 # did we now get the right result?
1867 $a = _acmp($c,$trial,$x_org);
1868
1869 if ($a == 0) # yes, exactly
1870 {
1871 _inc($c, $x);
1872 return ($x,1);
1873 }
1874 return ($x,0) if $a > 0;
1875
1876 # Result still too small (we should come here only if the estimate above
1877 # was very off base):
1878
1879 # Now let the normal trial run obtain the real result
1880 # Simple loop that increments $x by 2 in each step, possible overstepping
1881 # the real result
091c87b1 1882
2ebb273f 1883 my $base_mul = _mul($c, _copy($c,$base), $base); # $base * $base
8df1e0a2 1884
9b924220 1885 while (($a = _acmp($c,$trial,$x_org)) < 0)
091c87b1 1886 {
2ebb273f 1887# print STDERR "# small _log_int at ", _str($c,$x), "\n";
8df1e0a2 1888 _mul($c,$trial,$base_mul); _add($c, $x, [2]);
091c87b1 1889 }
8df1e0a2
T
1890
1891 my $exact = 1;
1892 if ($a > 0)
091c87b1 1893 {
8df1e0a2
T
1894 # overstepped the result
1895 _dec($c, $x);
1896 _div($c,$trial,$base);
9b924220 1897 $a = _acmp($c,$trial,$x_org);
8df1e0a2 1898 if ($a > 0)
091c87b1 1899 {
8df1e0a2 1900 _dec($c, $x);
091c87b1 1901 }
2ebb273f 1902 $exact = 0 if $a != 0; # a = -1 => not exact result, a = 0 => exact
091c87b1
T
1903 }
1904
8df1e0a2 1905 ($x,$exact); # return result
b3abae2a
JH
1906 }
1907
56d9de68
T
1908# for debugging:
1909 use constant DEBUG => 0;
1910 my $steps = 0;
1911 sub steps { $steps };
b3abae2a
JH
1912
1913sub _sqrt
0716bf9b 1914 {
56d9de68 1915 # square-root of $x in place
990fb837 1916 # Compute a guess of the result (by rule of thumb), then improve it via
56d9de68 1917 # Newton's method.
394e6ffb 1918 my ($c,$x) = @_;
0716bf9b 1919
394e6ffb
JH
1920 if (scalar @$x == 1)
1921 {
50109ad0 1922 # fits into one Perl scalar, so result can be computed directly
394e6ffb
JH
1923 $x->[0] = int(sqrt($x->[0]));
1924 return $x;
1925 }
1926 my $y = _copy($c,$x);
b3abae2a
JH
1927 # hopefully _len/2 is < $BASE, the -1 is to always undershot the guess
1928 # since our guess will "grow"
1929 my $l = int((_len($c,$x)-1) / 2);
1930
56d9de68 1931 my $lastelem = $x->[-1]; # for guess
b3abae2a
JH
1932 my $elems = scalar @$x - 1;
1933 # not enough digits, but could have more?
56d9de68 1934 if ((length($lastelem) <= 3) && ($elems > 1))
b3abae2a
JH
1935 {
1936 # right-align with zero pad
1937 my $len = length($lastelem) & 1;
1938 print "$lastelem => " if DEBUG;
1939 $lastelem .= substr($x->[-2] . '0' x $BASE_LEN,0,$BASE_LEN);
1940 # former odd => make odd again, or former even to even again
56d9de68 1941 $lastelem = $lastelem / 10 if (length($lastelem) & 1) != $len;
b3abae2a
JH
1942 print "$lastelem\n" if DEBUG;
1943 }
0716bf9b 1944
61f5c3f5
T
1945 # construct $x (instead of _lsft($c,$x,$l,10)
1946 my $r = $l % $BASE_LEN; # 10000 00000 00000 00000 ($BASE_LEN=5)
1947 $l = int($l / $BASE_LEN);
b3abae2a 1948 print "l = $l " if DEBUG;
56d9de68
T
1949
1950 splice @$x,$l; # keep ref($x), but modify it
1951
b3abae2a
JH
1952 # we make the first part of the guess not '1000...0' but int(sqrt($lastelem))
1953 # that gives us:
56d9de68
T
1954 # 14400 00000 => sqrt(14400) => guess first digits to be 120
1955 # 144000 000000 => sqrt(144000) => guess 379
b3abae2a 1956
b3abae2a
JH
1957 print "$lastelem (elems $elems) => " if DEBUG;
1958 $lastelem = $lastelem / 10 if ($elems & 1 == 1); # odd or even?
1959 my $g = sqrt($lastelem); $g =~ s/\.//; # 2.345 => 2345
1960 $r -= 1 if $elems & 1 == 0; # 70 => 7
1961
1962 # padd with zeros if result is too short
1963 $x->[$l--] = int(substr($g . '0' x $r,0,$r+1));
1964 print "now ",$x->[-1] if DEBUG;
1965 print " would have been ", int('1' . '0' x $r),"\n" if DEBUG;
56d9de68 1966
b3abae2a 1967 # If @$x > 1, we could compute the second elem of the guess, too, to create
56d9de68 1968 # an even better guess. Not implemented yet. Does it improve performance?
b3abae2a 1969 $x->[$l--] = 0 while ($l >= 0); # all other digits of guess are zero
56d9de68 1970
9b924220 1971 print "start x= ",_str($c,$x),"\n" if DEBUG;
394e6ffb
JH
1972 my $two = _two();
1973 my $last = _zero();
1974 my $lastlast = _zero();
b3abae2a 1975 $steps = 0 if DEBUG;
394e6ffb
JH
1976 while (_acmp($c,$last,$x) != 0 && _acmp($c,$lastlast,$x) != 0)
1977 {
b3abae2a 1978 $steps++ if DEBUG;
394e6ffb
JH
1979 $lastlast = _copy($c,$last);
1980 $last = _copy($c,$x);
1981 _add($c,$x, _div($c,_copy($c,$y),$x));
1982 _div($c,$x, $two );
9b924220 1983 print " x= ",_str($c,$x),"\n" if DEBUG;
394e6ffb 1984 }
b3abae2a 1985 print "\nsteps in sqrt: $steps, " if DEBUG;
394e6ffb 1986 _dec($c,$x) if _acmp($c,$y,_mul($c,_copy($c,$x),$x)) < 0; # overshot?
b3abae2a 1987 print " final ",$x->[-1],"\n" if DEBUG;
394e6ffb 1988 $x;
0716bf9b
JH
1989 }
1990
990fb837
RGS
1991sub _root
1992 {
1993 # take n'th root of $x in place (n >= 3)
990fb837
RGS
1994 my ($c,$x,$n) = @_;
1995
1996 if (scalar @$x == 1)
1997 {
1998 if (scalar @$n > 1)
1999 {
2000 # result will always be smaller than 2 so trunc to 1 at once
2001 $x->[0] = 1;
2002 }
2003 else
2004 {
50109ad0 2005 # fits into one Perl scalar, so result can be computed directly
091c87b1
T
2006 # cannot use int() here, because it rounds wrongly (try
2007 # (81 ** 3) ** (1/3) to see what I mean)
2008 #$x->[0] = int( $x->[0] ** (1 / $n->[0]) );
2009 # round to 8 digits, then truncate result to integer
2010 $x->[0] = int ( sprintf ("%.8f", $x->[0] ** (1 / $n->[0]) ) );
990fb837
RGS
2011 }
2012 return $x;
2013 }
2014
3a427a11
RGS
2015 # we know now that X is more than one element long
2016
c38b2de2
JH
2017 # if $n is a power of two, we can repeatedly take sqrt($X) and find the
2018 # proper result, because sqrt(sqrt($x)) == root($x,4)
2019 my $b = _as_bin($c,$n);
9b924220 2020 if ($b =~ /0b1(0+)$/)
c38b2de2
JH
2021 {
2022 my $count = CORE::length($1); # 0b100 => len('00') => 2
2023 my $cnt = $count; # counter for loop
2024 unshift (@$x, 0); # add one element, together with one
2025 # more below in the loop this makes 2
2026 while ($cnt-- > 0)
2027 {
2028 # 'inflate' $X by adding one element, basically computing
2029 # $x * $BASE * $BASE. This gives us more $BASE_LEN digits for result
2030 # since len(sqrt($X)) approx == len($x) / 2.
2031 unshift (@$x, 0);
2032 # calculate sqrt($x), $x is now one element to big, again. In the next
2033 # round we make that two, again.
2034 _sqrt($c,$x);
2035 }
2036 # $x is now one element to big, so truncate result by removing it
2037 splice (@$x,0,1);
2038 }
2039 else
2040 {
091c87b1 2041 # trial computation by starting with 2,4,8,16 etc until we overstep
3a427a11 2042 my $step;
091c87b1
T
2043 my $trial = _two();
2044
3a427a11
RGS
2045 # while still to do more than X steps
2046 do
091c87b1 2047 {
3a427a11
RGS
2048 $step = _two();
2049 while (_acmp($c, _pow($c, _copy($c, $trial), $n), $x) < 0)
2050 {
2051 _mul ($c, $step, [2]);
2052 _add ($c, $trial, $step);
2053 }
2054
2055 # hit exactly?
2056 if (_acmp($c, _pow($c, _copy($c, $trial), $n), $x) == 0)
2057 {
2058 @$x = @$trial; # make copy while preserving ref to $x
2059 return $x;
2060 }
2061 # overstepped, so go back on step
2062 _sub($c, $trial, $step);
2063 } while (scalar @$step > 1 || $step->[0] > 128);
091c87b1 2064
3a427a11
RGS
2065 # reset step to 2
2066 $step = _two();
091c87b1 2067 # add two, because $trial cannot be exactly the result (otherwise we would
9681bfa6 2068 # already have found it)
091c87b1
T
2069 _add($c, $trial, $step);
2070
3a427a11
RGS
2071 # and now add more and more (2,4,6,8,10 etc)
2072 while (_acmp($c, _pow($c, _copy($c, $trial), $n), $x) < 0)
2073 {
2074 _add ($c, $trial, $step);
2075 }
091c87b1
T
2076
2077 # hit not exactly? (overstepped)
091c87b1
T
2078 if (_acmp($c, _pow($c, _copy($c, $trial), $n), $x) > 0)
2079 {
2080 _dec($c,$trial);
2081 }
3a427a11
RGS
2082
2083 # hit not exactly? (overstepped)
2084 # 80 too small, 81 slightly too big, 82 too big
091c87b1
T
2085 if (_acmp($c, _pow($c, _copy($c, $trial), $n), $x) > 0)
2086 {
3a427a11 2087 _dec ($c, $trial);
091c87b1 2088 }
3a427a11 2089
091c87b1
T
2090 @$x = @$trial; # make copy while preserving ref to $x
2091 return $x;
c38b2de2 2092 }
990fb837
RGS
2093 $x;
2094 }
2095
394e6ffb
JH
2096##############################################################################
2097# binary stuff
0716bf9b 2098
394e6ffb
JH
2099sub _and
2100 {
2101 my ($c,$x,$y) = @_;
0716bf9b 2102
394e6ffb
JH
2103 # the shortcut makes equal, large numbers _really_ fast, and makes only a
2104 # very small performance drop for small numbers (e.g. something with less
2105 # than 32 bit) Since we optimize for large numbers, this is enabled.
2106 return $x if _acmp($c,$x,$y) == 0; # shortcut
0716bf9b 2107
394e6ffb
JH
2108 my $m = _one(); my ($xr,$yr);
2109 my $mask = $AND_MASK;
2110
2111 my $x1 = $x;
2112 my $y1 = _copy($c,$y); # make copy
2113 $x = _zero();
2114 my ($b,$xrr,$yrr);
2115 use integer;
2116 while (!_is_zero($c,$x1) && !_is_zero($c,$y1))
2117 {
2118 ($x1, $xr) = _div($c,$x1,$mask);
2119 ($y1, $yr) = _div($c,$y1,$mask);
2120
2121 # make ints() from $xr, $yr
9b924220 2122 # this is when the AND_BITS are greater than $BASE and is slower for
394e6ffb
JH
2123 # small (<256 bits) numbers, but faster for large numbers. Disabled
2124 # due to KISS principle
2125
2126# $b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
2127# $b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
9b924220 2128# _add($c,$x, _mul($c, _new( $c, ($xrr & $yrr) ), $m) );
394e6ffb 2129
61f5c3f5
T
2130 # 0+ due to '&' doesn't work in strings
2131 _add($c,$x, _mul($c, [ 0+$xr->[0] & 0+$yr->[0] ], $m) );
394e6ffb
JH
2132 _mul($c,$m,$mask);
2133 }
2134 $x;
0716bf9b
JH
2135 }
2136
394e6ffb 2137sub _xor
0716bf9b 2138 {
394e6ffb
JH
2139 my ($c,$x,$y) = @_;
2140
2141 return _zero() if _acmp($c,$x,$y) == 0; # shortcut (see -and)
2142
2143 my $m = _one(); my ($xr,$yr);
2144 my $mask = $XOR_MASK;
2145
2146 my $x1 = $x;
2147 my $y1 = _copy($c,$y); # make copy
2148 $x = _zero();
2149 my ($b,$xrr,$yrr);
2150 use integer;
2151 while (!_is_zero($c,$x1) && !_is_zero($c,$y1))
0716bf9b 2152 {
394e6ffb
JH
2153 ($x1, $xr) = _div($c,$x1,$mask);
2154 ($y1, $yr) = _div($c,$y1,$mask);
2155 # make ints() from $xr, $yr (see _and())
2156 #$b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
2157 #$b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
9b924220 2158 #_add($c,$x, _mul($c, _new( $c, ($xrr ^ $yrr) ), $m) );
61f5c3f5
T
2159
2160 # 0+ due to '^' doesn't work in strings
2161 _add($c,$x, _mul($c, [ 0+$xr->[0] ^ 0+$yr->[0] ], $m) );
394e6ffb 2162 _mul($c,$m,$mask);
0716bf9b 2163 }
394e6ffb
JH
2164 # the loop stops when the shorter of the two numbers is exhausted
2165 # the remainder of the longer one will survive bit-by-bit, so we simple
2166 # multiply-add it in
2167 _add($c,$x, _mul($c, $x1, $m) ) if !_is_zero($c,$x1);
2168 _add($c,$x, _mul($c, $y1, $m) ) if !_is_zero($c,$y1);
2169
2170 $x;
0716bf9b
JH
2171 }
2172
394e6ffb 2173sub _or
0716bf9b 2174 {
394e6ffb 2175 my ($c,$x,$y) = @_;
0716bf9b 2176
394e6ffb 2177 return $x if _acmp($c,$x,$y) == 0; # shortcut (see _and)
0716bf9b 2178
394e6ffb
JH
2179 my $m = _one(); my ($xr,$yr);
2180 my $mask = $OR_MASK;
0716bf9b 2181
394e6ffb
JH
2182 my $x1 = $x;
2183 my $y1 = _copy($c,$y); # make copy
2184 $x = _zero();
2185 my ($b,$xrr,$yrr);
2186 use integer;
2187 while (!_is_zero($c,$x1) && !_is_zero($c,$y1))
2188 {
2189 ($x1, $xr) = _div($c,$x1,$mask);
2190 ($y1, $yr) = _div($c,$y1,$mask);
2191 # make ints() from $xr, $yr (see _and())
2192# $b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
2193# $b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
9b924220 2194# _add($c,$x, _mul($c, _new( $c, ($xrr | $yrr) ), $m) );
394e6ffb 2195
61f5c3f5
T
2196 # 0+ due to '|' doesn't work in strings
2197 _add($c,$x, _mul($c, [ 0+$xr->[0] | 0+$yr->[0] ], $m) );
394e6ffb
JH
2198 _mul($c,$m,$mask);
2199 }
2200 # the loop stops when the shorter of the two numbers is exhausted
2201 # the remainder of the longer one will survive bit-by-bit, so we simple
2202 # multiply-add it in
2203 _add($c,$x, _mul($c, $x1, $m) ) if !_is_zero($c,$x1);
2204 _add($c,$x, _mul($c, $y1, $m) ) if !_is_zero($c,$y1);
2205
2206 $x;
0716bf9b
JH
2207 }
2208
61f5c3f5
T
2209sub _as_hex
2210 {
2211 # convert a decimal number to hex (ref to array, return ref to string)
2212 my ($c,$x) = @_;
2213
50109ad0 2214 # fits into one element (handle also 0x0 case)
03874afe 2215 return sprintf("0x%x",$x->[0]) if @$x == 1;
990fb837 2216
61f5c3f5
T
2217 my $x1 = _copy($c,$x);
2218
2219 my $es = '';
1ddff52a
JH
2220 my ($xr, $h, $x10000);
2221 if ($] >= 5.006)
2222 {
2223 $x10000 = [ 0x10000 ]; $h = 'h4';
2224 }
2225 else
2226 {
2227 $x10000 = [ 0x1000 ]; $h = 'h3';
2228 }
091c87b1 2229 while (@$x1 != 1 || $x1->[0] != 0) # _is_zero()
61f5c3f5
T
2230 {
2231 ($x1, $xr) = _div($c,$x1,$x10000);
7b29e1e6 2232 $es .= unpack($h,pack('V',$xr->[0]));
61f5c3f5
T
2233 }
2234 $es = reverse $es;
2235 $es =~ s/^[0]+//; # strip leading zeros
03874afe 2236 '0x' . $es; # return result prepended with 0x
61f5c3f5
T
2237 }
2238
2239sub _as_bin
2240 {
2241 # convert a decimal number to bin (ref to array, return ref to string)
2242 my ($c,$x) = @_;
2243
50109ad0 2244 # fits into one element (and Perl recent enough), handle also 0b0 case
091c87b1
T
2245 # handle zero case for older Perls
2246 if ($] <= 5.005 && @$x == 1 && $x->[0] == 0)
2247 {
9b924220 2248 my $t = '0b0'; return $t;
091c87b1
T
2249 }
2250 if (@$x == 1 && $] >= 5.006)
990fb837 2251 {
091c87b1 2252 my $t = sprintf("0b%b",$x->[0]);
9b924220 2253 return $t;
990fb837 2254 }
61f5c3f5
T
2255 my $x1 = _copy($c,$x);
2256
2257 my $es = '';
1ddff52a
JH
2258 my ($xr, $b, $x10000);
2259 if ($] >= 5.006)
2260 {
2261 $x10000 = [ 0x10000 ]; $b = 'b16';
2262 }
2263 else
2264 {
2265 $x10000 = [ 0x1000 ]; $b = 'b12';
2266 }
091c87b1 2267 while (!(@$x1 == 1 && $x1->[0] == 0)) # _is_zero()
61f5c3f5
T
2268 {
2269 ($x1, $xr) = _div($c,$x1,$x10000);
7b29e1e6 2270 $es .= unpack($b,pack('v',$xr->[0]));
61f5c3f5
T
2271 }
2272 $es = reverse $es;
2273 $es =~ s/^[0]+//; # strip leading zeros
03874afe 2274 '0b' . $es; # return result prepended with 0b
61f5c3f5
T
2275 }
2276
7b29e1e6
T
2277sub _as_oct
2278 {
2279 # convert a decimal number to octal (ref to array, return ref to string)
2280 my ($c,$x) = @_;
2281
50109ad0 2282 # fits into one element (handle also 0 case)
7b29e1e6
T
2283 return sprintf("0%o",$x->[0]) if @$x == 1;
2284
2285 my $x1 = _copy($c,$x);
2286
2287 my $es = '';
2288 my $xr;
2289 my $x1000 = [ 0100000 ];
2290 while (@$x1 != 1 || $x1->[0] != 0) # _is_zero()
2291 {
2292 ($x1, $xr) = _div($c,$x1,$x1000);
2293 $es .= reverse sprintf("%05o", $xr->[0]);
2294 }
2295 $es = reverse $es;
2296 $es =~ s/^[0]+//; # strip leading zeros
2297 '0' . $es; # return result prepended with 0
2298 }
2299
2300sub _from_oct
2301 {
50109ad0 2302 # convert a octal number to decimal (string, return ref to array)
7b29e1e6
T
2303 my ($c,$os) = @_;
2304
2305 # for older Perls, play safe
2306 my $m = [ 0100000 ];
2307 my $d = 5; # 5 digits at a time
2308
2309 my $mul = _one();
2310 my $x = _zero();
2311
2312 my $len = int( (length($os)-1)/$d ); # $d digit parts, w/o the '0'
2313 my $val; my $i = -$d;
2314 while ($len >= 0)
2315 {
2316 $val = substr($os,$i,$d); # get oct digits
0617f807 2317 $val = CORE::oct($val);
7b29e1e6
T
2318 $i -= $d; $len --;
2319 my $adder = [ $val ];
2320 _add ($c, $x, _mul ($c, $adder, $mul ) ) if $val != 0;
2321 _mul ($c, $mul, $m ) if $len >= 0; # skip last mul
2322 }
2323 $x;
2324 }
2325
394e6ffb 2326sub _from_hex
0716bf9b 2327 {
50109ad0 2328 # convert a hex number to decimal (string, return ref to array)
394e6ffb 2329 my ($c,$hs) = @_;
0716bf9b 2330
2d2b2744 2331 my $m = _new($c, 0x10000000); # 28 bit at a time (<32 bit!)
03874afe
T
2332 my $d = 7; # 7 digits at a time
2333 if ($] <= 5.006)
2334 {
2335 # for older Perls, play safe
2336 $m = [ 0x10000 ]; # 16 bit at a time (<32 bit!)
2337 $d = 4; # 4 digits at a time
2338 }
2339
394e6ffb 2340 my $mul = _one();
394e6ffb 2341 my $x = _zero();
0716bf9b 2342
03874afe
T
2343 my $len = int( (length($hs)-2)/$d ); # $d digit parts, w/o the '0x'
2344 my $val; my $i = -$d;
394e6ffb
JH
2345 while ($len >= 0)
2346 {
03874afe 2347 $val = substr($hs,$i,$d); # get hex digits
7b29e1e6 2348 $val =~ s/^0x// if $len == 0; # for last part only because
0617f807 2349 $val = CORE::hex($val); # hex does not like wrong chars
03874afe 2350 $i -= $d; $len --;
2d2b2744
T
2351 my $adder = [ $val ];
2352 # if the resulting number was to big to fit into one element, create a
2353 # two-element version (bug found by Mark Lakata - Thanx!)
2354 if (CORE::length($val) > $BASE_LEN)
2355 {
2356 $adder = _new($c,$val);
2357 }
2358 _add ($c, $x, _mul ($c, $adder, $mul ) ) if $val != 0;
394e6ffb
JH
2359 _mul ($c, $mul, $m ) if $len >= 0; # skip last mul
2360 }
2361 $x;
2362 }
2363
2364sub _from_bin
0716bf9b 2365 {
50109ad0 2366 # convert a hex number to decimal (string, return ref to array)
394e6ffb 2367 my ($c,$bs) = @_;
0716bf9b 2368
091c87b1 2369 # instead of converting X (8) bit at a time, it is faster to "convert" the
13a12e00
JH
2370 # number to hex, and then call _from_hex.
2371
9b924220 2372 my $hs = $bs;
13a12e00
JH
2373 $hs =~ s/^[+-]?0b//; # remove sign and 0b
2374 my $l = length($hs); # bits
2375 $hs = '0' x (8-($l % 8)) . $hs if ($l % 8) != 0; # padd left side w/ 0
03874afe 2376 my $h = '0x' . unpack('H*', pack ('B*', $hs)); # repack as hex
091c87b1 2377
03874afe 2378 $c->_from_hex($h);
0716bf9b
JH
2379 }
2380
07d34614
T
2381##############################################################################
2382# special modulus functions
2383
56d9de68 2384sub _modinv
d614cd8b 2385 {
de1ac46b 2386 # modular multiplicative inverse
56d9de68 2387 my ($c,$x,$y) = @_;
1ddff52a 2388
de1ac46b
PJA
2389 # modulo zero
2390 if (_is_zero($c, $y)) {
2391 return (undef, undef);
2392 }
2393
2394 # modulo one
2395 if (_is_one($c, $y)) {
2396 return (_zero($c), '+');
2397 }
2398
2399 my $u = _zero($c);
2400 my $v = _one($c);
2401 my $a = _copy($c,$y);
2402 my $b = _copy($c,$x);
1ddff52a 2403
de1ac46b
PJA
2404 # Euclid's Algorithm for bgcd(), only that we calc bgcd() ($a) and the result
2405 # ($u) at the same time. See comments in BigInt for why this works.
56d9de68 2406 my $q;
56d9de68 2407 my $sign = 1;
de1ac46b
PJA
2408 {
2409 ($a, $q, $b) = ($b, _div($c, $a, $b)); # step 1
2410 last if _is_zero($c, $b);
2411
2412 my $t = _add($c, # step 2:
2413 _mul($c, _copy($c, $v), $q) , # t = v * q
2414 $u ); # + u
2415 $u = $v; # u = v
2416 $v = $t; # v = t
2417 $sign = -$sign;
2418 redo;
2419 }
1ddff52a
JH
2420
2421 # if the gcd is not 1, then return NaN
de1ac46b
PJA
2422 return (undef, undef) unless _is_one($c, $a);
2423
2424 ($v, $sign == 1 ? '+' : '-');
d614cd8b
JH
2425 }
2426
2427sub _modpow
2428 {
2429 # modulus of power ($x ** $y) % $z
07d34614
T
2430 my ($c,$num,$exp,$mod) = @_;
2431
ddb27680 2432 # a^b (mod 1) = 0 for all a and b
07d34614
T
2433 if (_is_one($c,$mod))
2434 {
ddb27680
PJA
2435 @$num = 0;
2436 return $num;
07d34614 2437 }
3a7a056a
PJA
2438
2439 # 0^a (mod m) = 0 if m != 0, a != 0
2440 # 0^0 (mod m) = 1 if m != 0
d7ccd3f3 2441 if (_is_zero($c, $num)) {
3a7a056a
PJA
2442 if (_is_zero($c, $exp)) {
2443 @$num = 1;
2444 } else {
2445 @$num = 0;
2446 }
2447 return $num;
2448 }
1ddff52a
JH
2449
2450# $num = _mod($c,$num,$mod); # this does not make it faster
07d34614
T
2451
2452 my $acc = _copy($c,$num); my $t = _one();
2453
9b924220 2454 my $expbin = _as_bin($c,$exp); $expbin =~ s/^0b//;
1ddff52a
JH
2455 my $len = length($expbin);
2456 while (--$len >= 0)
07d34614 2457 {
1ddff52a 2458 if ( substr($expbin,$len,1) eq '1') # is_odd
07d34614
T
2459 {
2460 _mul($c,$t,$acc);
2461 $t = _mod($c,$t,$mod);
2462 }
2463 _mul($c,$acc,$acc);
2464 $acc = _mod($c,$acc,$mod);
07d34614
T
2465 }
2466 @$num = @$t;
2467 $num;
d614cd8b
JH
2468 }
2469
de1ac46b
PJA
2470sub _gcd {
2471 # Greatest common divisor.
9b924220 2472
de1ac46b
PJA
2473 my ($c, $x, $y) = @_;
2474
2475 # gcd(0,0) = 0
2476 # gcd(0,a) = a, if a != 0
2477
2478 if (@$x == 1 && $x->[0] == 0) {
2479 if (@$y == 1 && $y->[0] == 0) {
2480 @$x = 0;
2481 } else {
2482 @$x = @$y;
2483 }
2484 return $x;
9b924220 2485 }
de1ac46b
PJA
2486
2487 # Until $y is zero ...
2488
2489 until (@$y == 1 && $y->[0] == 0) {
2490
2491 # Compute remainder.
2492
2493 _mod($c, $x, $y);
2494
2495 # Swap $x and $y.
2496
2497 my $tmp = [ @$x ];
2498 @$x = @$y;
2499 $y = $tmp; # no deref here; that would modify input $y
2500 }
2501
2502 return $x;
2503}
9b924220 2504
394e6ffb
JH
2505##############################################################################
2506##############################################################################
2507
0716bf9b
JH
25081;
2509__END__
2510
db2e1fb3 2511=pod
2512
0716bf9b
JH
2513=head1 NAME
2514
2515Math::BigInt::Calc - Pure Perl module to support Math::BigInt
2516
2517=head1 SYNOPSIS
2518
db2e1fb3 2519This library provides support for big integer calculations. It is not
2520intended to be used by other modules. Other modules which support the same
2521API (see below) can also be used to support Math::BigInt, like
2522Math::BigInt::GMP and Math::BigInt::Pari.
0716bf9b
JH
2523
2524=head1 DESCRIPTION
2525
db2e1fb3 2526In this library, the numbers are represented in base B = 10**N, where N is
2527the largest possible value that does not cause overflow in the intermediate
2528computations. The base B elements are stored in an array, with the least
2529significant element stored in array element zero. There are no leading zero
2530elements, except a single zero element when the number is zero.
2531
2532For instance, if B = 10000, the number 1234567890 is represented internally
2533as [3456, 7890, 12].
2534
2535=head1 THE Math::BigInt API
2536
027dc388 2537In order to allow for multiple big integer libraries, Math::BigInt was
db2e1fb3 2538rewritten to use a plug-in library for core math routines. Any module which
2539conforms to the API can be used by Math::BigInt by using this in your program:
0716bf9b 2540
ee15d750 2541 use Math::BigInt lib => 'libname';
0716bf9b 2542
db2e1fb3 2543'libname' is either the long name, like 'Math::BigInt::Pari', or only the short
2544version, like 'Pari'.
2545
2546=head2 General Notes
2547
2548A library only needs to deal with unsigned big integers. Testing of input
2549parameter validity is done by the caller, so there is no need to worry about
2550underflow (e.g., in C<_sub()> and C<_dec()>) nor about division by zero (e.g.,
2551in C<_div()>) or similar cases.
2552
2553For some methods, the first parameter can be modified. That includes the
2554possibility that you return a reference to a completely different object
2555instead. Although keeping the reference and just changing its contents is
2556preferred over creating and returning a different reference.
2557
2558Return values are always objects, strings, Perl scalars, or true/false for
2559comparison routines.
2560
2561=head2 API version 1
2562
2563The following methods must be defined in order to support the use by
2564Math::BigInt v1.70 or later.
2565
2566=head3 API version
2567
2568=over 4
2569
2570=item I<api_version()>
2571
2572Return API version as a Perl scalar, 1 for Math::BigInt v1.70, 2 for
2573Math::BigInt v1.83.
2574
2575=back
2576
2577=head3 Constructors
2578
2579=over 4
2580
2581=item I<_new(STR)>
2582
2583Convert a string representing an unsigned decimal number to an object
2584representing the same number. The input is normalize, i.e., it matches
2585C<^(0|[1-9]\d*)$>.
2586
2587=item I<_zero()>
2588
2589Return an object representing the number zero.
2590
2591=item I<_one()>
2592
2593Return an object representing the number one.
2594
2595=item I<_two()>
2596
2597Return an object representing the number two.
2598
2599=item I<_ten()>
2600
2601Return an object representing the number ten.
2602
2603=item I<_from_bin(STR)>
2604
2605Return an object given a string representing a binary number. The input has a
2606'0b' prefix and matches the regular expression C<^0[bB](0|1[01]*)$>.
2607
2608=item I<_from_oct(STR)>
2609
2610Return an object given a string representing an octal number. The input has a
2611'0' prefix and matches the regular expression C<^0[1-7]*$>.
2612
2613=item I<_from_hex(STR)>
2614
2615Return an object given a string representing a hexadecimal number. The input
2616has a '0x' prefix and matches the regular expression
2617C<^0x(0|[1-9a-fA-F][\da-fA-F]*)$>.
2618
2619=back
2620
2621=head3 Mathematical functions
2622
2623Each of these methods may modify the first input argument, except I<_bgcd()>,
2624which shall not modify any input argument, and I<_sub()> which may modify the
2625second input argument.
2626
2627=over 4
2628
2629=item I<_add(OBJ1, OBJ2)>
2630
2631Returns the result of adding OBJ2 to OBJ1.
2632
2633=item I<_mul(OBJ1, OBJ2)>
2634
2635Returns the result of multiplying OBJ2 and OBJ1.
2636
2637=item I<_div(OBJ1, OBJ2)>
2638
2639Returns the result of dividing OBJ1 by OBJ2 and truncating the result to an
2640integer.
2641
2642=item I<_sub(OBJ1, OBJ2, FLAG)>
2643
2644=item I<_sub(OBJ1, OBJ2)>
2645
2646Returns the result of subtracting OBJ2 by OBJ1. If C<flag> is false or omitted,
2647OBJ1 might be modified. If C<flag> is true, OBJ2 might be modified.
2648
2649=item I<_dec(OBJ)>
2650
2651Decrement OBJ by one.
2652
2653=item I<_inc(OBJ)>
2654
2655Increment OBJ by one.
2656
2657=item I<_mod(OBJ1, OBJ2)>
2658
2659Return OBJ1 modulo OBJ2, i.e., the remainder after dividing OBJ1 by OBJ2.
2660
2661=item I<_sqrt(OBJ)>
2662
2663Return the square root of the object, truncated to integer.
2664
2665=item I<_root(OBJ, N)>
2666
2667Return Nth root of the object, truncated to int. N is E<gt>= 3.
2668
2669=item I<_fac(OBJ)>
2670
2671Return factorial of object (1*2*3*4*...).
2672
2673=item I<_pow(OBJ1, OBJ2)>
2674
2675Return OBJ1 to the power of OBJ2. By convention, 0**0 = 1.
2676
2677=item I<_modinv(OBJ1, OBJ2)>
2678
2679Return modular multiplicative inverse, i.e., return OBJ3 so that
2680
2681 (OBJ3 * OBJ1) % OBJ2 = 1 % OBJ2
2682
2683The result is returned as two arguments. If the modular multiplicative
2684inverse does not exist, both arguments are undefined. Otherwise, the
2685arguments are a number (object) and its sign ("+" or "-").
2686
2687The output value, with its sign, must either be a positive value in the
2688range 1,2,...,OBJ2-1 or the same value subtracted OBJ2. For instance, if the
2689input arguments are objects representing the numbers 7 and 5, the method
2690must either return an object representing the number 3 and a "+" sign, since
2691(3*7) % 5 = 1 % 5, or an object representing the number 2 and "-" sign,
2692since (-2*7) % 5 = 1 % 5.
2693
2694=item I<_modpow(OBJ1, OBJ2, OBJ3)>
2695
2696Return modular exponentiation, (OBJ1 ** OBJ2) % OBJ3.
2697
2698=item I<_rsft(OBJ, N, B)>
2699
2700Shift object N digits right in base B and return the resulting object. This is
2701equivalent to performing integer division by B**N and discarding the remainder,
2702except that it might be much faster, depending on how the number is represented
2703internally.
2704
2705For instance, if the object $obj represents the hexadecimal number 0xabcde,
2706then C<_rsft($obj, 2, 16)> returns an object representing the number 0xabc. The
2707"remainer", 0xde, is discarded and not returned.
2708
2709=item I<_lsft(OBJ, N, B)>
2710
2711Shift the object N digits left in base B. This is equivalent to multiplying by
2712B**N, except that it might be much faster, depending on how the number is
2713represented internally.
2714
2715=item I<_log_int(OBJ, B)>
2716
2717Return integer log of OBJ to base BASE. This method has two output arguments,
2718the OBJECT and a STATUS. The STATUS is Perl scalar; it is 1 if OBJ is the exact
2719result, 0 if the result was truncted to give OBJ, and undef if it is unknown
2720whether OBJ is the exact result.
2721
2722=item I<_gcd(OBJ1, OBJ2)>
2723
2724Return the greatest common divisor of OBJ1 and OBJ2.
2725
2726=back
2727
2728=head3 Bitwise operators
2729
2730Each of these methods may modify the first input argument.
2731
2732=over 4
2733
2734=item I<_and(OBJ1, OBJ2)>
2735
2736Return bitwise and. If necessary, the smallest number is padded with leading
2737zeros.
2738
2739=item I<_or(OBJ1, OBJ2)>
2740
2741Return bitwise or. If necessary, the smallest number is padded with leading
2742zeros.
2743
2744=item I<_xor(OBJ1, OBJ2)>
2745
2746Return bitwise exclusive or. If necessary, the smallest number is padded
2747with leading zeros.
2748
2749=back
2750
2751=head3 Boolean operators
2752
2753=over 4
2754
2755=item I<_is_zero(OBJ)>
2756
2757Returns a true value if OBJ is zero, and false value otherwise.
2758
2759=item I<_is_one(OBJ)>
2760
2761Returns a true value if OBJ is one, and false value otherwise.
2762
2763=item I<_is_two(OBJ)>
2764
2765Returns a true value if OBJ is two, and false value otherwise.
2766
2767=item I<_is_ten(OBJ)>
2768
2769Returns a true value if OBJ is ten, and false value otherwise.
2770
2771=item I<_is_even(OBJ)>
2772
2773Return a true value if OBJ is an even integer, and a false value otherwise.
2774
2775=item I<_is_odd(OBJ)>
2776
2777Return a true value if OBJ is an even integer, and a false value otherwise.
2778
2779=item I<_acmp(OBJ1, OBJ2)>
2780
2781Compare OBJ1 and OBJ2 and return -1, 0, or 1, if OBJ1 is less than, equal
2782to, or larger than OBJ2, respectively.
2783
2784=back
2785
2786=head3 String conversion
2787
2788=over 4
2789
2790=item I<_str(OBJ)>
2791
2792Return a string representing the object. The returned string should have no
2793leading zeros, i.e., it should match C<^(0|[1-9]\d*)$>.
2794
2795=item I<_as_bin(OBJ)>
2796
2797Return the binary string representation of the number. The string must have a
2798'0b' prefix.
2799
2800=item I<_as_oct(OBJ)>
2801
2802Return the octal string representation of the number. The string must have
2803a '0x' prefix.
2804
2805Note: This method was required from Math::BigInt version 1.78, but the required
2806API version number was not incremented, so there are older libraries that
2807support API version 1, but do not support C<_as_oct()>.
2808
2809=item I<_as_hex(OBJ)>
2810
2811Return the hexadecimal string representation of the number. The string must
2812have a '0x' prefix.
2813
2814=back
2815
2816=head3 Numeric conversion
2817
2818=over 4
2819
2820=item I<_num(OBJ)>
2821
2822Given an object, return a Perl scalar number (int/float) representing this
2823number.
2824
2825=back
2826
2827=head3 Miscellaneous
2828
2829=over 4
2830
2831=item I<_copy(OBJ)>
2832
2833Return a true copy of the object.
2834
2835=item I<_len(OBJ)>
2836
2837Returns the number of the decimal digits in the number. The output is a
2838Perl scalar.
2839
2840=item I<_zeros(OBJ)>
2841
2842Return the number of trailing decimal zeros. The output is a Perl scalar.
2843
2844=item I<_digit(OBJ, N)>
2845
2846Return the Nth digit as a Perl scalar. N is a Perl scalar, where zero refers to
2847the rightmost (least significant) digit, and negative values count from the
2848left (most significant digit). If $obj represents the number 123, then
2849I<_digit($obj, 0)> is 3 and I<_digit(123, -1)> is 1.
2850
2851=item I<_check(OBJ)>
2852
2853Return a true value if the object is OK, and a false value otherwise. This is a
2854check routine to test the internal state of the object for corruption.
2855
2856=back
2857
2858=head2 API version 2
2859
2860The following methods are required for an API version of 2 or greater.
2861
2862=head3 Constructors
2863
2864=over 4
2865
2866=item I<_1ex(N)>
2867
2868Return an object representing the number 10**N where N E<gt>= 0 is a Perl
2869scalar.
2870
2871=back
2872
2873=head3 Mathematical functions
2874
2875=over 4
2876
2877=item I<_nok(OBJ1, OBJ2)>
2878
2879Return the binomial coefficient OBJ1 over OBJ1.
2880
2881=back
2882
2883=head3 Miscellaneous
2884
2885=over 4
2886
2887=item I<_alen(OBJ)>
2888
2889Return the approximate number of decimal digits of the object. The
2890output is one Perl scalar. This estimate must be greater than or equal
2891to what C<_len()> returns.
2892
2893=back
2894
2895=head2 API optional methods
2896
2897The following methods are optional, and can be defined if the underlying lib
9b924220
RGS
2898has a fast way to do them. If undefined, Math::BigInt will use pure Perl (hence
2899slow) fallback routines to emulate these:
53ea20b1 2900
db2e1fb3 2901=head3 Signed bitwise operators.
9b924220 2902
db2e1fb3 2903Each of these methods may modify the first input argument.
0716bf9b 2904
db2e1fb3 2905=over 4
574bacfe 2906
db2e1fb3 2907=item I<_signed_or(OBJ1, OBJ2, SIGN1, SIGN2)>
574bacfe 2908
db2e1fb3 2909Return the signed bitwise or.
2910
2911=item I<_signed_and(OBJ1, OBJ2, SIGN1, SIGN2)>
2912
2913Return the signed bitwise and.
2914
2915=item I<_signed_xor(OBJ1, OBJ2, SIGN1, SIGN2)>
2916
2917Return the signed bitwise exclusive or.
2918
2919=back
990fb837 2920
574bacfe
JH
2921=head1 WRAP YOUR OWN
2922
2923If you want to port your own favourite c-lib for big numbers to the
2924Math::BigInt interface, you can take any of the already existing modules as
2925a rough guideline. You should really wrap up the latest BigInt and BigFloat
bd05a461 2926testsuites with your module, and replace in them any of the following:
574bacfe
JH
2927
2928 use Math::BigInt;
2929
bd05a461 2930by this:
574bacfe
JH
2931
2932 use Math::BigInt lib => 'yourlib';
2933
2934This way you ensure that your library really works 100% within Math::BigInt.
0716bf9b 2935
945313f0
PJA
2936=head1 BUGS
2937
2938Please report any bugs or feature requests to
2939C<bug-math-bigint at rt.cpan.org>, or through the web interface at
2940L<https://rt.cpan.org/Ticket/Create.html?Queue=Math-BigInt>
2941(requires login).
2942We will be notified, and then you'll automatically be notified of progress on
2943your bug as I make changes.
2944
2945=head1 SUPPORT
2946
2947You can find documentation for this module with the perldoc command.
2948
2949 perldoc Math::BigInt::Calc
2950
2951You can also look for information at:
2952
2953=over 4
2954
2955=item * RT: CPAN's request tracker
2956
2957L<https://rt.cpan.org/Public/Dist/Display.html?Name=Math-BigInt>
2958
2959=item * AnnoCPAN: Annotated CPAN documentation
2960
2961L<http://annocpan.org/dist/Math-BigInt>
2962
2963=item * CPAN Ratings
2964
2965L<http://cpanratings.perl.org/dist/Math-BigInt>
2966
2967=item * Search CPAN
2968
2969L<http://search.cpan.org/dist/Math-BigInt/>
2970
2971=item * CPAN Testers Matrix
2972
2973L<http://matrix.cpantesters.org/?dist=Math-BigInt>
2974
2975=item * The Bignum mailing list
2976
2977=over 4
2978
2979=item * Post to mailing list
2980
2981C<bignum at lists.scsys.co.uk>
2982
2983=item * View mailing list
2984
2985L<http://lists.scsys.co.uk/pipermail/bignum/>
2986
2987=item * Subscribe/Unsubscribe
2988
2989L<http://lists.scsys.co.uk/cgi-bin/mailman/listinfo/bignum>
2990
2991=back
2992
2993=back
2994
0716bf9b 2995=head1 LICENSE
53ea20b1 2996
0716bf9b
JH
2997This program is free software; you may redistribute it and/or modify it under
2998the same terms as Perl itself.
2999
3000=head1 AUTHORS
3001
db2e1fb3 3002=over 4
3003
3004=item *
3005
0716bf9b 3006Original math code by Mark Biggar, rewritten by Tels L<http://bloodgate.com/>
990fb837 3007in late 2000.
db2e1fb3 3008
3009=item *
3010
df0693ed 3011Separated from BigInt and shaped API with the help of John Peacock.
b68b7ab1 3012
db2e1fb3 3013=item *
3014
7b29e1e6 3015Fixed, speed-up, streamlined and enhanced by Tels 2001 - 2007.
0716bf9b 3016
db2e1fb3 3017=item *
3018
3019API documentation corrected and extended by Peter John Acklam,
3020E<lt>pjacklam@online.noE<gt>
3021
3022=back
3023
0716bf9b
JH
3024=head1 SEE ALSO
3025
50109ad0 3026L<Math::BigInt>, L<Math::BigFloat>,
990fb837 3027L<Math::BigInt::GMP>, L<Math::BigInt::FastCalc> and L<Math::BigInt::Pari>.
0716bf9b
JH
3028
3029=cut