Bumped Math-BigInt, Math-BigInt-FastCalc and Math-BigRat versions for release per...
[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; # dont use warnings for older Perls
6
7 our $VERSION = '1.99_03';
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 whoknows 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 clobbers up 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 to big
1359   if (scalar @$yo > 1)
1360     {
1361     my ($xo,$rem) = _div($c,$x,$yo);
1362     return $rem;
1363     }
1364
1365   my $y = $yo->[0];
1366   # both are single element arrays
1367   if (scalar @$x == 1)
1368     {
1369     $x->[0] %= $y;
1370     return $x;
1371     }
1372
1373   # @y is a single element, but @x has more than one element
1374   my $b = $BASE % $y;
1375   if ($b == 0)
1376     {
1377     # when BASE % Y == 0 then (B * BASE) % Y == 0
1378     # (B * BASE) % $y + A % Y => A % Y
1379     # so need to consider only last element: O(1)
1380     $x->[0] %= $y;
1381     }
1382   elsif ($b == 1)
1383     {
1384     # else need to go through all elements: O(N), but loop is a bit simplified
1385     my $r = 0;
1386     foreach (@$x)
1387       {
1388       $r = ($r + $_) % $y;              # not much faster, but heh...
1389       #$r += $_ % $y; $r %= $y;
1390       }
1391     $r = 0 if $r == $y;
1392     $x->[0] = $r;
1393     }
1394   else
1395     {
1396     # else need to go through all elements: O(N)
1397     my $r = 0; my $bm = 1;
1398     foreach (@$x)
1399       {
1400       $r = ($_ * $bm + $r) % $y;
1401       $bm = ($bm * $b) % $y;
1402
1403       #$r += ($_ % $y) * $bm;
1404       #$bm *= $b;
1405       #$bm %= $y;
1406       #$r %= $y;
1407       }
1408     $r = 0 if $r == $y;
1409     $x->[0] = $r;
1410     }
1411   splice (@$x,1);               # keep one element of $x
1412   $x;
1413   }
1414
1415 ##############################################################################
1416 # shifts
1417
1418 sub _rsft
1419   {
1420   my ($c,$x,$y,$n) = @_;
1421
1422   if ($n != 10)
1423     {
1424     $n = _new($c,$n); return _div($c,$x, _pow($c,$n,$y));
1425     }
1426
1427   # shortcut (faster) for shifting by 10)
1428   # multiples of $BASE_LEN
1429   my $dst = 0;                          # destination
1430   my $src = _num($c,$y);                # as normal int
1431   my $xlen = (@$x-1)*$BASE_LEN+length(int($x->[-1]));  # len of x in digits
1432   if ($src >= $xlen or ($src == $xlen and ! defined $x->[1]))
1433     {
1434     # 12345 67890 shifted right by more than 10 digits => 0
1435     splice (@$x,1);                    # leave only one element
1436     $x->[0] = 0;                       # set to zero
1437     return $x;
1438     }
1439   my $rem = $src % $BASE_LEN;           # remainder to shift
1440   $src = int($src / $BASE_LEN);         # source
1441   if ($rem == 0)
1442     {
1443     splice (@$x,0,$src);                # even faster, 38.4 => 39.3
1444     }
1445   else
1446     {
1447     my $len = scalar @$x - $src;        # elems to go
1448     my $vd; my $z = '0'x $BASE_LEN;
1449     $x->[scalar @$x] = 0;               # avoid || 0 test inside loop
1450     while ($dst < $len)
1451       {
1452       $vd = $z.$x->[$src];
1453       $vd = substr($vd,-$BASE_LEN,$BASE_LEN-$rem);
1454       $src++;
1455       $vd = substr($z.$x->[$src],-$rem,$rem) . $vd;
1456       $vd = substr($vd,-$BASE_LEN,$BASE_LEN) if length($vd) > $BASE_LEN;
1457       $x->[$dst] = int($vd);
1458       $dst++;
1459       }
1460     splice (@$x,$dst) if $dst > 0;              # kill left-over array elems
1461     pop @$x if $x->[-1] == 0 && @$x > 1;        # kill last element if 0
1462     } # else rem == 0
1463   $x;
1464   }
1465
1466 sub _lsft
1467   {
1468   my ($c,$x,$y,$n) = @_;
1469
1470   if ($n != 10)
1471     {
1472     $n = _new($c,$n); return _mul($c,$x, _pow($c,$n,$y));
1473     }
1474
1475   # shortcut (faster) for shifting by 10) since we are in base 10eX
1476   # multiples of $BASE_LEN:
1477   my $src = scalar @$x;                 # source
1478   my $len = _num($c,$y);                # shift-len as normal int
1479   my $rem = $len % $BASE_LEN;           # remainder to shift
1480   my $dst = $src + int($len/$BASE_LEN); # destination
1481   my $vd;                               # further speedup
1482   $x->[$src] = 0;                       # avoid first ||0 for speed
1483   my $z = '0' x $BASE_LEN;
1484   while ($src >= 0)
1485     {
1486     $vd = $x->[$src]; $vd = $z.$vd;
1487     $vd = substr($vd,-$BASE_LEN+$rem,$BASE_LEN-$rem);
1488     $vd .= $src > 0 ? substr($z.$x->[$src-1],-$BASE_LEN,$rem) : '0' x $rem;
1489     $vd = substr($vd,-$BASE_LEN,$BASE_LEN) if length($vd) > $BASE_LEN;
1490     $x->[$dst] = int($vd);
1491     $dst--; $src--;
1492     }
1493   # set lowest parts to 0
1494   while ($dst >= 0) { $x->[$dst--] = 0; }
1495   # fix spurious last zero element
1496   splice @$x,-1 if $x->[-1] == 0;
1497   $x;
1498   }
1499
1500 sub _pow
1501   {
1502   # power of $x to $y
1503   # ref to array, ref to array, return ref to array
1504   my ($c,$cx,$cy) = @_;
1505
1506   if (scalar @$cy == 1 && $cy->[0] == 0)
1507     {
1508     splice (@$cx,1); $cx->[0] = 1;              # y == 0 => x => 1
1509     return $cx;
1510     }
1511   if ((scalar @$cx == 1 && $cx->[0] == 1) ||    #    x == 1
1512       (scalar @$cy == 1 && $cy->[0] == 1))      # or y == 1
1513     {
1514     return $cx;
1515     }
1516   if (scalar @$cx == 1 && $cx->[0] == 0)
1517     {
1518     splice (@$cx,1); $cx->[0] = 0;              # 0 ** y => 0 (if not y <= 0)
1519     return $cx;
1520     }
1521
1522   my $pow2 = _one();
1523
1524   my $y_bin = _as_bin($c,$cy); $y_bin =~ s/^0b//;
1525   my $len = length($y_bin);
1526   while (--$len > 0)
1527     {
1528     _mul($c,$pow2,$cx) if substr($y_bin,$len,1) eq '1';         # is odd?
1529     _mul($c,$cx,$cx);
1530     }
1531
1532   _mul($c,$cx,$pow2);
1533   $cx;
1534   }
1535
1536 sub _nok
1537   {
1538   # n over k
1539   # ref to array, return ref to array
1540   my ($c,$n,$k) = @_;
1541
1542   # ( 7 )       7!       1*2*3*4 * 5*6*7   5 * 6 * 7       6   7
1543   # ( - ) = --------- =  --------------- = --------- = 5 * - * -
1544   # ( 3 )   (7-3)! 3!    1*2*3*4 * 1*2*3   1 * 2 * 3       2   3
1545
1546   if (!_is_zero($c,$k))
1547     {
1548     my $x = _copy($c,$n);
1549     _sub($c,$n,$k);
1550     _inc($c,$n);
1551     my $f = _copy($c,$n); _inc($c,$f);          # n = 5, f = 6, d = 2
1552     my $d = _two($c);
1553     while (_acmp($c,$f,$x) <= 0)                # f <= n ?
1554       {
1555       # n = (n * f / d) == 5 * 6 / 2
1556       $n = _mul($c,$n,$f); $n = _div($c,$n,$d);
1557       # f = 7, d = 3
1558       _inc($c,$f); _inc($c,$d);
1559       }
1560     }
1561   else
1562     {
1563     # keep ref to $n and set it to 1
1564     splice (@$n,1); $n->[0] = 1;
1565     }
1566   $n;
1567   }
1568
1569 my @factorials = (
1570   1,
1571   1,
1572   2,
1573   2*3,
1574   2*3*4,
1575   2*3*4*5,
1576   2*3*4*5*6,
1577   2*3*4*5*6*7,
1578 );
1579
1580 sub _fac
1581   {
1582   # factorial of $x
1583   # ref to array, return ref to array
1584   my ($c,$cx) = @_;
1585
1586   if ((@$cx == 1) && ($cx->[0] <= 7))
1587     {
1588     $cx->[0] = $factorials[$cx->[0]];           # 0 => 1, 1 => 1, 2 => 2 etc.
1589     return $cx;
1590     }
1591
1592   if ((@$cx == 1) &&            # we do this only if $x >= 12 and $x <= 7000
1593       ($cx->[0] >= 12 && $cx->[0] < 7000))
1594     {
1595
1596   # Calculate (k-j) * (k-j+1) ... k .. (k+j-1) * (k + j)
1597   # See http://blogten.blogspot.com/2007/01/calculating-n.html
1598   # The above series can be expressed as factors:
1599   #   k * k - (j - i) * 2
1600   # We cache k*k, and calculate (j * j) as the sum of the first j odd integers
1601
1602   # This will not work when N exceeds the storage of a Perl scalar, however,
1603   # in this case the algorithm would be way to slow to terminate, anyway.
1604
1605   # As soon as the last element of $cx is 0, we split it up and remember
1606   # how many zeors we got so far. The reason is that n! will accumulate
1607   # zeros at the end rather fast.
1608   my $zero_elements = 0;
1609
1610   # If n is even, set n = n -1
1611   my $k = _num($c,$cx); my $even = 1;
1612   if (($k & 1) == 0)
1613     {
1614     $even = $k; $k --;
1615     }
1616   # set k to the center point
1617   $k = ($k + 1) / 2;
1618 #  print "k $k even: $even\n";
1619   # now calculate k * k
1620   my $k2 = $k * $k;
1621   my $odd = 1; my $sum = 1;
1622   my $i = $k - 1;
1623   # keep reference to x
1624   my $new_x = _new($c, $k * $even);
1625   @$cx = @$new_x;
1626   if ($cx->[0] == 0)
1627     {
1628     $zero_elements ++; shift @$cx;
1629     }
1630 #  print STDERR "x = ", _str($c,$cx),"\n";
1631   my $BASE2 = int(sqrt($BASE))-1;
1632   my $j = 1; 
1633   while ($j <= $i)
1634     {
1635     my $m = ($k2 - $sum); $odd += 2; $sum += $odd; $j++;
1636     while ($j <= $i && ($m < $BASE2) && (($k2 - $sum) < $BASE2))
1637       {
1638       $m *= ($k2 - $sum);
1639       $odd += 2; $sum += $odd; $j++;
1640 #      print STDERR "\n k2 $k2 m $m sum $sum odd $odd\n"; sleep(1);
1641       }
1642     if ($m < $BASE)
1643       {
1644       _mul($c,$cx,[$m]);
1645       }
1646     else
1647       {
1648       _mul($c,$cx,$c->_new($m));
1649       }
1650     if ($cx->[0] == 0)
1651       {
1652       $zero_elements ++; shift @$cx;
1653       }
1654 #    print STDERR "Calculate $k2 - $sum = $m (x = ", _str($c,$cx),")\n";
1655     }
1656   # multiply in the zeros again
1657   unshift @$cx, (0) x $zero_elements; 
1658   return $cx;
1659   }
1660
1661   # go forward until $base is exceeded
1662   # limit is either $x steps (steps == 100 means a result always too high) or
1663   # $base.
1664   my $steps = 100; $steps = $cx->[0] if @$cx == 1;
1665   my $r = 2; my $cf = 3; my $step = 2; my $last = $r;
1666   while ($r*$cf < $BASE && $step < $steps)
1667     {
1668     $last = $r; $r *= $cf++; $step++;
1669     }
1670   if ((@$cx == 1) && $step == $cx->[0])
1671     {
1672     # completely done, so keep reference to $x and return
1673     $cx->[0] = $r;
1674     return $cx;
1675     }
1676   
1677   # now we must do the left over steps
1678   my $n;                                        # steps still to do
1679   if (scalar @$cx == 1)
1680     {
1681     $n = $cx->[0];
1682     }
1683   else
1684     {
1685     $n = _copy($c,$cx);
1686     }
1687
1688   # Set $cx to the last result below $BASE (but keep ref to $x)
1689   $cx->[0] = $last; splice (@$cx,1);
1690   # As soon as the last element of $cx is 0, we split it up and remember
1691   # how many zeors we got so far. The reason is that n! will accumulate
1692   # zeros at the end rather fast.
1693   my $zero_elements = 0;
1694
1695   # do left-over steps fit into a scalar?
1696   if (ref $n eq 'ARRAY')
1697     {
1698     # No, so use slower inc() & cmp()
1699     # ($n is at least $BASE here)
1700     my $base_2 = int(sqrt($BASE)) - 1;
1701     #print STDERR "base_2: $base_2\n"; 
1702     while ($step < $base_2)
1703       {
1704       if ($cx->[0] == 0)
1705         {
1706         $zero_elements ++; shift @$cx;
1707         }
1708       my $b = $step * ($step + 1); $step += 2;
1709       _mul($c,$cx,[$b]);
1710       }
1711     $step = [$step];
1712     while (_acmp($c,$step,$n) <= 0)
1713       {
1714       if ($cx->[0] == 0)
1715         {
1716         $zero_elements ++; shift @$cx;
1717         }
1718       _mul($c,$cx,$step); _inc($c,$step);
1719       }
1720     }
1721   else
1722     {
1723     # Yes, so we can speed it up slightly
1724   
1725 #    print "# left over steps $n\n";
1726
1727     my $base_4 = int(sqrt(sqrt($BASE))) - 2;
1728     #print STDERR "base_4: $base_4\n";
1729     my $n4 = $n - 4; 
1730     while ($step < $n4 && $step < $base_4)
1731       {
1732       if ($cx->[0] == 0)
1733         {
1734         $zero_elements ++; shift @$cx;
1735         }
1736       my $b = $step * ($step + 1); $step += 2; $b *= $step * ($step + 1); $step += 2;
1737       _mul($c,$cx,[$b]);
1738       }
1739     my $base_2 = int(sqrt($BASE)) - 1;
1740     my $n2 = $n - 2; 
1741     #print STDERR "base_2: $base_2\n"; 
1742     while ($step < $n2 && $step < $base_2)
1743       {
1744       if ($cx->[0] == 0)
1745         {
1746         $zero_elements ++; shift @$cx;
1747         }
1748       my $b = $step * ($step + 1); $step += 2;
1749       _mul($c,$cx,[$b]);
1750       }
1751     # do what's left over
1752     while ($step <= $n)
1753       {
1754       _mul($c,$cx,[$step]); $step++;
1755       if ($cx->[0] == 0)
1756         {
1757         $zero_elements ++; shift @$cx;
1758         }
1759       }
1760     }
1761   # multiply in the zeros again
1762   unshift @$cx, (0) x $zero_elements;
1763   $cx;                  # return result
1764   }
1765
1766 #############################################################################
1767
1768 sub _log_int
1769   {
1770   # calculate integer log of $x to base $base
1771   # ref to array, ref to array - return ref to array
1772   my ($c,$x,$base) = @_;
1773
1774   # X == 0 => NaN
1775   return if (scalar @$x == 1 && $x->[0] == 0);
1776   # BASE 0 or 1 => NaN
1777   return if (scalar @$base == 1 && $base->[0] < 2);
1778   my $cmp = _acmp($c,$x,$base); # X == BASE => 1
1779   if ($cmp == 0)
1780     {
1781     splice (@$x,1); $x->[0] = 1;
1782     return ($x,1)
1783     }
1784   # X < BASE
1785   if ($cmp < 0)
1786     {
1787     splice (@$x,1); $x->[0] = 0;
1788     return ($x,undef);
1789     }
1790
1791   my $x_org = _copy($c,$x);             # preserve x
1792   splice(@$x,1); $x->[0] = 1;           # keep ref to $x
1793
1794   # Compute a guess for the result based on:
1795   # $guess = int ( length_in_base_10(X) / ( log(base) / log(10) ) )
1796   my $len = _len($c,$x_org);
1797   my $log = log($base->[-1]) / log(10);
1798
1799   # for each additional element in $base, we add $BASE_LEN to the result,
1800   # based on the observation that log($BASE,10) is BASE_LEN and
1801   # log(x*y) == log(x) + log(y):
1802   $log += ((scalar @$base)-1) * $BASE_LEN;
1803
1804   # calculate now a guess based on the values obtained above:
1805   my $res = int($len / $log);
1806
1807   $x->[0] = $res;
1808   my $trial = _pow ($c, _copy($c, $base), $x);
1809   my $a = _acmp($c,$trial,$x_org);
1810
1811 #  print STDERR "# trial ", _str($c,$x)," was: $a (0 = exact, -1 too small, +1 too big)\n";
1812
1813   # found an exact result?
1814   return ($x,1) if $a == 0;
1815
1816   if ($a > 0)
1817     {
1818     # or too big
1819     _div($c,$trial,$base); _dec($c, $x);
1820     while (($a = _acmp($c,$trial,$x_org)) > 0)
1821       {
1822 #      print STDERR "# big _log_int at ", _str($c,$x), "\n"; 
1823       _div($c,$trial,$base); _dec($c, $x);
1824       }
1825     # result is now exact (a == 0), or too small (a < 0)
1826     return ($x, $a == 0 ? 1 : 0);
1827     }
1828
1829   # else: result was to small
1830   _mul($c,$trial,$base);
1831
1832   # did we now get the right result?
1833   $a = _acmp($c,$trial,$x_org);
1834
1835   if ($a == 0)                          # yes, exactly
1836     {
1837     _inc($c, $x);
1838     return ($x,1); 
1839     }
1840   return ($x,0) if $a > 0;  
1841
1842   # Result still too small (we should come here only if the estimate above
1843   # was very off base):
1844  
1845   # Now let the normal trial run obtain the real result
1846   # Simple loop that increments $x by 2 in each step, possible overstepping
1847   # the real result
1848
1849   my $base_mul = _mul($c, _copy($c,$base), $base);      # $base * $base
1850
1851   while (($a = _acmp($c,$trial,$x_org)) < 0)
1852     {
1853 #    print STDERR "# small _log_int at ", _str($c,$x), "\n"; 
1854     _mul($c,$trial,$base_mul); _add($c, $x, [2]);
1855     }
1856
1857   my $exact = 1;
1858   if ($a > 0)
1859     {
1860     # overstepped the result
1861     _dec($c, $x);
1862     _div($c,$trial,$base);
1863     $a = _acmp($c,$trial,$x_org);
1864     if ($a > 0)
1865       {
1866       _dec($c, $x);
1867       }
1868     $exact = 0 if $a != 0;              # a = -1 => not exact result, a = 0 => exact
1869     }
1870   
1871   ($x,$exact);                          # return result
1872   }
1873
1874 # for debugging:
1875   use constant DEBUG => 0;
1876   my $steps = 0;
1877   sub steps { $steps };
1878
1879 sub _sqrt
1880   {
1881   # square-root of $x in place
1882   # Compute a guess of the result (by rule of thumb), then improve it via
1883   # Newton's method.
1884   my ($c,$x) = @_;
1885
1886   if (scalar @$x == 1)
1887     {
1888     # fits into one Perl scalar, so result can be computed directly
1889     $x->[0] = int(sqrt($x->[0]));
1890     return $x;
1891     } 
1892   my $y = _copy($c,$x);
1893   # hopefully _len/2 is < $BASE, the -1 is to always undershot the guess
1894   # since our guess will "grow"
1895   my $l = int((_len($c,$x)-1) / 2);     
1896
1897   my $lastelem = $x->[-1];                                      # for guess
1898   my $elems = scalar @$x - 1;
1899   # not enough digits, but could have more?
1900   if ((length($lastelem) <= 3) && ($elems > 1))
1901     {
1902     # right-align with zero pad
1903     my $len = length($lastelem) & 1;
1904     print "$lastelem => " if DEBUG;
1905     $lastelem .= substr($x->[-2] . '0' x $BASE_LEN,0,$BASE_LEN);
1906     # former odd => make odd again, or former even to even again
1907     $lastelem = $lastelem / 10 if (length($lastelem) & 1) != $len;
1908     print "$lastelem\n" if DEBUG;
1909     }
1910
1911   # construct $x (instead of _lsft($c,$x,$l,10)
1912   my $r = $l % $BASE_LEN;       # 10000 00000 00000 00000 ($BASE_LEN=5)
1913   $l = int($l / $BASE_LEN);
1914   print "l =  $l " if DEBUG;
1915
1916   splice @$x,$l;                # keep ref($x), but modify it
1917
1918   # we make the first part of the guess not '1000...0' but int(sqrt($lastelem))
1919   # that gives us:
1920   # 14400 00000 => sqrt(14400) => guess first digits to be 120
1921   # 144000 000000 => sqrt(144000) => guess 379
1922
1923   print "$lastelem (elems $elems) => " if DEBUG;
1924   $lastelem = $lastelem / 10 if ($elems & 1 == 1);              # odd or even?
1925   my $g = sqrt($lastelem); $g =~ s/\.//;                        # 2.345 => 2345
1926   $r -= 1 if $elems & 1 == 0;                                   # 70 => 7
1927
1928   # padd with zeros if result is too short
1929   $x->[$l--] = int(substr($g . '0' x $r,0,$r+1));
1930   print "now ",$x->[-1] if DEBUG;
1931   print " would have been ", int('1' . '0' x $r),"\n" if DEBUG;
1932
1933   # If @$x > 1, we could compute the second elem of the guess, too, to create
1934   # an even better guess. Not implemented yet. Does it improve performance?
1935   $x->[$l--] = 0 while ($l >= 0);       # all other digits of guess are zero
1936
1937   print "start x= ",_str($c,$x),"\n" if DEBUG;
1938   my $two = _two();
1939   my $last = _zero();
1940   my $lastlast = _zero();
1941   $steps = 0 if DEBUG;
1942   while (_acmp($c,$last,$x) != 0 && _acmp($c,$lastlast,$x) != 0)
1943     {
1944     $steps++ if DEBUG;
1945     $lastlast = _copy($c,$last);
1946     $last = _copy($c,$x);
1947     _add($c,$x, _div($c,_copy($c,$y),$x));
1948     _div($c,$x, $two );
1949     print " x= ",_str($c,$x),"\n" if DEBUG;
1950     }
1951   print "\nsteps in sqrt: $steps, " if DEBUG;
1952   _dec($c,$x) if _acmp($c,$y,_mul($c,_copy($c,$x),$x)) < 0;     # overshot? 
1953   print " final ",$x->[-1],"\n" if DEBUG;
1954   $x;
1955   }
1956
1957 sub _root
1958   {
1959   # take n'th root of $x in place (n >= 3)
1960   my ($c,$x,$n) = @_;
1961  
1962   if (scalar @$x == 1)
1963     {
1964     if (scalar @$n > 1)
1965       {
1966       # result will always be smaller than 2 so trunc to 1 at once
1967       $x->[0] = 1;
1968       }
1969     else
1970       {
1971       # fits into one Perl scalar, so result can be computed directly
1972       # cannot use int() here, because it rounds wrongly (try 
1973       # (81 ** 3) ** (1/3) to see what I mean)
1974       #$x->[0] = int( $x->[0] ** (1 / $n->[0]) );
1975       # round to 8 digits, then truncate result to integer
1976       $x->[0] = int ( sprintf ("%.8f", $x->[0] ** (1 / $n->[0]) ) );
1977       }
1978     return $x;
1979     } 
1980
1981   # we know now that X is more than one element long
1982
1983   # if $n is a power of two, we can repeatedly take sqrt($X) and find the
1984   # proper result, because sqrt(sqrt($x)) == root($x,4)
1985   my $b = _as_bin($c,$n);
1986   if ($b =~ /0b1(0+)$/)
1987     {
1988     my $count = CORE::length($1);       # 0b100 => len('00') => 2
1989     my $cnt = $count;                   # counter for loop
1990     unshift (@$x, 0);                   # add one element, together with one
1991                                         # more below in the loop this makes 2
1992     while ($cnt-- > 0)
1993       {
1994       # 'inflate' $X by adding one element, basically computing
1995       # $x * $BASE * $BASE. This gives us more $BASE_LEN digits for result
1996       # since len(sqrt($X)) approx == len($x) / 2.
1997       unshift (@$x, 0);
1998       # calculate sqrt($x), $x is now one element to big, again. In the next
1999       # round we make that two, again.
2000       _sqrt($c,$x);
2001       }
2002     # $x is now one element to big, so truncate result by removing it
2003     splice (@$x,0,1);
2004     } 
2005   else
2006     {
2007     # trial computation by starting with 2,4,8,16 etc until we overstep
2008     my $step;
2009     my $trial = _two();
2010
2011     # while still to do more than X steps
2012     do
2013       {
2014       $step = _two();
2015       while (_acmp($c, _pow($c, _copy($c, $trial), $n), $x) < 0)
2016         {
2017         _mul ($c, $step, [2]);
2018         _add ($c, $trial, $step);
2019         }
2020
2021       # hit exactly?
2022       if (_acmp($c, _pow($c, _copy($c, $trial), $n), $x) == 0)
2023         {
2024         @$x = @$trial;                  # make copy while preserving ref to $x
2025         return $x;
2026         }
2027       # overstepped, so go back on step
2028       _sub($c, $trial, $step);
2029       } while (scalar @$step > 1 || $step->[0] > 128);
2030
2031     # reset step to 2
2032     $step = _two();
2033     # add two, because $trial cannot be exactly the result (otherwise we would
2034     # already have found it)
2035     _add($c, $trial, $step);
2036  
2037     # and now add more and more (2,4,6,8,10 etc)
2038     while (_acmp($c, _pow($c, _copy($c, $trial), $n), $x) < 0)
2039       {
2040       _add ($c, $trial, $step);
2041       }
2042
2043     # hit not exactly? (overstepped)
2044     if (_acmp($c, _pow($c, _copy($c, $trial), $n), $x) > 0)
2045       {
2046       _dec($c,$trial);
2047       }
2048
2049     # hit not exactly? (overstepped)
2050     # 80 too small, 81 slightly too big, 82 too big
2051     if (_acmp($c, _pow($c, _copy($c, $trial), $n), $x) > 0)
2052       {
2053       _dec ($c, $trial); 
2054       }
2055
2056     @$x = @$trial;                      # make copy while preserving ref to $x
2057     return $x;
2058     }
2059   $x; 
2060   }
2061
2062 ##############################################################################
2063 # binary stuff
2064
2065 sub _and
2066   {
2067   my ($c,$x,$y) = @_;
2068
2069   # the shortcut makes equal, large numbers _really_ fast, and makes only a
2070   # very small performance drop for small numbers (e.g. something with less
2071   # than 32 bit) Since we optimize for large numbers, this is enabled.
2072   return $x if _acmp($c,$x,$y) == 0;            # shortcut
2073   
2074   my $m = _one(); my ($xr,$yr);
2075   my $mask = $AND_MASK;
2076
2077   my $x1 = $x;
2078   my $y1 = _copy($c,$y);                        # make copy
2079   $x = _zero();
2080   my ($b,$xrr,$yrr);
2081   use integer;
2082   while (!_is_zero($c,$x1) && !_is_zero($c,$y1))
2083     {
2084     ($x1, $xr) = _div($c,$x1,$mask);
2085     ($y1, $yr) = _div($c,$y1,$mask);
2086
2087     # make ints() from $xr, $yr
2088     # this is when the AND_BITS are greater than $BASE and is slower for
2089     # small (<256 bits) numbers, but faster for large numbers. Disabled
2090     # due to KISS principle
2091
2092 #    $b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
2093 #    $b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
2094 #    _add($c,$x, _mul($c, _new( $c, ($xrr & $yrr) ), $m) );
2095     
2096     # 0+ due to '&' doesn't work in strings
2097     _add($c,$x, _mul($c, [ 0+$xr->[0] & 0+$yr->[0] ], $m) );
2098     _mul($c,$m,$mask);
2099     }
2100   $x;
2101   }
2102
2103 sub _xor
2104   {
2105   my ($c,$x,$y) = @_;
2106
2107   return _zero() if _acmp($c,$x,$y) == 0;       # shortcut (see -and)
2108
2109   my $m = _one(); my ($xr,$yr);
2110   my $mask = $XOR_MASK;
2111
2112   my $x1 = $x;
2113   my $y1 = _copy($c,$y);                        # make copy
2114   $x = _zero();
2115   my ($b,$xrr,$yrr);
2116   use integer;
2117   while (!_is_zero($c,$x1) && !_is_zero($c,$y1))
2118     {
2119     ($x1, $xr) = _div($c,$x1,$mask);
2120     ($y1, $yr) = _div($c,$y1,$mask);
2121     # make ints() from $xr, $yr (see _and())
2122     #$b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
2123     #$b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
2124     #_add($c,$x, _mul($c, _new( $c, ($xrr ^ $yrr) ), $m) );
2125
2126     # 0+ due to '^' doesn't work in strings
2127     _add($c,$x, _mul($c, [ 0+$xr->[0] ^ 0+$yr->[0] ], $m) );
2128     _mul($c,$m,$mask);
2129     }
2130   # the loop stops when the shorter of the two numbers is exhausted
2131   # the remainder of the longer one will survive bit-by-bit, so we simple
2132   # multiply-add it in
2133   _add($c,$x, _mul($c, $x1, $m) ) if !_is_zero($c,$x1);
2134   _add($c,$x, _mul($c, $y1, $m) ) if !_is_zero($c,$y1);
2135   
2136   $x;
2137   }
2138
2139 sub _or
2140   {
2141   my ($c,$x,$y) = @_;
2142
2143   return $x if _acmp($c,$x,$y) == 0;            # shortcut (see _and)
2144
2145   my $m = _one(); my ($xr,$yr);
2146   my $mask = $OR_MASK;
2147
2148   my $x1 = $x;
2149   my $y1 = _copy($c,$y);                        # make copy
2150   $x = _zero();
2151   my ($b,$xrr,$yrr);
2152   use integer;
2153   while (!_is_zero($c,$x1) && !_is_zero($c,$y1))
2154     {
2155     ($x1, $xr) = _div($c,$x1,$mask);
2156     ($y1, $yr) = _div($c,$y1,$mask);
2157     # make ints() from $xr, $yr (see _and())
2158 #    $b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
2159 #    $b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
2160 #    _add($c,$x, _mul($c, _new( $c, ($xrr | $yrr) ), $m) );
2161     
2162     # 0+ due to '|' doesn't work in strings
2163     _add($c,$x, _mul($c, [ 0+$xr->[0] | 0+$yr->[0] ], $m) );
2164     _mul($c,$m,$mask);
2165     }
2166   # the loop stops when the shorter of the two numbers is exhausted
2167   # the remainder of the longer one will survive bit-by-bit, so we simple
2168   # multiply-add it in
2169   _add($c,$x, _mul($c, $x1, $m) ) if !_is_zero($c,$x1);
2170   _add($c,$x, _mul($c, $y1, $m) ) if !_is_zero($c,$y1);
2171   
2172   $x;
2173   }
2174
2175 sub _as_hex
2176   {
2177   # convert a decimal number to hex (ref to array, return ref to string)
2178   my ($c,$x) = @_;
2179
2180   # fits into one element (handle also 0x0 case)
2181   return sprintf("0x%x",$x->[0]) if @$x == 1;
2182
2183   my $x1 = _copy($c,$x);
2184
2185   my $es = '';
2186   my ($xr, $h, $x10000);
2187   if ($] >= 5.006)
2188     {
2189     $x10000 = [ 0x10000 ]; $h = 'h4';
2190     }
2191   else
2192     {
2193     $x10000 = [ 0x1000 ]; $h = 'h3';
2194     }
2195   while (@$x1 != 1 || $x1->[0] != 0)            # _is_zero()
2196     {
2197     ($x1, $xr) = _div($c,$x1,$x10000);
2198     $es .= unpack($h,pack('V',$xr->[0]));
2199     }
2200   $es = reverse $es;
2201   $es =~ s/^[0]+//;   # strip leading zeros
2202   '0x' . $es;                                   # return result prepended with 0x
2203   }
2204
2205 sub _as_bin
2206   {
2207   # convert a decimal number to bin (ref to array, return ref to string)
2208   my ($c,$x) = @_;
2209
2210   # fits into one element (and Perl recent enough), handle also 0b0 case
2211   # handle zero case for older Perls
2212   if ($] <= 5.005 && @$x == 1 && $x->[0] == 0)
2213     {
2214     my $t = '0b0'; return $t;
2215     }
2216   if (@$x == 1 && $] >= 5.006)
2217     {
2218     my $t = sprintf("0b%b",$x->[0]);
2219     return $t;
2220     }
2221   my $x1 = _copy($c,$x);
2222
2223   my $es = '';
2224   my ($xr, $b, $x10000);
2225   if ($] >= 5.006)
2226     {
2227     $x10000 = [ 0x10000 ]; $b = 'b16';
2228     }
2229   else
2230     {
2231     $x10000 = [ 0x1000 ]; $b = 'b12';
2232     }
2233   while (!(@$x1 == 1 && $x1->[0] == 0))         # _is_zero()
2234     {
2235     ($x1, $xr) = _div($c,$x1,$x10000);
2236     $es .= unpack($b,pack('v',$xr->[0]));
2237     }
2238   $es = reverse $es;
2239   $es =~ s/^[0]+//;   # strip leading zeros
2240   '0b' . $es;                                   # return result prepended with 0b
2241   }
2242
2243 sub _as_oct
2244   {
2245   # convert a decimal number to octal (ref to array, return ref to string)
2246   my ($c,$x) = @_;
2247
2248   # fits into one element (handle also 0 case)
2249   return sprintf("0%o",$x->[0]) if @$x == 1;
2250
2251   my $x1 = _copy($c,$x);
2252
2253   my $es = '';
2254   my $xr;
2255   my $x1000 = [ 0100000 ];
2256   while (@$x1 != 1 || $x1->[0] != 0)            # _is_zero()
2257     {
2258     ($x1, $xr) = _div($c,$x1,$x1000);
2259     $es .= reverse sprintf("%05o", $xr->[0]);
2260     }
2261   $es = reverse $es;
2262   $es =~ s/^[0]+//;   # strip leading zeros
2263   '0' . $es;                                    # return result prepended with 0
2264   }
2265
2266 sub _from_oct
2267   {
2268   # convert a octal number to decimal (string, return ref to array)
2269   my ($c,$os) = @_;
2270
2271   # for older Perls, play safe
2272   my $m = [ 0100000 ];
2273   my $d = 5;                                    # 5 digits at a time
2274
2275   my $mul = _one();
2276   my $x = _zero();
2277
2278   my $len = int( (length($os)-1)/$d );          # $d digit parts, w/o the '0'
2279   my $val; my $i = -$d;
2280   while ($len >= 0)
2281     {
2282     $val = substr($os,$i,$d);                   # get oct digits
2283     $val = CORE::oct($val);
2284     $i -= $d; $len --;
2285     my $adder = [ $val ];
2286     _add ($c, $x, _mul ($c, $adder, $mul ) ) if $val != 0;
2287     _mul ($c, $mul, $m ) if $len >= 0;          # skip last mul
2288     }
2289   $x;
2290   }
2291
2292 sub _from_hex
2293   {
2294   # convert a hex number to decimal (string, return ref to array)
2295   my ($c,$hs) = @_;
2296
2297   my $m = _new($c, 0x10000000);                 # 28 bit at a time (<32 bit!)
2298   my $d = 7;                                    # 7 digits at a time
2299   if ($] <= 5.006)
2300     {
2301     # for older Perls, play safe
2302     $m = [ 0x10000 ];                           # 16 bit at a time (<32 bit!)
2303     $d = 4;                                     # 4 digits at a time
2304     }
2305
2306   my $mul = _one();
2307   my $x = _zero();
2308
2309   my $len = int( (length($hs)-2)/$d );          # $d digit parts, w/o the '0x'
2310   my $val; my $i = -$d;
2311   while ($len >= 0)
2312     {
2313     $val = substr($hs,$i,$d);                   # get hex digits
2314     $val =~ s/^0x// if $len == 0;               # for last part only because
2315     $val = CORE::hex($val);                     # hex does not like wrong chars
2316     $i -= $d; $len --;
2317     my $adder = [ $val ];
2318     # if the resulting number was to big to fit into one element, create a
2319     # two-element version (bug found by Mark Lakata - Thanx!)
2320     if (CORE::length($val) > $BASE_LEN)
2321       {
2322       $adder = _new($c,$val);
2323       }
2324     _add ($c, $x, _mul ($c, $adder, $mul ) ) if $val != 0;
2325     _mul ($c, $mul, $m ) if $len >= 0;          # skip last mul
2326     }
2327   $x;
2328   }
2329
2330 sub _from_bin
2331   {
2332   # convert a hex number to decimal (string, return ref to array)
2333   my ($c,$bs) = @_;
2334
2335   # instead of converting X (8) bit at a time, it is faster to "convert" the
2336   # number to hex, and then call _from_hex.
2337
2338   my $hs = $bs;
2339   $hs =~ s/^[+-]?0b//;                                  # remove sign and 0b
2340   my $l = length($hs);                                  # bits
2341   $hs = '0' x (8-($l % 8)) . $hs if ($l % 8) != 0;      # padd left side w/ 0
2342   my $h = '0x' . unpack('H*', pack ('B*', $hs));        # repack as hex
2343   
2344   $c->_from_hex($h);
2345   }
2346
2347 ##############################################################################
2348 # special modulus functions
2349
2350 sub _modinv
2351   {
2352   # modular inverse
2353   my ($c,$x,$y) = @_;
2354
2355   my $u = _zero($c); my $u1 = _one($c);
2356   my $a = _copy($c,$y); my $b = _copy($c,$x);
2357
2358   # Euclid's Algorithm for bgcd(), only that we calc bgcd() ($a) and the
2359   # result ($u) at the same time. See comments in BigInt for why this works.
2360   my $q;
2361   ($a, $q, $b) = ($b, _div($c,$a,$b));          # step 1
2362   my $sign = 1;
2363   while (!_is_zero($c,$b))
2364     {
2365     my $t = _add($c,                            # step 2:
2366        _mul($c,_copy($c,$u1), $q) ,             #  t =  u1 * q
2367        $u );                                    #     + u
2368     $u = $u1;                                   #  u = u1, u1 = t
2369     $u1 = $t;
2370     $sign = -$sign;
2371     ($a, $q, $b) = ($b, _div($c,$a,$b));        # step 1
2372     }
2373
2374   # if the gcd is not 1, then return NaN
2375   return (undef,undef) unless _is_one($c,$a);
2376  
2377   ($u1, $sign == 1 ? '+' : '-');
2378   }
2379
2380 sub _modpow
2381   {
2382   # modulus of power ($x ** $y) % $z
2383   my ($c,$num,$exp,$mod) = @_;
2384
2385   # a^b (mod 1) = 0 for all a and b
2386   if (_is_one($c,$mod))
2387     {
2388         @$num = 0;
2389         return $num;
2390     }
2391
2392   # 0^a (mod m) = 0 if m != 0, a != 0
2393   # 0^0 (mod m) = 1 if m != 0
2394   if (_is_zero($c, $num)) {
2395       if (_is_zero($c, $exp)) {
2396           @$num = 1;
2397       } else {
2398           @$num = 0;
2399       }
2400       return $num;
2401   }
2402
2403 #  $num = _mod($c,$num,$mod);   # this does not make it faster
2404
2405   my $acc = _copy($c,$num); my $t = _one();
2406
2407   my $expbin = _as_bin($c,$exp); $expbin =~ s/^0b//;
2408   my $len = length($expbin);
2409   while (--$len >= 0)
2410     {
2411     if ( substr($expbin,$len,1) eq '1')                 # is_odd
2412       {
2413       _mul($c,$t,$acc);
2414       $t = _mod($c,$t,$mod);
2415       }
2416     _mul($c,$acc,$acc);
2417     $acc = _mod($c,$acc,$mod);
2418     }
2419   @$num = @$t;
2420   $num;
2421   }
2422
2423 sub _gcd
2424   {
2425   # greatest common divisor
2426   my ($c,$x,$y) = @_;
2427
2428   while ( (scalar @$y != 1) || ($y->[0] != 0) )         # while ($y != 0)
2429     {
2430     my $t = _copy($c,$y);
2431     $y = _mod($c, $x, $y);
2432     $x = $t;
2433     }
2434   $x;
2435   }
2436
2437 ##############################################################################
2438 ##############################################################################
2439
2440 1;
2441 __END__
2442
2443 =head1 NAME
2444
2445 Math::BigInt::Calc - Pure Perl module to support Math::BigInt
2446
2447 =head1 SYNOPSIS
2448
2449 Provides support for big integer calculations. Not intended to be used by other
2450 modules. Other modules which sport the same functions can also be used to support
2451 Math::BigInt, like Math::BigInt::GMP or Math::BigInt::Pari.
2452
2453 =head1 DESCRIPTION
2454
2455 In order to allow for multiple big integer libraries, Math::BigInt was
2456 rewritten to use library modules for core math routines. Any module which
2457 follows the same API as this can be used instead by using the following:
2458
2459         use Math::BigInt lib => 'libname';
2460
2461 'libname' is either the long name ('Math::BigInt::Pari'), or only the short
2462 version like 'Pari'.
2463
2464 =head1 STORAGE
2465
2466 =head1 METHODS
2467
2468 The following functions MUST be defined in order to support the use by
2469 Math::BigInt v1.70 or later:
2470
2471         api_version()   return API version, 1 for v1.70, 2 for v1.83
2472         _new(string)    return ref to new object from ref to decimal string
2473         _zero()         return a new object with value 0
2474         _one()          return a new object with value 1
2475         _two()          return a new object with value 2
2476         _ten()          return a new object with value 10
2477
2478         _str(obj)       return ref to a string representing the object
2479         _num(obj)       returns a Perl integer/floating point number
2480                         NOTE: because of Perl numeric notation defaults,
2481                         the _num'ified obj may lose accuracy due to 
2482                         machine-dependent floating point size limitations
2483                     
2484         _add(obj,obj)   Simple addition of two objects
2485         _mul(obj,obj)   Multiplication of two objects
2486         _div(obj,obj)   Division of the 1st object by the 2nd
2487                         In list context, returns (result,remainder).
2488                         NOTE: this is integer math, so no
2489                         fractional part will be returned.
2490                         The second operand will be not be 0, so no need to
2491                         check for that.
2492         _sub(obj,obj)   Simple subtraction of 1 object from another
2493                         a third, optional parameter indicates that the params
2494                         are swapped. In this case, the first param needs to
2495                         be preserved, while you can destroy the second.
2496                         sub (x,y,1) => return x - y and keep x intact!
2497         _dec(obj)       decrement object by one (input is guaranteed to be > 0)
2498         _inc(obj)       increment object by one
2499
2500
2501         _acmp(obj,obj)  <=> operator for objects (return -1, 0 or 1)
2502
2503         _len(obj)       returns count of the decimal digits of the object
2504         _digit(obj,n)   returns the n'th decimal digit of object
2505
2506         _is_one(obj)    return true if argument is 1
2507         _is_two(obj)    return true if argument is 2
2508         _is_ten(obj)    return true if argument is 10
2509         _is_zero(obj)   return true if argument is 0
2510         _is_even(obj)   return true if argument is even (0,2,4,6..)
2511         _is_odd(obj)    return true if argument is odd (1,3,5,7..)
2512
2513         _copy           return a ref to a true copy of the object
2514
2515         _check(obj)     check whether internal representation is still intact
2516                         return 0 for ok, otherwise error message as string
2517
2518         _from_hex(str)  return new object from a hexadecimal string
2519         _from_bin(str)  return new object from a binary string
2520         _from_oct(str)  return new object from an octal string
2521         
2522         _as_hex(str)    return string containing the value as
2523                         unsigned hex string, with the '0x' prepended.
2524                         Leading zeros must be stripped.
2525         _as_bin(str)    Like as_hex, only as binary string containing only
2526                         zeros and ones. Leading zeros must be stripped and a
2527                         '0b' must be prepended.
2528         
2529         _rsft(obj,N,B)  shift object in base B by N 'digits' right
2530         _lsft(obj,N,B)  shift object in base B by N 'digits' left
2531         
2532         _xor(obj1,obj2) XOR (bit-wise) object 1 with object 2
2533                         Note: XOR, AND and OR pad with zeros if size mismatches
2534         _and(obj1,obj2) AND (bit-wise) object 1 with object 2
2535         _or(obj1,obj2)  OR (bit-wise) object 1 with object 2
2536
2537         _mod(obj1,obj2) Return remainder of div of the 1st by the 2nd object
2538         _sqrt(obj)      return the square root of object (truncated to int)
2539         _root(obj)      return the n'th (n >= 3) root of obj (truncated to int)
2540         _fac(obj)       return factorial of object 1 (1*2*3*4..)
2541         _pow(obj1,obj2) return object 1 to the power of object 2
2542                         return undef for NaN
2543         _zeros(obj)     return number of trailing decimal zeros
2544         _modinv         return inverse modulus
2545         _modpow         return modulus of power ($x ** $y) % $z
2546         _log_int(X,N)   calculate integer log() of X in base N
2547                         X >= 0, N >= 0 (return undef for NaN)
2548                         returns (RESULT, EXACT) where EXACT is:
2549                          1     : result is exactly RESULT
2550                          0     : result was truncated to RESULT
2551                          undef : unknown whether result is exactly RESULT
2552         _gcd(obj,obj)   return Greatest Common Divisor of two objects
2553
2554 The following functions are REQUIRED for an api_version of 2 or greater:
2555
2556         _1ex($x)        create the number 1Ex where x >= 0
2557         _alen(obj)      returns approximate count of the decimal digits of the
2558                         object. This estimate MUST always be greater or equal
2559                         to what _len() returns.
2560         _nok(n,k)       calculate n over k (binomial coefficient)
2561
2562 The following functions are optional, and can be defined if the underlying lib
2563 has a fast way to do them. If undefined, Math::BigInt will use pure Perl (hence
2564 slow) fallback routines to emulate these:
2565         
2566         _signed_or
2567         _signed_and
2568         _signed_xor
2569
2570 Input strings come in as unsigned but with prefix (i.e. as '123', '0xabc'
2571 or '0b1101').
2572
2573 So the library needs only to deal with unsigned big integers. Testing of input
2574 parameter validity is done by the caller, so you need not worry about
2575 underflow (f.i. in C<_sub()>, C<_dec()>) nor about division by zero or similar
2576 cases.
2577
2578 The first parameter can be modified, that includes the possibility that you
2579 return a reference to a completely different object instead. Although keeping
2580 the reference and just changing its contents is preferred over creating and
2581 returning a different reference.
2582
2583 Return values are always references to objects, strings, or true/false for
2584 comparison routines.
2585
2586 =head1 WRAP YOUR OWN
2587
2588 If you want to port your own favourite c-lib for big numbers to the
2589 Math::BigInt interface, you can take any of the already existing modules as
2590 a rough guideline. You should really wrap up the latest BigInt and BigFloat
2591 testsuites with your module, and replace in them any of the following:
2592
2593         use Math::BigInt;
2594
2595 by this:
2596
2597         use Math::BigInt lib => 'yourlib';
2598
2599 This way you ensure that your library really works 100% within Math::BigInt.
2600
2601 =head1 LICENSE
2602  
2603 This program is free software; you may redistribute it and/or modify it under
2604 the same terms as Perl itself. 
2605
2606 =head1 AUTHORS
2607
2608 Original math code by Mark Biggar, rewritten by Tels L<http://bloodgate.com/>
2609 in late 2000.
2610 Separated from BigInt and shaped API with the help of John Peacock.
2611
2612 Fixed, speed-up, streamlined and enhanced by Tels 2001 - 2007.
2613
2614 =head1 SEE ALSO
2615
2616 L<Math::BigInt>, L<Math::BigFloat>,
2617 L<Math::BigInt::GMP>, L<Math::BigInt::FastCalc> and L<Math::BigInt::Pari>.
2618
2619 =cut