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