Upgrade to Math::BigInt 1.48.
[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.17';
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 ($BASE,$RBASE,$BASE_LEN,$MAX_VAL,$BASE_LEN2);
34 my ($AND_BITS,$XOR_BITS,$OR_BITS);
35 my ($AND_MASK,$XOR_MASK,$OR_MASK);
36
37 sub _base_len 
38   {
39   # set/get the BASE_LEN and assorted other, connected values
40   # used only be the testsuite, set is used only by the BEGIN block below
41   shift;
42
43   my $b = shift;
44   if (defined $b)
45     {
46     $b = 5 if $^O =~ /^uts/;    # UTS needs 5, because 6 and 7 break
47     $BASE_LEN = $b+1;
48     my $caught;
49     while (--$BASE_LEN > 5)
50       {
51       $BASE = int("1e".$BASE_LEN);
52       $RBASE = abs('1e-'.$BASE_LEN);                    # see USE_MUL
53       $caught = 0;
54       $caught += 1 if (int($BASE * $RBASE) != 1);       # should be 1
55       $caught += 2 if (int($BASE / $BASE) != 1);        # should be 1
56       # print "caught $caught\n";
57       last if $caught != 3;
58       }
59     $BASE = int("1e".$BASE_LEN);
60     $RBASE = abs('1e-'.$BASE_LEN);                      # see USE_MUL
61     $MAX_VAL = $BASE-1;
62     $BASE_LEN2 = int($BASE_LEN / 2);                    # for mul shortcut
63     # print "BASE_LEN: $BASE_LEN MAX_VAL: $MAX_VAL BASE: $BASE RBASE: $RBASE\n";
64     
65     if ($caught & 1 != 0)
66       {
67       # must USE_MUL
68       *{_mul} = \&_mul_use_mul;
69       *{_div} = \&_div_use_mul;
70       }
71     else                # $caught must be 2, since it can't be 1 nor 3
72       {
73       # can USE_DIV instead
74       *{_mul} = \&_mul_use_div;
75       *{_div} = \&_div_use_div;
76       }
77     }
78   if (wantarray)
79     {
80     return ($BASE_LEN, $AND_BITS, $XOR_BITS, $OR_BITS);
81     }
82   $BASE_LEN;
83   }
84
85 BEGIN
86   {
87   # from Daniel Pfeiffer: determine largest group of digits that is precisely
88   # multipliable with itself plus carry
89   # Test now changed to expect the proper pattern, not a result off by 1 or 2
90   my ($e, $num) = 3;    # lowest value we will use is 3+1-1 = 3
91   do 
92     {
93     $num = ('9' x ++$e) + 0;
94     $num *= $num + 1.0;
95     # print "$num $e\n";
96     } while ("$num" =~ /9{$e}0{$e}/);   # must be a certain pattern
97   $e--;                                 # last test failed, so retract one step
98   # the limits below brush the problems with the test above under the rug:
99   # the test should be able to find the proper $e automatically
100   $e = 5 if $^O =~ /^uts/;      # UTS get's some special treatment
101   $e = 5 if $^O =~ /^unicos/;   # unicos is also problematic (6 seems to work
102                                 # there, but we play safe)
103   $e = 8 if $e > 8;             # cap, for VMS, OS/390 and other 64 bit systems
104
105   __PACKAGE__->_base_len($e);   # set and store
106
107   # find out how many bits _and, _or and _xor can take (old default = 16)
108   # I don't think anybody has yet 128 bit scalars, so let's play safe.
109   use integer;
110   local $^W = 0;        # don't warn about 'nonportable number'
111   $AND_BITS = 15; $XOR_BITS = 15; $OR_BITS  = 15;
112
113   # find max bits, we will not go higher than numberofbits that fit into $BASE
114   # to make _and etc simpler (and faster for smaller, slower for large numbers)
115   my $max = 16;
116   while (2 ** $max < $BASE) { $max++; }
117   my ($x,$y,$z);
118   do {
119     $AND_BITS++;
120     $x = oct('0b' . '1' x $AND_BITS); $y = $x & $x;
121     $z = (2 ** $AND_BITS) - 1;
122     } while ($AND_BITS < $max && $x == $z && $y == $x);
123   $AND_BITS --;                                         # retreat one step
124   do {
125     $XOR_BITS++;
126     $x = oct('0b' . '1' x $XOR_BITS); $y = $x ^ 0;
127     $z = (2 ** $XOR_BITS) - 1;
128     } while ($XOR_BITS < $max && $x == $z && $y == $x);
129   $XOR_BITS --;                                         # retreat one step
130   do {
131     $OR_BITS++;
132     $x = oct('0b' . '1' x $OR_BITS); $y = $x | $x;
133     $z = (2 ** $OR_BITS) - 1;
134     } while ($OR_BITS < $max && $x == $z && $y == $x);
135   $OR_BITS --;                                          # retreat one step
136   
137   # print "AND $AND_BITS XOR $XOR_BITS OR $OR_BITS\n";
138   }
139
140 ##############################################################################
141 # create objects from various representations
142
143 sub _new
144   {
145   # (ref to string) return ref to num_array
146   # Convert a number from string format to internal base 100000 format.
147   # Assumes normalized value as input.
148   my $d = $_[1];
149   my $il = CORE::length($$d)-1;
150   # these leaves '00000' instead of int 0 and will be corrected after any op
151   return [ reverse(unpack("a" . ($il % $BASE_LEN+1) 
152     . ("a$BASE_LEN" x ($il / $BASE_LEN)), $$d)) ];
153   }                                                                             
154   
155 BEGIN
156   {
157   $AND_MASK = __PACKAGE__->_new( \( 2 ** $AND_BITS ));
158   $XOR_MASK = __PACKAGE__->_new( \( 2 ** $XOR_BITS ));
159   $OR_MASK = __PACKAGE__->_new( \( 2 ** $OR_BITS ));
160   }
161
162 sub _zero
163   {
164   # create a zero
165   return [ 0 ];
166   }
167
168 sub _one
169   {
170   # create a one
171   return [ 1 ];
172   }
173
174 sub _two
175   {
176   # create a two (for _pow)
177   return [ 2 ];
178   }
179
180 sub _copy
181   {
182   return [ @{$_[1]} ];
183   }
184
185 # catch and throw away
186 sub import { }
187
188 ##############################################################################
189 # convert back to string and number
190
191 sub _str
192   {
193   # (ref to BINT) return num_str
194   # Convert number from internal base 100000 format to string format.
195   # internal format is always normalized (no leading zeros, "-0" => "+0")
196   my $ar = $_[1];
197   my $ret = "";
198   my $l = scalar @$ar;         # number of parts
199   return $nan if $l < 1;       # should not happen
200   # handle first one different to strip leading zeros from it (there are no
201   # leading zero parts in internal representation)
202   $l --; $ret .= $ar->[$l]; $l--;
203   # Interestingly, the pre-padd method uses more time
204   # the old grep variant takes longer (14 to 10 sec)
205   my $z = '0' x ($BASE_LEN-1);                            
206   while ($l >= 0)
207     {
208     $ret .= substr($z.$ar->[$l],-$BASE_LEN); # fastest way I could think of
209     $l--;
210     }
211   return \$ret;
212   }                                                                             
213
214 sub _num
215   {
216   # Make a number (scalar int/float) from a BigInt object
217   my $x = $_[1];
218   return $x->[0] if scalar @$x == 1;  # below $BASE
219   my $fac = 1;
220   my $num = 0;
221   foreach (@$x)
222     {
223     $num += $fac*$_; $fac *= $BASE;
224     }
225   return $num; 
226   }
227
228 ##############################################################################
229 # actual math code
230
231 sub _add
232   {
233   # (ref to int_num_array, ref to int_num_array)
234   # routine to add two base 1eX numbers
235   # stolen from Knuth Vol 2 Algorithm A pg 231
236   # there are separate routines to add and sub as per Knuth pg 233
237   # This routine clobbers up array x, but not y.
238  
239   my ($c,$x,$y) = @_;
240  
241   # for each in Y, add Y to X and carry. If after that, something is left in
242   # X, foreach in X add carry to X and then return X, carry
243   # Trades one "$j++" for having to shift arrays, $j could be made integer
244   # but this would impose a limit to number-length of 2**32.
245   my $i; my $car = 0; my $j = 0;
246   for $i (@$y)
247     {
248     $x->[$j] -= $BASE if $car = (($x->[$j] += $i + $car) >= $BASE) ? 1 : 0;
249     $j++;
250     }
251   while ($car != 0)
252     {
253     $x->[$j] -= $BASE if $car = (($x->[$j] += $car) >= $BASE) ? 1 : 0; $j++;
254     }
255   return $x;
256   }                                                                             
257
258 sub _inc
259   {
260   # (ref to int_num_array, ref to int_num_array)
261   # routine to add 1 to a base 1eX numbers
262   # This routine clobbers up array x, but not y.
263   my ($c,$x) = @_;
264
265   for my $i (@$x)
266     {
267     return $x if (($i += 1) < $BASE);           # early out
268     $i -= $BASE;
269     }
270   if ($x->[-1] == 0)                            # last overflowed
271     {
272     push @$x,1;                                 # extend
273     }
274   return $x;
275   }                                                                             
276
277 sub _dec
278   {
279   # (ref to int_num_array, ref to int_num_array)
280   # routine to add 1 to a base 1eX numbers
281   # This routine clobbers up array x, but not y.
282   my ($c,$x) = @_;
283
284   for my $i (@$x)
285     {
286     last if (($i -= 1) >= 0);                   # early out
287     $i = $MAX_VAL;
288     }
289   pop @$x if $x->[-1] == 0 && @$x > 1;          # last overflowed (but leave 0)
290   return $x;
291   }                                                                             
292
293 sub _sub
294   {
295   # (ref to int_num_array, ref to int_num_array)
296   # subtract base 1eX numbers -- stolen from Knuth Vol 2 pg 232, $x > $y
297   # subtract Y from X (X is always greater/equal!) by modifying x in place
298   my ($c,$sx,$sy,$s) = @_;
299  
300   my $car = 0; my $i; my $j = 0;
301   if (!$s)
302     {
303     #print "case 2\n";
304     for $i (@$sx)
305       {
306       last unless defined $sy->[$j] || $car;
307       $i += $BASE if $car = (($i -= ($sy->[$j] || 0) + $car) < 0); $j++;
308       }
309     # might leave leading zeros, so fix that
310     return __strip_zeros($sx);
311     }
312   #print "case 1 (swap)\n";
313   for $i (@$sx)
314     {
315     last unless defined $sy->[$j] || $car;
316     $sy->[$j] += $BASE
317      if $car = (($sy->[$j] = $i-($sy->[$j]||0) - $car) < 0);
318     $j++;
319     }
320   # might leave leading zeros, so fix that
321   __strip_zeros($sy);
322   }                                                                             
323
324 sub _mul_use_mul
325   {
326   # (BINT, BINT) return nothing
327   # multiply two numbers in internal representation
328   # modifies first arg, second need not be different from first
329   my ($c,$xv,$yv) = @_;
330
331   # shortcut for two very short numbers
332   # +0 since part maybe string '00001' from new()
333   if ((@$xv == 1) && (@$yv == 1)
334    && (length($xv->[0]+0) <= $BASE_LEN2)
335    && (length($yv->[0]+0) <= $BASE_LEN2))
336    {
337    $xv->[0] *= $yv->[0];
338    return $xv;
339    }
340   
341   my @prod = (); my ($prod,$car,$cty,$xi,$yi);
342   # since multiplying $x with $x fails, make copy in this case
343   $yv = [@$xv] if "$xv" eq "$yv";       # same references?
344   for $xi (@$xv)
345     {
346     $car = 0; $cty = 0;
347
348     # slow variant
349 #    for $yi (@$yv)
350 #      {
351 #      $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
352 #      $prod[$cty++] =
353 #       $prod - ($car = int($prod * RBASE)) * $BASE;  # see USE_MUL
354 #      }
355 #    $prod[$cty] += $car if $car; # need really to check for 0?
356 #    $xi = shift @prod;
357
358     # faster variant
359     # looping through this if $xi == 0 is silly - so optimize it away!
360     $xi = (shift @prod || 0), next if $xi == 0;
361     for $yi (@$yv)
362       {
363       $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
364 ##     this is actually a tad slower
365 ##        $prod = $prod[$cty]; $prod += ($car + $xi * $yi);     # no ||0 here
366       $prod[$cty++] =
367        $prod - ($car = int($prod * $RBASE)) * $BASE;  # see USE_MUL
368       }
369     $prod[$cty] += $car if $car; # need really to check for 0?
370     $xi = shift @prod || 0;     # || 0 makes v5.005_3 happy
371     }
372   push @$xv, @prod;
373   __strip_zeros($xv);
374   }                                                                             
375
376 sub _mul_use_div
377   {
378   # (BINT, BINT) return nothing
379   # multiply two numbers in internal representation
380   # modifies first arg, second need not be different from first
381   my ($c,$xv,$yv) = @_;
382  
383   # shortcut for two very short numbers
384   # +0 since part maybe string '00001' from new()
385   if ((@$xv == 1) && (@$yv == 1)
386    && (length($xv->[0]+0) <= $BASE_LEN2)
387    && (length($yv->[0]+0) <= $BASE_LEN2))
388    {
389    $xv->[0] *= $yv->[0];
390    return $xv;
391    }
392   
393   my @prod = (); my ($prod,$car,$cty,$xi,$yi);
394   # since multiplying $x with $x fails, make copy in this case
395   $yv = [@$xv] if "$xv" eq "$yv";       # same references?
396   for $xi (@$xv)
397     {
398     $car = 0; $cty = 0;
399     # looping through this if $xi == 0 is silly - so optimize it away!
400     $xi = (shift @prod || 0), next if $xi == 0;
401     for $yi (@$yv)
402       {
403       $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
404       $prod[$cty++] =
405        $prod - ($car = int($prod / $BASE)) * $BASE;
406       }
407     $prod[$cty] += $car if $car; # need really to check for 0?
408     $xi = shift @prod || 0;     # || 0 makes v5.005_3 happy
409     }
410   push @$xv, @prod;
411   __strip_zeros($xv);
412   }                                                                             
413
414 sub _div_use_mul
415   {
416   # ref to array, ref to array, modify first array and return remainder if 
417   # in list context
418   my ($c,$x,$yorg) = @_;
419   my ($car,$bar,$prd,$dd,$xi,$yi,@q,$v2,$v1);
420
421   my (@d,$tmp,$q,$u2,$u1,$u0);
422
423   $car = $bar = $prd = 0;
424   
425   my $y = [ @$yorg ];
426   if (($dd = int($BASE/($y->[-1]+1))) != 1) 
427     {
428     for $xi (@$x) 
429       {
430       $xi = $xi * $dd + $car;
431       $xi -= ($car = int($xi * $RBASE)) * $BASE;        # see USE_MUL
432       }
433     push(@$x, $car); $car = 0;
434     for $yi (@$y) 
435       {
436       $yi = $yi * $dd + $car;
437       $yi -= ($car = int($yi * $RBASE)) * $BASE;        # see USE_MUL
438       }
439     }
440   else 
441     {
442     push(@$x, 0);
443     }
444   @q = (); ($v2,$v1) = @$y[-2,-1];
445   $v2 = 0 unless $v2;
446   while ($#$x > $#$y) 
447     {
448     ($u2,$u1,$u0) = @$x[-3..-1];
449     $u2 = 0 unless $u2;
450     #warn "oups v1 is 0, u0: $u0 $y->[-2] $y->[-1] l ",scalar @$y,"\n"
451     # if $v1 == 0;
452     # $q = (($u0 == $v1) ? 99999 : int(($u0*$BASE+$u1)/$v1));
453      $q = (($u0 == $v1) ? $MAX_VAL : int(($u0*$BASE+$u1)/$v1));
454     --$q while ($v2*$q > ($u0*$BASE+$u1-$q*$v1)*$BASE+$u2);
455     if ($q)
456       {
457       ($car, $bar) = (0,0);
458       for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi) 
459         {
460         $prd = $q * $y->[$yi] + $car;
461         $prd -= ($car = int($prd * $RBASE)) * $BASE;    # see USE_MUL
462         $x->[$xi] += $BASE if ($bar = (($x->[$xi] -= $prd + $bar) < 0));
463         }
464       if ($x->[-1] < $car + $bar) 
465         {
466         $car = 0; --$q;
467         for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi) 
468           {
469           $x->[$xi] -= $BASE
470            if ($car = (($x->[$xi] += $y->[$yi] + $car) > $BASE));
471           }
472         }   
473       }
474       pop(@$x); unshift(@q, $q);
475     }
476   if (wantarray) 
477     {
478     @d = ();
479     if ($dd != 1)  
480       {
481       $car = 0; 
482       for $xi (reverse @$x) 
483         {
484         $prd = $car * $BASE + $xi;
485         $car = $prd - ($tmp = int($prd / $dd)) * $dd; # see USE_MUL
486         unshift(@d, $tmp);
487         }
488       }
489     else 
490       {
491       @d = @$x;
492       }
493     @$x = @q;
494     __strip_zeros($x); 
495     __strip_zeros(\@d);
496     _check('',$x);
497     _check('',\@d);
498     return ($x,\@d);
499     }
500   @$x = @q;
501   __strip_zeros($x); 
502     _check('',$x);
503   }
504
505 sub _div_use_div
506   {
507   # ref to array, ref to array, modify first array and return remainder if 
508   # in list context
509   my ($c,$x,$yorg) = @_;
510   my ($car,$bar,$prd,$dd,$xi,$yi,@q,$v2,$v1);
511
512   my (@d,$tmp,$q,$u2,$u1,$u0);
513
514   $car = $bar = $prd = 0;
515   
516   my $y = [ @$yorg ];
517   if (($dd = int($BASE/($y->[-1]+1))) != 1) 
518     {
519     for $xi (@$x) 
520       {
521       $xi = $xi * $dd + $car;
522       $xi -= ($car = int($xi / $BASE)) * $BASE;
523       }
524     push(@$x, $car); $car = 0;
525     for $yi (@$y) 
526       {
527       $yi = $yi * $dd + $car;
528       $yi -= ($car = int($yi / $BASE)) * $BASE;
529       }
530     }
531   else 
532     {
533     push(@$x, 0);
534     }
535   @q = (); ($v2,$v1) = @$y[-2,-1];
536   $v2 = 0 unless $v2;
537   while ($#$x > $#$y) 
538     {
539     ($u2,$u1,$u0) = @$x[-3..-1];
540     $u2 = 0 unless $u2;
541     #warn "oups v1 is 0, u0: $u0 $y->[-2] $y->[-1] l ",scalar @$y,"\n"
542     # if $v1 == 0;
543     # $q = (($u0 == $v1) ? 99999 : int(($u0*$BASE+$u1)/$v1));
544      $q = (($u0 == $v1) ? $MAX_VAL : int(($u0*$BASE+$u1)/$v1));
545     --$q while ($v2*$q > ($u0*$BASE+$u1-$q*$v1)*$BASE+$u2);
546     if ($q)
547       {
548       ($car, $bar) = (0,0);
549       for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi) 
550         {
551         $prd = $q * $y->[$yi] + $car;
552         $prd -= ($car = int($prd / $BASE)) * $BASE;
553         $x->[$xi] += $BASE if ($bar = (($x->[$xi] -= $prd + $bar) < 0));
554         }
555       if ($x->[-1] < $car + $bar) 
556         {
557         $car = 0; --$q;
558         for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi) 
559           {
560           $x->[$xi] -= $BASE
561            if ($car = (($x->[$xi] += $y->[$yi] + $car) > $BASE));
562           }
563         }   
564       }
565       pop(@$x); unshift(@q, $q);
566     }
567   if (wantarray) 
568     {
569     @d = ();
570     if ($dd != 1)  
571       {
572       $car = 0; 
573       for $xi (reverse @$x) 
574         {
575         $prd = $car * $BASE + $xi;
576         $car = $prd - ($tmp = int($prd / $dd)) * $dd;
577         unshift(@d, $tmp);
578         }
579       }
580     else 
581       {
582       @d = @$x;
583       }
584     @$x = @q;
585     __strip_zeros($x); 
586     __strip_zeros(\@d);
587     return ($x,\@d);
588     }
589   @$x = @q;
590   __strip_zeros($x); 
591   }
592
593 ##############################################################################
594 # testing
595
596 sub _acmp
597   {
598   # internal absolute post-normalized compare (ignore signs)
599   # ref to array, ref to array, return <0, 0, >0
600   # arrays must have at least one entry; this is not checked for
601
602   my ($c,$cx,$cy) = @_;
603
604   # fat comp based on array elements
605   my $lxy = scalar @$cx - scalar @$cy;
606   return -1 if $lxy < 0;                                # already differs, ret
607   return 1 if $lxy > 0;                                 # ditto
608   
609   # now calculate length based on digits, not parts
610   $lxy = _len($c,$cx) - _len($c,$cy);                   # difference
611   return -1 if $lxy < 0;
612   return 1 if $lxy > 0;
613
614   # hm, same lengths,  but same contents?
615   my $i = 0; my $a;
616   # first way takes 5.49 sec instead of 4.87, but has the early out advantage
617   # so grep is slightly faster, but more inflexible. hm. $_ instead of $k
618   # yields 5.6 instead of 5.5 sec huh?
619   # manual way (abort if unequal, good for early ne)
620   my $j = scalar @$cx - 1;
621   while ($j >= 0)
622    {
623    last if ($a = $cx->[$j] - $cy->[$j]); $j--;
624    }
625   return 1 if $a > 0;
626   return -1 if $a < 0;
627   return 0;                                     # equal
628   # while it early aborts, it is even slower than the manual variant
629   #grep { return $a if ($a = $_ - $cy->[$i++]); } @$cx;
630   # grep way, go trough all (bad for early ne)
631   #grep { $a = $_ - $cy->[$i++]; } @$cx;
632   #return $a;
633   }
634
635 sub _len
636   {
637   # compute number of digits in bigint, minus the sign
638
639   # int() because add/sub sometimes leaves strings (like '00005') instead of
640   # '5' in this place, thus causing length() to report wrong length
641   my $cx = $_[1];
642
643   return (@$cx-1)*$BASE_LEN+length(int($cx->[-1]));
644   }
645
646 sub _digit
647   {
648   # return the nth digit, negative values count backward
649   # zero is rightmost, so _digit(123,0) will give 3
650   my ($c,$x,$n) = @_;
651
652   my $len = _len('',$x);
653
654   $n = $len+$n if $n < 0;               # -1 last, -2 second-to-last
655   $n = abs($n);                         # if negative was too big
656   $len--; $n = $len if $n > $len;       # n to big?
657   
658   my $elem = int($n / $BASE_LEN);       # which array element
659   my $digit = $n % $BASE_LEN;           # which digit in this element
660   $elem = '0000'.@$x[$elem];            # get element padded with 0's
661   return substr($elem,-$digit-1,1);
662   }
663
664 sub _zeros
665   {
666   # return amount of trailing zeros in decimal
667   # check each array elem in _m for having 0 at end as long as elem == 0
668   # Upon finding a elem != 0, stop
669   my $x = $_[1];
670   my $zeros = 0; my $elem;
671   foreach my $e (@$x)
672     {
673     if ($e != 0)
674       {
675       $elem = "$e";                             # preserve x
676       $elem =~ s/.*?(0*$)/$1/;                  # strip anything not zero
677       $zeros *= $BASE_LEN;                      # elems * 5
678       $zeros += CORE::length($elem);            # count trailing zeros
679       last;                                     # early out
680       }
681     $zeros ++;                                  # real else branch: 50% slower!
682     }
683   return $zeros;
684   }
685
686 ##############################################################################
687 # _is_* routines
688
689 sub _is_zero
690   {
691   # return true if arg (BINT or num_str) is zero (array '+', '0')
692   my $x = $_[1];
693   return (((scalar @$x == 1) && ($x->[0] == 0))) <=> 0;
694   }
695
696 sub _is_even
697   {
698   # return true if arg (BINT or num_str) is even
699   my $x = $_[1];
700   return (!($x->[0] & 1)) <=> 0; 
701   }
702
703 sub _is_odd
704   {
705   # return true if arg (BINT or num_str) is even
706   my $x = $_[1];
707   return (($x->[0] & 1)) <=> 0; 
708   }
709
710 sub _is_one
711   {
712   # return true if arg (BINT or num_str) is one (array '+', '1')
713   my $x = $_[1];
714   return (scalar @$x == 1) && ($x->[0] == 1) <=> 0; 
715   }
716
717 sub __strip_zeros
718   {
719   # internal normalization function that strips leading zeros from the array
720   # args: ref to array
721   my $s = shift;
722  
723   my $cnt = scalar @$s; # get count of parts
724   my $i = $cnt-1;
725   push @$s,0 if $i < 0;         # div might return empty results, so fix it
726
727   #print "strip: cnt $cnt i $i\n";
728   # '0', '3', '4', '0', '0',
729   #  0    1    2    3    4
730   # cnt = 5, i = 4
731   # i = 4
732   # i = 3
733   # => fcnt = cnt - i (5-2 => 3, cnt => 5-1 = 4, throw away from 4th pos)
734   # >= 1: skip first part (this can be zero)
735   while ($i > 0) { last if $s->[$i] != 0; $i--; }
736   $i++; splice @$s,$i if ($i < $cnt); # $i cant be 0
737   $s;                                                                    
738   }                                                                             
739
740 ###############################################################################
741 # check routine to test internal state of corruptions
742
743 sub _check
744   {
745   # used by the test suite
746   my $x = $_[1];
747
748   return "$x is not a reference" if !ref($x);
749
750   # are all parts are valid?
751   my $i = 0; my $j = scalar @$x; my ($e,$try);
752   while ($i < $j)
753     {
754     $e = $x->[$i]; $e = 'undef' unless defined $e;
755     $try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e)";
756     last if $e !~ /^[+]?[0-9]+$/;
757     $try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e) (stringify)";
758     last if "$e" !~ /^[+]?[0-9]+$/;
759     $try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e) (cat-stringify)";
760     last if '' . "$e" !~ /^[+]?[0-9]+$/;
761     $try = ' < 0 || >= $BASE; '."($x, $e)";
762     last if $e <0 || $e >= $BASE;
763     # this test is disabled, since new/bnorm and certain ops (like early out
764     # in add/sub) are allowed/expected to leave '00000' in some elements
765     #$try = '=~ /^00+/; '."($x, $e)";
766     #last if $e =~ /^00+/;
767     $i++;
768     }
769   return "Illegal part '$e' at pos $i (tested: $try)" if $i < $j;
770   return 0;
771   }
772
773
774 ###############################################################################
775 ###############################################################################
776 # some optional routines to make BigInt faster
777
778 sub _mod
779   {
780   # if possible, use mod shortcut
781   my ($c,$x,$yo) = @_;
782
783   # slow way since $y to big
784   if (scalar @$yo > 1)
785     {
786     my ($xo,$rem) = _div($c,$x,$yo);
787     return $rem;
788     }
789   my $y = $yo->[0];
790   # both are single element arrays
791   if (scalar @$x == 1)
792     {
793     $x->[0] %= $y;
794     return $x;
795     }
796
797   # @y is single element, but  @x has more than one
798   my $b = $BASE % $y;
799   if ($b == 0)
800     {
801     # when BASE % Y == 0 then (B * BASE) % Y == 0
802     # (B * BASE) % $y + A % Y => A % Y
803     # so need to consider only last element: O(1)
804     $x->[0] %= $y;
805     }
806   elsif ($b == 1)
807     {
808     # else need to go trough all elements: O(N),  but loop is a bit simplified
809     my $r = 0;
810     foreach (@$x)
811       {
812       $r += $_ % $y;
813       $r %= $y;
814       }
815     $r = 0 if $r == $y;
816     $x->[0] = $r;
817     }
818   else
819     {
820     # else need to go trough all elements: O(N)
821     my $r = 0; my $bm = 1;
822     foreach (@$x)
823       {
824       $r += ($_ % $y) * $bm;
825       $bm *= $b;
826       $bm %= $y;
827       $r %= $y;
828       }
829     $r = 0 if $r == $y;
830     $x->[0] = $r;
831     }
832   splice (@$x,1);
833   return $x;
834   }
835
836 ##############################################################################
837 # shifts
838
839 sub _rsft
840   {
841   my ($c,$x,$y,$n) = @_;
842
843   if ($n != 10)
844     {
845     return;     # we cant do this here, due to now _pow, so signal failure
846     }
847   else
848     {
849     # shortcut (faster) for shifting by 10)
850     # multiples of $BASE_LEN
851     my $dst = 0;                                # destination
852     my $src = _num($c,$y);                      # as normal int
853     my $rem = $src % $BASE_LEN;                 # remainder to shift
854     $src = int($src / $BASE_LEN);               # source
855     if ($rem == 0)
856       {
857       splice (@$x,0,$src);                      # even faster, 38.4 => 39.3
858       }
859     else
860       {
861       my $len = scalar @$x - $src;              # elems to go
862       my $vd; my $z = '0'x $BASE_LEN;
863       $x->[scalar @$x] = 0;                     # avoid || 0 test inside loop
864       while ($dst < $len)
865         {
866         $vd = $z.$x->[$src];
867         $vd = substr($vd,-$BASE_LEN,$BASE_LEN-$rem);
868         $src++;
869         $vd = substr($z.$x->[$src],-$rem,$rem) . $vd;
870         $vd = substr($vd,-$BASE_LEN,$BASE_LEN) if length($vd) > $BASE_LEN;
871         $x->[$dst] = int($vd);
872         $dst++;
873         }
874       splice (@$x,$dst) if $dst > 0;            # kill left-over array elems
875       pop @$x if $x->[-1] == 0;                 # kill last element if 0
876       } # else rem == 0
877     }
878   $x;
879   }
880
881 sub _lsft
882   {
883   my ($c,$x,$y,$n) = @_;
884
885   if ($n != 10)
886     {
887     return;     # we cant do this here, due to now _pow, so signal failure
888     }
889   else
890     {
891     # shortcut (faster) for shifting by 10) since we are in base 10eX
892     # multiples of $BASE_LEN:
893     my $src = scalar @$x;                       # source
894     my $len = _num($c,$y);                      # shift-len as normal int
895     my $rem = $len % $BASE_LEN;                 # remainder to shift
896     my $dst = $src + int($len/$BASE_LEN);       # destination
897     my $vd;                                     # further speedup
898     $x->[$src] = 0;                             # avoid first ||0 for speed
899     my $z = '0' x $BASE_LEN;
900     while ($src >= 0)
901       {
902       $vd = $x->[$src]; $vd = $z.$vd;
903       $vd = substr($vd,-$BASE_LEN+$rem,$BASE_LEN-$rem);
904       $vd .= $src > 0 ? substr($z.$x->[$src-1],-$BASE_LEN,$rem) : '0' x $rem;
905       $vd = substr($vd,-$BASE_LEN,$BASE_LEN) if length($vd) > $BASE_LEN;
906       $x->[$dst] = int($vd);
907       $dst--; $src--;
908       }
909     # set lowest parts to 0
910     while ($dst >= 0) { $x->[$dst--] = 0; }
911     # fix spurios last zero element
912     splice @$x,-1 if $x->[-1] == 0;
913     }
914   $x;
915   }
916
917 sub _pow
918   {
919   # power of $x to $y
920   # ref to array, ref to array, return ref to array
921   my ($c,$cx,$cy) = @_;
922
923   my $pow2 = _one();
924   my $two = _two();
925   my $y1 = _copy($c,$cy);
926   while (!_is_one($c,$y1))
927     {
928     _mul($c,$pow2,$cx) if _is_odd($c,$y1);
929     _div($c,$y1,$two);
930     _mul($c,$cx,$cx);
931     }
932   _mul($c,$cx,$pow2) unless _is_one($c,$pow2);
933   return $cx;
934   }
935
936 sub _sqrt
937   {
938   # square-root of $x
939   # ref to array, return ref to array
940   my ($c,$x) = @_;
941
942   if (scalar @$x == 1)
943     {
944     # fit's into one Perl scalar
945     $x->[0] = int(sqrt($x->[0]));
946     return $x;
947     } 
948   my $y = _copy($c,$x);
949   my $l = [ _len($c,$x) / 2 ];
950
951   splice @$x,0; $x->[0] = 1;    # keep ref($x), but modify it
952
953   _lsft($c,$x,$l,10);
954
955   my $two = _two();
956   my $last = _zero();
957   my $lastlast = _zero();
958   while (_acmp($c,$last,$x) != 0 && _acmp($c,$lastlast,$x) != 0)
959     {
960     $lastlast = _copy($c,$last);
961     $last = _copy($c,$x);
962     _add($c,$x, _div($c,_copy($c,$y),$x));
963     _div($c,$x, $two );
964     }
965   _dec($c,$x) if _acmp($c,$y,_mul($c,_copy($c,$x),$x)) < 0;     # overshot? 
966   $x;
967   }
968
969 ##############################################################################
970 # binary stuff
971
972 sub _and
973   {
974   my ($c,$x,$y) = @_;
975
976   # the shortcut makes equal, large numbers _really_ fast, and makes only a
977   # very small performance drop for small numbers (e.g. something with less
978   # than 32 bit) Since we optimize for large numbers, this is enabled.
979   return $x if _acmp($c,$x,$y) == 0;            # shortcut
980   
981   my $m = _one(); my ($xr,$yr);
982   my $mask = $AND_MASK;
983
984   my $x1 = $x;
985   my $y1 = _copy($c,$y);                        # make copy
986   $x = _zero();
987   my ($b,$xrr,$yrr);
988   use integer;
989   while (!_is_zero($c,$x1) && !_is_zero($c,$y1))
990     {
991     ($x1, $xr) = _div($c,$x1,$mask);
992     ($y1, $yr) = _div($c,$y1,$mask);
993
994     # make ints() from $xr, $yr
995     # this is when the AND_BITS are greater tahn $BASE and is slower for
996     # small (<256 bits) numbers, but faster for large numbers. Disabled
997     # due to KISS principle
998
999 #    $b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
1000 #    $b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
1001 #    _add($c,$x, _mul($c, _new( $c, \($xrr & $yrr) ), $m) );
1002     
1003     _add($c,$x, _mul($c, [ $xr->[0] & $yr->[0] ], $m) );
1004     _mul($c,$m,$mask);
1005     }
1006   $x;
1007   }
1008
1009 sub _xor
1010   {
1011   my ($c,$x,$y) = @_;
1012
1013   return _zero() if _acmp($c,$x,$y) == 0;       # shortcut (see -and)
1014
1015   my $m = _one(); my ($xr,$yr);
1016   my $mask = $XOR_MASK;
1017
1018   my $x1 = $x;
1019   my $y1 = _copy($c,$y);                        # make copy
1020   $x = _zero();
1021   my ($b,$xrr,$yrr);
1022   use integer;
1023   while (!_is_zero($c,$x1) && !_is_zero($c,$y1))
1024     {
1025     ($x1, $xr) = _div($c,$x1,$mask);
1026     ($y1, $yr) = _div($c,$y1,$mask);
1027     # make ints() from $xr, $yr (see _and())
1028     #$b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
1029     #$b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
1030     #_add($c,$x, _mul($c, _new( $c, \($xrr ^ $yrr) ), $m) );
1031     
1032     _add($c,$x, _mul($c, [ $xr->[0] ^ $yr->[0] ], $m) );
1033     _mul($c,$m,$mask);
1034     }
1035   # the loop stops when the shorter of the two numbers is exhausted
1036   # the remainder of the longer one will survive bit-by-bit, so we simple
1037   # multiply-add it in
1038   _add($c,$x, _mul($c, $x1, $m) ) if !_is_zero($c,$x1);
1039   _add($c,$x, _mul($c, $y1, $m) ) if !_is_zero($c,$y1);
1040   
1041   $x;
1042   }
1043
1044 sub _or
1045   {
1046   my ($c,$x,$y) = @_;
1047
1048   return $x if _acmp($c,$x,$y) == 0;            # shortcut (see _and)
1049
1050   my $m = _one(); my ($xr,$yr);
1051   my $mask = $OR_MASK;
1052
1053   my $x1 = $x;
1054   my $y1 = _copy($c,$y);                        # make copy
1055   $x = _zero();
1056   my ($b,$xrr,$yrr);
1057   use integer;
1058   while (!_is_zero($c,$x1) && !_is_zero($c,$y1))
1059     {
1060     ($x1, $xr) = _div($c,$x1,$mask);
1061     ($y1, $yr) = _div($c,$y1,$mask);
1062     # make ints() from $xr, $yr (see _and())
1063 #    $b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
1064 #    $b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
1065 #    _add($c,$x, _mul($c, _new( $c, \($xrr | $yrr) ), $m) );
1066     
1067     _add($c,$x, _mul($c, [ $xr->[0] | $yr->[0] ], $m) );
1068     _mul($c,$m,$mask);
1069     }
1070   # the loop stops when the shorter of the two numbers is exhausted
1071   # the remainder of the longer one will survive bit-by-bit, so we simple
1072   # multiply-add it in
1073   _add($c,$x, _mul($c, $x1, $m) ) if !_is_zero($c,$x1);
1074   _add($c,$x, _mul($c, $y1, $m) ) if !_is_zero($c,$y1);
1075   
1076   $x;
1077   }
1078
1079 sub _from_hex
1080   {
1081   # convert a hex number to decimal (ref to string, return ref to array)
1082   my ($c,$hs) = @_;
1083
1084   my $mul = _one();
1085   my $m = [ 0x10000 ];                          # 16 bit at a time
1086   my $x = _zero();
1087
1088   my $len = CORE::length($$hs)-2;
1089   $len = int($len/4);                           # 4-digit parts, w/o '0x'
1090   my $val; my $i = -4;
1091   while ($len >= 0)
1092     {
1093     $val = substr($$hs,$i,4);
1094     $val =~ s/^[+-]?0x// if $len == 0;          # for last part only because
1095     $val = hex($val);                           # hex does not like wrong chars
1096     $i -= 4; $len --;
1097     _add ($c, $x, _mul ($c, [ $val ], $mul ) ) if $val != 0;
1098     _mul ($c, $mul, $m ) if $len >= 0;          # skip last mul
1099     }
1100   $x;
1101   }
1102
1103 sub _from_bin
1104   {
1105   # convert a hex number to decimal (ref to string, return ref to array)
1106   my ($c,$bs) = @_;
1107
1108   my $mul = _one();
1109   my $m = [ 0x100 ];                            # 8 bit at a time
1110   my $x = _zero();
1111
1112   my $len = CORE::length($$bs)-2;
1113   $len = int($len/8);                           # 4-digit parts, w/o '0x'
1114   my $val; my $i = -8;
1115   while ($len >= 0)
1116     {
1117     $val = substr($$bs,$i,8);
1118     $val =~ s/^[+-]?0b// if $len == 0;          # for last part only
1119
1120     #$val = oct('0b'.$val);   # does not work on Perl prior to 5.6.0
1121     # $val = ('0' x (8-CORE::length($val))).$val if CORE::length($val) < 8;
1122     $val = ord(pack('B8',substr('00000000'.$val,-8,8))); 
1123
1124     $i -= 8; $len --;
1125     _add ($c, $x, _mul ($c, [ $val ], $mul ) ) if $val != 0;
1126     _mul ($c, $mul, $m ) if $len >= 0;          # skip last mul
1127     }
1128   $x;
1129   }
1130
1131 ##############################################################################
1132 ##############################################################################
1133
1134 1;
1135 __END__
1136
1137 =head1 NAME
1138
1139 Math::BigInt::Calc - Pure Perl module to support Math::BigInt
1140
1141 =head1 SYNOPSIS
1142
1143 Provides support for big integer calculations. Not intended to be used by other
1144 modules (except Math::BigInt::Cached). Other modules which sport the same
1145 functions can also be used to support Math::Bigint, like Math::BigInt::Pari.
1146
1147 =head1 DESCRIPTION
1148
1149 In order to allow for multiple big integer libraries, Math::BigInt was
1150 rewritten to use library modules for core math routines. Any module which
1151 follows the same API as this can be used instead by using the following:
1152
1153         use Math::BigInt lib => 'libname';
1154
1155 'libname' is either the long name ('Math::BigInt::Pari'), or only the short
1156 version like 'Pari'.
1157
1158 =head1 EXPORT
1159
1160 The following functions MUST be defined in order to support the use by
1161 Math::BigInt:
1162
1163         _new(string)    return ref to new object from ref to decimal string
1164         _zero()         return a new object with value 0
1165         _one()          return a new object with value 1
1166
1167         _str(obj)       return ref to a string representing the object
1168         _num(obj)       returns a Perl integer/floating point number
1169                         NOTE: because of Perl numeric notation defaults,
1170                         the _num'ified obj may lose accuracy due to 
1171                         machine-dependend floating point size limitations
1172                     
1173         _add(obj,obj)   Simple addition of two objects
1174         _mul(obj,obj)   Multiplication of two objects
1175         _div(obj,obj)   Division of the 1st object by the 2nd
1176                         In list context, returns (result,remainder).
1177                         NOTE: this is integer math, so no
1178                         fractional part will be returned.
1179         _sub(obj,obj)   Simple subtraction of 1 object from another
1180                         a third, optional parameter indicates that the params
1181                         are swapped. In this case, the first param needs to
1182                         be preserved, while you can destroy the second.
1183                         sub (x,y,1) => return x - y and keep x intact!
1184         _dec(obj)       decrement object by one (input is garant. to be > 0)
1185         _inc(obj)       increment object by one
1186
1187
1188         _acmp(obj,obj)  <=> operator for objects (return -1, 0 or 1)
1189
1190         _len(obj)       returns count of the decimal digits of the object
1191         _digit(obj,n)   returns the n'th decimal digit of object
1192
1193         _is_one(obj)    return true if argument is +1
1194         _is_zero(obj)   return true if argument is 0
1195         _is_even(obj)   return true if argument is even (0,2,4,6..)
1196         _is_odd(obj)    return true if argument is odd (1,3,5,7..)
1197
1198         _copy           return a ref to a true copy of the object
1199
1200         _check(obj)     check whether internal representation is still intact
1201                         return 0 for ok, otherwise error message as string
1202
1203 The following functions are optional, and can be defined if the underlying lib
1204 has a fast way to do them. If undefined, Math::BigInt will use pure Perl (hence
1205 slow) fallback routines to emulate these:
1206
1207         _from_hex(str)  return ref to new object from ref to hexadecimal string
1208         _from_bin(str)  return ref to new object from ref to binary string
1209         
1210         _as_hex(str)    return ref to scalar string containing the value as
1211                         unsigned hex string, with the '0x' prepended.
1212                         Leading zeros must be stripped.
1213         _as_bin(str)    Like as_hex, only as binary string containing only
1214                         zeros and ones. Leading zeros must be stripped and a
1215                         '0b' must be prepended.
1216         
1217         _rsft(obj,N,B)  shift object in base B by N 'digits' right
1218                         For unsupported bases B, return undef to signal failure
1219         _lsft(obj,N,B)  shift object in base B by N 'digits' left
1220                         For unsupported bases B, return undef to signal failure
1221         
1222         _xor(obj1,obj2) XOR (bit-wise) object 1 with object 2
1223                         Note: XOR, AND and OR pad with zeros if size mismatches
1224         _and(obj1,obj2) AND (bit-wise) object 1 with object 2
1225         _or(obj1,obj2)  OR (bit-wise) object 1 with object 2
1226
1227         _mod(obj,obj)   Return remainder of div of the 1st by the 2nd object
1228         _sqrt(obj)      return the square root of object (truncate to int)
1229         _pow(obj,obj)   return object 1 to the power of object 2
1230         _gcd(obj,obj)   return Greatest Common Divisor of two objects
1231         
1232         _zeros(obj)     return number of trailing decimal zeros
1233
1234 Input strings come in as unsigned but with prefix (i.e. as '123', '0xabc'
1235 or '0b1101').
1236
1237 Testing of input parameter validity is done by the caller, so you need not
1238 worry about underflow (f.i. in C<_sub()>, C<_dec()>) nor about division by
1239 zero or similar cases.
1240
1241 The first parameter can be modified, that includes the possibility that you
1242 return a reference to a completely different object instead. Although keeping
1243 the reference and just changing it's contents is prefered over creating and
1244 returning a different reference.
1245
1246 Return values are always references to objects or strings. Exceptions are
1247 C<_lsft()> and C<_rsft()>, which return undef if they can not shift the
1248 argument. This is used to delegate shifting of bases different than the one
1249 you can support back to Math::BigInt, which will use some generic code to
1250 calculate the result.
1251
1252 =head1 WRAP YOUR OWN
1253
1254 If you want to port your own favourite c-lib for big numbers to the
1255 Math::BigInt interface, you can take any of the already existing modules as
1256 a rough guideline. You should really wrap up the latest BigInt and BigFloat
1257 testsuites with your module, and replace in them any of the following:
1258
1259         use Math::BigInt;
1260
1261 by this:
1262
1263         use Math::BigInt lib => 'yourlib';
1264
1265 This way you ensure that your library really works 100% within Math::BigInt.
1266
1267 =head1 LICENSE
1268  
1269 This program is free software; you may redistribute it and/or modify it under
1270 the same terms as Perl itself. 
1271
1272 =head1 AUTHORS
1273
1274 Original math code by Mark Biggar, rewritten by Tels L<http://bloodgate.com/>
1275 in late 2000, 2001.
1276 Seperated from BigInt and shaped API with the help of John Peacock.
1277
1278 =head1 SEE ALSO
1279
1280 L<Math::BigInt>, L<Math::BigFloat>, L<Math::BigInt::BitVect>,
1281 L<Math::BigInt::GMP>, L<Math::BigInt::Cached> and L<Math::BigInt::Pari>.
1282
1283 =cut