This is a live mirror of the Perl 5 development currently hosted at https://github.com/perl/perl5
Upgrade to Math::BigInt pre-rel 1.66 as of
[perl5.git] / lib / Math / BigInt / Calc.pm
1 package Math::BigInt::Calc;
2
3 use 5.005;
4 use strict;
5 # use warnings; # dont use warnings for older Perls
6
7 require Exporter;
8 use vars qw/@ISA $VERSION/;
9 @ISA = qw(Exporter);
10
11 $VERSION = '0.36';
12
13 # Package to store unsigned big integers in decimal and do math with them
14
15 # Internally the numbers are stored in an array with at least 1 element, no
16 # leading zero parts (except the first) and in base 1eX where X is determined
17 # automatically at loading time to be the maximum possible value
18
19 # todo:
20 # - fully remove funky $# stuff (maybe)
21
22 # USE_MUL: due to problems on certain os (os390, posix-bc) "* 1e-5" is used
23 # instead of "/ 1e5" at some places, (marked with USE_MUL). Other platforms
24 # BS2000, some Crays need USE_DIV instead.
25 # The BEGIN block is used to determine which of the two variants gives the
26 # correct result.
27
28 # Beware of things like:
29 # $i = $i * $y + $car; $car = int($i / $MBASE); $i = $i % $MBASE;
30 # This works on x86, but fails on ARM (SA1100, iPAQ) due to whoknows what
31 # reasons. So, use this instead (slower, but correct):
32 # $i = $i * $y + $car; $car = int($i / $MBASE); $i -= $MBASE * $car;
33
34 ##############################################################################
35 # global constants, flags and accessory
36  
37 # constants for easier life
38 my $nan = 'NaN';
39 my ($MBASE,$BASE,$RBASE,$BASE_LEN,$MAX_VAL,$BASE_LEN2,$BASE_LEN_SMALL);
40 my ($AND_BITS,$XOR_BITS,$OR_BITS);
41 my ($AND_MASK,$XOR_MASK,$OR_MASK);
42
43 sub _base_len 
44   {
45   # set/get the BASE_LEN and assorted other, connected values
46   # used only be the testsuite, set is used only by the BEGIN block below
47   shift;
48
49   my $b = shift;
50   if (defined $b)
51     {
52     # find whether we can use mul or div or none in mul()/div()
53     # (in last case reduce BASE_LEN_SMALL)
54     $BASE_LEN_SMALL = $b+1;
55     my $caught = 0;
56     while (--$BASE_LEN_SMALL > 5)
57       {
58       $MBASE = int("1e".$BASE_LEN_SMALL);
59       $RBASE = abs('1e-'.$BASE_LEN_SMALL);              # see USE_MUL
60       $caught = 0;
61       $caught += 1 if (int($MBASE * $RBASE) != 1);      # should be 1
62       $caught += 2 if (int($MBASE / $MBASE) != 1);      # should be 1
63       last if $caught != 3;
64       }
65     # BASE_LEN is used for anything else than mul()/div()
66     $BASE_LEN = $BASE_LEN_SMALL;
67     $BASE_LEN = shift if (defined $_[0]);               # one more arg?
68     $BASE = int("1e".$BASE_LEN);
69
70     $BASE_LEN2 = int($BASE_LEN_SMALL / 2);              # for mul shortcut
71     $MBASE = int("1e".$BASE_LEN_SMALL);
72     $RBASE = abs('1e-'.$BASE_LEN_SMALL);                # see USE_MUL
73     $MAX_VAL = $MBASE-1;
74     
75     #print "BASE_LEN: $BASE_LEN MAX_VAL: $MAX_VAL BASE: $BASE RBASE: $RBASE ";
76     #print "BASE_LEN_SMALL: $BASE_LEN_SMALL MBASE: $MBASE\n";
77
78     undef &_mul;
79     undef &_div;
80
81     # $caught & 1 != 0 => cannot use MUL
82     # $caught & 2 != 0 => cannot use DIV
83     # The parens around ($caught & 1) were important, indeed, if we would use
84     # & here.
85     if ($caught == 2)                           # 2
86       {
87       # print "# use mul\n";
88       # must USE_MUL since we cannot use DIV
89       *{_mul} = \&_mul_use_mul;
90       *{_div} = \&_div_use_mul;
91       }
92     else                                        # 0 or 1
93       {
94       # print "# use div\n";
95       # can USE_DIV instead
96       *{_mul} = \&_mul_use_div;
97       *{_div} = \&_div_use_div;
98       }
99     }
100   return $BASE_LEN unless wantarray;
101   return ($BASE_LEN, $AND_BITS, $XOR_BITS, $OR_BITS, $BASE_LEN_SMALL, $MAX_VAL);
102   }
103
104 BEGIN
105   {
106   # from Daniel Pfeiffer: determine largest group of digits that is precisely
107   # multipliable with itself plus carry
108   # Test now changed to expect the proper pattern, not a result off by 1 or 2
109   my ($e, $num) = 3;    # lowest value we will use is 3+1-1 = 3
110   do 
111     {
112     $num = ('9' x ++$e) + 0;
113     $num *= $num + 1.0;
114     } while ("$num" =~ /9{$e}0{$e}/);   # must be a certain pattern
115   $e--;                                 # last test failed, so retract one step
116   # the limits below brush the problems with the test above under the rug:
117   # the test should be able to find the proper $e automatically
118   $e = 5 if $^O =~ /^uts/;      # UTS get's some special treatment
119   $e = 5 if $^O =~ /^unicos/;   # unicos is also problematic (6 seems to work
120                                 # there, but we play safe)
121   $e = 5 if $] < 5.006;         # cap, for older Perls
122   $e = 7 if $e > 7;             # cap, for VMS, OS/390 and other 64 bit systems
123                                 # 8 fails inside random testsuite, so take 7
124
125   # determine how many digits fit into an integer and can be safely added 
126   # together plus carry w/o causing an overflow
127
128   # this below detects 15 on a 64 bit system, because after that it becomes
129   # 1e16  and not 1000000 :/ I can make it detect 18, but then I get a lot of
130   # test failures. Ugh! (Tomake detect 18: uncomment lines marked with *)
131   use integer;
132   my $bi = 5;                   # approx. 16 bit
133   $num = int('9' x $bi);
134   # $num = 99999; # *
135   # while ( ($num+$num+1) eq '1' . '9' x $bi)   # *
136   while ( int($num+$num+1) eq '1' . '9' x $bi)
137     {
138     $bi++; $num = int('9' x $bi);
139     # $bi++; $num *= 10; $num += 9;     # *
140     }
141   $bi--;                                # back off one step
142   # by setting them equal, we ignore the findings and use the default
143   # one-size-fits-all approach from former versions
144   $bi = $e;                             # XXX, this should work always
145
146   __PACKAGE__->_base_len($e,$bi);       # set and store
147
148   # find out how many bits _and, _or and _xor can take (old default = 16)
149   # I don't think anybody has yet 128 bit scalars, so let's play safe.
150   local $^W = 0;        # don't warn about 'nonportable number'
151   $AND_BITS = 15; $XOR_BITS = 15; $OR_BITS  = 15;
152
153   # find max bits, we will not go higher than numberofbits that fit into $BASE
154   # to make _and etc simpler (and faster for smaller, slower for large numbers)
155   my $max = 16;
156   while (2 ** $max < $BASE) { $max++; }
157   {
158     no integer;
159     $max = 16 if $] < 5.006;    # older Perls might not take >16 too well
160   }
161   my ($x,$y,$z);
162   do {
163     $AND_BITS++;
164     $x = oct('0b' . '1' x $AND_BITS); $y = $x & $x;
165     $z = (2 ** $AND_BITS) - 1;
166     } while ($AND_BITS < $max && $x == $z && $y == $x);
167   $AND_BITS --;                                         # retreat one step
168   do {
169     $XOR_BITS++;
170     $x = oct('0b' . '1' x $XOR_BITS); $y = $x ^ 0;
171     $z = (2 ** $XOR_BITS) - 1;
172     } while ($XOR_BITS < $max && $x == $z && $y == $x);
173   $XOR_BITS --;                                         # retreat one step
174   do {
175     $OR_BITS++;
176     $x = oct('0b' . '1' x $OR_BITS); $y = $x | $x;
177     $z = (2 ** $OR_BITS) - 1;
178     } while ($OR_BITS < $max && $x == $z && $y == $x);
179   $OR_BITS --;                                          # retreat one step
180   
181   }
182
183 ###############################################################################
184
185 sub _new
186   {
187   # (ref to string) return ref to num_array
188   # Convert a number from string format (without sign) to internal base
189   # 1ex format. Assumes normalized value as input.
190   my $d = $_[1];
191   my $il = length($$d)-1;
192   # this leaves '00000' instead of int 0 and will be corrected after any op
193   [ reverse(unpack("a" . ($il % $BASE_LEN+1) 
194     . ("a$BASE_LEN" x ($il / $BASE_LEN)), $$d)) ];
195   }                                                                             
196   
197 BEGIN
198   {
199   $AND_MASK = __PACKAGE__->_new( \( 2 ** $AND_BITS ));
200   $XOR_MASK = __PACKAGE__->_new( \( 2 ** $XOR_BITS ));
201   $OR_MASK = __PACKAGE__->_new( \( 2 ** $OR_BITS ));
202   }
203
204 sub _zero
205   {
206   # create a zero
207   [ 0 ];
208   }
209
210 sub _one
211   {
212   # create a one
213   [ 1 ];
214   }
215
216 sub _two
217   {
218   # create a two (used internally for shifting)
219   [ 2 ];
220   }
221
222 sub _copy
223   {
224   [ @{$_[1]} ];
225   }
226
227 # catch and throw away
228 sub import { }
229
230 ##############################################################################
231 # convert back to string and number
232
233 sub _str
234   {
235   # (ref to BINT) return num_str
236   # Convert number from internal base 100000 format to string format.
237   # internal format is always normalized (no leading zeros, "-0" => "+0")
238   my $ar = $_[1];
239   my $ret = "";
240
241   my $l = scalar @$ar;          # number of parts
242   return $nan if $l < 1;        # should not happen
243
244   # handle first one different to strip leading zeros from it (there are no
245   # leading zero parts in internal representation)
246   $l --; $ret .= int($ar->[$l]); $l--;
247   # Interestingly, the pre-padd method uses more time
248   # the old grep variant takes longer (14 to 10 sec)
249   my $z = '0' x ($BASE_LEN-1);                            
250   while ($l >= 0)
251     {
252     $ret .= substr($z.$ar->[$l],-$BASE_LEN); # fastest way I could think of
253     $l--;
254     }
255   \$ret;
256   }                                                                             
257
258 sub _num
259   {
260   # Make a number (scalar int/float) from a BigInt object
261   my $x = $_[1];
262   return $x->[0] if scalar @$x == 1;  # below $BASE
263   my $fac = 1;
264   my $num = 0;
265   foreach (@$x)
266     {
267     $num += $fac*$_; $fac *= $BASE;
268     }
269   $num; 
270   }
271
272 ##############################################################################
273 # actual math code
274
275 sub _add
276   {
277   # (ref to int_num_array, ref to int_num_array)
278   # routine to add two base 1eX numbers
279   # stolen from Knuth Vol 2 Algorithm A pg 231
280   # there are separate routines to add and sub as per Knuth pg 233
281   # This routine clobbers up array x, but not y.
282  
283   my ($c,$x,$y) = @_;
284
285   return $x if (@$y == 1) && $y->[0] == 0;              # $x + 0 => $x
286   if ((@$x == 1) && $x->[0] == 0)                       # 0 + $y => $y->copy
287     {
288     # twice as slow as $x = [ @$y ], but necc. to retain $x as ref :(
289     @$x = @$y; return $x;               
290     }
291  
292   # for each in Y, add Y to X and carry. If after that, something is left in
293   # X, foreach in X add carry to X and then return X, carry
294   # Trades one "$j++" for having to shift arrays, $j could be made integer
295   # but this would impose a limit to number-length of 2**32.
296   my $i; my $car = 0; my $j = 0;
297   for $i (@$y)
298     {
299     $x->[$j] -= $BASE if $car = (($x->[$j] += $i + $car) >= $BASE) ? 1 : 0;
300     $j++;
301     }
302   while ($car != 0)
303     {
304     $x->[$j] -= $BASE if $car = (($x->[$j] += $car) >= $BASE) ? 1 : 0; $j++;
305     }
306   $x;
307   }                                                                             
308
309 sub _inc
310   {
311   # (ref to int_num_array, ref to int_num_array)
312   # routine to add 1 to a base 1eX numbers
313   # This routine modifies array x
314   my ($c,$x) = @_;
315
316   for my $i (@$x)
317     {
318     return $x if (($i += 1) < $BASE);           # early out
319     $i = 0;                                     # overflow, next
320     }
321   push @$x,1 if ($x->[-1] == 0);                # last overflowed, so extend
322   $x;
323   }                                                                             
324
325 sub _dec
326   {
327   # (ref to int_num_array, ref to int_num_array)
328   # routine to add 1 to a base 1eX numbers
329   # This routine modifies array x
330   my ($c,$x) = @_;
331
332   my $MAX = $BASE-1;                            # since MAX_VAL based on MBASE
333   for my $i (@$x)
334     {
335     last if (($i -= 1) >= 0);                   # early out
336     $i = $MAX;                                  # overflow, next
337     }
338   pop @$x if $x->[-1] == 0 && @$x > 1;          # last overflowed (but leave 0)
339   $x;
340   }                                                                             
341
342 sub _sub
343   {
344   # (ref to int_num_array, ref to int_num_array, swap)
345   # subtract base 1eX numbers -- stolen from Knuth Vol 2 pg 232, $x > $y
346   # subtract Y from X by modifying x in place
347   my ($c,$sx,$sy,$s) = @_;
348  
349   my $car = 0; my $i; my $j = 0;
350   if (!$s)
351     {
352     #print "case 2\n";
353     for $i (@$sx)
354       {
355       last unless defined $sy->[$j] || $car;
356       $i += $BASE if $car = (($i -= ($sy->[$j] || 0) + $car) < 0); $j++;
357       }
358     # might leave leading zeros, so fix that
359     return __strip_zeros($sx);
360     }
361   #print "case 1 (swap)\n";
362   for $i (@$sx)
363     {
364     # we can't do an early out if $x is < than $y, since we
365     # need to copy the high chunks from $y. Found by Bob Mathews.
366     #last unless defined $sy->[$j] || $car;
367     $sy->[$j] += $BASE
368      if $car = (($sy->[$j] = $i-($sy->[$j]||0) - $car) < 0);
369     $j++;
370     }
371   # might leave leading zeros, so fix that
372   __strip_zeros($sy);
373   }                                                                             
374
375 sub _mul_use_mul
376   {
377   # (ref to int_num_array, ref to int_num_array)
378   # multiply two numbers in internal representation
379   # modifies first arg, second need not be different from first
380   my ($c,$xv,$yv) = @_;
381
382   if (@$yv == 1)
383     {
384     # shortcut for two very short numbers (improved by Nathan Zook)
385     # works also if xv and yv are the same reference, and handles also $x == 0
386     if (@$xv == 1)
387       {
388       if (($xv->[0] *= $yv->[0]) >= $MBASE)
389          {
390          $xv->[0] = $xv->[0] - ($xv->[1] = int($xv->[0] * $RBASE)) * $MBASE;
391          };
392       return $xv;
393       }
394     # $x * 0 => 0
395     if ($yv->[0] == 0)
396       {
397       @$xv = (0);
398       return $xv;
399       }
400     # multiply a large number a by a single element one, so speed up
401     my $y = $yv->[0]; my $car = 0;
402     foreach my $i (@$xv)
403       {
404       $i = $i * $y + $car; $car = int($i * $RBASE); $i -= $car * $MBASE;
405       }
406     push @$xv, $car if $car != 0;
407     return $xv;
408     }
409   # shortcut for result $x == 0 => result = 0
410   return $xv if ( ((@$xv == 1) && ($xv->[0] == 0)) ); 
411
412   # since multiplying $x with $x fails, make copy in this case
413   $yv = [@$xv] if $xv == $yv;   # same references?
414
415   my @prod = (); my ($prod,$car,$cty,$xi,$yi);
416
417   for $xi (@$xv)
418     {
419     $car = 0; $cty = 0;
420
421     # slow variant
422 #    for $yi (@$yv)
423 #      {
424 #      $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
425 #      $prod[$cty++] =
426 #       $prod - ($car = int($prod * RBASE)) * $MBASE;  # see USE_MUL
427 #      }
428 #    $prod[$cty] += $car if $car; # need really to check for 0?
429 #    $xi = shift @prod;
430
431     # faster variant
432     # looping through this if $xi == 0 is silly - so optimize it away!
433     $xi = (shift @prod || 0), next if $xi == 0;
434     for $yi (@$yv)
435       {
436       $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
437 ##     this is actually a tad slower
438 ##        $prod = $prod[$cty]; $prod += ($car + $xi * $yi);     # no ||0 here
439       $prod[$cty++] =
440        $prod - ($car = int($prod * $RBASE)) * $MBASE;  # see USE_MUL
441       }
442     $prod[$cty] += $car if $car; # need really to check for 0?
443     $xi = shift @prod || 0;     # || 0 makes v5.005_3 happy
444     }
445   push @$xv, @prod;
446   __strip_zeros($xv);
447   $xv;
448   }                                                                             
449
450 sub _mul_use_div
451   {
452   # (ref to int_num_array, ref to int_num_array)
453   # multiply two numbers in internal representation
454   # modifies first arg, second need not be different from first
455   my ($c,$xv,$yv) = @_;
456  
457   if (@$yv == 1)
458     {
459     # shortcut for two small numbers, also handles $x == 0
460     if (@$xv == 1)
461       {
462       # shortcut for two very short numbers (improved by Nathan Zook)
463       # works also if xv and yv are the same reference, and handles also $x == 0
464       if (($xv->[0] *= $yv->[0]) >= $MBASE)
465           {
466           $xv->[0] =
467               $xv->[0] - ($xv->[1] = int($xv->[0] / $MBASE)) * $MBASE;
468           };
469       return $xv;
470       }
471     # $x * 0 => 0
472     if ($yv->[0] == 0)
473       {
474       @$xv = (0);
475       return $xv;
476       }
477     # multiply a large number a by a single element one, so speed up
478     my $y = $yv->[0]; my $car = 0;
479     foreach my $i (@$xv)
480       {
481       $i = $i * $y + $car; $car = int($i / $MBASE); $i -= $car * $MBASE;
482       }
483     push @$xv, $car if $car != 0;
484     return $xv;
485     }
486   # shortcut for result $x == 0 => result = 0
487   return $xv if ( ((@$xv == 1) && ($xv->[0] == 0)) ); 
488
489   # since multiplying $x with $x fails, make copy in this case
490   $yv = [@$xv] if $xv == $yv;   # same references?
491
492   my @prod = (); my ($prod,$car,$cty,$xi,$yi);
493   for $xi (@$xv)
494     {
495     $car = 0; $cty = 0;
496     # looping through this if $xi == 0 is silly - so optimize it away!
497     $xi = (shift @prod || 0), next if $xi == 0;
498     for $yi (@$yv)
499       {
500       $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
501       $prod[$cty++] =
502        $prod - ($car = int($prod / $MBASE)) * $MBASE;
503       }
504     $prod[$cty] += $car if $car; # need really to check for 0?
505     $xi = shift @prod || 0;     # || 0 makes v5.005_3 happy
506     }
507   push @$xv, @prod;
508   __strip_zeros($xv);
509   $xv;
510   }                                                                             
511
512 sub _div_use_mul
513   {
514   # ref to array, ref to array, modify first array and return remainder if 
515   # in list context
516
517   # see comments in _div_use_div() for more explanations
518
519   my ($c,$x,$yorg) = @_;
520   
521   # the general div algorithmn here is about O(N*N) and thus quite slow, so
522   # we first check for some special cases and use shortcuts to handle them.
523
524   # This works, because we store the numbers in a chunked format where each
525   # element contains 5..7 digits (depending on system).
526
527   # if both numbers have only one element:
528   if (@$x == 1 && @$yorg == 1)
529     {
530     # shortcut, $yorg and $x are two small numbers
531     if (wantarray)
532       {
533       my $r = [ $x->[0] % $yorg->[0] ];
534       $x->[0] = int($x->[0] / $yorg->[0]);
535       return ($x,$r); 
536       }
537     else
538       {
539       $x->[0] = int($x->[0] / $yorg->[0]);
540       return $x; 
541       }
542     }
543
544   # if x has more than one, but y has only one element:
545   if (@$yorg == 1)
546     {
547     my $rem;
548     $rem = _mod($c,[ @$x ],$yorg) if wantarray;
549
550     # shortcut, $y is < $BASE
551     my $j = scalar @$x; my $r = 0; 
552     my $y = $yorg->[0]; my $b;
553     while ($j-- > 0)
554       {
555       $b = $r * $MBASE + $x->[$j];
556       $x->[$j] = int($b/$y);
557       $r = $b % $y;
558       }
559     pop @$x if @$x > 1 && $x->[-1] == 0;        # splice up a leading zero 
560     return ($x,$rem) if wantarray;
561     return $x;
562     }
563
564   # now x and y have more than one element
565
566   # check whether y has more elements than x, if yet, the result will be 0
567   if (@$yorg > @$x)
568     {
569     my $rem;
570     $rem = [@$x] if wantarray;                  # make copy
571     splice (@$x,1);                             # keep ref to original array
572     $x->[0] = 0;                                # set to 0
573     return ($x,$rem) if wantarray;              # including remainder?
574     return $x;                                  # only x, which is [0] now
575     }
576   # check whether the numbers have the same number of elements, in that case
577   # the result will fit into one element and can be computed efficiently
578   if (@$yorg == @$x)
579     {
580     my $rem;
581     # if $yorg has more digits than $x (it's leading element is longer than
582     # the one from $x), the result will also be 0:
583     if (length(int($yorg->[-1])) > length(int($x->[-1])))
584       {
585       $rem = [@$x] if wantarray;                # make copy
586       splice (@$x,1);                           # keep ref to org array
587       $x->[0] = 0;                              # set to 0
588       return ($x,$rem) if wantarray;            # including remainder?
589       return $x;
590       }
591     # now calculate $x / $yorg
592     if (length(int($yorg->[-1])) == length(int($x->[-1])))
593       {
594       # same length, so make full compare, and if equal, return 1
595       # hm, same lengths, but same contents? So we need to check all parts:
596       my $a = 0; my $j = scalar @$x - 1;
597       # manual way (abort if unequal, good for early ne)
598       while ($j >= 0)
599         {
600         last if ($a = $x->[$j] - $yorg->[$j]); $j--;
601         }
602       # $a contains the result of the compare between X and Y
603       # a < 0: x < y, a == 0 => x == y, a > 0: x > y
604       if ($a <= 0)
605         {
606         if (wantarray)
607           {
608           $rem = [ 0 ];                 # a = 0 => x == y => rem 1
609           $rem = [@$x] if $a != 0;      # a < 0 => x < y => rem = x
610           }
611         splice(@$x,1);                  # keep single element
612         $x->[0] = 0;                    # if $a < 0
613         if ($a == 0)
614           {
615           # $x == $y
616           $x->[0] = 1;
617           }
618         return ($x,$rem) if wantarray;
619         return $x;
620         }
621       # $x >= $y, proceed normally
622       }
623     }
624
625   # all other cases:
626
627   my $y = [ @$yorg ];                           # always make copy to preserve
628
629   my ($car,$bar,$prd,$dd,$xi,$yi,@q,$v2,$v1,@d,$tmp,$q,$u2,$u1,$u0);
630
631   $car = $bar = $prd = 0;
632   if (($dd = int($MBASE/($y->[-1]+1))) != 1) 
633     {
634     for $xi (@$x) 
635       {
636       $xi = $xi * $dd + $car;
637       $xi -= ($car = int($xi * $RBASE)) * $MBASE;       # see USE_MUL
638       }
639     push(@$x, $car); $car = 0;
640     for $yi (@$y) 
641       {
642       $yi = $yi * $dd + $car;
643       $yi -= ($car = int($yi * $RBASE)) * $MBASE;       # see USE_MUL
644       }
645     }
646   else 
647     {
648     push(@$x, 0);
649     }
650   @q = (); ($v2,$v1) = @$y[-2,-1];
651   $v2 = 0 unless $v2;
652   while ($#$x > $#$y) 
653     {
654     ($u2,$u1,$u0) = @$x[-3..-1];
655     $u2 = 0 unless $u2;
656     #warn "oups v1 is 0, u0: $u0 $y->[-2] $y->[-1] l ",scalar @$y,"\n"
657     # if $v1 == 0;
658     $q = (($u0 == $v1) ? $MAX_VAL : int(($u0*$MBASE+$u1)/$v1));
659     --$q while ($v2*$q > ($u0*$MBASE+$u1-$q*$v1)*$MBASE+$u2);
660     if ($q)
661       {
662       ($car, $bar) = (0,0);
663       for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi) 
664         {
665         $prd = $q * $y->[$yi] + $car;
666         $prd -= ($car = int($prd * $RBASE)) * $MBASE;   # see USE_MUL
667         $x->[$xi] += $MBASE if ($bar = (($x->[$xi] -= $prd + $bar) < 0));
668         }
669       if ($x->[-1] < $car + $bar) 
670         {
671         $car = 0; --$q;
672         for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi) 
673           {
674           $x->[$xi] -= $MBASE
675            if ($car = (($x->[$xi] += $y->[$yi] + $car) >= $MBASE));
676           }
677         }   
678       }
679     pop(@$x);
680     unshift(@q, $q);
681     }
682   if (wantarray) 
683     {
684     @d = ();
685     if ($dd != 1)  
686       {
687       $car = 0; 
688       for $xi (reverse @$x) 
689         {
690         $prd = $car * $MBASE + $xi;
691         $car = $prd - ($tmp = int($prd / $dd)) * $dd; # see USE_MUL
692         unshift(@d, $tmp);
693         }
694       }
695     else 
696       {
697       @d = @$x;
698       }
699     @$x = @q;
700     my $d = \@d; 
701     __strip_zeros($x);
702     __strip_zeros($d);
703     return ($x,$d);
704     }
705   @$x = @q;
706   __strip_zeros($x);
707   $x;
708   }
709
710 sub _div_use_div
711   {
712   # ref to array, ref to array, modify first array and return remainder if 
713   # in list context
714   my ($c,$x,$yorg) = @_;
715
716   # the general div algorithmn here is about O(N*N) and thus quite slow, so
717   # we first check for some special cases and use shortcuts to handle them.
718
719   # This works, because we store the numbers in a chunked format where each
720   # element contains 5..7 digits (depending on system).
721
722   # if both numbers have only one element:
723   if (@$x == 1 && @$yorg == 1)
724     {
725     # shortcut, $yorg and $x are two small numbers
726     if (wantarray)
727       {
728       my $r = [ $x->[0] % $yorg->[0] ];
729       $x->[0] = int($x->[0] / $yorg->[0]);
730       return ($x,$r); 
731       }
732     else
733       {
734       $x->[0] = int($x->[0] / $yorg->[0]);
735       return $x; 
736       }
737     }
738   # if x has more than one, but y has only one element:
739   if (@$yorg == 1)
740     {
741     my $rem;
742     $rem = _mod($c,[ @$x ],$yorg) if wantarray;
743
744     # shortcut, $y is < $BASE
745     my $j = scalar @$x; my $r = 0; 
746     my $y = $yorg->[0]; my $b;
747     while ($j-- > 0)
748       {
749       $b = $r * $MBASE + $x->[$j];
750       $x->[$j] = int($b/$y);
751       $r = $b % $y;
752       }
753     pop @$x if @$x > 1 && $x->[-1] == 0;        # splice up a leading zero 
754     return ($x,$rem) if wantarray;
755     return $x;
756     }
757   # now x and y have more than one element
758
759   # check whether y has more elements than x, if yet, the result will be 0
760   if (@$yorg > @$x)
761     {
762     my $rem;
763     $rem = [@$x] if wantarray;                  # make copy
764     splice (@$x,1);                             # keep ref to original array
765     $x->[0] = 0;                                # set to 0
766     return ($x,$rem) if wantarray;              # including remainder?
767     return $x;                                  # only x, which is [0] now
768     }
769   # check whether the numbers have the same number of elements, in that case
770   # the result will fit into one element and can be computed efficiently
771   if (@$yorg == @$x)
772     {
773     my $rem;
774     # if $yorg has more digits than $x (it's leading element is longer than
775     # the one from $x), the result will also be 0:
776     if (length(int($yorg->[-1])) > length(int($x->[-1])))
777       {
778       $rem = [@$x] if wantarray;                # make copy
779       splice (@$x,1);                           # keep ref to org array
780       $x->[0] = 0;                              # set to 0
781       return ($x,$rem) if wantarray;            # including remainder?
782       return $x;
783       }
784     # now calculate $x / $yorg
785     if (length(int($yorg->[-1])) == length(int($x->[-1])))
786       {
787       # same length, so make full compare, and if equal, return 1
788       # hm, same lengths, but same contents? So we need to check all parts:
789       my $a = 0; my $j = scalar @$x - 1;
790       # manual way (abort if unequal, good for early ne)
791       while ($j >= 0)
792         {
793         last if ($a = $x->[$j] - $yorg->[$j]); $j--;
794         }
795       # $a contains the result of the compare between X and Y
796       # a < 0: x < y, a == 0 => x == y, a > 0: x > y
797       if ($a <= 0)
798         {
799         if (wantarray)
800           {
801           $rem = [ 0 ];                 # a = 0 => x == y => rem 1
802           $rem = [@$x] if $a != 0;      # a < 0 => x < y => rem = x
803           }
804         splice(@$x,1);                  # keep single element
805         $x->[0] = 0;                    # if $a < 0
806         if ($a == 0)
807           {
808           # $x == $y
809           $x->[0] = 1;
810           }
811         return ($x,$rem) if wantarray;
812         return $x;
813         }
814       # $x >= $y, so proceed normally
815       }
816     }
817
818   # all other cases:
819
820   my $y = [ @$yorg ];                           # always make copy to preserve
821  
822   my ($car,$bar,$prd,$dd,$xi,$yi,@q,$v2,$v1,@d,$tmp,$q,$u2,$u1,$u0);
823
824   $car = $bar = $prd = 0;
825   if (($dd = int($MBASE/($y->[-1]+1))) != 1) 
826     {
827     for $xi (@$x) 
828       {
829       $xi = $xi * $dd + $car;
830       $xi -= ($car = int($xi / $MBASE)) * $MBASE;
831       }
832     push(@$x, $car); $car = 0;
833     for $yi (@$y) 
834       {
835       $yi = $yi * $dd + $car;
836       $yi -= ($car = int($yi / $MBASE)) * $MBASE;
837       }
838     }
839   else 
840     {
841     push(@$x, 0);
842     }
843
844   # @q will accumulate the final result, $q contains the current computed
845   # part of the final result
846
847   @q = (); ($v2,$v1) = @$y[-2,-1];
848   $v2 = 0 unless $v2;
849   while ($#$x > $#$y) 
850     {
851     ($u2,$u1,$u0) = @$x[-3..-1];
852     $u2 = 0 unless $u2;
853     #warn "oups v1 is 0, u0: $u0 $y->[-2] $y->[-1] l ",scalar @$y,"\n"
854     # if $v1 == 0;
855     $q = (($u0 == $v1) ? $MAX_VAL : int(($u0*$MBASE+$u1)/$v1));
856     --$q while ($v2*$q > ($u0*$MBASE+$u1-$q*$v1)*$MBASE+$u2);
857     if ($q)
858       {
859       ($car, $bar) = (0,0);
860       for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi) 
861         {
862         $prd = $q * $y->[$yi] + $car;
863         $prd -= ($car = int($prd / $MBASE)) * $MBASE;
864         $x->[$xi] += $MBASE if ($bar = (($x->[$xi] -= $prd + $bar) < 0));
865         }
866       if ($x->[-1] < $car + $bar) 
867         {
868         $car = 0; --$q;
869         for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi) 
870           {
871           $x->[$xi] -= $MBASE
872            if ($car = (($x->[$xi] += $y->[$yi] + $car) >= $MBASE));
873           }
874         }   
875       }
876     pop(@$x); unshift(@q, $q);
877     }
878   if (wantarray) 
879     {
880     @d = ();
881     if ($dd != 1)  
882       {
883       $car = 0; 
884       for $xi (reverse @$x) 
885         {
886         $prd = $car * $MBASE + $xi;
887         $car = $prd - ($tmp = int($prd / $dd)) * $dd;
888         unshift(@d, $tmp);
889         }
890       }
891     else 
892       {
893       @d = @$x;
894       }
895     @$x = @q;
896     my $d = \@d; 
897     __strip_zeros($x);
898     __strip_zeros($d);
899     return ($x,$d);
900     }
901   @$x = @q;
902   __strip_zeros($x);
903   $x;
904   }
905
906 ##############################################################################
907 # testing
908
909 sub _acmp
910   {
911   # internal absolute post-normalized compare (ignore signs)
912   # ref to array, ref to array, return <0, 0, >0
913   # arrays must have at least one entry; this is not checked for
914
915   my ($c,$cx,$cy) = @_;
916
917   # fast comp based on number of array elements (aka pseudo-length)
918   my $lxy = scalar @$cx - scalar @$cy;
919   return -1 if $lxy < 0;                                # already differs, ret
920   return 1 if $lxy > 0;                                 # ditto
921
922   # now calculate length based on digits, not parts
923   # we need only the length of the last element, since both array have the
924   # same number of parts
925   $lxy = length(int($cx->[-1])) - length(int($cy->[-1]));
926   return -1 if $lxy < 0;
927   return 1 if $lxy > 0;
928
929   # hm, same lengths,  but same contents? So we need to check all parts:
930   my $a; my $j = scalar @$cx - 1;
931   # manual way (abort if unequal, good for early ne)
932   while ($j >= 0)
933     {
934     last if ($a = $cx->[$j] - $cy->[$j]); $j--;
935     }
936   return 1 if $a > 0;
937   return -1 if $a < 0;
938   0;                                            # numbers are equal
939   }
940
941 sub _len
942   {
943   # compute number of digits
944
945   # int() because add/sub sometimes leaves strings (like '00005') instead of
946   # '5' in this place, thus causing length() to report wrong length
947   my $cx = $_[1];
948
949   (@$cx-1)*$BASE_LEN+length(int($cx->[-1]));
950   }
951
952 sub _digit
953   {
954   # return the nth digit, negative values count backward
955   # zero is rightmost, so _digit(123,0) will give 3
956   my ($c,$x,$n) = @_;
957
958   my $len = _len('',$x);
959
960   $n = $len+$n if $n < 0;               # -1 last, -2 second-to-last
961   $n = abs($n);                         # if negative was too big
962   $len--; $n = $len if $n > $len;       # n to big?
963   
964   my $elem = int($n / $BASE_LEN);       # which array element
965   my $digit = $n % $BASE_LEN;           # which digit in this element
966   $elem = '0000'.@$x[$elem];            # get element padded with 0's
967   substr($elem,-$digit-1,1);
968   }
969
970 sub _zeros
971   {
972   # return amount of trailing zeros in decimal
973   # check each array elem in _m for having 0 at end as long as elem == 0
974   # Upon finding a elem != 0, stop
975   my $x = $_[1];
976   my $zeros = 0; my $elem;
977   foreach my $e (@$x)
978     {
979     if ($e != 0)
980       {
981       $elem = "$e";                             # preserve x
982       $elem =~ s/.*?(0*$)/$1/;                  # strip anything not zero
983       $zeros *= $BASE_LEN;                      # elems * 5
984       $zeros += length($elem);                  # count trailing zeros
985       last;                                     # early out
986       }
987     $zeros ++;                                  # real else branch: 50% slower!
988     }
989   $zeros;
990   }
991
992 ##############################################################################
993 # _is_* routines
994
995 sub _is_zero
996   {
997   # return true if arg (BINT or num_str) is zero (array '+', '0')
998   my $x = $_[1];
999
1000   (((scalar @$x == 1) && ($x->[0] == 0))) <=> 0;
1001   }
1002
1003 sub _is_even
1004   {
1005   # return true if arg (BINT or num_str) is even
1006   my $x = $_[1];
1007   (!($x->[0] & 1)) <=> 0; 
1008   }
1009
1010 sub _is_odd
1011   {
1012   # return true if arg (BINT or num_str) is even
1013   my $x = $_[1];
1014
1015   (($x->[0] & 1)) <=> 0; 
1016   }
1017
1018 sub _is_one
1019   {
1020   # return true if arg (BINT or num_str) is one (array '+', '1')
1021   my $x = $_[1];
1022
1023   (scalar @$x == 1) && ($x->[0] == 1) <=> 0; 
1024   }
1025
1026 sub __strip_zeros
1027   {
1028   # internal normalization function that strips leading zeros from the array
1029   # args: ref to array
1030   my $s = shift;
1031  
1032   my $cnt = scalar @$s; # get count of parts
1033   my $i = $cnt-1;
1034   push @$s,0 if $i < 0;         # div might return empty results, so fix it
1035
1036   return $s if @$s == 1;                # early out
1037
1038   #print "strip: cnt $cnt i $i\n";
1039   # '0', '3', '4', '0', '0',
1040   #  0    1    2    3    4
1041   # cnt = 5, i = 4
1042   # i = 4
1043   # i = 3
1044   # => fcnt = cnt - i (5-2 => 3, cnt => 5-1 = 4, throw away from 4th pos)
1045   # >= 1: skip first part (this can be zero)
1046   while ($i > 0) { last if $s->[$i] != 0; $i--; }
1047   $i++; splice @$s,$i if ($i < $cnt); # $i cant be 0
1048   $s;                                                                    
1049   }                                                                             
1050
1051 ###############################################################################
1052 # check routine to test internal state of corruptions
1053
1054 sub _check
1055   {
1056   # used by the test suite
1057   my $x = $_[1];
1058
1059   return "$x is not a reference" if !ref($x);
1060
1061   # are all parts are valid?
1062   my $i = 0; my $j = scalar @$x; my ($e,$try);
1063   while ($i < $j)
1064     {
1065     $e = $x->[$i]; $e = 'undef' unless defined $e;
1066     $try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e)";
1067     last if $e !~ /^[+]?[0-9]+$/;
1068     $try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e) (stringify)";
1069     last if "$e" !~ /^[+]?[0-9]+$/;
1070     $try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e) (cat-stringify)";
1071     last if '' . "$e" !~ /^[+]?[0-9]+$/;
1072     $try = ' < 0 || >= $BASE; '."($x, $e)";
1073     last if $e <0 || $e >= $BASE;
1074     # this test is disabled, since new/bnorm and certain ops (like early out
1075     # in add/sub) are allowed/expected to leave '00000' in some elements
1076     #$try = '=~ /^00+/; '."($x, $e)";
1077     #last if $e =~ /^00+/;
1078     $i++;
1079     }
1080   return "Illegal part '$e' at pos $i (tested: $try)" if $i < $j;
1081   return 0;
1082   }
1083
1084
1085 ###############################################################################
1086 ###############################################################################
1087 # some optional routines to make BigInt faster
1088
1089 sub _mod
1090   {
1091   # if possible, use mod shortcut
1092   my ($c,$x,$yo) = @_;
1093
1094   # slow way since $y to big
1095   if (scalar @$yo > 1)
1096     {
1097     my ($xo,$rem) = _div($c,$x,$yo);
1098     return $rem;
1099     }
1100
1101   my $y = $yo->[0];
1102   # both are single element arrays
1103   if (scalar @$x == 1)
1104     {
1105     $x->[0] %= $y;
1106     return $x;
1107     }
1108
1109   # @y is a single element, but @x has more than one element
1110   my $b = $BASE % $y;
1111   if ($b == 0)
1112     {
1113     # when BASE % Y == 0 then (B * BASE) % Y == 0
1114     # (B * BASE) % $y + A % Y => A % Y
1115     # so need to consider only last element: O(1)
1116     $x->[0] %= $y;
1117     }
1118   elsif ($b == 1)
1119     {
1120     # else need to go trough all elements: O(N), but loop is a bit simplified
1121     my $r = 0;
1122     foreach (@$x)
1123       {
1124       $r = ($r + $_) % $y;              # not much faster, but heh...
1125       #$r += $_ % $y; $r %= $y;
1126       }
1127     $r = 0 if $r == $y;
1128     $x->[0] = $r;
1129     }
1130   else
1131     {
1132     # else need to go trough all elements: O(N)
1133     my $r = 0; my $bm = 1;
1134     foreach (@$x)
1135       {
1136       $r = ($_ * $bm + $r) % $y;
1137       $bm = ($bm * $b) % $y;
1138
1139       #$r += ($_ % $y) * $bm;
1140       #$bm *= $b;
1141       #$bm %= $y;
1142       #$r %= $y;
1143       }
1144     $r = 0 if $r == $y;
1145     $x->[0] = $r;
1146     }
1147   splice (@$x,1);
1148   $x;
1149   }
1150
1151 ##############################################################################
1152 # shifts
1153
1154 sub _rsft
1155   {
1156   my ($c,$x,$y,$n) = @_;
1157
1158   if ($n != 10)
1159     {
1160     $n = _new($c,\$n); return _div($c,$x, _pow($c,$n,$y));
1161     }
1162
1163   # shortcut (faster) for shifting by 10)
1164   # multiples of $BASE_LEN
1165   my $dst = 0;                          # destination
1166   my $src = _num($c,$y);                # as normal int
1167   my $xlen = (@$x-1)*$BASE_LEN+length(int($x->[-1]));  # len of x in digits
1168   if ($src > $xlen or ($src == $xlen and ! defined $x->[1]))
1169     {
1170     # 12345 67890 shifted right by more than 10 digits => 0
1171     splice (@$x,1);                    # leave only one element
1172     $x->[0] = 0;                       # set to zero
1173     return $x;
1174     }
1175   my $rem = $src % $BASE_LEN;           # remainder to shift
1176   $src = int($src / $BASE_LEN);         # source
1177   if ($rem == 0)
1178     {
1179     splice (@$x,0,$src);                # even faster, 38.4 => 39.3
1180     }
1181   else
1182     {
1183     my $len = scalar @$x - $src;        # elems to go
1184     my $vd; my $z = '0'x $BASE_LEN;
1185     $x->[scalar @$x] = 0;               # avoid || 0 test inside loop
1186     while ($dst < $len)
1187       {
1188       $vd = $z.$x->[$src];
1189       $vd = substr($vd,-$BASE_LEN,$BASE_LEN-$rem);
1190       $src++;
1191       $vd = substr($z.$x->[$src],-$rem,$rem) . $vd;
1192       $vd = substr($vd,-$BASE_LEN,$BASE_LEN) if length($vd) > $BASE_LEN;
1193       $x->[$dst] = int($vd);
1194       $dst++;
1195       }
1196     splice (@$x,$dst) if $dst > 0;              # kill left-over array elems
1197     pop @$x if $x->[-1] == 0 && @$x > 1;        # kill last element if 0
1198     } # else rem == 0
1199   $x;
1200   }
1201
1202 sub _lsft
1203   {
1204   my ($c,$x,$y,$n) = @_;
1205
1206   if ($n != 10)
1207     {
1208     $n = _new($c,\$n); return _mul($c,$x, _pow($c,$n,$y));
1209     }
1210
1211   # shortcut (faster) for shifting by 10) since we are in base 10eX
1212   # multiples of $BASE_LEN:
1213   my $src = scalar @$x;                 # source
1214   my $len = _num($c,$y);                # shift-len as normal int
1215   my $rem = $len % $BASE_LEN;           # remainder to shift
1216   my $dst = $src + int($len/$BASE_LEN); # destination
1217   my $vd;                               # further speedup
1218   $x->[$src] = 0;                       # avoid first ||0 for speed
1219   my $z = '0' x $BASE_LEN;
1220   while ($src >= 0)
1221     {
1222     $vd = $x->[$src]; $vd = $z.$vd;
1223     $vd = substr($vd,-$BASE_LEN+$rem,$BASE_LEN-$rem);
1224     $vd .= $src > 0 ? substr($z.$x->[$src-1],-$BASE_LEN,$rem) : '0' x $rem;
1225     $vd = substr($vd,-$BASE_LEN,$BASE_LEN) if length($vd) > $BASE_LEN;
1226     $x->[$dst] = int($vd);
1227     $dst--; $src--;
1228     }
1229   # set lowest parts to 0
1230   while ($dst >= 0) { $x->[$dst--] = 0; }
1231   # fix spurios last zero element
1232   splice @$x,-1 if $x->[-1] == 0;
1233   $x;
1234   }
1235
1236 sub _pow
1237   {
1238   # power of $x to $y
1239   # ref to array, ref to array, return ref to array
1240   my ($c,$cx,$cy) = @_;
1241
1242   my $pow2 = _one();
1243
1244   my $y_bin = ${_as_bin($c,$cy)}; $y_bin =~ s/^0b//;
1245   my $len = length($y_bin);
1246   while (--$len > 0)
1247     {
1248     _mul($c,$pow2,$cx) if substr($y_bin,$len,1) eq '1';         # is odd?
1249     _mul($c,$cx,$cx);
1250     }
1251
1252   _mul($c,$cx,$pow2);
1253   $cx;
1254   }
1255
1256 sub _fac
1257   {
1258   # factorial of $x
1259   # ref to array, return ref to array
1260   my ($c,$cx) = @_;
1261
1262   if ((@$cx == 1) && ($cx->[0] <= 2))
1263     {
1264     $cx->[0] = 1 * ($cx->[0]||1); # 0,1 => 1, 2 => 2
1265     return $cx;
1266     }
1267
1268   # go forward until $base is exceeded
1269   # limit is either $x or $base (x == 100 means as result too high)
1270   my $steps = 100; $steps = $cx->[0] if @$cx == 1;
1271   my $r = 2; my $cf = 3; my $step = 1; my $last = $r;
1272   while ($r < $BASE && $step < $steps)
1273     {
1274     $last = $r; $r *= $cf++; $step++;
1275     }
1276   if ((@$cx == 1) && ($step == $cx->[0]))
1277     {
1278     # completely done
1279     $cx = [$last];
1280     return $cx;
1281     }
1282   # now we must do the left over steps
1283
1284   # do so as long as n has more than one element
1285   my $n = $cx->[0];
1286   # as soon as the last element of $cx is 0, we split it up and remember how
1287   # many zeors we got so far. The reason is that n! will accumulate zeros at
1288   # the end rather fast.
1289   my $zero_elements = 0;
1290   $cx = [$last];
1291   if (scalar @$cx == 1)
1292     {
1293     my $n = _copy($c,$cx);
1294     # no need to test for $steps, since $steps is a scalar and we stop before
1295     while (scalar @$n != 1)
1296       {
1297       if ($cx->[0] == 0)
1298         {
1299         $zero_elements ++; shift @$cx;
1300         }
1301       _mul($c,$cx,$n); _dec($c,$n);
1302       }
1303     $n = $n->[0];               # "convert" to scalar
1304     }
1305   
1306   # the left over steps will fit into a scalar, so we can speed it up
1307   while ($n != $step)
1308     {
1309     if ($cx->[0] == 0)
1310       {
1311       $zero_elements ++; shift @$cx;
1312       }
1313     _mul($c,$cx,[$n]); $n--;
1314     }
1315   # multiply in the zeros again
1316   while ($zero_elements-- > 0)
1317     {
1318     unshift @$cx, 0; 
1319     }
1320   $cx;
1321   }
1322
1323 # for debugging:
1324   use constant DEBUG => 0;
1325   my $steps = 0;
1326   sub steps { $steps };
1327
1328 sub _sqrt
1329   {
1330   # square-root of $x in place
1331   # Compute a guess of the result (by rule of thumb), then improve it via
1332   # Newton's method.
1333   my ($c,$x) = @_;
1334
1335   if (scalar @$x == 1)
1336     {
1337     # fit's into one Perl scalar, so result can be computed directly
1338     $x->[0] = int(sqrt($x->[0]));
1339     return $x;
1340     } 
1341   my $y = _copy($c,$x);
1342   # hopefully _len/2 is < $BASE, the -1 is to always undershot the guess
1343   # since our guess will "grow"
1344   my $l = int((_len($c,$x)-1) / 2);     
1345
1346   my $lastelem = $x->[-1];                                      # for guess
1347   my $elems = scalar @$x - 1;
1348   # not enough digits, but could have more?
1349   if ((length($lastelem) <= 3) && ($elems > 1))
1350     {
1351     # right-align with zero pad
1352     my $len = length($lastelem) & 1;
1353     print "$lastelem => " if DEBUG;
1354     $lastelem .= substr($x->[-2] . '0' x $BASE_LEN,0,$BASE_LEN);
1355     # former odd => make odd again, or former even to even again
1356     $lastelem = $lastelem / 10 if (length($lastelem) & 1) != $len;
1357     print "$lastelem\n" if DEBUG;
1358     }
1359
1360   # construct $x (instead of _lsft($c,$x,$l,10)
1361   my $r = $l % $BASE_LEN;       # 10000 00000 00000 00000 ($BASE_LEN=5)
1362   $l = int($l / $BASE_LEN);
1363   print "l =  $l " if DEBUG;
1364
1365   splice @$x,$l;                # keep ref($x), but modify it
1366
1367   # we make the first part of the guess not '1000...0' but int(sqrt($lastelem))
1368   # that gives us:
1369   # 14400 00000 => sqrt(14400) => guess first digits to be 120
1370   # 144000 000000 => sqrt(144000) => guess 379
1371
1372   print "$lastelem (elems $elems) => " if DEBUG;
1373   $lastelem = $lastelem / 10 if ($elems & 1 == 1);              # odd or even?
1374   my $g = sqrt($lastelem); $g =~ s/\.//;                        # 2.345 => 2345
1375   $r -= 1 if $elems & 1 == 0;                                   # 70 => 7
1376
1377   # padd with zeros if result is too short
1378   $x->[$l--] = int(substr($g . '0' x $r,0,$r+1));
1379   print "now ",$x->[-1] if DEBUG;
1380   print " would have been ", int('1' . '0' x $r),"\n" if DEBUG;
1381
1382   # If @$x > 1, we could compute the second elem of the guess, too, to create
1383   # an even better guess. Not implemented yet. Does it improve performance?
1384   $x->[$l--] = 0 while ($l >= 0);       # all other digits of guess are zero
1385
1386   print "start x= ",${_str($c,$x)},"\n" if DEBUG;
1387   my $two = _two();
1388   my $last = _zero();
1389   my $lastlast = _zero();
1390   $steps = 0 if DEBUG;
1391   while (_acmp($c,$last,$x) != 0 && _acmp($c,$lastlast,$x) != 0)
1392     {
1393     $steps++ if DEBUG;
1394     $lastlast = _copy($c,$last);
1395     $last = _copy($c,$x);
1396     _add($c,$x, _div($c,_copy($c,$y),$x));
1397     _div($c,$x, $two );
1398     print " x= ",${_str($c,$x)},"\n" if DEBUG;
1399     }
1400   print "\nsteps in sqrt: $steps, " if DEBUG;
1401   _dec($c,$x) if _acmp($c,$y,_mul($c,_copy($c,$x),$x)) < 0;     # overshot? 
1402   print " final ",$x->[-1],"\n" if DEBUG;
1403   $x;
1404   }
1405
1406 sub _root
1407   {
1408   # take n'th root of $x in place (n >= 3)
1409   # Compute a guess of the result (by rule of thumb), then improve it via
1410   # Newton's method.
1411   my ($c,$x,$n) = @_;
1412  
1413   if (scalar @$x == 1)
1414     {
1415     if (scalar @$n > 1)
1416       {
1417       # result will always be smaller than 2 so trunc to 1 at once
1418       $x->[0] = 1;
1419       }
1420     else
1421       {
1422       # fit's into one Perl scalar, so result can be computed directly
1423       $x->[0] = int( $x->[0] ** (1 / $n->[0]) );
1424       }
1425     return $x;
1426     } 
1427
1428   # XXX TODO
1429
1430   $x; 
1431   }
1432
1433 ##############################################################################
1434 # binary stuff
1435
1436 sub _and
1437   {
1438   my ($c,$x,$y) = @_;
1439
1440   # the shortcut makes equal, large numbers _really_ fast, and makes only a
1441   # very small performance drop for small numbers (e.g. something with less
1442   # than 32 bit) Since we optimize for large numbers, this is enabled.
1443   return $x if _acmp($c,$x,$y) == 0;            # shortcut
1444   
1445   my $m = _one(); my ($xr,$yr);
1446   my $mask = $AND_MASK;
1447
1448   my $x1 = $x;
1449   my $y1 = _copy($c,$y);                        # make copy
1450   $x = _zero();
1451   my ($b,$xrr,$yrr);
1452   use integer;
1453   while (!_is_zero($c,$x1) && !_is_zero($c,$y1))
1454     {
1455     ($x1, $xr) = _div($c,$x1,$mask);
1456     ($y1, $yr) = _div($c,$y1,$mask);
1457
1458     # make ints() from $xr, $yr
1459     # this is when the AND_BITS are greater tahn $BASE and is slower for
1460     # small (<256 bits) numbers, but faster for large numbers. Disabled
1461     # due to KISS principle
1462
1463 #    $b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
1464 #    $b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
1465 #    _add($c,$x, _mul($c, _new( $c, \($xrr & $yrr) ), $m) );
1466     
1467     # 0+ due to '&' doesn't work in strings
1468     _add($c,$x, _mul($c, [ 0+$xr->[0] & 0+$yr->[0] ], $m) );
1469     _mul($c,$m,$mask);
1470     }
1471   $x;
1472   }
1473
1474 sub _xor
1475   {
1476   my ($c,$x,$y) = @_;
1477
1478   return _zero() if _acmp($c,$x,$y) == 0;       # shortcut (see -and)
1479
1480   my $m = _one(); my ($xr,$yr);
1481   my $mask = $XOR_MASK;
1482
1483   my $x1 = $x;
1484   my $y1 = _copy($c,$y);                        # make copy
1485   $x = _zero();
1486   my ($b,$xrr,$yrr);
1487   use integer;
1488   while (!_is_zero($c,$x1) && !_is_zero($c,$y1))
1489     {
1490     ($x1, $xr) = _div($c,$x1,$mask);
1491     ($y1, $yr) = _div($c,$y1,$mask);
1492     # make ints() from $xr, $yr (see _and())
1493     #$b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
1494     #$b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
1495     #_add($c,$x, _mul($c, _new( $c, \($xrr ^ $yrr) ), $m) );
1496
1497     # 0+ due to '^' doesn't work in strings
1498     _add($c,$x, _mul($c, [ 0+$xr->[0] ^ 0+$yr->[0] ], $m) );
1499     _mul($c,$m,$mask);
1500     }
1501   # the loop stops when the shorter of the two numbers is exhausted
1502   # the remainder of the longer one will survive bit-by-bit, so we simple
1503   # multiply-add it in
1504   _add($c,$x, _mul($c, $x1, $m) ) if !_is_zero($c,$x1);
1505   _add($c,$x, _mul($c, $y1, $m) ) if !_is_zero($c,$y1);
1506   
1507   $x;
1508   }
1509
1510 sub _or
1511   {
1512   my ($c,$x,$y) = @_;
1513
1514   return $x if _acmp($c,$x,$y) == 0;            # shortcut (see _and)
1515
1516   my $m = _one(); my ($xr,$yr);
1517   my $mask = $OR_MASK;
1518
1519   my $x1 = $x;
1520   my $y1 = _copy($c,$y);                        # make copy
1521   $x = _zero();
1522   my ($b,$xrr,$yrr);
1523   use integer;
1524   while (!_is_zero($c,$x1) && !_is_zero($c,$y1))
1525     {
1526     ($x1, $xr) = _div($c,$x1,$mask);
1527     ($y1, $yr) = _div($c,$y1,$mask);
1528     # make ints() from $xr, $yr (see _and())
1529 #    $b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
1530 #    $b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
1531 #    _add($c,$x, _mul($c, _new( $c, \($xrr | $yrr) ), $m) );
1532     
1533     # 0+ due to '|' doesn't work in strings
1534     _add($c,$x, _mul($c, [ 0+$xr->[0] | 0+$yr->[0] ], $m) );
1535     _mul($c,$m,$mask);
1536     }
1537   # the loop stops when the shorter of the two numbers is exhausted
1538   # the remainder of the longer one will survive bit-by-bit, so we simple
1539   # multiply-add it in
1540   _add($c,$x, _mul($c, $x1, $m) ) if !_is_zero($c,$x1);
1541   _add($c,$x, _mul($c, $y1, $m) ) if !_is_zero($c,$y1);
1542   
1543   $x;
1544   }
1545
1546 sub _as_hex
1547   {
1548   # convert a decimal number to hex (ref to array, return ref to string)
1549   my ($c,$x) = @_;
1550
1551   # fit's into one element
1552   if (@$x == 1)
1553     {
1554     my $t = '0x' . sprintf("%x",$x->[0]);
1555     return \$t;
1556     }
1557
1558   my $x1 = _copy($c,$x);
1559
1560   my $es = '';
1561   my ($xr, $h, $x10000);
1562   if ($] >= 5.006)
1563     {
1564     $x10000 = [ 0x10000 ]; $h = 'h4';
1565     }
1566   else
1567     {
1568     $x10000 = [ 0x1000 ]; $h = 'h3';
1569     }
1570   while (! _is_zero($c,$x1))
1571     {
1572     ($x1, $xr) = _div($c,$x1,$x10000);
1573     $es .= unpack($h,pack('v',$xr->[0]));
1574     }
1575   $es = reverse $es;
1576   $es =~ s/^[0]+//;   # strip leading zeros
1577   $es = '0x' . $es;
1578   \$es;
1579   }
1580
1581 sub _as_bin
1582   {
1583   # convert a decimal number to bin (ref to array, return ref to string)
1584   my ($c,$x) = @_;
1585
1586   # fit's into one element
1587   if (@$x == 1)
1588     {
1589     my $t = '0b' . sprintf("%b",$x->[0]);
1590     return \$t;
1591     }
1592   my $x1 = _copy($c,$x);
1593
1594   my $es = '';
1595   my ($xr, $b, $x10000);
1596   if ($] >= 5.006)
1597     {
1598     $x10000 = [ 0x10000 ]; $b = 'b16';
1599     }
1600   else
1601     {
1602     $x10000 = [ 0x1000 ]; $b = 'b12';
1603     }
1604   while (! _is_zero($c,$x1))
1605     {
1606     ($x1, $xr) = _div($c,$x1,$x10000);
1607     $es .= unpack($b,pack('v',$xr->[0]));
1608     }
1609   $es = reverse $es;
1610   $es =~ s/^[0]+//;   # strip leading zeros
1611   $es = '0b' . $es;
1612   \$es;
1613   }
1614
1615 sub _from_hex
1616   {
1617   # convert a hex number to decimal (ref to string, return ref to array)
1618   my ($c,$hs) = @_;
1619
1620   my $mul = _one();
1621   my $m = [ 0x10000 ];                          # 16 bit at a time
1622   my $x = _zero();
1623
1624   my $len = length($$hs)-2;
1625   $len = int($len/4);                           # 4-digit parts, w/o '0x'
1626   my $val; my $i = -4;
1627   while ($len >= 0)
1628     {
1629     $val = substr($$hs,$i,4);
1630     $val =~ s/^[+-]?0x// if $len == 0;          # for last part only because
1631     $val = hex($val);                           # hex does not like wrong chars
1632     $i -= 4; $len --;
1633     _add ($c, $x, _mul ($c, [ $val ], $mul ) ) if $val != 0;
1634     _mul ($c, $mul, $m ) if $len >= 0;          # skip last mul
1635     }
1636   $x;
1637   }
1638
1639 sub _from_bin
1640   {
1641   # convert a hex number to decimal (ref to string, return ref to array)
1642   my ($c,$bs) = @_;
1643
1644   # instead of converting 8 bit at a time, it is faster to convert the
1645   # number to hex, and then call _from_hex.
1646
1647   my $hs = $$bs;
1648   $hs =~ s/^[+-]?0b//;                                  # remove sign and 0b
1649   my $l = length($hs);                                  # bits
1650   $hs = '0' x (8-($l % 8)) . $hs if ($l % 8) != 0;      # padd left side w/ 0
1651   my $h = unpack('H*', pack ('B*', $hs));               # repack as hex
1652   return $c->_from_hex(\('0x'.$h));
1653  
1654   my $mul = _one();
1655   my $m = [ 0x100 ];                            # 8 bit at a time
1656   my $x = _zero();
1657
1658   my $len = length($$bs)-2;
1659   $len = int($len/8);                           # 4-digit parts, w/o '0x'
1660   my $val; my $i = -8;
1661   while ($len >= 0)
1662     {
1663     $val = substr($$bs,$i,8);
1664     $val =~ s/^[+-]?0b// if $len == 0;          # for last part only
1665
1666     $val = ord(pack('B8',substr('00000000'.$val,-8,8))); 
1667
1668     $i -= 8; $len --;
1669     _add ($c, $x, _mul ($c, [ $val ], $mul ) ) if $val != 0;
1670     _mul ($c, $mul, $m ) if $len >= 0;          # skip last mul
1671     }
1672   $x;
1673   }
1674
1675 ##############################################################################
1676 # special modulus functions
1677
1678 sub _modinv
1679   {
1680   # modular inverse
1681   my ($c,$x,$y) = @_;
1682
1683   my $u = _zero($c); my $u1 = _one($c);
1684   my $a = _copy($c,$y); my $b = _copy($c,$x);
1685
1686   # Euclid's Algorithm for bgcd(), only that we calc bgcd() ($a) and the
1687   # result ($u) at the same time. See comments in BigInt for why this works.
1688   my $q;
1689   ($a, $q, $b) = ($b, _div($c,$a,$b));          # step 1
1690   my $sign = 1;
1691   while (!_is_zero($c,$b))
1692     {
1693     my $t = _add($c,                            # step 2:
1694        _mul($c,_copy($c,$u1), $q) ,             #  t =  u1 * q
1695        $u );                                    #     + u
1696     $u = $u1;                                   #  u = u1, u1 = t
1697     $u1 = $t;
1698     $sign = -$sign;
1699     ($a, $q, $b) = ($b, _div($c,$a,$b));        # step 1
1700     }
1701
1702   # if the gcd is not 1, then return NaN
1703   return (undef,undef) unless _is_one($c,$a);
1704  
1705   $sign = $sign == 1 ? '+' : '-';
1706   ($u1,$sign);
1707   }
1708
1709 sub _modpow
1710   {
1711   # modulus of power ($x ** $y) % $z
1712   my ($c,$num,$exp,$mod) = @_;
1713
1714   # in the trivial case,
1715   if (_is_one($c,$mod))
1716     {
1717     splice @$num,0,1; $num->[0] = 0;
1718     return $num;
1719     }
1720   if ((scalar @$num == 1) && (($num->[0] == 0) || ($num->[0] == 1)))
1721     {
1722     $num->[0] = 1;
1723     return $num;
1724     }
1725
1726 #  $num = _mod($c,$num,$mod);   # this does not make it faster
1727
1728   my $acc = _copy($c,$num); my $t = _one();
1729
1730   my $expbin = ${_as_bin($c,$exp)}; $expbin =~ s/^0b//;
1731   my $len = length($expbin);
1732   while (--$len >= 0)
1733     {
1734     if ( substr($expbin,$len,1) eq '1')                 # is_odd
1735       {
1736       _mul($c,$t,$acc);
1737       $t = _mod($c,$t,$mod);
1738       }
1739     _mul($c,$acc,$acc);
1740     $acc = _mod($c,$acc,$mod);
1741     }
1742   @$num = @$t;
1743   $num;
1744   }
1745
1746 ##############################################################################
1747 ##############################################################################
1748
1749 1;
1750 __END__
1751
1752 =head1 NAME
1753
1754 Math::BigInt::Calc - Pure Perl module to support Math::BigInt
1755
1756 =head1 SYNOPSIS
1757
1758 Provides support for big integer calculations. Not intended to be used by other
1759 modules (except Math::BigInt::Cached). Other modules which sport the same
1760 functions can also be used to support Math::BigInt, like Math::BigInt::Pari.
1761
1762 =head1 DESCRIPTION
1763
1764 In order to allow for multiple big integer libraries, Math::BigInt was
1765 rewritten to use library modules for core math routines. Any module which
1766 follows the same API as this can be used instead by using the following:
1767
1768         use Math::BigInt lib => 'libname';
1769
1770 'libname' is either the long name ('Math::BigInt::Pari'), or only the short
1771 version like 'Pari'.
1772
1773 =head1 STORAGE
1774
1775 =head1 METHODS
1776
1777 The following functions MUST be defined in order to support the use by
1778 Math::BigInt:
1779
1780         _new(string)    return ref to new object from ref to decimal string
1781         _zero()         return a new object with value 0
1782         _one()          return a new object with value 1
1783
1784         _str(obj)       return ref to a string representing the object
1785         _num(obj)       returns a Perl integer/floating point number
1786                         NOTE: because of Perl numeric notation defaults,
1787                         the _num'ified obj may lose accuracy due to 
1788                         machine-dependend floating point size limitations
1789                     
1790         _add(obj,obj)   Simple addition of two objects
1791         _mul(obj,obj)   Multiplication of two objects
1792         _div(obj,obj)   Division of the 1st object by the 2nd
1793                         In list context, returns (result,remainder).
1794                         NOTE: this is integer math, so no
1795                         fractional part will be returned.
1796                         The second operand will be not be 0, so no need to
1797                         check for that.
1798         _sub(obj,obj)   Simple subtraction of 1 object from another
1799                         a third, optional parameter indicates that the params
1800                         are swapped. In this case, the first param needs to
1801                         be preserved, while you can destroy the second.
1802                         sub (x,y,1) => return x - y and keep x intact!
1803         _dec(obj)       decrement object by one (input is garant. to be > 0)
1804         _inc(obj)       increment object by one
1805
1806
1807         _acmp(obj,obj)  <=> operator for objects (return -1, 0 or 1)
1808
1809         _len(obj)       returns count of the decimal digits of the object
1810         _digit(obj,n)   returns the n'th decimal digit of object
1811
1812         _is_one(obj)    return true if argument is +1
1813         _is_zero(obj)   return true if argument is 0
1814         _is_even(obj)   return true if argument is even (0,2,4,6..)
1815         _is_odd(obj)    return true if argument is odd (1,3,5,7..)
1816
1817         _copy           return a ref to a true copy of the object
1818
1819         _check(obj)     check whether internal representation is still intact
1820                         return 0 for ok, otherwise error message as string
1821
1822 The following functions are optional, and can be defined if the underlying lib
1823 has a fast way to do them. If undefined, Math::BigInt will use pure Perl (hence
1824 slow) fallback routines to emulate these:
1825
1826         _from_hex(str)  return ref to new object from ref to hexadecimal string
1827         _from_bin(str)  return ref to new object from ref to binary string
1828         
1829         _as_hex(str)    return ref to scalar string containing the value as
1830                         unsigned hex string, with the '0x' prepended.
1831                         Leading zeros must be stripped.
1832         _as_bin(str)    Like as_hex, only as binary string containing only
1833                         zeros and ones. Leading zeros must be stripped and a
1834                         '0b' must be prepended.
1835         
1836         _rsft(obj,N,B)  shift object in base B by N 'digits' right
1837                         For unsupported bases B, return undef to signal failure
1838         _lsft(obj,N,B)  shift object in base B by N 'digits' left
1839                         For unsupported bases B, return undef to signal failure
1840         
1841         _xor(obj1,obj2) XOR (bit-wise) object 1 with object 2
1842                         Note: XOR, AND and OR pad with zeros if size mismatches
1843         _and(obj1,obj2) AND (bit-wise) object 1 with object 2
1844         _or(obj1,obj2)  OR (bit-wise) object 1 with object 2
1845
1846         _mod(obj,obj)   Return remainder of div of the 1st by the 2nd object
1847         _sqrt(obj)      return the square root of object (truncated to int)
1848         _root(obj)      return the n'th (n >= 3) root of obj (truncated to int)
1849         _fac(obj)       return factorial of object 1 (1*2*3*4..)
1850         _pow(obj,obj)   return object 1 to the power of object 2
1851         _gcd(obj,obj)   return Greatest Common Divisor of two objects
1852         
1853         _zeros(obj)     return number of trailing decimal zeros
1854         _modinv         return inverse modulus
1855         _modpow         return modulus of power ($x ** $y) % $z
1856
1857 Input strings come in as unsigned but with prefix (i.e. as '123', '0xabc'
1858 or '0b1101').
1859
1860 So the library needs only to deal with unsigned big integers. Testing of input
1861 parameter validity is done by the caller, so you need not worry about
1862 underflow (f.i. in C<_sub()>, C<_dec()>) nor about division by zero or similar
1863 cases.
1864
1865 The first parameter can be modified, that includes the possibility that you
1866 return a reference to a completely different object instead. Although keeping
1867 the reference and just changing it's contents is prefered over creating and
1868 returning a different reference.
1869
1870 Return values are always references to objects, strings, or true/false for
1871 comparisation routines.
1872
1873 Exceptions are C<_lsft()> and C<_rsft()>, which return undef if they can not
1874 shift the argument. This is used to delegate shifting of bases different than
1875 the one you can support back to Math::BigInt, which will use some generic code
1876 to calculate the result.
1877
1878 =head1 WRAP YOUR OWN
1879
1880 If you want to port your own favourite c-lib for big numbers to the
1881 Math::BigInt interface, you can take any of the already existing modules as
1882 a rough guideline. You should really wrap up the latest BigInt and BigFloat
1883 testsuites with your module, and replace in them any of the following:
1884
1885         use Math::BigInt;
1886
1887 by this:
1888
1889         use Math::BigInt lib => 'yourlib';
1890
1891 This way you ensure that your library really works 100% within Math::BigInt.
1892
1893 =head1 LICENSE
1894  
1895 This program is free software; you may redistribute it and/or modify it under
1896 the same terms as Perl itself. 
1897
1898 =head1 AUTHORS
1899
1900 Original math code by Mark Biggar, rewritten by Tels L<http://bloodgate.com/>
1901 in late 2000.
1902 Seperated from BigInt and shaped API with the help of John Peacock.
1903 Fixed/enhanced by Tels 2001-2002.
1904
1905 =head1 SEE ALSO
1906
1907 L<Math::BigInt>, L<Math::BigFloat>, L<Math::BigInt::BitVect>,
1908 L<Math::BigInt::GMP>, L<Math::BigInt::FastCalc> and L<Math::BigInt::Pari>.
1909
1910 =cut