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