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