This is a live mirror of the Perl 5 development currently hosted at https://github.com/perl/perl5
e7b2102524d845a79d8bbaf4378a98cdcfb85ec0
[perl5.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 = '0.57';
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 number (scalar int/float) from a BigInt object 
276   my $x = $_[1];
277
278   return 0+$x->[0] if scalar @$x == 1;  # below $BASE
279   my $fac = 1;
280   my $num = 0;
281   foreach (@$x)
282     {
283     $num += $fac*$_; $fac *= $BASE;
284     }
285   $num; 
286   }
287
288 ##############################################################################
289 # actual math code
290
291 sub _add
292   {
293   # (ref to int_num_array, ref to int_num_array)
294   # routine to add two base 1eX numbers
295   # stolen from Knuth Vol 2 Algorithm A pg 231
296   # there are separate routines to add and sub as per Knuth pg 233
297   # This routine clobbers up array x, but not y.
298  
299   my ($c,$x,$y) = @_;
300
301   return $x if (@$y == 1) && $y->[0] == 0;              # $x + 0 => $x
302   if ((@$x == 1) && $x->[0] == 0)                       # 0 + $y => $y->copy
303     {
304     # twice as slow as $x = [ @$y ], but nec. to retain $x as ref :(
305     @$x = @$y; return $x;               
306     }
307  
308   # for each in Y, add Y to X and carry. If after that, something is left in
309   # X, foreach in X add carry to X and then return X, carry
310   # Trades one "$j++" for having to shift arrays
311   my $i; my $car = 0; my $j = 0;
312   for $i (@$y)
313     {
314     $x->[$j] -= $BASE if $car = (($x->[$j] += $i + $car) >= $BASE) ? 1 : 0;
315     $j++;
316     }
317   while ($car != 0)
318     {
319     $x->[$j] -= $BASE if $car = (($x->[$j] += $car) >= $BASE) ? 1 : 0; $j++;
320     }
321   $x;
322   }                                                                             
323
324 sub _inc
325   {
326   # (ref to int_num_array, ref to int_num_array)
327   # Add 1 to $x, modify $x in place
328   my ($c,$x) = @_;
329
330   for my $i (@$x)
331     {
332     return $x if (($i += 1) < $BASE);           # early out
333     $i = 0;                                     # overflow, next
334     }
335   push @$x,1 if (($x->[-1] || 0) == 0);         # last overflowed, so extend
336   $x;
337   }                                                                             
338
339 sub _dec
340   {
341   # (ref to int_num_array, ref to int_num_array)
342   # Sub 1 from $x, modify $x in place
343   my ($c,$x) = @_;
344
345   my $MAX = $BASE-1;                            # since MAX_VAL based on BASE
346   for my $i (@$x)
347     {
348     last if (($i -= 1) >= 0);                   # early out
349     $i = $MAX;                                  # underflow, next
350     }
351   pop @$x if $x->[-1] == 0 && @$x > 1;          # last underflowed (but leave 0)
352   $x;
353   }                                                                             
354
355 sub _sub
356   {
357   # (ref to int_num_array, ref to int_num_array, swap)
358   # subtract base 1eX numbers -- stolen from Knuth Vol 2 pg 232, $x > $y
359   # subtract Y from X by modifying x in place
360   my ($c,$sx,$sy,$s) = @_;
361  
362   my $car = 0; my $i; my $j = 0;
363   if (!$s)
364     {
365     for $i (@$sx)
366       {
367       last unless defined $sy->[$j] || $car;
368       $i += $BASE if $car = (($i -= ($sy->[$j] || 0) + $car) < 0); $j++;
369       }
370     # might leave leading zeros, so fix that
371     return __strip_zeros($sx);
372     }
373   for $i (@$sx)
374     {
375     # we can't do an early out if $x is < than $y, since we
376     # need to copy the high chunks from $y. Found by Bob Mathews.
377     #last unless defined $sy->[$j] || $car;
378     $sy->[$j] += $BASE
379      if $car = (($sy->[$j] = $i-($sy->[$j]||0) - $car) < 0);
380     $j++;
381     }
382   # might leave leading zeros, so fix that
383   __strip_zeros($sy);
384   }                                                                             
385
386 sub _mul_use_mul
387   {
388   # (ref to int_num_array, ref to int_num_array)
389   # multiply two numbers in internal representation
390   # modifies first arg, second need not be different from first
391   my ($c,$xv,$yv) = @_;
392
393   if (@$yv == 1)
394     {
395     # shortcut for two very short numbers (improved by Nathan Zook)
396     # works also if xv and yv are the same reference, and handles also $x == 0
397     if (@$xv == 1)
398       {
399       if (($xv->[0] *= $yv->[0]) >= $BASE)
400          {
401          $xv->[0] = $xv->[0] - ($xv->[1] = int($xv->[0] * $RBASE)) * $BASE;
402          };
403       return $xv;
404       }
405     # $x * 0 => 0
406     if ($yv->[0] == 0)
407       {
408       @$xv = (0);
409       return $xv;
410       }
411     # multiply a large number a by a single element one, so speed up
412     my $y = $yv->[0]; my $car = 0;
413     foreach my $i (@$xv)
414       {
415       $i = $i * $y + $car; $car = int($i * $RBASE); $i -= $car * $BASE;
416       }
417     push @$xv, $car if $car != 0;
418     return $xv;
419     }
420   # shortcut for result $x == 0 => result = 0
421   return $xv if ( ((@$xv == 1) && ($xv->[0] == 0)) ); 
422
423   # since multiplying $x with $x fails, make copy in this case
424   $yv = [@$xv] if $xv == $yv;   # same references?
425
426   my @prod = (); my ($prod,$car,$cty,$xi,$yi);
427
428   for $xi (@$xv)
429     {
430     $car = 0; $cty = 0;
431
432     # slow variant
433 #    for $yi (@$yv)
434 #      {
435 #      $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
436 #      $prod[$cty++] =
437 #       $prod - ($car = int($prod * RBASE)) * $BASE;  # see USE_MUL
438 #      }
439 #    $prod[$cty] += $car if $car; # need really to check for 0?
440 #    $xi = shift @prod;
441
442     # faster variant
443     # looping through this if $xi == 0 is silly - so optimize it away!
444     $xi = (shift @prod || 0), next if $xi == 0;
445     for $yi (@$yv)
446       {
447       $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
448 ##     this is actually a tad slower
449 ##        $prod = $prod[$cty]; $prod += ($car + $xi * $yi);     # no ||0 here
450       $prod[$cty++] =
451        $prod - ($car = int($prod * $RBASE)) * $BASE;  # see USE_MUL
452       }
453     $prod[$cty] += $car if $car; # need really to check for 0?
454     $xi = shift @prod || 0;     # || 0 makes v5.005_3 happy
455     }
456   push @$xv, @prod;
457   # can't have leading zeros
458 #  __strip_zeros($xv);
459   $xv;
460   }                                                                             
461
462 sub _mul_use_div_64
463   {
464   # (ref to int_num_array, ref to int_num_array)
465   # multiply two numbers in internal representation
466   # modifies first arg, second need not be different from first
467   # works for 64 bit integer with "use integer"
468   my ($c,$xv,$yv) = @_;
469
470   use integer;
471   if (@$yv == 1)
472     {
473     # shortcut for two small numbers, also handles $x == 0
474     if (@$xv == 1)
475       {
476       # shortcut for two very short numbers (improved by Nathan Zook)
477       # works also if xv and yv are the same reference, and handles also $x == 0
478       if (($xv->[0] *= $yv->[0]) >= $BASE)
479           {
480           $xv->[0] =
481               $xv->[0] - ($xv->[1] = $xv->[0] / $BASE) * $BASE;
482           };
483       return $xv;
484       }
485     # $x * 0 => 0
486     if ($yv->[0] == 0)
487       {
488       @$xv = (0);
489       return $xv;
490       }
491     # multiply a large number a by a single element one, so speed up
492     my $y = $yv->[0]; my $car = 0;
493     foreach my $i (@$xv)
494       {
495       #$i = $i * $y + $car; $car = $i / $BASE; $i -= $car * $BASE;
496       $i = $i * $y + $car; $i -= ($car = $i / $BASE) * $BASE;
497       }
498     push @$xv, $car if $car != 0;
499     return $xv;
500     }
501   # shortcut for result $x == 0 => result = 0
502   return $xv if ( ((@$xv == 1) && ($xv->[0] == 0)) ); 
503
504   # since multiplying $x with $x fails, make copy in this case
505   $yv = [@$xv] if $xv == $yv;   # same references?
506
507   my @prod = (); my ($prod,$car,$cty,$xi,$yi);
508   for $xi (@$xv)
509     {
510     $car = 0; $cty = 0;
511     # looping through this if $xi == 0 is silly - so optimize it away!
512     $xi = (shift @prod || 0), next if $xi == 0;
513     for $yi (@$yv)
514       {
515       $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
516       $prod[$cty++] = $prod - ($car = $prod / $BASE) * $BASE;
517       }
518     $prod[$cty] += $car if $car; # need really to check for 0?
519     $xi = shift @prod || 0;     # || 0 makes v5.005_3 happy
520     }
521   push @$xv, @prod;
522   $xv;
523   }                                                                             
524
525 sub _mul_use_div
526   {
527   # (ref to int_num_array, ref to int_num_array)
528   # multiply two numbers in internal representation
529   # modifies first arg, second need not be different from first
530   my ($c,$xv,$yv) = @_;
531
532   if (@$yv == 1)
533     {
534     # shortcut for two small numbers, also handles $x == 0
535     if (@$xv == 1)
536       {
537       # shortcut for two very short numbers (improved by Nathan Zook)
538       # works also if xv and yv are the same reference, and handles also $x == 0
539       if (($xv->[0] *= $yv->[0]) >= $BASE)
540           {
541           $xv->[0] =
542               $xv->[0] - ($xv->[1] = int($xv->[0] / $BASE)) * $BASE;
543           };
544       return $xv;
545       }
546     # $x * 0 => 0
547     if ($yv->[0] == 0)
548       {
549       @$xv = (0);
550       return $xv;
551       }
552     # multiply a large number a by a single element one, so speed up
553     my $y = $yv->[0]; my $car = 0;
554     foreach my $i (@$xv)
555       {
556       $i = $i * $y + $car; $car = int($i / $BASE); $i -= $car * $BASE;
557       # This (together with use integer;) does not work on 32-bit Perls
558       #$i = $i * $y + $car; $i -= ($car = $i / $BASE) * $BASE;
559       }
560     push @$xv, $car if $car != 0;
561     return $xv;
562     }
563   # shortcut for result $x == 0 => result = 0
564   return $xv if ( ((@$xv == 1) && ($xv->[0] == 0)) ); 
565
566   # since multiplying $x with $x fails, make copy in this case
567   $yv = [@$xv] if $xv == $yv;   # same references?
568
569   my @prod = (); my ($prod,$car,$cty,$xi,$yi);
570   for $xi (@$xv)
571     {
572     $car = 0; $cty = 0;
573     # looping through this if $xi == 0 is silly - so optimize it away!
574     $xi = (shift @prod || 0), next if $xi == 0;
575     for $yi (@$yv)
576       {
577       $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
578       $prod[$cty++] = $prod - ($car = int($prod / $BASE)) * $BASE;
579       }
580     $prod[$cty] += $car if $car; # need really to check for 0?
581     $xi = shift @prod || 0;     # || 0 makes v5.005_3 happy
582     }
583   push @$xv, @prod;
584   # can't have leading zeros
585 #  __strip_zeros($xv);
586   $xv;
587   }                                                                             
588
589 sub _div_use_mul
590   {
591   # ref to array, ref to array, modify first array and return remainder if 
592   # in list context
593
594   # see comments in _div_use_div() for more explanations
595
596   my ($c,$x,$yorg) = @_;
597   
598   # the general div algorithm here is about O(N*N) and thus quite slow, so
599   # we first check for some special cases and use shortcuts to handle them.
600
601   # This works, because we store the numbers in a chunked format where each
602   # element contains 5..7 digits (depending on system).
603
604   # if both numbers have only one element:
605   if (@$x == 1 && @$yorg == 1)
606     {
607     # shortcut, $yorg and $x are two small numbers
608     if (wantarray)
609       {
610       my $r = [ $x->[0] % $yorg->[0] ];
611       $x->[0] = int($x->[0] / $yorg->[0]);
612       return ($x,$r); 
613       }
614     else
615       {
616       $x->[0] = int($x->[0] / $yorg->[0]);
617       return $x; 
618       }
619     }
620
621   # if x has more than one, but y has only one element:
622   if (@$yorg == 1)
623     {
624     my $rem;
625     $rem = _mod($c,[ @$x ],$yorg) if wantarray;
626
627     # shortcut, $y is < $BASE
628     my $j = scalar @$x; my $r = 0; 
629     my $y = $yorg->[0]; my $b;
630     while ($j-- > 0)
631       {
632       $b = $r * $BASE + $x->[$j];
633       $x->[$j] = int($b/$y);
634       $r = $b % $y;
635       }
636     pop @$x if @$x > 1 && $x->[-1] == 0;        # splice up a leading zero 
637     return ($x,$rem) if wantarray;
638     return $x;
639     }
640
641   # now x and y have more than one element
642
643   # check whether y has more elements than x, if yet, the result will be 0
644   if (@$yorg > @$x)
645     {
646     my $rem;
647     $rem = [@$x] if wantarray;                  # make copy
648     splice (@$x,1);                             # keep ref to original array
649     $x->[0] = 0;                                # set to 0
650     return ($x,$rem) if wantarray;              # including remainder?
651     return $x;                                  # only x, which is [0] now
652     }
653   # check whether the numbers have the same number of elements, in that case
654   # the result will fit into one element and can be computed efficiently
655   if (@$yorg == @$x)
656     {
657     my $rem;
658     # if $yorg has more digits than $x (it's leading element is longer than
659     # the one from $x), the result will also be 0:
660     if (length(int($yorg->[-1])) > length(int($x->[-1])))
661       {
662       $rem = [@$x] if wantarray;                # make copy
663       splice (@$x,1);                           # keep ref to org array
664       $x->[0] = 0;                              # set to 0
665       return ($x,$rem) if wantarray;            # including remainder?
666       return $x;
667       }
668     # now calculate $x / $yorg
669     if (length(int($yorg->[-1])) == length(int($x->[-1])))
670       {
671       # same length, so make full compare
672
673       my $a = 0; my $j = scalar @$x - 1;
674       # manual way (abort if unequal, good for early ne)
675       while ($j >= 0)
676         {
677         last if ($a = $x->[$j] - $yorg->[$j]); $j--;
678         }
679       # $a contains the result of the compare between X and Y
680       # a < 0: x < y, a == 0: x == y, a > 0: x > y
681       if ($a <= 0)
682         {
683         $rem = [ 0 ];                   # a = 0 => x == y => rem 0
684         $rem = [@$x] if $a != 0;        # a < 0 => x < y => rem = x
685         splice(@$x,1);                  # keep single element
686         $x->[0] = 0;                    # if $a < 0
687         $x->[0] = 1 if $a == 0;         # $x == $y
688         return ($x,$rem) if wantarray;
689         return $x;
690         }
691       # $x >= $y, so proceed normally
692       }
693     }
694
695   # all other cases:
696
697   my $y = [ @$yorg ];                           # always make copy to preserve
698
699   my ($car,$bar,$prd,$dd,$xi,$yi,@q,$v2,$v1,@d,$tmp,$q,$u2,$u1,$u0);
700
701   $car = $bar = $prd = 0;
702   if (($dd = int($BASE/($y->[-1]+1))) != 1) 
703     {
704     for $xi (@$x) 
705       {
706       $xi = $xi * $dd + $car;
707       $xi -= ($car = int($xi * $RBASE)) * $BASE;        # see USE_MUL
708       }
709     push(@$x, $car); $car = 0;
710     for $yi (@$y) 
711       {
712       $yi = $yi * $dd + $car;
713       $yi -= ($car = int($yi * $RBASE)) * $BASE;        # see USE_MUL
714       }
715     }
716   else 
717     {
718     push(@$x, 0);
719     }
720   @q = (); ($v2,$v1) = @$y[-2,-1];
721   $v2 = 0 unless $v2;
722   while ($#$x > $#$y) 
723     {
724     ($u2,$u1,$u0) = @$x[-3..-1];
725     $u2 = 0 unless $u2;
726     #warn "oups v1 is 0, u0: $u0 $y->[-2] $y->[-1] l ",scalar @$y,"\n"
727     # if $v1 == 0;
728     $q = (($u0 == $v1) ? $MAX_VAL : int(($u0*$BASE+$u1)/$v1));
729     --$q while ($v2*$q > ($u0*$BASE+$u1-$q*$v1)*$BASE+$u2);
730     if ($q)
731       {
732       ($car, $bar) = (0,0);
733       for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi) 
734         {
735         $prd = $q * $y->[$yi] + $car;
736         $prd -= ($car = int($prd * $RBASE)) * $BASE;    # see USE_MUL
737         $x->[$xi] += $BASE if ($bar = (($x->[$xi] -= $prd + $bar) < 0));
738         }
739       if ($x->[-1] < $car + $bar) 
740         {
741         $car = 0; --$q;
742         for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi) 
743           {
744           $x->[$xi] -= $BASE
745            if ($car = (($x->[$xi] += $y->[$yi] + $car) >= $BASE));
746           }
747         }   
748       }
749     pop(@$x);
750     unshift(@q, $q);
751     }
752   if (wantarray) 
753     {
754     @d = ();
755     if ($dd != 1)  
756       {
757       $car = 0; 
758       for $xi (reverse @$x) 
759         {
760         $prd = $car * $BASE + $xi;
761         $car = $prd - ($tmp = int($prd / $dd)) * $dd; # see USE_MUL
762         unshift(@d, $tmp);
763         }
764       }
765     else 
766       {
767       @d = @$x;
768       }
769     @$x = @q;
770     my $d = \@d; 
771     __strip_zeros($x);
772     __strip_zeros($d);
773     return ($x,$d);
774     }
775   @$x = @q;
776   __strip_zeros($x);
777   $x;
778   }
779
780 sub _div_use_div_64
781   {
782   # ref to array, ref to array, modify first array and return remainder if 
783   # in list context
784   # This version works on 64 bit integers
785   my ($c,$x,$yorg) = @_;
786
787   use integer;
788   # the general div algorithm here is about O(N*N) and thus quite slow, so
789   # we first check for some special cases and use shortcuts to handle them.
790
791   # This works, because we store the numbers in a chunked format where each
792   # element contains 5..7 digits (depending on system).
793
794   # if both numbers have only one element:
795   if (@$x == 1 && @$yorg == 1)
796     {
797     # shortcut, $yorg and $x are two small numbers
798     if (wantarray)
799       {
800       my $r = [ $x->[0] % $yorg->[0] ];
801       $x->[0] = int($x->[0] / $yorg->[0]);
802       return ($x,$r); 
803       }
804     else
805       {
806       $x->[0] = int($x->[0] / $yorg->[0]);
807       return $x; 
808       }
809     }
810   # if x has more than one, but y has only one element:
811   if (@$yorg == 1)
812     {
813     my $rem;
814     $rem = _mod($c,[ @$x ],$yorg) if wantarray;
815
816     # shortcut, $y is < $BASE
817     my $j = scalar @$x; my $r = 0; 
818     my $y = $yorg->[0]; my $b;
819     while ($j-- > 0)
820       {
821       $b = $r * $BASE + $x->[$j];
822       $x->[$j] = int($b/$y);
823       $r = $b % $y;
824       }
825     pop @$x if @$x > 1 && $x->[-1] == 0;        # splice up a leading zero 
826     return ($x,$rem) if wantarray;
827     return $x;
828     }
829   # now x and y have more than one element
830
831   # check whether y has more elements than x, if yet, the result will be 0
832   if (@$yorg > @$x)
833     {
834     my $rem;
835     $rem = [@$x] if wantarray;                  # make copy
836     splice (@$x,1);                             # keep ref to original array
837     $x->[0] = 0;                                # set to 0
838     return ($x,$rem) if wantarray;              # including remainder?
839     return $x;                                  # only x, which is [0] now
840     }
841   # check whether the numbers have the same number of elements, in that case
842   # the result will fit into one element and can be computed efficiently
843   if (@$yorg == @$x)
844     {
845     my $rem;
846     # if $yorg has more digits than $x (it's leading element is longer than
847     # the one from $x), the result will also be 0:
848     if (length(int($yorg->[-1])) > length(int($x->[-1])))
849       {
850       $rem = [@$x] if wantarray;                # make copy
851       splice (@$x,1);                           # keep ref to org array
852       $x->[0] = 0;                              # set to 0
853       return ($x,$rem) if wantarray;            # including remainder?
854       return $x;
855       }
856     # now calculate $x / $yorg
857
858     if (length(int($yorg->[-1])) == length(int($x->[-1])))
859       {
860       # same length, so make full compare
861
862       my $a = 0; my $j = scalar @$x - 1;
863       # manual way (abort if unequal, good for early ne)
864       while ($j >= 0)
865         {
866         last if ($a = $x->[$j] - $yorg->[$j]); $j--;
867         }
868       # $a contains the result of the compare between X and Y
869       # a < 0: x < y, a == 0: x == y, a > 0: x > y
870       if ($a <= 0)
871         {
872         $rem = [ 0 ];                   # a = 0 => x == y => rem 0
873         $rem = [@$x] if $a != 0;        # a < 0 => x < y => rem = x
874         splice(@$x,1);                  # keep single element
875         $x->[0] = 0;                    # if $a < 0
876         $x->[0] = 1 if $a == 0;         # $x == $y
877         return ($x,$rem) if wantarray;  # including remainder?
878         return $x;
879         }
880       # $x >= $y, so proceed normally
881
882       }
883     }
884
885   # all other cases:
886
887   my $y = [ @$yorg ];                           # always make copy to preserve
888  
889   my ($car,$bar,$prd,$dd,$xi,$yi,@q,$v2,$v1,@d,$tmp,$q,$u2,$u1,$u0);
890
891   $car = $bar = $prd = 0;
892   if (($dd = int($BASE/($y->[-1]+1))) != 1) 
893     {
894     for $xi (@$x) 
895       {
896       $xi = $xi * $dd + $car;
897       $xi -= ($car = int($xi / $BASE)) * $BASE;
898       }
899     push(@$x, $car); $car = 0;
900     for $yi (@$y) 
901       {
902       $yi = $yi * $dd + $car;
903       $yi -= ($car = int($yi / $BASE)) * $BASE;
904       }
905     }
906   else 
907     {
908     push(@$x, 0);
909     }
910
911   # @q will accumulate the final result, $q contains the current computed
912   # part of the final result
913
914   @q = (); ($v2,$v1) = @$y[-2,-1];
915   $v2 = 0 unless $v2;
916   while ($#$x > $#$y) 
917     {
918     ($u2,$u1,$u0) = @$x[-3..-1];
919     $u2 = 0 unless $u2;
920     #warn "oups v1 is 0, u0: $u0 $y->[-2] $y->[-1] l ",scalar @$y,"\n"
921     # if $v1 == 0;
922     $q = (($u0 == $v1) ? $MAX_VAL : int(($u0*$BASE+$u1)/$v1));
923     --$q while ($v2*$q > ($u0*$BASE+$u1-$q*$v1)*$BASE+$u2);
924     if ($q)
925       {
926       ($car, $bar) = (0,0);
927       for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi) 
928         {
929         $prd = $q * $y->[$yi] + $car;
930         $prd -= ($car = int($prd / $BASE)) * $BASE;
931         $x->[$xi] += $BASE if ($bar = (($x->[$xi] -= $prd + $bar) < 0));
932         }
933       if ($x->[-1] < $car + $bar) 
934         {
935         $car = 0; --$q;
936         for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi) 
937           {
938           $x->[$xi] -= $BASE
939            if ($car = (($x->[$xi] += $y->[$yi] + $car) >= $BASE));
940           }
941         }   
942       }
943     pop(@$x); unshift(@q, $q);
944     }
945   if (wantarray) 
946     {
947     @d = ();
948     if ($dd != 1)  
949       {
950       $car = 0; 
951       for $xi (reverse @$x) 
952         {
953         $prd = $car * $BASE + $xi;
954         $car = $prd - ($tmp = int($prd / $dd)) * $dd;
955         unshift(@d, $tmp);
956         }
957       }
958     else 
959       {
960       @d = @$x;
961       }
962     @$x = @q;
963     my $d = \@d; 
964     __strip_zeros($x);
965     __strip_zeros($d);
966     return ($x,$d);
967     }
968   @$x = @q;
969   __strip_zeros($x);
970   $x;
971   }
972
973 sub _div_use_div
974   {
975   # ref to array, ref to array, modify first array and return remainder if 
976   # in list context
977   my ($c,$x,$yorg) = @_;
978
979   # the general div algorithm here is about O(N*N) and thus quite slow, so
980   # we first check for some special cases and use shortcuts to handle them.
981
982   # This works, because we store the numbers in a chunked format where each
983   # element contains 5..7 digits (depending on system).
984
985   # if both numbers have only one element:
986   if (@$x == 1 && @$yorg == 1)
987     {
988     # shortcut, $yorg and $x are two small numbers
989     if (wantarray)
990       {
991       my $r = [ $x->[0] % $yorg->[0] ];
992       $x->[0] = int($x->[0] / $yorg->[0]);
993       return ($x,$r); 
994       }
995     else
996       {
997       $x->[0] = int($x->[0] / $yorg->[0]);
998       return $x; 
999       }
1000     }
1001   # if x has more than one, but y has only one element:
1002   if (@$yorg == 1)
1003     {
1004     my $rem;
1005     $rem = _mod($c,[ @$x ],$yorg) if wantarray;
1006
1007     # shortcut, $y is < $BASE
1008     my $j = scalar @$x; my $r = 0; 
1009     my $y = $yorg->[0]; my $b;
1010     while ($j-- > 0)
1011       {
1012       $b = $r * $BASE + $x->[$j];
1013       $x->[$j] = int($b/$y);
1014       $r = $b % $y;
1015       }
1016     pop @$x if @$x > 1 && $x->[-1] == 0;        # splice up a leading zero 
1017     return ($x,$rem) if wantarray;
1018     return $x;
1019     }
1020   # now x and y have more than one element
1021
1022   # check whether y has more elements than x, if yet, the result will be 0
1023   if (@$yorg > @$x)
1024     {
1025     my $rem;
1026     $rem = [@$x] if wantarray;                  # make copy
1027     splice (@$x,1);                             # keep ref to original array
1028     $x->[0] = 0;                                # set to 0
1029     return ($x,$rem) if wantarray;              # including remainder?
1030     return $x;                                  # only x, which is [0] now
1031     }
1032   # check whether the numbers have the same number of elements, in that case
1033   # the result will fit into one element and can be computed efficiently
1034   if (@$yorg == @$x)
1035     {
1036     my $rem;
1037     # if $yorg has more digits than $x (it's leading element is longer than
1038     # the one from $x), the result will also be 0:
1039     if (length(int($yorg->[-1])) > length(int($x->[-1])))
1040       {
1041       $rem = [@$x] if wantarray;                # make copy
1042       splice (@$x,1);                           # keep ref to org array
1043       $x->[0] = 0;                              # set to 0
1044       return ($x,$rem) if wantarray;            # including remainder?
1045       return $x;
1046       }
1047     # now calculate $x / $yorg
1048
1049     if (length(int($yorg->[-1])) == length(int($x->[-1])))
1050       {
1051       # same length, so make full compare
1052
1053       my $a = 0; my $j = scalar @$x - 1;
1054       # manual way (abort if unequal, good for early ne)
1055       while ($j >= 0)
1056         {
1057         last if ($a = $x->[$j] - $yorg->[$j]); $j--;
1058         }
1059       # $a contains the result of the compare between X and Y
1060       # a < 0: x < y, a == 0: x == y, a > 0: x > y
1061       if ($a <= 0)
1062         {
1063         $rem = [ 0 ];                   # a = 0 => x == y => rem 0
1064         $rem = [@$x] if $a != 0;        # a < 0 => x < y => rem = x
1065         splice(@$x,1);                  # keep single element
1066         $x->[0] = 0;                    # if $a < 0
1067         $x->[0] = 1 if $a == 0;         # $x == $y
1068         return ($x,$rem) if wantarray;  # including remainder?
1069         return $x;
1070         }
1071       # $x >= $y, so proceed normally
1072
1073       }
1074     }
1075
1076   # all other cases:
1077
1078   my $y = [ @$yorg ];                           # always make copy to preserve
1079  
1080   my ($car,$bar,$prd,$dd,$xi,$yi,@q,$v2,$v1,@d,$tmp,$q,$u2,$u1,$u0);
1081
1082   $car = $bar = $prd = 0;
1083   if (($dd = int($BASE/($y->[-1]+1))) != 1) 
1084     {
1085     for $xi (@$x) 
1086       {
1087       $xi = $xi * $dd + $car;
1088       $xi -= ($car = int($xi / $BASE)) * $BASE;
1089       }
1090     push(@$x, $car); $car = 0;
1091     for $yi (@$y) 
1092       {
1093       $yi = $yi * $dd + $car;
1094       $yi -= ($car = int($yi / $BASE)) * $BASE;
1095       }
1096     }
1097   else 
1098     {
1099     push(@$x, 0);
1100     }
1101
1102   # @q will accumulate the final result, $q contains the current computed
1103   # part of the final result
1104
1105   @q = (); ($v2,$v1) = @$y[-2,-1];
1106   $v2 = 0 unless $v2;
1107   while ($#$x > $#$y) 
1108     {
1109     ($u2,$u1,$u0) = @$x[-3..-1];
1110     $u2 = 0 unless $u2;
1111     #warn "oups v1 is 0, u0: $u0 $y->[-2] $y->[-1] l ",scalar @$y,"\n"
1112     # if $v1 == 0;
1113     $q = (($u0 == $v1) ? $MAX_VAL : int(($u0*$BASE+$u1)/$v1));
1114     --$q while ($v2*$q > ($u0*$BASE+$u1-$q*$v1)*$BASE+$u2);
1115     if ($q)
1116       {
1117       ($car, $bar) = (0,0);
1118       for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi) 
1119         {
1120         $prd = $q * $y->[$yi] + $car;
1121         $prd -= ($car = int($prd / $BASE)) * $BASE;
1122         $x->[$xi] += $BASE if ($bar = (($x->[$xi] -= $prd + $bar) < 0));
1123         }
1124       if ($x->[-1] < $car + $bar) 
1125         {
1126         $car = 0; --$q;
1127         for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi) 
1128           {
1129           $x->[$xi] -= $BASE
1130            if ($car = (($x->[$xi] += $y->[$yi] + $car) >= $BASE));
1131           }
1132         }   
1133       }
1134     pop(@$x); unshift(@q, $q);
1135     }
1136   if (wantarray) 
1137     {
1138     @d = ();
1139     if ($dd != 1)  
1140       {
1141       $car = 0; 
1142       for $xi (reverse @$x) 
1143         {
1144         $prd = $car * $BASE + $xi;
1145         $car = $prd - ($tmp = int($prd / $dd)) * $dd;
1146         unshift(@d, $tmp);
1147         }
1148       }
1149     else 
1150       {
1151       @d = @$x;
1152       }
1153     @$x = @q;
1154     my $d = \@d; 
1155     __strip_zeros($x);
1156     __strip_zeros($d);
1157     return ($x,$d);
1158     }
1159   @$x = @q;
1160   __strip_zeros($x);
1161   $x;
1162   }
1163
1164 ##############################################################################
1165 # testing
1166
1167 sub _acmp
1168   {
1169   # internal absolute post-normalized compare (ignore signs)
1170   # ref to array, ref to array, return <0, 0, >0
1171   # arrays must have at least one entry; this is not checked for
1172   my ($c,$cx,$cy) = @_;
1173  
1174   # shortcut for short numbers 
1175   return (($cx->[0] <=> $cy->[0]) <=> 0) 
1176    if scalar @$cx == scalar @$cy && scalar @$cx == 1;
1177
1178   # fast comp based on number of array elements (aka pseudo-length)
1179   my $lxy = (scalar @$cx - scalar @$cy)
1180   # or length of first element if same number of elements (aka difference 0)
1181     ||
1182   # need int() here because sometimes the last element is '00018' vs '18'
1183    (length(int($cx->[-1])) - length(int($cy->[-1])));
1184   return -1 if $lxy < 0;                                # already differs, ret
1185   return 1 if $lxy > 0;                                 # ditto
1186
1187   # manual way (abort if unequal, good for early ne)
1188   my $a; my $j = scalar @$cx;
1189   while (--$j >= 0)
1190     {
1191     last if ($a = $cx->[$j] - $cy->[$j]);
1192     }
1193   $a <=> 0;
1194   }
1195
1196 sub _len
1197   {
1198   # compute number of digits in base 10
1199
1200   # int() because add/sub sometimes leaves strings (like '00005') instead of
1201   # '5' in this place, thus causing length() to report wrong length
1202   my $cx = $_[1];
1203
1204   (@$cx-1)*$BASE_LEN+length(int($cx->[-1]));
1205   }
1206
1207 sub _digit
1208   {
1209   # Return the nth digit. Zero is rightmost, so _digit(123,0) gives 3.
1210   # Negative values count from the left, so _digit(123, -1) gives 1.
1211   my ($c,$x,$n) = @_;
1212
1213   my $len = _len('',$x);
1214
1215   $n += $len if $n < 0;                 # -1 last, -2 second-to-last
1216   return "0" if $n < 0 || $n >= $len;   # return 0 for digits out of range
1217
1218   my $elem = int($n / $BASE_LEN);       # which array element
1219   my $digit = $n % $BASE_LEN;           # which digit in this element
1220   substr("$x->[$elem]", -$digit-1, 1);
1221   }
1222
1223 sub _zeros
1224   {
1225   # return amount of trailing zeros in decimal
1226   # check each array elem in _m for having 0 at end as long as elem == 0
1227   # Upon finding a elem != 0, stop
1228   my $x = $_[1];
1229
1230   return 0 if scalar @$x == 1 && $x->[0] == 0;
1231
1232   my $zeros = 0; my $elem;
1233   foreach my $e (@$x)
1234     {
1235     if ($e != 0)
1236       {
1237       $elem = "$e";                             # preserve x
1238       $elem =~ s/.*?(0*$)/$1/;                  # strip anything not zero
1239       $zeros *= $BASE_LEN;                      # elems * 5
1240       $zeros += length($elem);                  # count trailing zeros
1241       last;                                     # early out
1242       }
1243     $zeros ++;                                  # real else branch: 50% slower!
1244     }
1245   $zeros;
1246   }
1247
1248 ##############################################################################
1249 # _is_* routines
1250
1251 sub _is_zero
1252   {
1253   # return true if arg is zero 
1254   (((scalar @{$_[1]} == 1) && ($_[1]->[0] == 0))) <=> 0;
1255   }
1256
1257 sub _is_even
1258   {
1259   # return true if arg is even
1260   (!($_[1]->[0] & 1)) <=> 0; 
1261   }
1262
1263 sub _is_odd
1264   {
1265   # return true if arg is odd
1266   (($_[1]->[0] & 1)) <=> 0;
1267   }
1268
1269 sub _is_one
1270   {
1271   # return true if arg is one
1272   (scalar @{$_[1]} == 1) && ($_[1]->[0] == 1) <=> 0; 
1273   }
1274
1275 sub _is_two
1276   {
1277   # return true if arg is two 
1278   (scalar @{$_[1]} == 1) && ($_[1]->[0] == 2) <=> 0; 
1279   }
1280
1281 sub _is_ten
1282   {
1283   # return true if arg is ten 
1284   (scalar @{$_[1]} == 1) && ($_[1]->[0] == 10) <=> 0; 
1285   }
1286
1287 sub __strip_zeros
1288   {
1289   # internal normalization function that strips leading zeros from the array
1290   # args: ref to array
1291   my $s = shift;
1292  
1293   my $cnt = scalar @$s; # get count of parts
1294   my $i = $cnt-1;
1295   push @$s,0 if $i < 0;         # div might return empty results, so fix it
1296
1297   return $s if @$s == 1;                # early out
1298
1299   #print "strip: cnt $cnt i $i\n";
1300   # '0', '3', '4', '0', '0',
1301   #  0    1    2    3    4
1302   # cnt = 5, i = 4
1303   # i = 4
1304   # i = 3
1305   # => fcnt = cnt - i (5-2 => 3, cnt => 5-1 = 4, throw away from 4th pos)
1306   # >= 1: skip first part (this can be zero)
1307   while ($i > 0) { last if $s->[$i] != 0; $i--; }
1308   $i++; splice @$s,$i if ($i < $cnt); # $i cant be 0
1309   $s;                                                                    
1310   }                                                                             
1311
1312 ###############################################################################
1313 # check routine to test internal state for corruptions
1314
1315 sub _check
1316   {
1317   # used by the test suite
1318   my $x = $_[1];
1319
1320   return "$x is not a reference" if !ref($x);
1321
1322   # are all parts are valid?
1323   my $i = 0; my $j = scalar @$x; my ($e,$try);
1324   while ($i < $j)
1325     {
1326     $e = $x->[$i]; $e = 'undef' unless defined $e;
1327     $try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e)";
1328     last if $e !~ /^[+]?[0-9]+$/;
1329     $try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e) (stringify)";
1330     last if "$e" !~ /^[+]?[0-9]+$/;
1331     $try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e) (cat-stringify)";
1332     last if '' . "$e" !~ /^[+]?[0-9]+$/;
1333     $try = ' < 0 || >= $BASE; '."($x, $e)";
1334     last if $e <0 || $e >= $BASE;
1335     # this test is disabled, since new/bnorm and certain ops (like early out
1336     # in add/sub) are allowed/expected to leave '00000' in some elements
1337     #$try = '=~ /^00+/; '."($x, $e)";
1338     #last if $e =~ /^00+/;
1339     $i++;
1340     }
1341   return "Illegal part '$e' at pos $i (tested: $try)" if $i < $j;
1342   0;
1343   }
1344
1345
1346 ###############################################################################
1347
1348 sub _mod
1349   {
1350   # if possible, use mod shortcut
1351   my ($c,$x,$yo) = @_;
1352
1353   # slow way since $y to big
1354   if (scalar @$yo > 1)
1355     {
1356     my ($xo,$rem) = _div($c,$x,$yo);
1357     return $rem;
1358     }
1359
1360   my $y = $yo->[0];
1361   # both are single element arrays
1362   if (scalar @$x == 1)
1363     {
1364     $x->[0] %= $y;
1365     return $x;
1366     }
1367
1368   # @y is a single element, but @x has more than one element
1369   my $b = $BASE % $y;
1370   if ($b == 0)
1371     {
1372     # when BASE % Y == 0 then (B * BASE) % Y == 0
1373     # (B * BASE) % $y + A % Y => A % Y
1374     # so need to consider only last element: O(1)
1375     $x->[0] %= $y;
1376     }
1377   elsif ($b == 1)
1378     {
1379     # else need to go through all elements: O(N), but loop is a bit simplified
1380     my $r = 0;
1381     foreach (@$x)
1382       {
1383       $r = ($r + $_) % $y;              # not much faster, but heh...
1384       #$r += $_ % $y; $r %= $y;
1385       }
1386     $r = 0 if $r == $y;
1387     $x->[0] = $r;
1388     }
1389   else
1390     {
1391     # else need to go through all elements: O(N)
1392     my $r = 0; my $bm = 1;
1393     foreach (@$x)
1394       {
1395       $r = ($_ * $bm + $r) % $y;
1396       $bm = ($bm * $b) % $y;
1397
1398       #$r += ($_ % $y) * $bm;
1399       #$bm *= $b;
1400       #$bm %= $y;
1401       #$r %= $y;
1402       }
1403     $r = 0 if $r == $y;
1404     $x->[0] = $r;
1405     }
1406   splice (@$x,1);               # keep one element of $x
1407   $x;
1408   }
1409
1410 ##############################################################################
1411 # shifts
1412
1413 sub _rsft
1414   {
1415   my ($c,$x,$y,$n) = @_;
1416
1417   if ($n != 10)
1418     {
1419     $n = _new($c,$n); return _div($c,$x, _pow($c,$n,$y));
1420     }
1421
1422   # shortcut (faster) for shifting by 10)
1423   # multiples of $BASE_LEN
1424   my $dst = 0;                          # destination
1425   my $src = _num($c,$y);                # as normal int
1426   my $xlen = (@$x-1)*$BASE_LEN+length(int($x->[-1]));  # len of x in digits
1427   if ($src >= $xlen or ($src == $xlen and ! defined $x->[1]))
1428     {
1429     # 12345 67890 shifted right by more than 10 digits => 0
1430     splice (@$x,1);                    # leave only one element
1431     $x->[0] = 0;                       # set to zero
1432     return $x;
1433     }
1434   my $rem = $src % $BASE_LEN;           # remainder to shift
1435   $src = int($src / $BASE_LEN);         # source
1436   if ($rem == 0)
1437     {
1438     splice (@$x,0,$src);                # even faster, 38.4 => 39.3
1439     }
1440   else
1441     {
1442     my $len = scalar @$x - $src;        # elems to go
1443     my $vd; my $z = '0'x $BASE_LEN;
1444     $x->[scalar @$x] = 0;               # avoid || 0 test inside loop
1445     while ($dst < $len)
1446       {
1447       $vd = $z.$x->[$src];
1448       $vd = substr($vd,-$BASE_LEN,$BASE_LEN-$rem);
1449       $src++;
1450       $vd = substr($z.$x->[$src],-$rem,$rem) . $vd;
1451       $vd = substr($vd,-$BASE_LEN,$BASE_LEN) if length($vd) > $BASE_LEN;
1452       $x->[$dst] = int($vd);
1453       $dst++;
1454       }
1455     splice (@$x,$dst) if $dst > 0;              # kill left-over array elems
1456     pop @$x if $x->[-1] == 0 && @$x > 1;        # kill last element if 0
1457     } # else rem == 0
1458   $x;
1459   }
1460
1461 sub _lsft
1462   {
1463   my ($c,$x,$y,$n) = @_;
1464
1465   if ($n != 10)
1466     {
1467     $n = _new($c,$n); return _mul($c,$x, _pow($c,$n,$y));
1468     }
1469
1470   # shortcut (faster) for shifting by 10) since we are in base 10eX
1471   # multiples of $BASE_LEN:
1472   my $src = scalar @$x;                 # source
1473   my $len = _num($c,$y);                # shift-len as normal int
1474   my $rem = $len % $BASE_LEN;           # remainder to shift
1475   my $dst = $src + int($len/$BASE_LEN); # destination
1476   my $vd;                               # further speedup
1477   $x->[$src] = 0;                       # avoid first ||0 for speed
1478   my $z = '0' x $BASE_LEN;
1479   while ($src >= 0)
1480     {
1481     $vd = $x->[$src]; $vd = $z.$vd;
1482     $vd = substr($vd,-$BASE_LEN+$rem,$BASE_LEN-$rem);
1483     $vd .= $src > 0 ? substr($z.$x->[$src-1],-$BASE_LEN,$rem) : '0' x $rem;
1484     $vd = substr($vd,-$BASE_LEN,$BASE_LEN) if length($vd) > $BASE_LEN;
1485     $x->[$dst] = int($vd);
1486     $dst--; $src--;
1487     }
1488   # set lowest parts to 0
1489   while ($dst >= 0) { $x->[$dst--] = 0; }
1490   # fix spurios last zero element
1491   splice @$x,-1 if $x->[-1] == 0;
1492   $x;
1493   }
1494
1495 sub _pow
1496   {
1497   # power of $x to $y
1498   # ref to array, ref to array, return ref to array
1499   my ($c,$cx,$cy) = @_;
1500
1501   if (scalar @$cy == 1 && $cy->[0] == 0)
1502     {
1503     splice (@$cx,1); $cx->[0] = 1;              # y == 0 => x => 1
1504     return $cx;
1505     }
1506   if ((scalar @$cx == 1 && $cx->[0] == 1) ||    #    x == 1
1507       (scalar @$cy == 1 && $cy->[0] == 1))      # or y == 1
1508     {
1509     return $cx;
1510     }
1511   if (scalar @$cx == 1 && $cx->[0] == 0)
1512     {
1513     splice (@$cx,1); $cx->[0] = 0;              # 0 ** y => 0 (if not y <= 0)
1514     return $cx;
1515     }
1516
1517   my $pow2 = _one();
1518
1519   my $y_bin = _as_bin($c,$cy); $y_bin =~ s/^0b//;
1520   my $len = length($y_bin);
1521   while (--$len > 0)
1522     {
1523     _mul($c,$pow2,$cx) if substr($y_bin,$len,1) eq '1';         # is odd?
1524     _mul($c,$cx,$cx);
1525     }
1526
1527   _mul($c,$cx,$pow2);
1528   $cx;
1529   }
1530
1531 sub _nok
1532   {
1533   # n over k
1534   # ref to array, return ref to array
1535   my ($c,$n,$k) = @_;
1536
1537   # ( 7 )       7!       1*2*3*4 * 5*6*7   5 * 6 * 7       6   7
1538   # ( - ) = --------- =  --------------- = --------- = 5 * - * -
1539   # ( 3 )   (7-3)! 3!    1*2*3*4 * 1*2*3   1 * 2 * 3       2   3
1540
1541   if (!_is_zero($c,$k))
1542     {
1543     my $x = _copy($c,$n);
1544     _sub($c,$n,$k);
1545     _inc($c,$n);
1546     my $f = _copy($c,$n); _inc($c,$f);          # n = 5, f = 6, d = 2
1547     my $d = _two($c);
1548     while (_acmp($c,$f,$x) <= 0)                # f <= n ?
1549       {
1550       # n = (n * f / d) == 5 * 6 / 2
1551       $n = _mul($c,$n,$f); $n = _div($c,$n,$d);
1552       # f = 7, d = 3
1553       _inc($c,$f); _inc($c,$d);
1554       }
1555     }
1556   else
1557     {
1558     # keep ref to $n and set it to 1
1559     splice (@$n,1); $n->[0] = 1;
1560     }
1561   $n;
1562   }
1563
1564 my @factorials = (
1565   1,
1566   1,
1567   2,
1568   2*3,
1569   2*3*4,
1570   2*3*4*5,
1571   2*3*4*5*6,
1572   2*3*4*5*6*7,
1573 );
1574
1575 sub _fac
1576   {
1577   # factorial of $x
1578   # ref to array, return ref to array
1579   my ($c,$cx) = @_;
1580
1581   if ((@$cx == 1) && ($cx->[0] <= 7))
1582     {
1583     $cx->[0] = $factorials[$cx->[0]];           # 0 => 1, 1 => 1, 2 => 2 etc.
1584     return $cx;
1585     }
1586
1587   if ((@$cx == 1) &&            # we do this only if $x >= 12 and $x <= 7000
1588       ($cx->[0] >= 12 && $cx->[0] < 7000))
1589     {
1590
1591   # Calculate (k-j) * (k-j+1) ... k .. (k+j-1) * (k + j)
1592   # See http://blogten.blogspot.com/2007/01/calculating-n.html
1593   # The above series can be expressed as factors:
1594   #   k * k - (j - i) * 2
1595   # We cache k*k, and calculate (j * j) as the sum of the first j odd integers
1596
1597   # This will not work when N exceeds the storage of a Perl scalar, however,
1598   # in this case the algorithm would be way to slow to terminate, anyway.
1599
1600   # As soon as the last element of $cx is 0, we split it up and remember
1601   # how many zeors we got so far. The reason is that n! will accumulate
1602   # zeros at the end rather fast.
1603   my $zero_elements = 0;
1604
1605   # If n is even, set n = n -1
1606   my $k = _num($c,$cx); my $even = 1;
1607   if (($k & 1) == 0)
1608     {
1609     $even = $k; $k --;
1610     }
1611   # set k to the center point
1612   $k = ($k + 1) / 2;
1613 #  print "k $k even: $even\n";
1614   # now calculate k * k
1615   my $k2 = $k * $k;
1616   my $odd = 1; my $sum = 1;
1617   my $i = $k - 1;
1618   # keep reference to x
1619   my $new_x = _new($c, $k * $even);
1620   @$cx = @$new_x;
1621   if ($cx->[0] == 0)
1622     {
1623     $zero_elements ++; shift @$cx;
1624     }
1625 #  print STDERR "x = ", _str($c,$cx),"\n";
1626   my $BASE2 = int(sqrt($BASE))-1;
1627   my $j = 1; 
1628   while ($j <= $i)
1629     {
1630     my $m = ($k2 - $sum); $odd += 2; $sum += $odd; $j++;
1631     while ($j <= $i && ($m < $BASE2) && (($k2 - $sum) < $BASE2))
1632       {
1633       $m *= ($k2 - $sum);
1634       $odd += 2; $sum += $odd; $j++;
1635 #      print STDERR "\n k2 $k2 m $m sum $sum odd $odd\n"; sleep(1);
1636       }
1637     if ($m < $BASE)
1638       {
1639       _mul($c,$cx,[$m]);
1640       }
1641     else
1642       {
1643       _mul($c,$cx,$c->_new($m));
1644       }
1645     if ($cx->[0] == 0)
1646       {
1647       $zero_elements ++; shift @$cx;
1648       }
1649 #    print STDERR "Calculate $k2 - $sum = $m (x = ", _str($c,$cx),")\n";
1650     }
1651   # multiply in the zeros again
1652   unshift @$cx, (0) x $zero_elements; 
1653   return $cx;
1654   }
1655
1656   # go forward until $base is exceeded
1657   # limit is either $x steps (steps == 100 means a result always too high) or
1658   # $base.
1659   my $steps = 100; $steps = $cx->[0] if @$cx == 1;
1660   my $r = 2; my $cf = 3; my $step = 2; my $last = $r;
1661   while ($r*$cf < $BASE && $step < $steps)
1662     {
1663     $last = $r; $r *= $cf++; $step++;
1664     }
1665   if ((@$cx == 1) && $step == $cx->[0])
1666     {
1667     # completely done, so keep reference to $x and return
1668     $cx->[0] = $r;
1669     return $cx;
1670     }
1671   
1672   # now we must do the left over steps
1673   my $n;                                        # steps still to do
1674   if (scalar @$cx == 1)
1675     {
1676     $n = $cx->[0];
1677     }
1678   else
1679     {
1680     $n = _copy($c,$cx);
1681     }
1682
1683   # Set $cx to the last result below $BASE (but keep ref to $x)
1684   $cx->[0] = $last; splice (@$cx,1);
1685   # As soon as the last element of $cx is 0, we split it up and remember
1686   # how many zeors we got so far. The reason is that n! will accumulate
1687   # zeros at the end rather fast.
1688   my $zero_elements = 0;
1689
1690   # do left-over steps fit into a scalar?
1691   if (ref $n eq 'ARRAY')
1692     {
1693     # No, so use slower inc() & cmp()
1694     # ($n is at least $BASE here)
1695     my $base_2 = int(sqrt($BASE)) - 1;
1696     #print STDERR "base_2: $base_2\n"; 
1697     while ($step < $base_2)
1698       {
1699       if ($cx->[0] == 0)
1700         {
1701         $zero_elements ++; shift @$cx;
1702         }
1703       my $b = $step * ($step + 1); $step += 2;
1704       _mul($c,$cx,[$b]);
1705       }
1706     $step = [$step];
1707     while (_acmp($c,$step,$n) <= 0)
1708       {
1709       if ($cx->[0] == 0)
1710         {
1711         $zero_elements ++; shift @$cx;
1712         }
1713       _mul($c,$cx,$step); _inc($c,$step);
1714       }
1715     }
1716   else
1717     {
1718     # Yes, so we can speed it up slightly
1719   
1720 #    print "# left over steps $n\n";
1721
1722     my $base_4 = int(sqrt(sqrt($BASE))) - 2;
1723     #print STDERR "base_4: $base_4\n";
1724     my $n4 = $n - 4; 
1725     while ($step < $n4 && $step < $base_4)
1726       {
1727       if ($cx->[0] == 0)
1728         {
1729         $zero_elements ++; shift @$cx;
1730         }
1731       my $b = $step * ($step + 1); $step += 2; $b *= $step * ($step + 1); $step += 2;
1732       _mul($c,$cx,[$b]);
1733       }
1734     my $base_2 = int(sqrt($BASE)) - 1;
1735     my $n2 = $n - 2; 
1736     #print STDERR "base_2: $base_2\n"; 
1737     while ($step < $n2 && $step < $base_2)
1738       {
1739       if ($cx->[0] == 0)
1740         {
1741         $zero_elements ++; shift @$cx;
1742         }
1743       my $b = $step * ($step + 1); $step += 2;
1744       _mul($c,$cx,[$b]);
1745       }
1746     # do what's left over
1747     while ($step <= $n)
1748       {
1749       _mul($c,$cx,[$step]); $step++;
1750       if ($cx->[0] == 0)
1751         {
1752         $zero_elements ++; shift @$cx;
1753         }
1754       }
1755     }
1756   # multiply in the zeros again
1757   unshift @$cx, (0) x $zero_elements;
1758   $cx;                  # return result
1759   }
1760
1761 #############################################################################
1762
1763 sub _log_int
1764   {
1765   # calculate integer log of $x to base $base
1766   # ref to array, ref to array - return ref to array
1767   my ($c,$x,$base) = @_;
1768
1769   # X == 0 => NaN
1770   return if (scalar @$x == 1 && $x->[0] == 0);
1771   # BASE 0 or 1 => NaN
1772   return if (scalar @$base == 1 && $base->[0] < 2);
1773   my $cmp = _acmp($c,$x,$base); # X == BASE => 1
1774   if ($cmp == 0)
1775     {
1776     splice (@$x,1); $x->[0] = 1;
1777     return ($x,1)
1778     }
1779   # X < BASE
1780   if ($cmp < 0)
1781     {
1782     splice (@$x,1); $x->[0] = 0;
1783     return ($x,undef);
1784     }
1785
1786   my $x_org = _copy($c,$x);             # preserve x
1787   splice(@$x,1); $x->[0] = 1;           # keep ref to $x
1788
1789   # Compute a guess for the result based on:
1790   # $guess = int ( length_in_base_10(X) / ( log(base) / log(10) ) )
1791   my $len = _len($c,$x_org);
1792   my $log = log($base->[-1]) / log(10);
1793
1794   # for each additional element in $base, we add $BASE_LEN to the result,
1795   # based on the observation that log($BASE,10) is BASE_LEN and
1796   # log(x*y) == log(x) + log(y):
1797   $log += ((scalar @$base)-1) * $BASE_LEN;
1798
1799   # calculate now a guess based on the values obtained above:
1800   my $res = int($len / $log);
1801
1802   $x->[0] = $res;
1803   my $trial = _pow ($c, _copy($c, $base), $x);
1804   my $a = _acmp($c,$trial,$x_org);
1805
1806 #  print STDERR "# trial ", _str($c,$x)," was: $a (0 = exact, -1 too small, +1 too big)\n";
1807
1808   # found an exact result?
1809   return ($x,1) if $a == 0;
1810
1811   if ($a > 0)
1812     {
1813     # or too big
1814     _div($c,$trial,$base); _dec($c, $x);
1815     while (($a = _acmp($c,$trial,$x_org)) > 0)
1816       {
1817 #      print STDERR "# big _log_int at ", _str($c,$x), "\n"; 
1818       _div($c,$trial,$base); _dec($c, $x);
1819       }
1820     # result is now exact (a == 0), or too small (a < 0)
1821     return ($x, $a == 0 ? 1 : 0);
1822     }
1823
1824   # else: result was to small
1825   _mul($c,$trial,$base);
1826
1827   # did we now get the right result?
1828   $a = _acmp($c,$trial,$x_org);
1829
1830   if ($a == 0)                          # yes, exactly
1831     {
1832     _inc($c, $x);
1833     return ($x,1); 
1834     }
1835   return ($x,0) if $a > 0;  
1836
1837   # Result still too small (we should come here only if the estimate above
1838   # was very off base):
1839  
1840   # Now let the normal trial run obtain the real result
1841   # Simple loop that increments $x by 2 in each step, possible overstepping
1842   # the real result
1843
1844   my $base_mul = _mul($c, _copy($c,$base), $base);      # $base * $base
1845
1846   while (($a = _acmp($c,$trial,$x_org)) < 0)
1847     {
1848 #    print STDERR "# small _log_int at ", _str($c,$x), "\n"; 
1849     _mul($c,$trial,$base_mul); _add($c, $x, [2]);
1850     }
1851
1852   my $exact = 1;
1853   if ($a > 0)
1854     {
1855     # overstepped the result
1856     _dec($c, $x);
1857     _div($c,$trial,$base);
1858     $a = _acmp($c,$trial,$x_org);
1859     if ($a > 0)
1860       {
1861       _dec($c, $x);
1862       }
1863     $exact = 0 if $a != 0;              # a = -1 => not exact result, a = 0 => exact
1864     }
1865   
1866   ($x,$exact);                          # return result
1867   }
1868
1869 # for debugging:
1870   use constant DEBUG => 0;
1871   my $steps = 0;
1872   sub steps { $steps };
1873
1874 sub _sqrt
1875   {
1876   # square-root of $x in place
1877   # Compute a guess of the result (by rule of thumb), then improve it via
1878   # Newton's method.
1879   my ($c,$x) = @_;
1880
1881   if (scalar @$x == 1)
1882     {
1883     # fits into one Perl scalar, so result can be computed directly
1884     $x->[0] = int(sqrt($x->[0]));
1885     return $x;
1886     } 
1887   my $y = _copy($c,$x);
1888   # hopefully _len/2 is < $BASE, the -1 is to always undershot the guess
1889   # since our guess will "grow"
1890   my $l = int((_len($c,$x)-1) / 2);     
1891
1892   my $lastelem = $x->[-1];                                      # for guess
1893   my $elems = scalar @$x - 1;
1894   # not enough digits, but could have more?
1895   if ((length($lastelem) <= 3) && ($elems > 1))
1896     {
1897     # right-align with zero pad
1898     my $len = length($lastelem) & 1;
1899     print "$lastelem => " if DEBUG;
1900     $lastelem .= substr($x->[-2] . '0' x $BASE_LEN,0,$BASE_LEN);
1901     # former odd => make odd again, or former even to even again
1902     $lastelem = $lastelem / 10 if (length($lastelem) & 1) != $len;
1903     print "$lastelem\n" if DEBUG;
1904     }
1905
1906   # construct $x (instead of _lsft($c,$x,$l,10)
1907   my $r = $l % $BASE_LEN;       # 10000 00000 00000 00000 ($BASE_LEN=5)
1908   $l = int($l / $BASE_LEN);
1909   print "l =  $l " if DEBUG;
1910
1911   splice @$x,$l;                # keep ref($x), but modify it
1912
1913   # we make the first part of the guess not '1000...0' but int(sqrt($lastelem))
1914   # that gives us:
1915   # 14400 00000 => sqrt(14400) => guess first digits to be 120
1916   # 144000 000000 => sqrt(144000) => guess 379
1917
1918   print "$lastelem (elems $elems) => " if DEBUG;
1919   $lastelem = $lastelem / 10 if ($elems & 1 == 1);              # odd or even?
1920   my $g = sqrt($lastelem); $g =~ s/\.//;                        # 2.345 => 2345
1921   $r -= 1 if $elems & 1 == 0;                                   # 70 => 7
1922
1923   # padd with zeros if result is too short
1924   $x->[$l--] = int(substr($g . '0' x $r,0,$r+1));
1925   print "now ",$x->[-1] if DEBUG;
1926   print " would have been ", int('1' . '0' x $r),"\n" if DEBUG;
1927
1928   # If @$x > 1, we could compute the second elem of the guess, too, to create
1929   # an even better guess. Not implemented yet. Does it improve performance?
1930   $x->[$l--] = 0 while ($l >= 0);       # all other digits of guess are zero
1931
1932   print "start x= ",_str($c,$x),"\n" if DEBUG;
1933   my $two = _two();
1934   my $last = _zero();
1935   my $lastlast = _zero();
1936   $steps = 0 if DEBUG;
1937   while (_acmp($c,$last,$x) != 0 && _acmp($c,$lastlast,$x) != 0)
1938     {
1939     $steps++ if DEBUG;
1940     $lastlast = _copy($c,$last);
1941     $last = _copy($c,$x);
1942     _add($c,$x, _div($c,_copy($c,$y),$x));
1943     _div($c,$x, $two );
1944     print " x= ",_str($c,$x),"\n" if DEBUG;
1945     }
1946   print "\nsteps in sqrt: $steps, " if DEBUG;
1947   _dec($c,$x) if _acmp($c,$y,_mul($c,_copy($c,$x),$x)) < 0;     # overshot? 
1948   print " final ",$x->[-1],"\n" if DEBUG;
1949   $x;
1950   }
1951
1952 sub _root
1953   {
1954   # take n'th root of $x in place (n >= 3)
1955   my ($c,$x,$n) = @_;
1956  
1957   if (scalar @$x == 1)
1958     {
1959     if (scalar @$n > 1)
1960       {
1961       # result will always be smaller than 2 so trunc to 1 at once
1962       $x->[0] = 1;
1963       }
1964     else
1965       {
1966       # fits into one Perl scalar, so result can be computed directly
1967       # cannot use int() here, because it rounds wrongly (try 
1968       # (81 ** 3) ** (1/3) to see what I mean)
1969       #$x->[0] = int( $x->[0] ** (1 / $n->[0]) );
1970       # round to 8 digits, then truncate result to integer
1971       $x->[0] = int ( sprintf ("%.8f", $x->[0] ** (1 / $n->[0]) ) );
1972       }
1973     return $x;
1974     } 
1975
1976   # we know now that X is more than one element long
1977
1978   # if $n is a power of two, we can repeatedly take sqrt($X) and find the
1979   # proper result, because sqrt(sqrt($x)) == root($x,4)
1980   my $b = _as_bin($c,$n);
1981   if ($b =~ /0b1(0+)$/)
1982     {
1983     my $count = CORE::length($1);       # 0b100 => len('00') => 2
1984     my $cnt = $count;                   # counter for loop
1985     unshift (@$x, 0);                   # add one element, together with one
1986                                         # more below in the loop this makes 2
1987     while ($cnt-- > 0)
1988       {
1989       # 'inflate' $X by adding one element, basically computing
1990       # $x * $BASE * $BASE. This gives us more $BASE_LEN digits for result
1991       # since len(sqrt($X)) approx == len($x) / 2.
1992       unshift (@$x, 0);
1993       # calculate sqrt($x), $x is now one element to big, again. In the next
1994       # round we make that two, again.
1995       _sqrt($c,$x);
1996       }
1997     # $x is now one element to big, so truncate result by removing it
1998     splice (@$x,0,1);
1999     } 
2000   else
2001     {
2002     # trial computation by starting with 2,4,8,16 etc until we overstep
2003     my $step;
2004     my $trial = _two();
2005
2006     # while still to do more than X steps
2007     do
2008       {
2009       $step = _two();
2010       while (_acmp($c, _pow($c, _copy($c, $trial), $n), $x) < 0)
2011         {
2012         _mul ($c, $step, [2]);
2013         _add ($c, $trial, $step);
2014         }
2015
2016       # hit exactly?
2017       if (_acmp($c, _pow($c, _copy($c, $trial), $n), $x) == 0)
2018         {
2019         @$x = @$trial;                  # make copy while preserving ref to $x
2020         return $x;
2021         }
2022       # overstepped, so go back on step
2023       _sub($c, $trial, $step);
2024       } while (scalar @$step > 1 || $step->[0] > 128);
2025
2026     # reset step to 2
2027     $step = _two();
2028     # add two, because $trial cannot be exactly the result (otherwise we would
2029     # already have found it)
2030     _add($c, $trial, $step);
2031  
2032     # and now add more and more (2,4,6,8,10 etc)
2033     while (_acmp($c, _pow($c, _copy($c, $trial), $n), $x) < 0)
2034       {
2035       _add ($c, $trial, $step);
2036       }
2037
2038     # hit not exactly? (overstepped)
2039     if (_acmp($c, _pow($c, _copy($c, $trial), $n), $x) > 0)
2040       {
2041       _dec($c,$trial);
2042       }
2043
2044     # hit not exactly? (overstepped)
2045     # 80 too small, 81 slightly too big, 82 too big
2046     if (_acmp($c, _pow($c, _copy($c, $trial), $n), $x) > 0)
2047       {
2048       _dec ($c, $trial); 
2049       }
2050
2051     @$x = @$trial;                      # make copy while preserving ref to $x
2052     return $x;
2053     }
2054   $x; 
2055   }
2056
2057 ##############################################################################
2058 # binary stuff
2059
2060 sub _and
2061   {
2062   my ($c,$x,$y) = @_;
2063
2064   # the shortcut makes equal, large numbers _really_ fast, and makes only a
2065   # very small performance drop for small numbers (e.g. something with less
2066   # than 32 bit) Since we optimize for large numbers, this is enabled.
2067   return $x if _acmp($c,$x,$y) == 0;            # shortcut
2068   
2069   my $m = _one(); my ($xr,$yr);
2070   my $mask = $AND_MASK;
2071
2072   my $x1 = $x;
2073   my $y1 = _copy($c,$y);                        # make copy
2074   $x = _zero();
2075   my ($b,$xrr,$yrr);
2076   use integer;
2077   while (!_is_zero($c,$x1) && !_is_zero($c,$y1))
2078     {
2079     ($x1, $xr) = _div($c,$x1,$mask);
2080     ($y1, $yr) = _div($c,$y1,$mask);
2081
2082     # make ints() from $xr, $yr
2083     # this is when the AND_BITS are greater than $BASE and is slower for
2084     # small (<256 bits) numbers, but faster for large numbers. Disabled
2085     # due to KISS principle
2086
2087 #    $b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
2088 #    $b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
2089 #    _add($c,$x, _mul($c, _new( $c, ($xrr & $yrr) ), $m) );
2090     
2091     # 0+ due to '&' doesn't work in strings
2092     _add($c,$x, _mul($c, [ 0+$xr->[0] & 0+$yr->[0] ], $m) );
2093     _mul($c,$m,$mask);
2094     }
2095   $x;
2096   }
2097
2098 sub _xor
2099   {
2100   my ($c,$x,$y) = @_;
2101
2102   return _zero() if _acmp($c,$x,$y) == 0;       # shortcut (see -and)
2103
2104   my $m = _one(); my ($xr,$yr);
2105   my $mask = $XOR_MASK;
2106
2107   my $x1 = $x;
2108   my $y1 = _copy($c,$y);                        # make copy
2109   $x = _zero();
2110   my ($b,$xrr,$yrr);
2111   use integer;
2112   while (!_is_zero($c,$x1) && !_is_zero($c,$y1))
2113     {
2114     ($x1, $xr) = _div($c,$x1,$mask);
2115     ($y1, $yr) = _div($c,$y1,$mask);
2116     # make ints() from $xr, $yr (see _and())
2117     #$b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
2118     #$b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
2119     #_add($c,$x, _mul($c, _new( $c, ($xrr ^ $yrr) ), $m) );
2120
2121     # 0+ due to '^' doesn't work in strings
2122     _add($c,$x, _mul($c, [ 0+$xr->[0] ^ 0+$yr->[0] ], $m) );
2123     _mul($c,$m,$mask);
2124     }
2125   # the loop stops when the shorter of the two numbers is exhausted
2126   # the remainder of the longer one will survive bit-by-bit, so we simple
2127   # multiply-add it in
2128   _add($c,$x, _mul($c, $x1, $m) ) if !_is_zero($c,$x1);
2129   _add($c,$x, _mul($c, $y1, $m) ) if !_is_zero($c,$y1);
2130   
2131   $x;
2132   }
2133
2134 sub _or
2135   {
2136   my ($c,$x,$y) = @_;
2137
2138   return $x if _acmp($c,$x,$y) == 0;            # shortcut (see _and)
2139
2140   my $m = _one(); my ($xr,$yr);
2141   my $mask = $OR_MASK;
2142
2143   my $x1 = $x;
2144   my $y1 = _copy($c,$y);                        # make copy
2145   $x = _zero();
2146   my ($b,$xrr,$yrr);
2147   use integer;
2148   while (!_is_zero($c,$x1) && !_is_zero($c,$y1))
2149     {
2150     ($x1, $xr) = _div($c,$x1,$mask);
2151     ($y1, $yr) = _div($c,$y1,$mask);
2152     # make ints() from $xr, $yr (see _and())
2153 #    $b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
2154 #    $b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
2155 #    _add($c,$x, _mul($c, _new( $c, ($xrr | $yrr) ), $m) );
2156     
2157     # 0+ due to '|' doesn't work in strings
2158     _add($c,$x, _mul($c, [ 0+$xr->[0] | 0+$yr->[0] ], $m) );
2159     _mul($c,$m,$mask);
2160     }
2161   # the loop stops when the shorter of the two numbers is exhausted
2162   # the remainder of the longer one will survive bit-by-bit, so we simple
2163   # multiply-add it in
2164   _add($c,$x, _mul($c, $x1, $m) ) if !_is_zero($c,$x1);
2165   _add($c,$x, _mul($c, $y1, $m) ) if !_is_zero($c,$y1);
2166   
2167   $x;
2168   }
2169
2170 sub _as_hex
2171   {
2172   # convert a decimal number to hex (ref to array, return ref to string)
2173   my ($c,$x) = @_;
2174
2175   # fits into one element (handle also 0x0 case)
2176   return sprintf("0x%x",$x->[0]) if @$x == 1;
2177
2178   my $x1 = _copy($c,$x);
2179
2180   my $es = '';
2181   my ($xr, $h, $x10000);
2182   if ($] >= 5.006)
2183     {
2184     $x10000 = [ 0x10000 ]; $h = 'h4';
2185     }
2186   else
2187     {
2188     $x10000 = [ 0x1000 ]; $h = 'h3';
2189     }
2190   while (@$x1 != 1 || $x1->[0] != 0)            # _is_zero()
2191     {
2192     ($x1, $xr) = _div($c,$x1,$x10000);
2193     $es .= unpack($h,pack('V',$xr->[0]));
2194     }
2195   $es = reverse $es;
2196   $es =~ s/^[0]+//;   # strip leading zeros
2197   '0x' . $es;                                   # return result prepended with 0x
2198   }
2199
2200 sub _as_bin
2201   {
2202   # convert a decimal number to bin (ref to array, return ref to string)
2203   my ($c,$x) = @_;
2204
2205   # fits into one element (and Perl recent enough), handle also 0b0 case
2206   # handle zero case for older Perls
2207   if ($] <= 5.005 && @$x == 1 && $x->[0] == 0)
2208     {
2209     my $t = '0b0'; return $t;
2210     }
2211   if (@$x == 1 && $] >= 5.006)
2212     {
2213     my $t = sprintf("0b%b",$x->[0]);
2214     return $t;
2215     }
2216   my $x1 = _copy($c,$x);
2217
2218   my $es = '';
2219   my ($xr, $b, $x10000);
2220   if ($] >= 5.006)
2221     {
2222     $x10000 = [ 0x10000 ]; $b = 'b16';
2223     }
2224   else
2225     {
2226     $x10000 = [ 0x1000 ]; $b = 'b12';
2227     }
2228   while (!(@$x1 == 1 && $x1->[0] == 0))         # _is_zero()
2229     {
2230     ($x1, $xr) = _div($c,$x1,$x10000);
2231     $es .= unpack($b,pack('v',$xr->[0]));
2232     }
2233   $es = reverse $es;
2234   $es =~ s/^[0]+//;   # strip leading zeros
2235   '0b' . $es;                                   # return result prepended with 0b
2236   }
2237
2238 sub _as_oct
2239   {
2240   # convert a decimal number to octal (ref to array, return ref to string)
2241   my ($c,$x) = @_;
2242
2243   # fits into one element (handle also 0 case)
2244   return sprintf("0%o",$x->[0]) if @$x == 1;
2245
2246   my $x1 = _copy($c,$x);
2247
2248   my $es = '';
2249   my $xr;
2250   my $x1000 = [ 0100000 ];
2251   while (@$x1 != 1 || $x1->[0] != 0)            # _is_zero()
2252     {
2253     ($x1, $xr) = _div($c,$x1,$x1000);
2254     $es .= reverse sprintf("%05o", $xr->[0]);
2255     }
2256   $es = reverse $es;
2257   $es =~ s/^[0]+//;   # strip leading zeros
2258   '0' . $es;                                    # return result prepended with 0
2259   }
2260
2261 sub _from_oct
2262   {
2263   # convert a octal number to decimal (string, return ref to array)
2264   my ($c,$os) = @_;
2265
2266   # for older Perls, play safe
2267   my $m = [ 0100000 ];
2268   my $d = 5;                                    # 5 digits at a time
2269
2270   my $mul = _one();
2271   my $x = _zero();
2272
2273   my $len = int( (length($os)-1)/$d );          # $d digit parts, w/o the '0'
2274   my $val; my $i = -$d;
2275   while ($len >= 0)
2276     {
2277     $val = substr($os,$i,$d);                   # get oct digits
2278     $val = CORE::oct($val);
2279     $i -= $d; $len --;
2280     my $adder = [ $val ];
2281     _add ($c, $x, _mul ($c, $adder, $mul ) ) if $val != 0;
2282     _mul ($c, $mul, $m ) if $len >= 0;          # skip last mul
2283     }
2284   $x;
2285   }
2286
2287 sub _from_hex
2288   {
2289   # convert a hex number to decimal (string, return ref to array)
2290   my ($c,$hs) = @_;
2291
2292   my $m = _new($c, 0x10000000);                 # 28 bit at a time (<32 bit!)
2293   my $d = 7;                                    # 7 digits at a time
2294   if ($] <= 5.006)
2295     {
2296     # for older Perls, play safe
2297     $m = [ 0x10000 ];                           # 16 bit at a time (<32 bit!)
2298     $d = 4;                                     # 4 digits at a time
2299     }
2300
2301   my $mul = _one();
2302   my $x = _zero();
2303
2304   my $len = int( (length($hs)-2)/$d );          # $d digit parts, w/o the '0x'
2305   my $val; my $i = -$d;
2306   while ($len >= 0)
2307     {
2308     $val = substr($hs,$i,$d);                   # get hex digits
2309     $val =~ s/^0x// if $len == 0;               # for last part only because
2310     $val = CORE::hex($val);                     # hex does not like wrong chars
2311     $i -= $d; $len --;
2312     my $adder = [ $val ];
2313     # if the resulting number was to big to fit into one element, create a
2314     # two-element version (bug found by Mark Lakata - Thanx!)
2315     if (CORE::length($val) > $BASE_LEN)
2316       {
2317       $adder = _new($c,$val);
2318       }
2319     _add ($c, $x, _mul ($c, $adder, $mul ) ) if $val != 0;
2320     _mul ($c, $mul, $m ) if $len >= 0;          # skip last mul
2321     }
2322   $x;
2323   }
2324
2325 sub _from_bin
2326   {
2327   # convert a hex number to decimal (string, return ref to array)
2328   my ($c,$bs) = @_;
2329
2330   # instead of converting X (8) bit at a time, it is faster to "convert" the
2331   # number to hex, and then call _from_hex.
2332
2333   my $hs = $bs;
2334   $hs =~ s/^[+-]?0b//;                                  # remove sign and 0b
2335   my $l = length($hs);                                  # bits
2336   $hs = '0' x (8-($l % 8)) . $hs if ($l % 8) != 0;      # padd left side w/ 0
2337   my $h = '0x' . unpack('H*', pack ('B*', $hs));        # repack as hex
2338   
2339   $c->_from_hex($h);
2340   }
2341
2342 ##############################################################################
2343 # special modulus functions
2344
2345 sub _modinv
2346   {
2347   # modular inverse
2348   my ($c,$x,$y) = @_;
2349
2350   my $u = _zero($c); my $u1 = _one($c);
2351   my $a = _copy($c,$y); my $b = _copy($c,$x);
2352
2353   # Euclid's Algorithm for bgcd(), only that we calc bgcd() ($a) and the
2354   # result ($u) at the same time. See comments in BigInt for why this works.
2355   my $q;
2356   ($a, $q, $b) = ($b, _div($c,$a,$b));          # step 1
2357   my $sign = 1;
2358   while (!_is_zero($c,$b))
2359     {
2360     my $t = _add($c,                            # step 2:
2361        _mul($c,_copy($c,$u1), $q) ,             #  t =  u1 * q
2362        $u );                                    #     + u
2363     $u = $u1;                                   #  u = u1, u1 = t
2364     $u1 = $t;
2365     $sign = -$sign;
2366     ($a, $q, $b) = ($b, _div($c,$a,$b));        # step 1
2367     }
2368
2369   # if the gcd is not 1, then return NaN
2370   return (undef,undef) unless _is_one($c,$a);
2371  
2372   ($u1, $sign == 1 ? '+' : '-');
2373   }
2374
2375 sub _modpow
2376   {
2377   # modulus of power ($x ** $y) % $z
2378   my ($c,$num,$exp,$mod) = @_;
2379
2380   # a^b (mod 1) = 0 for all a and b
2381   if (_is_one($c,$mod))
2382     {
2383         @$num = 0;
2384         return $num;
2385     }
2386
2387   # 0^a (mod m) = 0 if m != 0, a != 0
2388   # 0^0 (mod m) = 1 if m != 0
2389   if (_is_one($c, $num)) {
2390       if (_is_zero($c, $exp)) {
2391           @$num = 1;
2392       } else {
2393           @$num = 0;
2394       }
2395       return $num;
2396   }
2397
2398 #  $num = _mod($c,$num,$mod);   # this does not make it faster
2399
2400   my $acc = _copy($c,$num); my $t = _one();
2401
2402   my $expbin = _as_bin($c,$exp); $expbin =~ s/^0b//;
2403   my $len = length($expbin);
2404   while (--$len >= 0)
2405     {
2406     if ( substr($expbin,$len,1) eq '1')                 # is_odd
2407       {
2408       _mul($c,$t,$acc);
2409       $t = _mod($c,$t,$mod);
2410       }
2411     _mul($c,$acc,$acc);
2412     $acc = _mod($c,$acc,$mod);
2413     }
2414   @$num = @$t;
2415   $num;
2416   }
2417
2418 sub _gcd
2419   {
2420   # greatest common divisor
2421   my ($c,$x,$y) = @_;
2422
2423   while ( (scalar @$y != 1) || ($y->[0] != 0) )         # while ($y != 0)
2424     {
2425     my $t = _copy($c,$y);
2426     $y = _mod($c, $x, $y);
2427     $x = $t;
2428     }
2429   $x;
2430   }
2431
2432 ##############################################################################
2433 ##############################################################################
2434
2435 1;
2436 __END__
2437
2438 =head1 NAME
2439
2440 Math::BigInt::Calc - Pure Perl module to support Math::BigInt
2441
2442 =head1 SYNOPSIS
2443
2444 Provides support for big integer calculations. Not intended to be used by other
2445 modules. Other modules which sport the same functions can also be used to support
2446 Math::BigInt, like Math::BigInt::GMP or Math::BigInt::Pari.
2447
2448 =head1 DESCRIPTION
2449
2450 In order to allow for multiple big integer libraries, Math::BigInt was
2451 rewritten to use library modules for core math routines. Any module which
2452 follows the same API as this can be used instead by using the following:
2453
2454         use Math::BigInt lib => 'libname';
2455
2456 'libname' is either the long name ('Math::BigInt::Pari'), or only the short
2457 version like 'Pari'.
2458
2459 =head1 STORAGE
2460
2461 =head1 METHODS
2462
2463 The following functions MUST be defined in order to support the use by
2464 Math::BigInt v1.70 or later:
2465
2466         api_version()   return API version, 1 for v1.70, 2 for v1.83
2467         _new(string)    return ref to new object from ref to decimal string
2468         _zero()         return a new object with value 0
2469         _one()          return a new object with value 1
2470         _two()          return a new object with value 2
2471         _ten()          return a new object with value 10
2472
2473         _str(obj)       return ref to a string representing the object
2474         _num(obj)       returns a Perl integer/floating point number
2475                         NOTE: because of Perl numeric notation defaults,
2476                         the _num'ified obj may lose accuracy due to 
2477                         machine-dependent floating point size limitations
2478                     
2479         _add(obj,obj)   Simple addition of two objects
2480         _mul(obj,obj)   Multiplication of two objects
2481         _div(obj,obj)   Division of the 1st object by the 2nd
2482                         In list context, returns (result,remainder).
2483                         NOTE: this is integer math, so no
2484                         fractional part will be returned.
2485                         The second operand will be not be 0, so no need to
2486                         check for that.
2487         _sub(obj,obj)   Simple subtraction of 1 object from another
2488                         a third, optional parameter indicates that the params
2489                         are swapped. In this case, the first param needs to
2490                         be preserved, while you can destroy the second.
2491                         sub (x,y,1) => return x - y and keep x intact!
2492         _dec(obj)       decrement object by one (input is guaranteed to be > 0)
2493         _inc(obj)       increment object by one
2494
2495
2496         _acmp(obj,obj)  <=> operator for objects (return -1, 0 or 1)
2497
2498         _len(obj)       returns count of the decimal digits of the object
2499         _digit(obj,n)   returns the n'th decimal digit of object
2500
2501         _is_one(obj)    return true if argument is 1
2502         _is_two(obj)    return true if argument is 2
2503         _is_ten(obj)    return true if argument is 10
2504         _is_zero(obj)   return true if argument is 0
2505         _is_even(obj)   return true if argument is even (0,2,4,6..)
2506         _is_odd(obj)    return true if argument is odd (1,3,5,7..)
2507
2508         _copy           return a ref to a true copy of the object
2509
2510         _check(obj)     check whether internal representation is still intact
2511                         return 0 for ok, otherwise error message as string
2512
2513         _from_hex(str)  return new object from a hexadecimal string
2514         _from_bin(str)  return new object from a binary string
2515         _from_oct(str)  return new object from an octal string
2516         
2517         _as_hex(str)    return string containing the value as
2518                         unsigned hex string, with the '0x' prepended.
2519                         Leading zeros must be stripped.
2520         _as_bin(str)    Like as_hex, only as binary string containing only
2521                         zeros and ones. Leading zeros must be stripped and a
2522                         '0b' must be prepended.
2523         
2524         _rsft(obj,N,B)  shift object in base B by N 'digits' right
2525         _lsft(obj,N,B)  shift object in base B by N 'digits' left
2526         
2527         _xor(obj1,obj2) XOR (bit-wise) object 1 with object 2
2528                         Note: XOR, AND and OR pad with zeros if size mismatches
2529         _and(obj1,obj2) AND (bit-wise) object 1 with object 2
2530         _or(obj1,obj2)  OR (bit-wise) object 1 with object 2
2531
2532         _mod(obj1,obj2) Return remainder of div of the 1st by the 2nd object
2533         _sqrt(obj)      return the square root of object (truncated to int)
2534         _root(obj)      return the n'th (n >= 3) root of obj (truncated to int)
2535         _fac(obj)       return factorial of object 1 (1*2*3*4..)
2536         _pow(obj1,obj2) return object 1 to the power of object 2
2537                         return undef for NaN
2538         _zeros(obj)     return number of trailing decimal zeros
2539         _modinv         return inverse modulus
2540         _modpow         return modulus of power ($x ** $y) % $z
2541         _log_int(X,N)   calculate integer log() of X in base N
2542                         X >= 0, N >= 0 (return undef for NaN)
2543                         returns (RESULT, EXACT) where EXACT is:
2544                          1     : result is exactly RESULT
2545                          0     : result was truncated to RESULT
2546                          undef : unknown whether result is exactly RESULT
2547         _gcd(obj,obj)   return Greatest Common Divisor of two objects
2548
2549 The following functions are REQUIRED for an api_version of 2 or greater:
2550
2551         _1ex($x)        create the number 1Ex where x >= 0
2552         _alen(obj)      returns approximate count of the decimal digits of the
2553                         object. This estimate MUST always be greater or equal
2554                         to what _len() returns.
2555         _nok(n,k)       calculate n over k (binomial coefficient)
2556
2557 The following functions are optional, and can be defined if the underlying lib
2558 has a fast way to do them. If undefined, Math::BigInt will use pure Perl (hence
2559 slow) fallback routines to emulate these:
2560         
2561         _signed_or
2562         _signed_and
2563         _signed_xor
2564
2565 Input strings come in as unsigned but with prefix (i.e. as '123', '0xabc'
2566 or '0b1101').
2567
2568 So the library needs only to deal with unsigned big integers. Testing of input
2569 parameter validity is done by the caller, so you need not worry about
2570 underflow (f.i. in C<_sub()>, C<_dec()>) nor about division by zero or similar
2571 cases.
2572
2573 The first parameter can be modified, that includes the possibility that you
2574 return a reference to a completely different object instead. Although keeping
2575 the reference and just changing its contents is preferred over creating and
2576 returning a different reference.
2577
2578 Return values are always references to objects, strings, or true/false for
2579 comparison routines.
2580
2581 =head1 WRAP YOUR OWN
2582
2583 If you want to port your own favourite c-lib for big numbers to the
2584 Math::BigInt interface, you can take any of the already existing modules as
2585 a rough guideline. You should really wrap up the latest BigInt and BigFloat
2586 testsuites with your module, and replace in them any of the following:
2587
2588         use Math::BigInt;
2589
2590 by this:
2591
2592         use Math::BigInt lib => 'yourlib';
2593
2594 This way you ensure that your library really works 100% within Math::BigInt.
2595
2596 =head1 LICENSE
2597  
2598 This program is free software; you may redistribute it and/or modify it under
2599 the same terms as Perl itself. 
2600
2601 =head1 AUTHORS
2602
2603 Original math code by Mark Biggar, rewritten by Tels L<http://bloodgate.com/>
2604 in late 2000.
2605 Separated from BigInt and shaped API with the help of John Peacock.
2606
2607 Fixed, speed-up, streamlined and enhanced by Tels 2001 - 2007.
2608
2609 =head1 SEE ALSO
2610
2611 L<Math::BigInt>, L<Math::BigFloat>,
2612 L<Math::BigInt::GMP>, L<Math::BigInt::FastCalc> and L<Math::BigInt::Pari>.
2613
2614 =cut