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