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