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