Bump all Math::BigInt $VERSIONs as per 945313f0ae following a0732aaa4b
[perl.git] / dist / Math-BigInt / lib / Math / BigInt / Calc.pm
1 package Math::BigInt::Calc;
2
3 use 5.006002;
4 use strict;
5 # use warnings; # do not use warnings for older Perls
6
7 our $VERSION = '1.9997';
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
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
14
15 # todo:
16 # - fully remove funky $# stuff in div() (maybe - that code scares me...)
17
18 # USE_MUL: due to problems on certain os (os390, posix-bc) "* 1e-5" is used
19 # instead of "/ 1e5" at some places, (marked with USE_MUL). Other platforms
20 # BS2000, some Crays need USE_DIV instead.
21 # The BEGIN block is used to determine which of the two variants gives the
22 # correct result.
23
24 # Beware of things like:
25 # $i = $i * $y + $car; $car = int($i / $BASE); $i = $i % $BASE;
26 # This works on x86, but fails on ARM (SA1100, iPAQ) due to who knows what
27 # reasons. So, use this instead (slower, but correct):
28 # $i = $i * $y + $car; $car = int($i / $BASE); $i -= $BASE * $car;
29
30 ##############################################################################
31 # global constants, flags and accessory
32
33 # announce that we are compatible with MBI v1.83 and up
34 sub api_version () { 2; }
35  
36 # constants for easier life
37 my ($BASE,$BASE_LEN,$RBASE,$MAX_VAL);
38 my ($AND_BITS,$XOR_BITS,$OR_BITS);
39 my ($AND_MASK,$XOR_MASK,$OR_MASK);
40
41 sub _base_len 
42   {
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:
46   shift;
47
48   my ($b, $int) = @_;
49   if (defined $b)
50     {
51     # avoid redefinitions
52     undef &_mul;
53     undef &_div;
54
55     if ($] >= 5.008 && $int && $b > 7)
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;
63       return ($BASE_LEN, $BASE, $AND_BITS, $XOR_BITS, $OR_BITS, $BASE_LEN, $MAX_VAL,);
64       }
65
66     # find whether we can use mul or div in mul()/div()
67     $BASE_LEN = $b+1;
68     my $caught = 0;
69     while (--$BASE_LEN > 5)
70       {
71       $BASE = int("1e".$BASE_LEN);
72       $RBASE = abs('1e-'.$BASE_LEN);                    # see USE_MUL
73       $caught = 0;
74       $caught += 1 if (int($BASE * $RBASE) != 1);       # should be 1
75       $caught += 2 if (int($BASE / $BASE) != 1);        # should be 1
76       last if $caught != 3;
77       }
78     $BASE = int("1e".$BASE_LEN);
79     $RBASE = abs('1e-'.$BASE_LEN);                      # see USE_MUL
80     $MAX_VAL = $BASE-1;
81    
82     # ($caught & 1) != 0 => cannot use MUL
83     # ($caught & 2) != 0 => cannot use DIV
84     if ($caught == 2)                           # 2
85       {
86       # must USE_MUL since we cannot use DIV
87       *_mul = \&_mul_use_mul;
88       *_div = \&_div_use_mul;
89       }
90     else                                        # 0 or 1
91       {
92       # can USE_DIV instead
93       *_mul = \&_mul_use_div;
94       *_div = \&_div_use_div;
95       }
96     }
97   return $BASE_LEN unless wantarray;
98   return ($BASE_LEN, $BASE, $AND_BITS, $XOR_BITS, $OR_BITS, $BASE_LEN, $MAX_VAL);
99   }
100
101 sub _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
116 BEGIN
117   {
118   # from Daniel Pfeiffer: determine largest group of digits that is precisely
119   # multipliable with itself plus carry
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
122   do 
123     {
124     $num = ('9' x ++$e) + 0;
125     $num *= $num + 1.0;
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)
133
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
153
154   use integer;
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.
157   local $^W = 0;        # don't warn about 'nonportable number'
158   $AND_BITS = 15; $XOR_BITS = 15; $OR_BITS = 15;
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++; }
164   {
165     no integer;
166     $max = 16 if $] < 5.006;    # older Perls might not take >16 too well
167   }
168   my ($x,$y,$z);
169   do {
170     $AND_BITS++;
171     $x = CORE::oct('0b' . '1' x $AND_BITS); $y = $x & $x;
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++;
177     $x = CORE::oct('0b' . '1' x $XOR_BITS); $y = $x ^ 0;
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++;
183     $x = CORE::oct('0b' . '1' x $OR_BITS); $y = $x | $x;
184     $z = (2 ** $OR_BITS) - 1;
185     } while ($OR_BITS < $max && $x == $z && $y == $x);
186   $OR_BITS --;                                          # retreat one step
187   
188   $AND_MASK = __PACKAGE__->_new( ( 2 ** $AND_BITS ));
189   $XOR_MASK = __PACKAGE__->_new( ( 2 ** $XOR_BITS ));
190   $OR_MASK = __PACKAGE__->_new( ( 2 ** $OR_BITS ));
191
192   # We can compute the approximate length no faster than the real length:
193   *_alen = \&_len;
194   }
195
196 ###############################################################################
197
198 sub _zero
199   {
200   # create a zero
201   [ 0 ];
202   }
203
204 sub _one
205   {
206   # create a one
207   [ 1 ];
208   }
209
210 sub _two
211   {
212   # create a two (used internally for shifting)
213   [ 2 ];
214   }
215
216 sub _ten
217   {
218   # create a 10 (used internally for shifting)
219   [ 10 ];
220   }
221
222 sub _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
232 sub _copy
233   {
234   # make a true copy
235   [ @{$_[1]} ];
236   }
237
238 # catch and throw away
239 sub import { }
240
241 ##############################################################################
242 # convert back to string and number
243
244 sub _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")
249   my $ar = $_[1];
250
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     }
257
258   my $ret = "";
259   # handle first one different to strip leading zeros from it (there are no
260   # leading zero parts in internal representation)
261   $l --; $ret .= int($ar->[$l]); $l--;
262   # Interestingly, the pre-padd method uses more time
263   # the old grep variant takes longer (14 vs. 10 sec)
264   my $z = '0' x ($BASE_LEN-1);                            
265   while ($l >= 0)
266     {
267     $ret .= substr($z.$ar->[$l],-$BASE_LEN); # fastest way I could think of
268     $l--;
269     }
270   $ret;
271   }                                                                             
272
273 sub _num
274   {
275     # Make a Perl scalar number (int/float) from a BigInt object.
276     my $x = $_[1];
277
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];
289     }
290     return $num;
291   }
292
293 ##############################################################################
294 # actual math code
295
296 sub _add
297   {
298   # (ref to int_num_array, ref to int_num_array)
299   # routine to add two base 1eX numbers
300   # stolen from Knuth Vol 2 Algorithm A pg 231
301   # there are separate routines to add and sub as per Knuth pg 233
302   # This routine modifies array x, but not y.
303  
304   my ($c,$x,$y) = @_;
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     {
309     # twice as slow as $x = [ @$y ], but nec. to retain $x as ref :(
310     @$x = @$y; return $x;               
311     }
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
315   # Trades one "$j++" for having to shift arrays
316   my $i; my $car = 0; my $j = 0;
317   for $i (@$y)
318     {
319     $x->[$j] -= $BASE if $car = (($x->[$j] += $i + $car) >= $BASE) ? 1 : 0;
320     $j++;
321     }
322   while ($car != 0)
323     {
324     $x->[$j] -= $BASE if $car = (($x->[$j] += $car) >= $BASE) ? 1 : 0; $j++;
325     }
326   $x;
327   }                                                                             
328
329 sub _inc
330   {
331   # (ref to int_num_array, ref to int_num_array)
332   # Add 1 to $x, modify $x in place
333   my ($c,$x) = @_;
334
335   for my $i (@$x)
336     {
337     return $x if (($i += 1) < $BASE);           # early out
338     $i = 0;                                     # overflow, next
339     }
340   push @$x,1 if (($x->[-1] || 0) == 0);         # last overflowed, so extend
341   $x;
342   }                                                                             
343
344 sub _dec
345   {
346   # (ref to int_num_array, ref to int_num_array)
347   # Sub 1 from $x, modify $x in place
348   my ($c,$x) = @_;
349
350   my $MAX = $BASE-1;                            # since MAX_VAL based on BASE
351   for my $i (@$x)
352     {
353     last if (($i -= 1) >= 0);                   # early out
354     $i = $MAX;                                  # underflow, next
355     }
356   pop @$x if $x->[-1] == 0 && @$x > 1;          # last underflowed (but leave 0)
357   $x;
358   }                                                                             
359
360 sub _sub
361   {
362   # (ref to int_num_array, ref to int_num_array, swap)
363   # subtract base 1eX numbers -- stolen from Knuth Vol 2 pg 232, $x > $y
364   # subtract Y from X by modifying x in place
365   my ($c,$sx,$sy,$s) = @_;
366  
367   my $car = 0; my $i; my $j = 0;
368   if (!$s)
369     {
370     for $i (@$sx)
371       {
372       last unless defined $sy->[$j] || $car;
373       $i += $BASE if $car = (($i -= ($sy->[$j] || 0) + $car) < 0); $j++;
374       }
375     # might leave leading zeros, so fix that
376     return __strip_zeros($sx);
377     }
378   for $i (@$sx)
379     {
380     # we can't do an early out if $x is < than $y, since we
381     # need to copy the high chunks from $y. Found by Bob Mathews.
382     #last unless defined $sy->[$j] || $car;
383     $sy->[$j] += $BASE
384      if $car = (($sy->[$j] = $i-($sy->[$j]||0) - $car) < 0);
385     $j++;
386     }
387   # might leave leading zeros, so fix that
388   __strip_zeros($sy);
389   }                                                                             
390
391 sub _mul_use_mul
392   {
393   # (ref to int_num_array, ref to int_num_array)
394   # multiply two numbers in internal representation
395   # modifies first arg, second need not be different from first
396   my ($c,$xv,$yv) = @_;
397
398   if (@$yv == 1)
399     {
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       {
404       if (($xv->[0] *= $yv->[0]) >= $BASE)
405          {
406          $xv->[0] = $xv->[0] - ($xv->[1] = int($xv->[0] * $RBASE)) * $BASE;
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       {
420       $i = $i * $y + $car; $car = int($i * $RBASE); $i -= $car * $BASE;
421       }
422     push @$xv, $car if $car != 0;
423     return $xv;
424     }
425   # shortcut for result $x == 0 => result = 0
426   return $xv if ( ((@$xv == 1) && ($xv->[0] == 0)) ); 
427
428   # since multiplying $x with $x fails, make copy in this case
429   $yv = [@$xv] if $xv == $yv;   # same references?
430
431   my @prod = (); my ($prod,$car,$cty,$xi,$yi);
432
433   for $xi (@$xv)
434     {
435     $car = 0; $cty = 0;
436
437     # slow variant
438 #    for $yi (@$yv)
439 #      {
440 #      $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
441 #      $prod[$cty++] =
442 #       $prod - ($car = int($prod * RBASE)) * $BASE;  # see USE_MUL
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;
450     for $yi (@$yv)
451       {
452       $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
453 ##     this is actually a tad slower
454 ##        $prod = $prod[$cty]; $prod += ($car + $xi * $yi);     # no ||0 here
455       $prod[$cty++] =
456        $prod - ($car = int($prod * $RBASE)) * $BASE;  # see USE_MUL
457       }
458     $prod[$cty] += $car if $car; # need really to check for 0?
459     $xi = shift @prod || 0;     # || 0 makes v5.005_3 happy
460     }
461   push @$xv, @prod;
462   # can't have leading zeros
463 #  __strip_zeros($xv);
464   $xv;
465   }                                                                             
466
467 sub _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
530 sub _mul_use_div
531   {
532   # (ref to int_num_array, ref to int_num_array)
533   # multiply two numbers in internal representation
534   # modifies first arg, second need not be different from first
535   my ($c,$xv,$yv) = @_;
536
537   if (@$yv == 1)
538     {
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
544       if (($xv->[0] *= $yv->[0]) >= $BASE)
545           {
546           $xv->[0] =
547               $xv->[0] - ($xv->[1] = int($xv->[0] / $BASE)) * $BASE;
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       {
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;
564       }
565     push @$xv, $car if $car != 0;
566     return $xv;
567     }
568   # shortcut for result $x == 0 => result = 0
569   return $xv if ( ((@$xv == 1) && ($xv->[0] == 0)) ); 
570
571   # since multiplying $x with $x fails, make copy in this case
572   $yv = [@$xv] if $xv == $yv;   # same references?
573
574   my @prod = (); my ($prod,$car,$cty,$xi,$yi);
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;
583       $prod[$cty++] = $prod - ($car = int($prod / $BASE)) * $BASE;
584       }
585     $prod[$cty] += $car if $car; # need really to check for 0?
586     $xi = shift @prod || 0;     # || 0 makes v5.005_3 happy
587     }
588   push @$xv, @prod;
589   # can't have leading zeros
590 #  __strip_zeros($xv);
591   $xv;
592   }                                                                             
593
594 sub _div_use_mul
595   {
596   # ref to array, ref to array, modify first array and return remainder if 
597   # in list context
598
599   # see comments in _div_use_div() for more explanations
600
601   my ($c,$x,$yorg) = @_;
602   
603   # the general div algorithm here is about O(N*N) and thus quite slow, so
604   # we first check for some special cases and use shortcuts to handle them.
605
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:
610   if (@$x == 1 && @$yorg == 1)
611     {
612     # shortcut, $yorg and $x are two small numbers
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     }
625
626   # if x has more than one, but y has only one element:
627   if (@$yorg == 1)
628     {
629     my $rem;
630     $rem = _mod($c,[ @$x ],$yorg) if wantarray;
631
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       {
637       $b = $r * $BASE + $x->[$j];
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     }
645
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       {
676       # same length, so make full compare
677
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
685       # a < 0: x < y, a == 0: x == y, a > 0: x > y
686       if ($a <= 0)
687         {
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
693         return ($x,$rem) if wantarray;
694         return $x;
695         }
696       # $x >= $y, so proceed normally
697       }
698     }
699
700   # all other cases:
701
702   my $y = [ @$yorg ];                           # always make copy to preserve
703
704   my ($car,$bar,$prd,$dd,$xi,$yi,@q,$v2,$v1,@d,$tmp,$q,$u2,$u1,$u0);
705
706   $car = $bar = $prd = 0;
707   if (($dd = int($BASE/($y->[-1]+1))) != 1) 
708     {
709     for $xi (@$x) 
710       {
711       $xi = $xi * $dd + $car;
712       $xi -= ($car = int($xi * $RBASE)) * $BASE;        # see USE_MUL
713       }
714     push(@$x, $car); $car = 0;
715     for $yi (@$y) 
716       {
717       $yi = $yi * $dd + $car;
718       $yi -= ($car = int($yi * $RBASE)) * $BASE;        # see USE_MUL
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;
733     $q = (($u0 == $v1) ? $MAX_VAL : int(($u0*$BASE+$u1)/$v1));
734     --$q while ($v2*$q > ($u0*$BASE+$u1-$q*$v1)*$BASE+$u2);
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;
741         $prd -= ($car = int($prd * $RBASE)) * $BASE;    # see USE_MUL
742         $x->[$xi] += $BASE if ($bar = (($x->[$xi] -= $prd + $bar) < 0));
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           {
749           $x->[$xi] -= $BASE
750            if ($car = (($x->[$xi] += $y->[$yi] + $car) >= $BASE));
751           }
752         }   
753       }
754     pop(@$x);
755     unshift(@q, $q);
756     }
757   if (wantarray) 
758     {
759     @d = ();
760     if ($dd != 1)  
761       {
762       $car = 0; 
763       for $xi (reverse @$x) 
764         {
765         $prd = $car * $BASE + $xi;
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;
775     my $d = \@d; 
776     __strip_zeros($x);
777     __strip_zeros($d);
778     return ($x,$d);
779     }
780   @$x = @q;
781   __strip_zeros($x);
782   $x;
783   }
784
785 sub _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;
793   # the general div algorithm here is about O(N*N) and thus quite slow, so
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
978 sub _div_use_div
979   {
980   # ref to array, ref to array, modify first array and return remainder if 
981   # in list context
982   my ($c,$x,$yorg) = @_;
983
984   # the general div algorithm here is about O(N*N) and thus quite slow, so
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:
991   if (@$x == 1 && @$yorg == 1)
992     {
993     # shortcut, $yorg and $x are two small numbers
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     }
1006   # if x has more than one, but y has only one element:
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       {
1017       $b = $r * $BASE + $x->[$j];
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     }
1025   # now x and y have more than one element
1026
1027   # check whether y has more elements than x, if yet, the result will be 0
1028   if (@$yorg > @$x)
1029     {
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?
1035     return $x;                                  # only x, which is [0] now
1036     }
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
1053
1054     if (length(int($yorg->[-1])) == length(int($x->[-1])))
1055       {
1056       # same length, so make full compare
1057
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         }
1064       # $a contains the result of the compare between X and Y
1065       # a < 0: x < y, a == 0: x == y, a > 0: x > y
1066       if ($a <= 0)
1067         {
1068         $rem = [ 0 ];                   # a = 0 => x == y => rem 0
1069         $rem = [@$x] if $a != 0;        # a < 0 => x < y => rem = x
1070         splice(@$x,1);                  # keep single element
1071         $x->[0] = 0;                    # if $a < 0
1072         $x->[0] = 1 if $a == 0;         # $x == $y
1073         return ($x,$rem) if wantarray;  # including remainder?
1074         return $x;
1075         }
1076       # $x >= $y, so proceed normally
1077
1078       }
1079     }
1080
1081   # all other cases:
1082
1083   my $y = [ @$yorg ];                           # always make copy to preserve
1084  
1085   my ($car,$bar,$prd,$dd,$xi,$yi,@q,$v2,$v1,@d,$tmp,$q,$u2,$u1,$u0);
1086
1087   $car = $bar = $prd = 0;
1088   if (($dd = int($BASE/($y->[-1]+1))) != 1) 
1089     {
1090     for $xi (@$x) 
1091       {
1092       $xi = $xi * $dd + $car;
1093       $xi -= ($car = int($xi / $BASE)) * $BASE;
1094       }
1095     push(@$x, $car); $car = 0;
1096     for $yi (@$y) 
1097       {
1098       $yi = $yi * $dd + $car;
1099       $yi -= ($car = int($yi / $BASE)) * $BASE;
1100       }
1101     }
1102   else 
1103     {
1104     push(@$x, 0);
1105     }
1106
1107   # @q will accumulate the final result, $q contains the current computed
1108   # part of the final result
1109
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;
1118     $q = (($u0 == $v1) ? $MAX_VAL : int(($u0*$BASE+$u1)/$v1));
1119     --$q while ($v2*$q > ($u0*$BASE+$u1-$q*$v1)*$BASE+$u2);
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;
1126         $prd -= ($car = int($prd / $BASE)) * $BASE;
1127         $x->[$xi] += $BASE if ($bar = (($x->[$xi] -= $prd + $bar) < 0));
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           {
1134           $x->[$xi] -= $BASE
1135            if ($car = (($x->[$xi] += $y->[$yi] + $car) >= $BASE));
1136           }
1137         }   
1138       }
1139     pop(@$x); unshift(@q, $q);
1140     }
1141   if (wantarray) 
1142     {
1143     @d = ();
1144     if ($dd != 1)  
1145       {
1146       $car = 0; 
1147       for $xi (reverse @$x) 
1148         {
1149         $prd = $car * $BASE + $xi;
1150         $car = $prd - ($tmp = int($prd / $dd)) * $dd;
1151         unshift(@d, $tmp);
1152         }
1153       }
1154     else 
1155       {
1156       @d = @$x;
1157       }
1158     @$x = @q;
1159     my $d = \@d; 
1160     __strip_zeros($x);
1161     __strip_zeros($d);
1162     return ($x,$d);
1163     }
1164   @$x = @q;
1165   __strip_zeros($x);
1166   $x;
1167   }
1168
1169 ##############################################################################
1170 # testing
1171
1172 sub _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
1177   my ($c,$cx,$cy) = @_;
1178  
1179   # shortcut for short numbers 
1180   return (($cx->[0] <=> $cy->[0]) <=> 0) 
1181    if scalar @$cx == scalar @$cy && scalar @$cx == 1;
1182
1183   # fast comp based on number of array elements (aka pseudo-length)
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])));
1189   return -1 if $lxy < 0;                                # already differs, ret
1190   return 1 if $lxy > 0;                                 # ditto
1191
1192   # manual way (abort if unequal, good for early ne)
1193   my $a; my $j = scalar @$cx;
1194   while (--$j >= 0)
1195     {
1196     last if ($a = $cx->[$j] - $cy->[$j]);
1197     }
1198   $a <=> 0;
1199   }
1200
1201 sub _len
1202   {
1203   # compute number of digits in base 10
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
1209   (@$cx-1)*$BASE_LEN+length(int($cx->[-1]));
1210   }
1211
1212 sub _digit
1213   {
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.
1216   my ($c,$x,$n) = @_;
1217
1218   my $len = _len('',$x);
1219
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);
1226   }
1227
1228 sub _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];
1234
1235   return 0 if scalar @$x == 1 && $x->[0] == 0;
1236
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
1245       $zeros += length($elem);                  # count trailing zeros
1246       last;                                     # early out
1247       }
1248     $zeros ++;                                  # real else branch: 50% slower!
1249     }
1250   $zeros;
1251   }
1252
1253 ##############################################################################
1254 # _is_* routines
1255
1256 sub _is_zero
1257   {
1258   # return true if arg is zero 
1259   (((scalar @{$_[1]} == 1) && ($_[1]->[0] == 0))) <=> 0;
1260   }
1261
1262 sub _is_even
1263   {
1264   # return true if arg is even
1265   (!($_[1]->[0] & 1)) <=> 0; 
1266   }
1267
1268 sub _is_odd
1269   {
1270   # return true if arg is odd
1271   (($_[1]->[0] & 1)) <=> 0;
1272   }
1273
1274 sub _is_one
1275   {
1276   # return true if arg is one
1277   (scalar @{$_[1]} == 1) && ($_[1]->[0] == 1) <=> 0; 
1278   }
1279
1280 sub _is_two
1281   {
1282   # return true if arg is two 
1283   (scalar @{$_[1]} == 1) && ($_[1]->[0] == 2) <=> 0; 
1284   }
1285
1286 sub _is_ten
1287   {
1288   # return true if arg is ten 
1289   (scalar @{$_[1]} == 1) && ($_[1]->[0] == 10) <=> 0; 
1290   }
1291
1292 sub __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
1302   return $s if @$s == 1;                # early out
1303
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 ###############################################################################
1318 # check routine to test internal state for corruptions
1319
1320 sub _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;
1347   0;
1348   }
1349
1350
1351 ###############################################################################
1352
1353 sub _mod
1354   {
1355   # if possible, use mod shortcut
1356   my ($c,$x,$yo) = @_;
1357
1358   # slow way since $y too big
1359   if (scalar @$yo > 1)
1360     {
1361     my ($xo,$rem) = _div($c,$x,$yo);
1362     @$x = @$rem;
1363     return $x;
1364     }
1365
1366   my $y = $yo->[0];
1367
1368   # if both are single element arrays
1369   if (scalar @$x == 1)
1370     {
1371     $x->[0] %= $y;
1372     return $x;
1373     }
1374
1375   # if @$x has more than one element, but @$y is a single element
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     }
1384   elsif ($b == 1)
1385     {
1386     # else need to go through all elements in @$x: O(N), but loop is a bit
1387     # simplified
1388     my $r = 0;
1389     foreach (@$x)
1390       {
1391       $r = ($r + $_) % $y;              # not much faster, but heh...
1392       #$r += $_ % $y; $r %= $y;
1393       }
1394     $r = 0 if $r == $y;
1395     $x->[0] = $r;
1396     }
1397   else
1398     {
1399     # else need to go through all elements in @$x: O(N)
1400     my $r = 0;
1401     my $bm = 1;
1402     foreach (@$x)
1403       {
1404       $r = ($_ * $bm + $r) % $y;
1405       $bm = ($bm * $b) % $y;
1406
1407       #$r += ($_ % $y) * $bm;
1408       #$bm *= $b;
1409       #$bm %= $y;
1410       #$r %= $y;
1411       }
1412     $r = 0 if $r == $y;
1413     $x->[0] = $r;
1414     }
1415   @$x = $x->[0];                # keep one element of @$x
1416   return $x;
1417   }
1418
1419 ##############################################################################
1420 # shifts
1421
1422 sub _rsft
1423   {
1424   my ($c,$x,$y,$n) = @_;
1425
1426   if ($n != 10)
1427     {
1428     $n = _new($c,$n); return _div($c,$x, _pow($c,$n,$y));
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
1435   my $xlen = (@$x-1)*$BASE_LEN+length(int($x->[-1]));  # len of x in digits
1436   if ($src >= $xlen or ($src == $xlen and ! defined $x->[1]))
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     }
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
1448     }
1449   else
1450     {
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)
1455       {
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++;
1463       }
1464     splice (@$x,$dst) if $dst > 0;              # kill left-over array elems
1465     pop @$x if $x->[-1] == 0 && @$x > 1;        # kill last element if 0
1466     } # else rem == 0
1467   $x;
1468   }
1469
1470 sub _lsft
1471   {
1472   my ($c,$x,$y,$n) = @_;
1473
1474   if ($n != 10)
1475     {
1476     $n = _new($c,$n); return _mul($c,$x, _pow($c,$n,$y));
1477     }
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)
1489     {
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--;
1496     }
1497   # set lowest parts to 0
1498   while ($dst >= 0) { $x->[$dst--] = 0; }
1499   # fix spurious last zero element
1500   splice @$x,-1 if $x->[-1] == 0;
1501   $x;
1502   }
1503
1504 sub _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
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
1526   my $pow2 = _one();
1527
1528   my $y_bin = _as_bin($c,$cy); $y_bin =~ s/^0b//;
1529   my $len = length($y_bin);
1530   while (--$len > 0)
1531     {
1532     _mul($c,$pow2,$cx) if substr($y_bin,$len,1) eq '1';         # is odd?
1533     _mul($c,$cx,$cx);
1534     }
1535
1536   _mul($c,$cx,$pow2);
1537   $cx;
1538   }
1539
1540 sub _nok {
1541     # Return binomial coefficient (n over k).
1542     # Given refs to arrays, return ref to array.
1543     # First input argument is modified.
1544
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.
1549
1550     {
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         }
1555     }
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;
1565     }
1566
1567     else {
1568
1569         # Make a copy of the original n, since we'll be modifying n in-place.
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 }
1602
1603 my @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
1614 sub _fac
1615   {
1616   # factorial of $x
1617   # ref to array, return ref to array
1618   my ($c,$cx) = @_;
1619
1620   if ((@$cx == 1) && ($cx->[0] <= 7))
1621     {
1622     $cx->[0] = $factorials[$cx->[0]];           # 0 => 1, 1 => 1, 2 => 2 etc.
1623     return $cx;
1624     }
1625
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
1695   # go forward until $base is exceeded
1696   # limit is either $x steps (steps == 100 means a result always too high) or
1697   # $base.
1698   my $steps = 100; $steps = $cx->[0] if @$cx == 1;
1699   my $r = 2; my $cf = 3; my $step = 2; my $last = $r;
1700   while ($r*$cf < $BASE && $step < $steps)
1701     {
1702     $last = $r; $r *= $cf++; $step++;
1703     }
1704   if ((@$cx == 1) && $step == $cx->[0])
1705     {
1706     # completely done, so keep reference to $x and return
1707     $cx->[0] = $r;
1708     return $cx;
1709     }
1710   
1711   # now we must do the left over steps
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     }
1721
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.
1727   my $zero_elements = 0;
1728
1729   # do left-over steps fit into a scalar?
1730   if (ref $n eq 'ARRAY')
1731     {
1732     # No, so use slower inc() & cmp()
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       }
1745     $step = [$step];
1746     while (_acmp($c,$step,$n) <= 0)
1747       {
1748       if ($cx->[0] == 0)
1749         {
1750         $zero_elements ++; shift @$cx;
1751         }
1752       _mul($c,$cx,$step); _inc($c,$step);
1753       }
1754     }
1755   else
1756     {
1757     # Yes, so we can speed it up slightly
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)
1765       {
1766       if ($cx->[0] == 0)
1767         {
1768         $zero_elements ++; shift @$cx;
1769         }
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       {
1788       _mul($c,$cx,[$step]); $step++;
1789       if ($cx->[0] == 0)
1790         {
1791         $zero_elements ++; shift @$cx;
1792         }
1793       }
1794     }
1795   # multiply in the zeros again
1796   unshift @$cx, (0) x $zero_elements;
1797   $cx;                  # return result
1798   }
1799
1800 #############################################################################
1801
1802 sub _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);
1812   my $cmp = _acmp($c,$x,$base); # X == BASE => 1
1813   if ($cmp == 0)
1814     {
1815     splice (@$x,1); $x->[0] = 1;
1816     return ($x,1)
1817     }
1818   # X < BASE
1819   if ($cmp < 0)
1820     {
1821     splice (@$x,1); $x->[0] = 0;
1822     return ($x,undef);
1823     }
1824
1825   my $x_org = _copy($c,$x);             # preserve x
1826   splice(@$x,1); $x->[0] = 1;           # keep ref to $x
1827
1828   # Compute a guess for the result based on:
1829   # $guess = int ( length_in_base_10(X) / ( log(base) / log(10) ) )
1830   my $len = _len($c,$x_org);
1831   my $log = log($base->[-1]) / log(10);
1832
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;
1837
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)
1855       {
1856 #      print STDERR "# big _log_int at ", _str($c,$x), "\n"; 
1857       _div($c,$trial,$base); _dec($c, $x);
1858       }
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
1882
1883   my $base_mul = _mul($c, _copy($c,$base), $base);      # $base * $base
1884
1885   while (($a = _acmp($c,$trial,$x_org)) < 0)
1886     {
1887 #    print STDERR "# small _log_int at ", _str($c,$x), "\n"; 
1888     _mul($c,$trial,$base_mul); _add($c, $x, [2]);
1889     }
1890
1891   my $exact = 1;
1892   if ($a > 0)
1893     {
1894     # overstepped the result
1895     _dec($c, $x);
1896     _div($c,$trial,$base);
1897     $a = _acmp($c,$trial,$x_org);
1898     if ($a > 0)
1899       {
1900       _dec($c, $x);
1901       }
1902     $exact = 0 if $a != 0;              # a = -1 => not exact result, a = 0 => exact
1903     }
1904   
1905   ($x,$exact);                          # return result
1906   }
1907
1908 # for debugging:
1909   use constant DEBUG => 0;
1910   my $steps = 0;
1911   sub steps { $steps };
1912
1913 sub _sqrt
1914   {
1915   # square-root of $x in place
1916   # Compute a guess of the result (by rule of thumb), then improve it via
1917   # Newton's method.
1918   my ($c,$x) = @_;
1919
1920   if (scalar @$x == 1)
1921     {
1922     # fits into one Perl scalar, so result can be computed directly
1923     $x->[0] = int(sqrt($x->[0]));
1924     return $x;
1925     } 
1926   my $y = _copy($c,$x);
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
1931   my $lastelem = $x->[-1];                                      # for guess
1932   my $elems = scalar @$x - 1;
1933   # not enough digits, but could have more?
1934   if ((length($lastelem) <= 3) && ($elems > 1))
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
1941     $lastelem = $lastelem / 10 if (length($lastelem) & 1) != $len;
1942     print "$lastelem\n" if DEBUG;
1943     }
1944
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);
1948   print "l =  $l " if DEBUG;
1949
1950   splice @$x,$l;                # keep ref($x), but modify it
1951
1952   # we make the first part of the guess not '1000...0' but int(sqrt($lastelem))
1953   # that gives us:
1954   # 14400 00000 => sqrt(14400) => guess first digits to be 120
1955   # 144000 000000 => sqrt(144000) => guess 379
1956
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;
1966
1967   # If @$x > 1, we could compute the second elem of the guess, too, to create
1968   # an even better guess. Not implemented yet. Does it improve performance?
1969   $x->[$l--] = 0 while ($l >= 0);       # all other digits of guess are zero
1970
1971   print "start x= ",_str($c,$x),"\n" if DEBUG;
1972   my $two = _two();
1973   my $last = _zero();
1974   my $lastlast = _zero();
1975   $steps = 0 if DEBUG;
1976   while (_acmp($c,$last,$x) != 0 && _acmp($c,$lastlast,$x) != 0)
1977     {
1978     $steps++ if DEBUG;
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 );
1983     print " x= ",_str($c,$x),"\n" if DEBUG;
1984     }
1985   print "\nsteps in sqrt: $steps, " if DEBUG;
1986   _dec($c,$x) if _acmp($c,$y,_mul($c,_copy($c,$x),$x)) < 0;     # overshot? 
1987   print " final ",$x->[-1],"\n" if DEBUG;
1988   $x;
1989   }
1990
1991 sub _root
1992   {
1993   # take n'th root of $x in place (n >= 3)
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       {
2005       # fits into one Perl scalar, so result can be computed directly
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]) ) );
2011       }
2012     return $x;
2013     } 
2014
2015   # we know now that X is more than one element long
2016
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);
2020   if ($b =~ /0b1(0+)$/)
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     {
2041     # trial computation by starting with 2,4,8,16 etc until we overstep
2042     my $step;
2043     my $trial = _two();
2044
2045     # while still to do more than X steps
2046     do
2047       {
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);
2064
2065     # reset step to 2
2066     $step = _two();
2067     # add two, because $trial cannot be exactly the result (otherwise we would
2068     # already have found it)
2069     _add($c, $trial, $step);
2070  
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       }
2076
2077     # hit not exactly? (overstepped)
2078     if (_acmp($c, _pow($c, _copy($c, $trial), $n), $x) > 0)
2079       {
2080       _dec($c,$trial);
2081       }
2082
2083     # hit not exactly? (overstepped)
2084     # 80 too small, 81 slightly too big, 82 too big
2085     if (_acmp($c, _pow($c, _copy($c, $trial), $n), $x) > 0)
2086       {
2087       _dec ($c, $trial); 
2088       }
2089
2090     @$x = @$trial;                      # make copy while preserving ref to $x
2091     return $x;
2092     }
2093   $x; 
2094   }
2095
2096 ##############################################################################
2097 # binary stuff
2098
2099 sub _and
2100   {
2101   my ($c,$x,$y) = @_;
2102
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
2107   
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
2122     # this is when the AND_BITS are greater than $BASE and is slower for
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; }
2128 #    _add($c,$x, _mul($c, _new( $c, ($xrr & $yrr) ), $m) );
2129     
2130     # 0+ due to '&' doesn't work in strings
2131     _add($c,$x, _mul($c, [ 0+$xr->[0] & 0+$yr->[0] ], $m) );
2132     _mul($c,$m,$mask);
2133     }
2134   $x;
2135   }
2136
2137 sub _xor
2138   {
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))
2152     {
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; }
2158     #_add($c,$x, _mul($c, _new( $c, ($xrr ^ $yrr) ), $m) );
2159
2160     # 0+ due to '^' doesn't work in strings
2161     _add($c,$x, _mul($c, [ 0+$xr->[0] ^ 0+$yr->[0] ], $m) );
2162     _mul($c,$m,$mask);
2163     }
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;
2171   }
2172
2173 sub _or
2174   {
2175   my ($c,$x,$y) = @_;
2176
2177   return $x if _acmp($c,$x,$y) == 0;            # shortcut (see _and)
2178
2179   my $m = _one(); my ($xr,$yr);
2180   my $mask = $OR_MASK;
2181
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; }
2194 #    _add($c,$x, _mul($c, _new( $c, ($xrr | $yrr) ), $m) );
2195     
2196     # 0+ due to '|' doesn't work in strings
2197     _add($c,$x, _mul($c, [ 0+$xr->[0] | 0+$yr->[0] ], $m) );
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;
2207   }
2208
2209 sub _as_hex
2210   {
2211   # convert a decimal number to hex (ref to array, return ref to string)
2212   my ($c,$x) = @_;
2213
2214   # fits into one element (handle also 0x0 case)
2215   return sprintf("0x%x",$x->[0]) if @$x == 1;
2216
2217   my $x1 = _copy($c,$x);
2218
2219   my $es = '';
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     }
2229   while (@$x1 != 1 || $x1->[0] != 0)            # _is_zero()
2230     {
2231     ($x1, $xr) = _div($c,$x1,$x10000);
2232     $es .= unpack($h,pack('V',$xr->[0]));
2233     }
2234   $es = reverse $es;
2235   $es =~ s/^[0]+//;   # strip leading zeros
2236   '0x' . $es;                                   # return result prepended with 0x
2237   }
2238
2239 sub _as_bin
2240   {
2241   # convert a decimal number to bin (ref to array, return ref to string)
2242   my ($c,$x) = @_;
2243
2244   # fits into one element (and Perl recent enough), handle also 0b0 case
2245   # handle zero case for older Perls
2246   if ($] <= 5.005 && @$x == 1 && $x->[0] == 0)
2247     {
2248     my $t = '0b0'; return $t;
2249     }
2250   if (@$x == 1 && $] >= 5.006)
2251     {
2252     my $t = sprintf("0b%b",$x->[0]);
2253     return $t;
2254     }
2255   my $x1 = _copy($c,$x);
2256
2257   my $es = '';
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     }
2267   while (!(@$x1 == 1 && $x1->[0] == 0))         # _is_zero()
2268     {
2269     ($x1, $xr) = _div($c,$x1,$x10000);
2270     $es .= unpack($b,pack('v',$xr->[0]));
2271     }
2272   $es = reverse $es;
2273   $es =~ s/^[0]+//;   # strip leading zeros
2274   '0b' . $es;                                   # return result prepended with 0b
2275   }
2276
2277 sub _as_oct
2278   {
2279   # convert a decimal number to octal (ref to array, return ref to string)
2280   my ($c,$x) = @_;
2281
2282   # fits into one element (handle also 0 case)
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
2300 sub _from_oct
2301   {
2302   # convert a octal number to decimal (string, return ref to array)
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
2317     $val = CORE::oct($val);
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
2326 sub _from_hex
2327   {
2328   # convert a hex number to decimal (string, return ref to array)
2329   my ($c,$hs) = @_;
2330
2331   my $m = _new($c, 0x10000000);                 # 28 bit at a time (<32 bit!)
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
2340   my $mul = _one();
2341   my $x = _zero();
2342
2343   my $len = int( (length($hs)-2)/$d );          # $d digit parts, w/o the '0x'
2344   my $val; my $i = -$d;
2345   while ($len >= 0)
2346     {
2347     $val = substr($hs,$i,$d);                   # get hex digits
2348     $val =~ s/^0x// if $len == 0;               # for last part only because
2349     $val = CORE::hex($val);                     # hex does not like wrong chars
2350     $i -= $d; $len --;
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;
2359     _mul ($c, $mul, $m ) if $len >= 0;          # skip last mul
2360     }
2361   $x;
2362   }
2363
2364 sub _from_bin
2365   {
2366   # convert a hex number to decimal (string, return ref to array)
2367   my ($c,$bs) = @_;
2368
2369   # instead of converting X (8) bit at a time, it is faster to "convert" the
2370   # number to hex, and then call _from_hex.
2371
2372   my $hs = $bs;
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
2376   my $h = '0x' . unpack('H*', pack ('B*', $hs));        # repack as hex
2377   
2378   $c->_from_hex($h);
2379   }
2380
2381 ##############################################################################
2382 # special modulus functions
2383
2384 sub _modinv
2385   {
2386   # modular multiplicative inverse
2387   my ($c,$x,$y) = @_;
2388
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);
2403
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.
2406   my $q;
2407   my $sign = 1;
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   }
2420
2421   # if the gcd is not 1, then return NaN
2422   return (undef, undef) unless _is_one($c, $a);
2423
2424   ($v, $sign == 1 ? '+' : '-');
2425   }
2426
2427 sub _modpow
2428   {
2429   # modulus of power ($x ** $y) % $z
2430   my ($c,$num,$exp,$mod) = @_;
2431
2432   # a^b (mod 1) = 0 for all a and b
2433   if (_is_one($c,$mod))
2434     {
2435         @$num = 0;
2436         return $num;
2437     }
2438
2439   # 0^a (mod m) = 0 if m != 0, a != 0
2440   # 0^0 (mod m) = 1 if m != 0
2441   if (_is_zero($c, $num)) {
2442       if (_is_zero($c, $exp)) {
2443           @$num = 1;
2444       } else {
2445           @$num = 0;
2446       }
2447       return $num;
2448   }
2449
2450 #  $num = _mod($c,$num,$mod);   # this does not make it faster
2451
2452   my $acc = _copy($c,$num); my $t = _one();
2453
2454   my $expbin = _as_bin($c,$exp); $expbin =~ s/^0b//;
2455   my $len = length($expbin);
2456   while (--$len >= 0)
2457     {
2458     if ( substr($expbin,$len,1) eq '1')                 # is_odd
2459       {
2460       _mul($c,$t,$acc);
2461       $t = _mod($c,$t,$mod);
2462       }
2463     _mul($c,$acc,$acc);
2464     $acc = _mod($c,$acc,$mod);
2465     }
2466   @$num = @$t;
2467   $num;
2468   }
2469
2470 sub _gcd {
2471     # Greatest common divisor.
2472
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;
2485     }
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 }
2504
2505 ##############################################################################
2506 ##############################################################################
2507
2508 1;
2509 __END__
2510
2511 =pod
2512
2513 =head1 NAME
2514
2515 Math::BigInt::Calc - Pure Perl module to support Math::BigInt
2516
2517 =head1 SYNOPSIS
2518
2519 This library provides support for big integer calculations. It is not
2520 intended to be used by other modules. Other modules which support the same
2521 API (see below) can also be used to support Math::BigInt, like
2522 Math::BigInt::GMP and Math::BigInt::Pari.
2523
2524 =head1 DESCRIPTION
2525
2526 In this library, the numbers are represented in base B = 10**N, where N is
2527 the largest possible value that does not cause overflow in the intermediate
2528 computations. The base B elements are stored in an array, with the least
2529 significant element stored in array element zero. There are no leading zero
2530 elements, except a single zero element when the number is zero.
2531
2532 For instance, if B = 10000, the number 1234567890 is represented internally
2533 as [3456, 7890, 12].
2534
2535 =head1 THE Math::BigInt API
2536
2537 In order to allow for multiple big integer libraries, Math::BigInt was
2538 rewritten to use a plug-in library for core math routines. Any module which
2539 conforms to the API can be used by Math::BigInt by using this in your program:
2540
2541         use Math::BigInt lib => 'libname';
2542
2543 'libname' is either the long name, like 'Math::BigInt::Pari', or only the short
2544 version, like 'Pari'.
2545
2546 =head2 General Notes
2547
2548 A library only needs to deal with unsigned big integers. Testing of input
2549 parameter validity is done by the caller, so there is no need to worry about
2550 underflow (e.g., in C<_sub()> and C<_dec()>) nor about division by zero (e.g.,
2551 in C<_div()>) or similar cases.
2552
2553 For some methods, the first parameter can be modified. That includes the
2554 possibility that you return a reference to a completely different object
2555 instead. Although keeping the reference and just changing its contents is
2556 preferred over creating and returning a different reference.
2557
2558 Return values are always objects, strings, Perl scalars, or true/false for
2559 comparison routines.
2560
2561 =head2 API version 1
2562
2563 The following methods must be defined in order to support the use by
2564 Math::BigInt v1.70 or later.
2565
2566 =head3 API version
2567
2568 =over 4
2569
2570 =item I<api_version()>
2571
2572 Return API version as a Perl scalar, 1 for Math::BigInt v1.70, 2 for
2573 Math::BigInt v1.83.
2574
2575 =back
2576
2577 =head3 Constructors
2578
2579 =over 4
2580
2581 =item I<_new(STR)>
2582
2583 Convert a string representing an unsigned decimal number to an object
2584 representing the same number. The input is normalize, i.e., it matches
2585 C<^(0|[1-9]\d*)$>.
2586
2587 =item I<_zero()>
2588
2589 Return an object representing the number zero.
2590
2591 =item I<_one()>
2592
2593 Return an object representing the number one.
2594
2595 =item I<_two()>
2596
2597 Return an object representing the number two.
2598
2599 =item I<_ten()>
2600
2601 Return an object representing the number ten.
2602
2603 =item I<_from_bin(STR)>
2604
2605 Return 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
2610 Return 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
2615 Return an object given a string representing a hexadecimal number. The input
2616 has a '0x' prefix and matches the regular expression
2617 C<^0x(0|[1-9a-fA-F][\da-fA-F]*)$>.
2618
2619 =back
2620
2621 =head3 Mathematical functions
2622
2623 Each of these methods may modify the first input argument, except I<_bgcd()>,
2624 which shall not modify any input argument, and I<_sub()> which may modify the
2625 second input argument.
2626
2627 =over 4
2628
2629 =item I<_add(OBJ1, OBJ2)>
2630
2631 Returns the result of adding OBJ2 to OBJ1.
2632
2633 =item I<_mul(OBJ1, OBJ2)>
2634
2635 Returns the result of multiplying OBJ2 and OBJ1.
2636
2637 =item I<_div(OBJ1, OBJ2)>
2638
2639 Returns the result of dividing OBJ1 by OBJ2 and truncating the result to an
2640 integer.
2641
2642 =item I<_sub(OBJ1, OBJ2, FLAG)>
2643
2644 =item I<_sub(OBJ1, OBJ2)>
2645
2646 Returns the result of subtracting OBJ2 by OBJ1. If C<flag> is false or omitted,
2647 OBJ1 might be modified. If C<flag> is true, OBJ2 might be modified.
2648
2649 =item I<_dec(OBJ)>
2650
2651 Decrement OBJ by one.
2652
2653 =item I<_inc(OBJ)>
2654
2655 Increment OBJ by one.
2656
2657 =item I<_mod(OBJ1, OBJ2)>
2658
2659 Return OBJ1 modulo OBJ2, i.e., the remainder after dividing OBJ1 by OBJ2.
2660
2661 =item I<_sqrt(OBJ)>
2662
2663 Return the square root of the object, truncated to integer.
2664
2665 =item I<_root(OBJ, N)>
2666
2667 Return Nth root of the object, truncated to int. N is E<gt>= 3.
2668
2669 =item I<_fac(OBJ)>
2670
2671 Return factorial of object (1*2*3*4*...).
2672
2673 =item I<_pow(OBJ1, OBJ2)>
2674
2675 Return OBJ1 to the power of OBJ2. By convention, 0**0 = 1.
2676
2677 =item I<_modinv(OBJ1, OBJ2)>
2678
2679 Return modular multiplicative inverse, i.e., return OBJ3 so that
2680
2681     (OBJ3 * OBJ1) % OBJ2 = 1 % OBJ2
2682
2683 The result is returned as two arguments. If the modular multiplicative
2684 inverse does not exist, both arguments are undefined. Otherwise, the
2685 arguments are a number (object) and its sign ("+" or "-").
2686
2687 The output value, with its sign, must either be a positive value in the
2688 range 1,2,...,OBJ2-1 or the same value subtracted OBJ2. For instance, if the
2689 input arguments are objects representing the numbers 7 and 5, the method
2690 must 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,
2692 since (-2*7) % 5 = 1 % 5.
2693
2694 =item I<_modpow(OBJ1, OBJ2, OBJ3)>
2695
2696 Return modular exponentiation, (OBJ1 ** OBJ2) % OBJ3.
2697
2698 =item I<_rsft(OBJ, N, B)>
2699
2700 Shift object N digits right in base B and return the resulting object. This is
2701 equivalent to performing integer division by B**N and discarding the remainder,
2702 except that it might be much faster, depending on how the number is represented
2703 internally.
2704
2705 For instance, if the object $obj represents the hexadecimal number 0xabcde,
2706 then 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
2711 Shift the object N digits left in base B. This is equivalent to multiplying by
2712 B**N, except that it might be much faster, depending on how the number is
2713 represented internally.
2714
2715 =item I<_log_int(OBJ, B)>
2716
2717 Return integer log of OBJ to base BASE. This method has two output arguments,
2718 the OBJECT and a STATUS. The STATUS is Perl scalar; it is 1 if OBJ is the exact
2719 result, 0 if the result was truncted to give OBJ, and undef if it is unknown
2720 whether OBJ is the exact result.
2721
2722 =item I<_gcd(OBJ1, OBJ2)>
2723
2724 Return the greatest common divisor of OBJ1 and OBJ2.
2725
2726 =back
2727
2728 =head3 Bitwise operators
2729
2730 Each of these methods may modify the first input argument.
2731
2732 =over 4
2733
2734 =item I<_and(OBJ1, OBJ2)>
2735
2736 Return bitwise and. If necessary, the smallest number is padded with leading
2737 zeros.
2738
2739 =item I<_or(OBJ1, OBJ2)>
2740
2741 Return bitwise or. If necessary, the smallest number is padded with leading
2742 zeros.
2743
2744 =item I<_xor(OBJ1, OBJ2)>
2745
2746 Return bitwise exclusive or. If necessary, the smallest number is padded
2747 with leading zeros.
2748
2749 =back
2750
2751 =head3 Boolean operators
2752
2753 =over 4
2754
2755 =item I<_is_zero(OBJ)>
2756
2757 Returns a true value if OBJ is zero, and false value otherwise.
2758
2759 =item I<_is_one(OBJ)>
2760
2761 Returns a true value if OBJ is one, and false value otherwise.
2762
2763 =item I<_is_two(OBJ)>
2764
2765 Returns a true value if OBJ is two, and false value otherwise.
2766
2767 =item I<_is_ten(OBJ)>
2768
2769 Returns a true value if OBJ is ten, and false value otherwise.
2770
2771 =item I<_is_even(OBJ)>
2772
2773 Return a true value if OBJ is an even integer, and a false value otherwise.
2774
2775 =item I<_is_odd(OBJ)>
2776
2777 Return a true value if OBJ is an even integer, and a false value otherwise.
2778
2779 =item I<_acmp(OBJ1, OBJ2)>
2780
2781 Compare OBJ1 and OBJ2 and return -1, 0, or 1, if OBJ1 is less than, equal
2782 to, or larger than OBJ2, respectively.
2783
2784 =back
2785
2786 =head3 String conversion
2787
2788 =over 4
2789
2790 =item I<_str(OBJ)>
2791
2792 Return a string representing the object. The returned string should have no
2793 leading zeros, i.e., it should match C<^(0|[1-9]\d*)$>.
2794
2795 =item I<_as_bin(OBJ)>
2796
2797 Return the binary string representation of the number. The string must have a
2798 '0b' prefix.
2799
2800 =item I<_as_oct(OBJ)>
2801
2802 Return the octal string representation of the number. The string must have
2803 a '0x' prefix.
2804
2805 Note: This method was required from Math::BigInt version 1.78, but the required
2806 API version number was not incremented, so there are older libraries that
2807 support API version 1, but do not support C<_as_oct()>.
2808
2809 =item I<_as_hex(OBJ)>
2810
2811 Return the hexadecimal string representation of the number. The string must
2812 have a '0x' prefix.
2813
2814 =back
2815
2816 =head3 Numeric conversion
2817
2818 =over 4
2819
2820 =item I<_num(OBJ)>
2821
2822 Given an object, return a Perl scalar number (int/float) representing this
2823 number.
2824
2825 =back
2826
2827 =head3 Miscellaneous
2828
2829 =over 4
2830
2831 =item I<_copy(OBJ)>
2832
2833 Return a true copy of the object.
2834
2835 =item I<_len(OBJ)>
2836
2837 Returns the number of the decimal digits in the number. The output is a
2838 Perl scalar.
2839
2840 =item I<_zeros(OBJ)>
2841
2842 Return the number of trailing decimal zeros. The output is a Perl scalar.
2843
2844 =item I<_digit(OBJ, N)>
2845
2846 Return the Nth digit as a Perl scalar. N is a Perl scalar, where zero refers to
2847 the rightmost (least significant) digit, and negative values count from the
2848 left (most significant digit). If $obj represents the number 123, then
2849 I<_digit($obj, 0)> is 3 and I<_digit(123, -1)> is 1.
2850
2851 =item I<_check(OBJ)>
2852
2853 Return a true value if the object is OK, and a false value otherwise. This is a
2854 check routine to test the internal state of the object for corruption.
2855
2856 =back
2857
2858 =head2 API version 2
2859
2860 The 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
2868 Return an object representing the number 10**N where N E<gt>= 0 is a Perl
2869 scalar.
2870
2871 =back
2872
2873 =head3 Mathematical functions
2874
2875 =over 4
2876
2877 =item I<_nok(OBJ1, OBJ2)>
2878
2879 Return the binomial coefficient OBJ1 over OBJ1.
2880
2881 =back
2882
2883 =head3 Miscellaneous
2884
2885 =over 4
2886
2887 =item I<_alen(OBJ)>
2888
2889 Return the approximate number of decimal digits of the object. The
2890 output is one Perl scalar. This estimate must be greater than or equal
2891 to what C<_len()> returns.
2892
2893 =back
2894
2895 =head2 API optional methods
2896
2897 The following methods are optional, and can be defined if the underlying lib
2898 has a fast way to do them. If undefined, Math::BigInt will use pure Perl (hence
2899 slow) fallback routines to emulate these:
2900
2901 =head3 Signed bitwise operators.
2902
2903 Each of these methods may modify the first input argument.
2904
2905 =over 4
2906
2907 =item I<_signed_or(OBJ1, OBJ2, SIGN1, SIGN2)>
2908
2909 Return the signed bitwise or.
2910
2911 =item I<_signed_and(OBJ1, OBJ2, SIGN1, SIGN2)>
2912
2913 Return the signed bitwise and.
2914
2915 =item I<_signed_xor(OBJ1, OBJ2, SIGN1, SIGN2)>
2916
2917 Return the signed bitwise exclusive or.
2918
2919 =back
2920
2921 =head1 WRAP YOUR OWN
2922
2923 If you want to port your own favourite c-lib for big numbers to the
2924 Math::BigInt interface, you can take any of the already existing modules as
2925 a rough guideline. You should really wrap up the latest BigInt and BigFloat
2926 testsuites with your module, and replace in them any of the following:
2927
2928         use Math::BigInt;
2929
2930 by this:
2931
2932         use Math::BigInt lib => 'yourlib';
2933
2934 This way you ensure that your library really works 100% within Math::BigInt.
2935
2936 =head1 BUGS
2937
2938 Please report any bugs or feature requests to
2939 C<bug-math-bigint at rt.cpan.org>, or through the web interface at
2940 L<https://rt.cpan.org/Ticket/Create.html?Queue=Math-BigInt>
2941 (requires login).
2942 We will be notified, and then you'll automatically be notified of progress on
2943 your bug as I make changes.
2944
2945 =head1 SUPPORT
2946
2947 You can find documentation for this module with the perldoc command.
2948
2949     perldoc Math::BigInt::Calc
2950
2951 You can also look for information at:
2952
2953 =over 4
2954
2955 =item * RT: CPAN's request tracker
2956
2957 L<https://rt.cpan.org/Public/Dist/Display.html?Name=Math-BigInt>
2958
2959 =item * AnnoCPAN: Annotated CPAN documentation
2960
2961 L<http://annocpan.org/dist/Math-BigInt>
2962
2963 =item * CPAN Ratings
2964
2965 L<http://cpanratings.perl.org/dist/Math-BigInt>
2966
2967 =item * Search CPAN
2968
2969 L<http://search.cpan.org/dist/Math-BigInt/>
2970
2971 =item * CPAN Testers Matrix
2972
2973 L<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
2981 C<bignum at lists.scsys.co.uk>
2982
2983 =item * View mailing list
2984
2985 L<http://lists.scsys.co.uk/pipermail/bignum/>
2986
2987 =item * Subscribe/Unsubscribe
2988
2989 L<http://lists.scsys.co.uk/cgi-bin/mailman/listinfo/bignum>
2990
2991 =back
2992
2993 =back
2994
2995 =head1 LICENSE
2996
2997 This program is free software; you may redistribute it and/or modify it under
2998 the same terms as Perl itself. 
2999
3000 =head1 AUTHORS
3001
3002 =over 4
3003
3004 =item *
3005
3006 Original math code by Mark Biggar, rewritten by Tels L<http://bloodgate.com/>
3007 in late 2000.
3008
3009 =item *
3010
3011 Separated from BigInt and shaped API with the help of John Peacock.
3012
3013 =item *
3014
3015 Fixed, speed-up, streamlined and enhanced by Tels 2001 - 2007.
3016
3017 =item *
3018
3019 API documentation corrected and extended by Peter John Acklam,
3020 E<lt>pjacklam@online.noE<gt>
3021
3022 =back
3023
3024 =head1 SEE ALSO
3025
3026 L<Math::BigInt>, L<Math::BigFloat>,
3027 L<Math::BigInt::GMP>, L<Math::BigInt::FastCalc> and L<Math::BigInt::Pari>.
3028
3029 =cut