[p5sagit/p5-mst-13.2.git] / lib / Math / BigInt / Calc.pm

package Math::BigInt::Calc;

use 5.005;
use strict;
use warnings;

require Exporter;

use vars qw/ @ISA @EXPORT $VERSION/;
@ISA = qw(Exporter);

@EXPORT = qw(
	_add _mul _div _mod _sub
	_new
	_str _num _acmp _len
	_digit
	_is_zero _is_one
	_is_even _is_odd
	_check _zero _one _copy _zeros
);
$VERSION = '0.06';

# Package to store unsigned big integers in decimal and do math with them

# Internally the numbers are stored in an array with at least 1 element, no
# leading zero parts (except the first) and in base 100000 

# todo:
# - fully remove funky $# stuff (maybe)
# - use integer; vs 1e7 as base 

# USE_MUL: due to problems on certain os (os390, posix-bc) "* 1e-5" is used
# instead of "/ 1e5" at some places, (marked with USE_MUL). But instead of
# using the reverse only on problematic machines, I used it everytime to avoid
# the costly comparisations. This _should_ work everywhere. Thanx Peter Prymmer

##############################################################################
# global constants, flags and accessory
 
# constants for easier life
my $nan = 'NaN';
my $BASE_LEN = 5;
my $BASE = int("1e".$BASE_LEN);		# var for trying to change it to 1e7
my $RBASE = 1e-5;           		# see USE_MUL
my $class = 'Math::BigInt::Calc';

##############################################################################
# create objects from various representations

sub _new
  {
  # (string) return ref to num_array
  # Convert a number from string format to internal base 100000 format.
  # Assumes normalized value as input.
  shift @_ if $_[0] eq $class;
  my $d = shift;
  # print "_new $d $$d\n";
  my $il = CORE::length($$d)-1;
  # these leaves '00000' instead of int 0 and will be corrected after any op
  return [ reverse(unpack("a" . ($il%5+1) . ("a5" x ($il/5)), $$d)) ];
  }                                                                             

sub _zero
  {
  # create a zero
  return [ 0 ];
  }

sub _one
  {
  # create a one
  return [ 1 ];
  }

sub _copy
  {
  shift @_ if $_[0] eq $class;
  my $x = shift;
  return [ @$x ];
  }

##############################################################################
# convert back to string and number

sub _str
  {
  # (ref to BINT) return num_str
  # Convert number from internal base 100000 format to string format.
  # internal format is always normalized (no leading zeros, "-0" => "+0")
  shift @_ if $_[0] eq $class;
  my $ar = shift;
  my $ret = "";
  my $l = scalar @$ar;         # number of parts
  return $nan if $l < 1;       # should not happen
  # handle first one different to strip leading zeros from it (there are no
  # leading zero parts in internal representation)
  $l --; $ret .= $ar->[$l]; $l--;
  # Interestingly, the pre-padd method uses more time
  # the old grep variant takes longer (14 to 10 sec)                            
  while ($l >= 0)
    {
    $ret .= substr('0000'.$ar->[$l],-5);   # fastest way I could think of
    $l--;
    }
  return \$ret;
  }                                                                             

sub _num
  {
  # Make a number (scalar int/float) from a BigInt object
  shift @_ if $_[0] eq $class;
  my $x = shift;
  return $x->[0] if scalar @$x == 1;  # below $BASE
  my $fac = 1;
  my $num = 0;
  foreach (@$x)
    {
    $num += $fac*$_; $fac *= $BASE;
    }
  return $num; 
  }

##############################################################################
# actual math code

sub _add
  {
  # (ref to int_num_array, ref to int_num_array)
  # routine to add two base 1e5 numbers
  # stolen from Knuth Vol 2 Algorithm A pg 231
  # there are separate routines to add and sub as per Kunth pg 233
  # This routine clobbers up array x, but not y.
 
  shift @_ if $_[0] eq $class;
  my ($x,$y) = @_;
 
  # for each in Y, add Y to X and carry. If after that, something is left in
  # X, foreach in X add carry to X and then return X, carry
  # Trades one "$j++" for having to shift arrays, $j could be made integer
  # but this would impose a limit to number-length to 2**32.
  my $i; my $car = 0; my $j = 0;
  for $i (@$y)
    {
    $x->[$j] -= $BASE
      if $car = (($x->[$j] += $i + $car) >= $BASE) ? 1 : 0;
    $j++;
    }
  while ($car != 0)
    {
    $x->[$j] -= $BASE if $car = (($x->[$j] += $car) >= $BASE) ? 1 : 0; $j++;
    }
    return $x;
  }                                                                             

sub _sub
  {
  # (ref to int_num_array, ref to int_num_array)
  # subtract base 1e5 numbers -- stolen from Knuth Vol 2 pg 232, $x > $y
  # subtract Y from X (X is always greater/equal!) by modifiyng x in place
  shift @_ if $_[0] eq $class;
  my ($sx,$sy,$s) = @_;
 
  my $car = 0; my $i; my $j = 0;
  if (!$s)
    {
    #print "case 2\n";
    for $i (@$sx)
      {
      last unless defined $sy->[$j] || $car;
      #print "x: $i y: $sy->[$j] c: $car\n";
      $i += $BASE if $car = (($i -= ($sy->[$j] || 0) + $car) < 0); $j++;
      #print "x: $i y: $sy->[$j-1] c: $car\n";
      }
    # might leave leading zeros, so fix that
    __strip_zeros($sx);
    return $sx;                                                                 
    }
  else
    {
    #print "case 1 (swap)\n";
    for $i (@$sx)
      {
      last unless defined $sy->[$j] || $car;
      #print "$sy->[$j] $i $car => $sx->[$j]\n";
      $sy->[$j] += $BASE
       if $car = (($sy->[$j] = $i-($sy->[$j]||0) - $car) < 0);
      #print "$sy->[$j] $i $car => $sy->[$j]\n";
      $j++;
      }
    # might leave leading zeros, so fix that
    __strip_zeros($sy);
    return $sy;
    }
  }                                                                             

sub _mul
  {
  # (BINT, BINT) return nothing
  # multiply two numbers in internal representation
  # modifies first arg, second needs not be different from first
  shift @_ if $_[0] eq $class;
  my ($xv,$yv) = @_;
 
  my @prod = (); my ($prod,$car,$cty,$xi,$yi);
  # since multiplying $x with $x fails, make copy in this case
  $yv = [@$xv] if "$xv" eq "$yv";
    # looping through @$y if $xi == 0 is silly! optimize it!
  for $xi (@$xv)
    {
    $car = 0; $cty = 0;
    for $yi (@$yv)
      {
      $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
      $prod[$cty++] =
       $prod - ($car = int($prod * 1e-5)) * $BASE;  # see USE_MUL
      }
    $prod[$cty] += $car if $car; # need really to check for 0?
    $xi = shift @prod;
    }
#  for $xi (@$xv)
#    {
#    $car = 0; $cty = 0;
#    # looping through this if $xi == 0 is silly! optimize it!
#    if (($xi||0) != 0)
#      {
#      for $yi (@$yv)
#        {
#        $prod = $prod[$cty]; $prod += ($car + $xi * $yi);	# no ||0 here
#        $prod[$cty++] =
#         $prod - ($car = int($prod * 1e-5)) * $BASE;  # see USE_MUL
#        }
#      }
#    $prod[$cty] += $car if $car; # need really to check for 0?
#    $xi = shift @prod;
#    }
  push @$xv, @prod;
  __strip_zeros($xv);
  # normalize (handled last to save check for $y->is_zero()
  return $xv;
  }                                                                             

sub _div
  {
  # ref to array, ref to array, modify first array and return reminder if 
  # in list context
  # does no longer handle sign
  shift @_ if $_[0] eq $class;
  my ($x,$yorg) = @_;
  my ($car,$bar,$prd,$dd,$xi,$yi,@q,$v2,$v1);

  my (@d,$tmp,$q,$u2,$u1,$u0);

  $car = $bar = $prd = 0;
  
  my $y = [ @$yorg ];
  if (($dd = int($BASE/($y->[-1]+1))) != 1) 
    {
    for $xi (@$x) 
      {
      $xi = $xi * $dd + $car;
      $xi -= ($car = int($xi * $RBASE)) * $BASE;	# see USE_MUL
      }
    push(@$x, $car); $car = 0;
    for $yi (@$y) 
      {
      $yi = $yi * $dd + $car;
      $yi -= ($car = int($yi * $RBASE)) * $BASE;	# see USE_MUL
      }
    }
  else 
    {
    push(@$x, 0);
    }
  @q = (); ($v2,$v1) = @$y[-2,-1];
  $v2 = 0 unless $v2;
  while ($#$x > $#$y) 
    {
    ($u2,$u1,$u0) = @$x[-3..-1];
    $u2 = 0 unless $u2;
    #warn "oups v1 is 0, u0: $u0 $y->[-2] $y->[-1] l ",scalar @$y,"\n"
    # if $v1 == 0;
    $q = (($u0 == $v1) ? 99999 : int(($u0*$BASE+$u1)/$v1));
    --$q while ($v2*$q > ($u0*1e5+$u1-$q*$v1)*$BASE+$u2);
    if ($q)
      {
      ($car, $bar) = (0,0);
      for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi) 
        {
        $prd = $q * $y->[$yi] + $car;
        $prd -= ($car = int($prd * $RBASE)) * $BASE;	# see USE_MUL
	$x->[$xi] += 1e5 if ($bar = (($x->[$xi] -= $prd + $bar) < 0));
	}
      if ($x->[-1] < $car + $bar) 
        {
        $car = 0; --$q;
	for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi) 
          {
	  $x->[$xi] -= 1e5
	   if ($car = (($x->[$xi] += $y->[$yi] + $car) > $BASE));
	  }
	}   
      }
      pop(@$x); unshift(@q, $q);
    }
  if (wantarray) 
    {
    @d = ();
    if ($dd != 1)  
      {
      $car = 0; 
      for $xi (reverse @$x) 
        {
        $prd = $car * $BASE + $xi;
        $car = $prd - ($tmp = int($prd / $dd)) * $dd; # see USE_MUL
        unshift(@d, $tmp);
        }
      }
    else 
      {
      @d = @$x;
      }
    @$x = @q;
    __strip_zeros($x); 
    __strip_zeros(\@d);
    return ($x,\@d);
    }
  @$x = @q;
  __strip_zeros($x); 
  return $x;
  }

##############################################################################
# testing

sub _acmp
  {
  # internal absolute post-normalized compare (ignore signs)
  # ref to array, ref to array, return <0, 0, >0
  # arrays must have at least on entry, this is not checked for

  shift @_ if $_[0] eq $class;
  my ($cx, $cy) = @_;

  #print "$cx $cy\n"; 
  my ($i,$a,$x,$y,$k);
  # calculate length based on digits, not parts
  $x = _len($cx); $y = _len($cy);
  # print "length: ",($x-$y),"\n";
  return $x-$y if ($x - $y);              # if different in length
  #print "full compare\n";
  $i = 0; $a = 0;
  # first way takes 5.49 sec instead of 4.87, but has the early out advantage
  # so grep is slightly faster, but more unflexible. hm. $_ instead if $k
  # yields 5.6 instead of 5.5 sec huh?
  # manual way (abort if unequal, good for early ne)
  my $j = scalar @$cx - 1;
  while ($j >= 0)
   {
   # print "$cx->[$j] $cy->[$j] $a",$cx->[$j]-$cy->[$j],"\n";
   last if ($a = $cx->[$j] - $cy->[$j]); $j--;
   }
  return $a;
  # while it early aborts, it is even slower than the manual variant
  #grep { return $a if ($a = $_ - $cy->[$i++]); } @$cx;
  # grep way, go trough all (bad for early ne)
  #grep { $a = $_ - $cy->[$i++]; } @$cx;
  #return $a;
  }

sub _len
  {
  # computer number of digits in bigint, minus the sign
  # int() because add/sub leaves sometimes strings (like '00005') instead of
  # int ('5') in this place, causing length to fail
  shift @_ if $_[0] eq $class;
  my $cx = shift;

  return (@$cx-1)*5+length(int($cx->[-1]));
  }

sub _digit
  {
  # return the nth digit, negative values count backward
  # zero is rightmost, so _digit(123,0) will give 3
  shift @_ if $_[0] eq $class;
  my $x = shift;
  my $n = shift || 0; 

  my $len = _len($x);

  $n = $len+$n if $n < 0;		# -1 last, -2 second-to-last
  $n = abs($n);				# if negative was too big
  $len--; $n = $len if $n > $len;	# n to big?
  
  my $elem = int($n / 5);		# which array element
  my $digit = $n % 5;			# which digit in this element
  $elem = '0000'.@$x[$elem];		# get element padded with 0's
  return substr($elem,-$digit-1,1);
  }

sub _zeros
  {
  # return amount of trailing zeros in decimal
  # check each array elem in _m for having 0 at end as long as elem == 0
  # Upon finding a elem != 0, stop
  shift @_ if $_[0] eq $class;
  my $x = shift;
  my $zeros = 0; my $elem;
  foreach my $e (@$x)
    {
    if ($e != 0)
      {
      $elem = "$e";                            # preserve x
      $elem =~ s/.*?(0*$)/$1/;                 # strip anything not zero
      $zeros *= 5;                             # elems * 5
      $zeros += CORE::length($elem);           # count trailing zeros
      last;                                    # early out
      }
    $zeros ++;                                 # real else branch: 50% slower!
    }
  return $zeros;
  }

##############################################################################
# _is_* routines

sub _is_zero
  {
  # return true if arg (BINT or num_str) is zero (array '+', '0')
  shift @_ if $_[0] eq $class;
  my ($x) = shift;
  return (((scalar @$x == 1) && ($x->[0] == 0))) <=> 0;
  }

sub _is_even
  {
  # return true if arg (BINT or num_str) is even
  shift @_ if $_[0] eq $class;
  my ($x) = shift;
  return (!($x->[0] & 1)) <=> 0; 
  }

sub _is_odd
  {
  # return true if arg (BINT or num_str) is even
  shift @_ if $_[0] eq $class;
  my ($x) = shift;
  return (($x->[0] & 1)) <=> 0; 
  }

sub _is_one
  {
  # return true if arg (BINT or num_str) is one (array '+', '1')
  shift @_ if $_[0] eq $class;
  my ($x) = shift;
  return (scalar @$x == 1) && ($x->[0] == 1) <=> 0; 
  }

sub __strip_zeros
  {
  # internal normalization function that strips leading zeros from the array
  # args: ref to array
  #trace(@_);
  shift @_ if $_[0] eq $class;
  my $s = shift;
 
  my $cnt = scalar @$s; # get count of parts
  my $i = $cnt-1;
  #print "strip: cnt $cnt i $i\n";
  # '0', '3', '4', '0', '0',
  #  0    1    2    3    4
  # cnt = 5, i = 4
  # i = 4
  # i = 3
  # => fcnt = cnt - i (5-2 => 3, cnt => 5-1 = 4, throw away from 4th pos)
  # >= 1: skip first part (this can be zero)
  while ($i > 0) { last if $s->[$i] != 0; $i--; }
  $i++; splice @$s,$i if ($i < $cnt); # $i cant be 0
  return $s;                                                                    
  }                                                                             

###############################################################################
# check routine to test internal state of corruptions

sub _check
  {
  # no checks yet, pull it out from the test suite
  shift @_ if $_[0] eq $class;

  my ($x) = shift;
  return "$x is not a reference" if !ref($x);

  # are all parts are valid?
  my $i = 0; my $j = scalar @$x; my ($e,$try);
  while ($i < $j)
    {
    $e = $x->[$i]; $e = 'undef' unless defined $e;
    $try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e)";
    last if $e !~ /^[+]?[0-9]+$/;
    $try = ' < 0 || >= $BASE; '."($x, $e)";
    last if $e <0 || $e >= $BASE;
    # this test is disabled, since new/bnorm and certain ops (like early out
    # in add/sub) are allowed/expected to leave '00000' in some elements
    #$try = '=~ /^00+/; '."($x, $e)";
    #last if $e =~ /^00+/;
    $i++;
    }
  return "Illegal part '$e' at pos $i (tested: $try)" if $i < $j;
  return 0;
  }

1;
__END__

=head1 NAME

Math::BigInt::Calc - Pure Perl module to support Math::BigInt

=head1 SYNOPSIS

Provides support for big integer calculations. Not intended
to be used by other modules. Other modules which export the
same functions can also be used to support Math::Bigint

=head1 DESCRIPTION

In order to allow for multiple big integer libraries, Math::BigInt
was rewritten to use library modules for core math routines. Any
module which follows the same API as this can be used instead by
using the following call:

	use Math::BigInt Calc => BigNum;

=head1 EXPORT

The following functions MUST be exported in order to support
the use by Math::BigInt:

	_new(string)	return ref to new object from ref to decimal string
	_zero()		return a new object with value 0
	_one()		return a new object with value 1

	_str(obj)	return ref to a string representing the object
	_num(obj)	returns a Perl integer/floating point number
			NOTE: because of Perl numeric notation defaults,
			the _num'ified obj may lose accuracy due to 
			machine-dependend floating point size limitations
                    
	_add(obj,obj)	Simple addition of two objects
	_mul(obj,obj)	Multiplication of two objects
	_div(obj,obj)	Division of the 1st object by the 2nd
			In list context, returns (result,reminder).
			NOTE: this is integer math, so there will no
			fractional part be returned.
	_sub(obj,obj)	Simple substraction of 1 object from another
			a third, optional parameter indicates that the params
			are swapped. In this case, the first param needs to
			be preserved, while you can destroy the second.
			sub (x,y,1) => return x - y and keep x intact!

	_acmp(obj,obj)	<=> operator for objects (return -1, 0 or 1)

	_len(obj)	returns count of the decimal digits of the object
	_digit(obj,n)	returns the n'th decimal digit of object

	_is_one(obj)	return true if argument is +1
	_is_zero(obj)	return true if argument is 0
	_is_even(obj)	return true if argument is even (0,2,4,6..)
	_is_odd(obj)	return true if argument is odd (1,3,5,7..)

	_copy		return a ref to a true copy of the object

	_check(obj)	check whether internal representation is still intact
			return 0 for ok, otherwise error message as string

The following functions are optional, and can be exported if the underlying lib
has a fast way to do them. If not defined, Math::BigInt will use a pure, but
slow Perl function as fallback to emulate these:

	_from_hex(str)	return ref to new object from ref to hexadecimal string
	_from_bin(str)	return ref to new object from ref to binary string
	
	_rsft(obj,N,B)	shift object in base B by N 'digits' right
	_lsft(obj,N,B)	shift object in base B by N 'digits' left
	
	_xor(obj1,obj2)	XOR (bit-wise) object 1 with object 2
			Mote: XOR, AND and OR pad with zeros if size mismatches
	_and(obj1,obj2)	AND (bit-wise) object 1 with object 2
	_or(obj1,obj2)	OR (bit-wise) object 1 with object 2

	_sqrt(obj)	return the square root of object
	_pow(obj,obj)	return object 1 to the power of object 2
	_gcd(obj,obj)	return Greatest Common Divisor of two objects
	
	_zeros(obj)	return amount of trailing decimal zeros

	_dec(obj)	decrement object by one (input is >= 1)
	_inc(obj)	increment object by one

Input strings come in as unsigned but with prefix (aka as '123', '0xabc'
or '0b1101').

Testing of input parameter validity is done by the caller, so you need not to
worry about underflow (C<_sub()>, C<_dec()>) nor division by zero or similiar
cases.

=head1 LICENSE
 
This program is free software; you may redistribute it and/or modify it under
the same terms as Perl itself. 

=head1 AUTHORS

Original math code by Mark Biggar, rewritten by Tels L<http://bloodgate.com/>
in late 2000, 2001.
Seperated from BigInt and shaped API with the help of John Peacock.

=head1 SEE ALSO

L<Math::BigInt>, L<Math::BigFloat>.

=cut
Commit	Line	Data
0716bf9b	1	package Math::BigInt::Calc;
	2
	3	use 5.005;
	4	use strict;
	5	use warnings;
	6
	7	require Exporter;
	8
	9	use vars qw/ @ISA @EXPORT $VERSION/;
	10	@ISA = qw(Exporter);
	11
	12	@EXPORT = qw(
	13	_add _mul _div _mod _sub
	14	_new
	15	_str _num _acmp _len
	16	_digit
	17	_is_zero _is_one
	18	_is_even _is_odd
	19	_check _zero _one _copy _zeros
	20	);
	21	$VERSION = '0.06';
	22
	23	# Package to store unsigned big integers in decimal and do math with them
	24
	25	# Internally the numbers are stored in an array with at least 1 element, no
	26	# leading zero parts (except the first) and in base 100000
	27
	28	# todo:
	29	# - fully remove funky $# stuff (maybe)
	30	# - use integer; vs 1e7 as base
	31
	32	# USE_MUL: due to problems on certain os (os390, posix-bc) "* 1e-5" is used
	33	# instead of "/ 1e5" at some places, (marked with USE_MUL). But instead of
	34	# using the reverse only on problematic machines, I used it everytime to avoid
	35	# the costly comparisations. This _should_ work everywhere. Thanx Peter Prymmer
	36
	37	##############################################################################
	38	# global constants, flags and accessory
	39
	40	# constants for easier life
	41	my $nan = 'NaN';
	42	my $BASE_LEN = 5;
	43	my $BASE = int("1e".$BASE_LEN); # var for trying to change it to 1e7
	44	my $RBASE = 1e-5; # see USE_MUL
	45	my $class = 'Math::BigInt::Calc';
	46
	47	##############################################################################
	48	# create objects from various representations
	49
	50	sub _new
	51	{
	52	# (string) return ref to num_array
	53	# Convert a number from string format to internal base 100000 format.
	54	# Assumes normalized value as input.
	55	shift @_ if $_[0] eq $class;
	56	my $d = shift;
	57	# print "_new $d $$d\n";
	58	my $il = CORE::length($$d)-1;
	59	# these leaves '00000' instead of int 0 and will be corrected after any op
	60	return [ reverse(unpack("a" . ($il%5+1) . ("a5" x ($il/5)), $$d)) ];
	61	}
	62
	63	sub _zero
	64	{
65	# create a zero
66	return [ 0 ];
67	}
68
69	sub _one
70	{
71	# create a one
72	return [ 1 ];
73	}
74
75	sub _copy
76	{
77	shift @_ if $_[0] eq $class;
78	my $x = shift;
79	return [ @$x ];
80	}
81
82	##############################################################################
83	# convert back to string and number
84
85	sub _str
86	{
87	# (ref to BINT) return num_str
88	# Convert number from internal base 100000 format to string format.
89	# internal format is always normalized (no leading zeros, "-0" => "+0")
90	shift @_ if $_[0] eq $class;
91	my $ar = shift;
92	my $ret = "";
93	my $l = scalar @$ar; # number of parts
94	return $nan if $l < 1; # should not happen
95	# handle first one different to strip leading zeros from it (there are no
96	# leading zero parts in internal representation)
97	$l --; $ret .= $ar->[$l]; $l--;
98	# Interestingly, the pre-padd method uses more time
99	# the old grep variant takes longer (14 to 10 sec)
100	while ($l >= 0)
101	{
102	$ret .= substr('0000'.$ar->[$l],-5); # fastest way I could think of
103	$l--;
104	}
105	return \$ret;
106	}
107
108	sub _num
109	{
110	# Make a number (scalar int/float) from a BigInt object
111	shift @_ if $_[0] eq $class;
112	my $x = shift;
113	return $x->[0] if scalar @$x == 1; # below $BASE
114	my $fac = 1;
115	my $num = 0;
116	foreach (@$x)
117	{
118	$num += $fac$_; $fac = $BASE;
119	}
120	return $num;
121	}
122
123	##############################################################################
124	# actual math code
125
126	sub _add
127	{
128	# (ref to int_num_array, ref to int_num_array)
129	# routine to add two base 1e5 numbers
130	# stolen from Knuth Vol 2 Algorithm A pg 231
131	# there are separate routines to add and sub as per Kunth pg 233
132	# This routine clobbers up array x, but not y.
133
134	shift @_ if $_[0] eq $class;
135	my ($x,$y) = @_;
136
137	# for each in Y, add Y to X and carry. If after that, something is left in
138	# X, foreach in X add carry to X and then return X, carry
139	# Trades one "$j++" for having to shift arrays, $j could be made integer
140	# but this would impose a limit to number-length to 2**32.
141	my $i; my $car = 0; my $j = 0;
142	for $i (@$y)
143	{
144	$x->[$j] -= $BASE
145	if $car = (($x->[$j] += $i + $car) >= $BASE) ? 1 : 0;
146	$j++;
147	}
148	while ($car != 0)
149	{
150	$x->[$j] -= $BASE if $car = (($x->[$j] += $car) >= $BASE) ? 1 : 0; $j++;
151	}
152	return $x;
153	}
154
155	sub _sub
156	{
157	# (ref to int_num_array, ref to int_num_array)
158	# subtract base 1e5 numbers -- stolen from Knuth Vol 2 pg 232, $x > $y
159	# subtract Y from X (X is always greater/equal!) by modifiyng x in place
160	shift @_ if $_[0] eq $class;
161	my ($sx,$sy,$s) = @_;
162
163	my $car = 0; my $i; my $j = 0;
164	if (!$s)
165	{
166	#print "case 2\n";
167	for $i (@$sx)
168	{
169	last unless defined $sy->[$j] \|\| $car;
170	#print "x: $i y: $sy->[$j] c: $car\n";
171	$i += $BASE if $car = (($i -= ($sy->[$j] \|\| 0) + $car) < 0); $j++;
172	#print "x: $i y: $sy->[$j-1] c: $car\n";
173	}
174	# might leave leading zeros, so fix that
175	__strip_zeros($sx);
176	return $sx;
177	}
178	else
179	{
180	#print "case 1 (swap)\n";
181	for $i (@$sx)
182	{
183	last unless defined $sy->[$j] \|\| $car;
184	#print "$sy->[$j] $i $car => $sx->[$j]\n";
185	$sy->[$j] += $BASE
186	if $car = (($sy->[$j] = $i-($sy->[$j]\|\|0) - $car) < 0);
187	#print "$sy->[$j] $i $car => $sy->[$j]\n";
188	$j++;
189	}
190	# might leave leading zeros, so fix that
191	__strip_zeros($sy);
192	return $sy;
193	}
194	}
195
196	sub _mul
197	{
198	# (BINT, BINT) return nothing
199	# multiply two numbers in internal representation
200	# modifies first arg, second needs not be different from first
201	shift @_ if $_[0] eq $class;
202	my ($xv,$yv) = @_;
203
204	my @prod = (); my ($prod,$car,$cty,$xi,$yi);
205	# since multiplying $x with $x fails, make copy in this case
206	$yv = [@$xv] if "$xv" eq "$yv";
207	# looping through @$y if $xi == 0 is silly! optimize it!
208	for $xi (@$xv)
209	{
210	$car = 0; $cty = 0;
211	for $yi (@$yv)
212	{
213	$prod = $xi * $yi + ($prod[$cty] \|\| 0) + $car;
214	$prod[$cty++] =
215	$prod - ($car = int($prod * 1e-5)) * $BASE; # see USE_MUL
216	}
217	$prod[$cty] += $car if $car; # need really to check for 0?
218	$xi = shift @prod;
219	}
220	# for $xi (@$xv)
221	# {
222	# $car = 0; $cty = 0;
223	# # looping through this if $xi == 0 is silly! optimize it!
224	# if (($xi\|\|0) != 0)
225	# {
226	# for $yi (@$yv)
227	# {
228	# $prod = $prod[$cty]; $prod += ($car + $xi * $yi); # no \|\|0 here
229	# $prod[$cty++] =
230	# $prod - ($car = int($prod * 1e-5)) * $BASE; # see USE_MUL
231	# }
232	# }
233	# $prod[$cty] += $car if $car; # need really to check for 0?
234	# $xi = shift @prod;
235	# }
236	push @$xv, @prod;
237	__strip_zeros($xv);
238	# normalize (handled last to save check for $y->is_zero()
239	return $xv;
240	}
241
242	sub _div
243	{
244	# ref to array, ref to array, modify first array and return reminder if
245	# in list context
246	# does no longer handle sign
247	shift @_ if $_[0] eq $class;
248	my ($x,$yorg) = @_;
249	my ($car,$bar,$prd,$dd,$xi,$yi,@q,$v2,$v1);
250
251	my (@d,$tmp,$q,$u2,$u1,$u0);
252
253	$car = $bar = $prd = 0;
254
255	my $y = [ @$yorg ];
256	if (($dd = int($BASE/($y->[-1]+1))) != 1)
257	{
258	for $xi (@$x)
259	{
260	$xi = $xi * $dd + $car;
261	$xi -= ($car = int($xi * $RBASE)) * $BASE; # see USE_MUL
262	}
263	push(@$x, $car); $car = 0;
264	for $yi (@$y)
265	{
266	$yi = $yi * $dd + $car;
267	$yi -= ($car = int($yi * $RBASE)) * $BASE; # see USE_MUL
268	}
269	}
270	else
271	{
272	push(@$x, 0);
273	}
274	@q = (); ($v2,$v1) = @$y[-2,-1];
275	$v2 = 0 unless $v2;
276	while ($#$x > $#$y)
277	{
278	($u2,$u1,$u0) = @$x[-3..-1];
279	$u2 = 0 unless $u2;
280	#warn "oups v1 is 0, u0: $u0 $y->[-2] $y->[-1] l ",scalar @$y,"\n"
281	# if $v1 == 0;
282	$q = (($u0 == $v1) ? 99999 : int(($u0*$BASE+$u1)/$v1));
283	--$q while ($v2$q > ($u01e5+$u1-$q$v1)$BASE+$u2);
284	if ($q)
285	{
286	($car, $bar) = (0,0);
287	for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
288	{
289	$prd = $q * $y->[$yi] + $car;
290	$prd -= ($car = int($prd * $RBASE)) * $BASE; # see USE_MUL
291	$x->[$xi] += 1e5 if ($bar = (($x->[$xi] -= $prd + $bar) < 0));
292	}
293	if ($x->[-1] < $car + $bar)
294	{
295	$car = 0; --$q;
296	for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
297	{
298	$x->[$xi] -= 1e5
299	if ($car = (($x->[$xi] += $y->[$yi] + $car) > $BASE));
300	}
301	}
302	}
303	pop(@$x); unshift(@q, $q);
304	}
305	if (wantarray)
306	{
307	@d = ();
308	if ($dd != 1)
309	{
310	$car = 0;
311	for $xi (reverse @$x)
312	{
313	$prd = $car * $BASE + $xi;
314	$car = $prd - ($tmp = int($prd / $dd)) * $dd; # see USE_MUL
315	unshift(@d, $tmp);
316	}
317	}
318	else
319	{
320	@d = @$x;
321	}
322	@$x = @q;
323	__strip_zeros($x);
324	__strip_zeros(\@d);
325	return ($x,\@d);
326	}
327	@$x = @q;
328	__strip_zeros($x);
329	return $x;
330	}
331
332	##############################################################################
333	# testing
334
335	sub _acmp
336	{
337	# internal absolute post-normalized compare (ignore signs)
338	# ref to array, ref to array, return <0, 0, >0
339	# arrays must have at least on entry, this is not checked for
340
341	shift @_ if $_[0] eq $class;
342	my ($cx, $cy) = @_;
343
344	#print "$cx $cy\n";
345	my ($i,$a,$x,$y,$k);
346	# calculate length based on digits, not parts
347	$x = _len($cx); $y = _len($cy);
348	# print "length: ",($x-$y),"\n";
349	return $x-$y if ($x - $y); # if different in length
350	#print "full compare\n";
351	$i = 0; $a = 0;
352	# first way takes 5.49 sec instead of 4.87, but has the early out advantage
353	# so grep is slightly faster, but more unflexible. hm. $_ instead if $k
354	# yields 5.6 instead of 5.5 sec huh?
355	# manual way (abort if unequal, good for early ne)
356	my $j = scalar @$cx - 1;
357	while ($j >= 0)
358	{
359	# print "$cx->[$j] $cy->[$j] $a",$cx->[$j]-$cy->[$j],"\n";
360	last if ($a = $cx->[$j] - $cy->[$j]); $j--;
361	}
362	return $a;
363	# while it early aborts, it is even slower than the manual variant
364	#grep { return $a if ($a = $_ - $cy->[$i++]); } @$cx;
365	# grep way, go trough all (bad for early ne)
366	#grep { $a = $_ - $cy->[$i++]; } @$cx;
367	#return $a;
368	}
369
370	sub _len
371	{
372	# computer number of digits in bigint, minus the sign
373	# int() because add/sub leaves sometimes strings (like '00005') instead of
374	# int ('5') in this place, causing length to fail
375	shift @_ if $_[0] eq $class;
376	my $cx = shift;
377
378	return (@$cx-1)*5+length(int($cx->[-1]));
379	}
380
381	sub _digit
382	{
383	# return the nth digit, negative values count backward
384	# zero is rightmost, so _digit(123,0) will give 3
385	shift @_ if $_[0] eq $class;
386	my $x = shift;
387	my $n = shift \|\| 0;
388
389	my $len = _len($x);
390
391	$n = $len+$n if $n < 0; # -1 last, -2 second-to-last
392	$n = abs($n); # if negative was too big
393	$len--; $n = $len if $n > $len; # n to big?
394
395	my $elem = int($n / 5); # which array element
396	my $digit = $n % 5; # which digit in this element
397	$elem = '0000'.@$x[$elem]; # get element padded with 0's
398	return substr($elem,-$digit-1,1);
399	}
400
401	sub _zeros
402	{
403	# return amount of trailing zeros in decimal
404	# check each array elem in _m for having 0 at end as long as elem == 0
405	# Upon finding a elem != 0, stop
406	shift @_ if $_[0] eq $class;
407	my $x = shift;
408	my $zeros = 0; my $elem;
409	foreach my $e (@$x)
410	{
411	if ($e != 0)
412	{
413	$elem = "$e"; # preserve x
414	$elem =~ s/.?(0$)/$1/; # strip anything not zero
415	$zeros = 5; # elems 5
416	$zeros += CORE::length($elem); # count trailing zeros
417	last; # early out
418	}
419	$zeros ++; # real else branch: 50% slower!
420	}
421	return $zeros;
422	}
423
424	##############################################################################
425	# _is_* routines
426
427	sub _is_zero
428	{
429	# return true if arg (BINT or num_str) is zero (array '+', '0')
430	shift @_ if $_[0] eq $class;
431	my ($x) = shift;
432	return (((scalar @$x == 1) && ($x->[0] == 0))) <=> 0;
433	}
434
435	sub _is_even
436	{
437	# return true if arg (BINT or num_str) is even
438	shift @_ if $_[0] eq $class;
439	my ($x) = shift;
440	return (!($x->[0] & 1)) <=> 0;
441	}
442
443	sub _is_odd
444	{
445	# return true if arg (BINT or num_str) is even
446	shift @_ if $_[0] eq $class;
447	my ($x) = shift;
448	return (($x->[0] & 1)) <=> 0;
449	}
450
451	sub _is_one
452	{
453	# return true if arg (BINT or num_str) is one (array '+', '1')
454	shift @_ if $_[0] eq $class;
455	my ($x) = shift;
456	return (scalar @$x == 1) && ($x->[0] == 1) <=> 0;
457	}
458
459	sub __strip_zeros
460	{
461	# internal normalization function that strips leading zeros from the array
462	# args: ref to array
463	#trace(@_);
464	shift @_ if $_[0] eq $class;
465	my $s = shift;
466
467	my $cnt = scalar @$s; # get count of parts
468	my $i = $cnt-1;
469	#print "strip: cnt $cnt i $i\n";
470	# '0', '3', '4', '0', '0',
471	# 0 1 2 3 4
472	# cnt = 5, i = 4
473	# i = 4
474	# i = 3
475	# => fcnt = cnt - i (5-2 => 3, cnt => 5-1 = 4, throw away from 4th pos)
476	# >= 1: skip first part (this can be zero)
477	while ($i > 0) { last if $s->[$i] != 0; $i--; }
478	$i++; splice @$s,$i if ($i < $cnt); # $i cant be 0
479	return $s;
480	}
481
482	###############################################################################
483	# check routine to test internal state of corruptions
484
485	sub _check
486	{
487	# no checks yet, pull it out from the test suite
488	shift @_ if $_[0] eq $class;
489
490	my ($x) = shift;
491	return "$x is not a reference" if !ref($x);
492
493	# are all parts are valid?
494	my $i = 0; my $j = scalar @$x; my ($e,$try);
495	while ($i < $j)
496	{
497	$e = $x->[$i]; $e = 'undef' unless defined $e;
498	$try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e)";
499	last if $e !~ /^[+]?[0-9]+$/;
500	$try = ' < 0 \|\| >= $BASE; '."($x, $e)";
501	last if $e <0 \|\| $e >= $BASE;
502	# this test is disabled, since new/bnorm and certain ops (like early out
503	# in add/sub) are allowed/expected to leave '00000' in some elements
504	#$try = '=~ /^00+/; '."($x, $e)";
505	#last if $e =~ /^00+/;
506	$i++;
507	}
508	return "Illegal part '$e' at pos $i (tested: $try)" if $i < $j;
509	return 0;
510	}
511
512	1;
513	__END__
514
515	=head1 NAME
516
517	Math::BigInt::Calc - Pure Perl module to support Math::BigInt
518
519	=head1 SYNOPSIS
520
521	Provides support for big integer calculations. Not intended
522	to be used by other modules. Other modules which export the
523	same functions can also be used to support Math::Bigint
524
525	=head1 DESCRIPTION
526
527	In order to allow for multiple big integer libraries, Math::BigInt
528	was rewritten to use library modules for core math routines. Any
529	module which follows the same API as this can be used instead by
530	using the following call:
531
532	use Math::BigInt Calc => BigNum;
533
534	=head1 EXPORT
535
536	The following functions MUST be exported in order to support
537	the use by Math::BigInt:
538
539	_new(string) return ref to new object from ref to decimal string
540	_zero() return a new object with value 0
541	_one() return a new object with value 1
542
543	_str(obj) return ref to a string representing the object
544	_num(obj) returns a Perl integer/floating point number
545	NOTE: because of Perl numeric notation defaults,
546	the _num'ified obj may lose accuracy due to
547	machine-dependend floating point size limitations
548
549	_add(obj,obj) Simple addition of two objects
550	_mul(obj,obj) Multiplication of two objects
551	_div(obj,obj) Division of the 1st object by the 2nd
552	In list context, returns (result,reminder).
553	NOTE: this is integer math, so there will no
554	fractional part be returned.
555	_sub(obj,obj) Simple substraction of 1 object from another
556	a third, optional parameter indicates that the params
557	are swapped. In this case, the first param needs to
558	be preserved, while you can destroy the second.
559	sub (x,y,1) => return x - y and keep x intact!
560
561	_acmp(obj,obj) <=> operator for objects (return -1, 0 or 1)
562
563	_len(obj) returns count of the decimal digits of the object
564	_digit(obj,n) returns the n'th decimal digit of object
565
566	_is_one(obj) return true if argument is +1
567	_is_zero(obj) return true if argument is 0
568	_is_even(obj) return true if argument is even (0,2,4,6..)
569	_is_odd(obj) return true if argument is odd (1,3,5,7..)
570
571	_copy return a ref to a true copy of the object
572
573	_check(obj) check whether internal representation is still intact
574	return 0 for ok, otherwise error message as string
575
576	The following functions are optional, and can be exported if the underlying lib
577	has a fast way to do them. If not defined, Math::BigInt will use a pure, but
578	slow Perl function as fallback to emulate these:
579
580	_from_hex(str) return ref to new object from ref to hexadecimal string
581	_from_bin(str) return ref to new object from ref to binary string
582
583	_rsft(obj,N,B) shift object in base B by N 'digits' right
584	_lsft(obj,N,B) shift object in base B by N 'digits' left
585
586	_xor(obj1,obj2) XOR (bit-wise) object 1 with object 2
587	Mote: XOR, AND and OR pad with zeros if size mismatches
588	_and(obj1,obj2) AND (bit-wise) object 1 with object 2
589	_or(obj1,obj2) OR (bit-wise) object 1 with object 2
590
591	_sqrt(obj) return the square root of object
592	_pow(obj,obj) return object 1 to the power of object 2
593	_gcd(obj,obj) return Greatest Common Divisor of two objects
594
595	_zeros(obj) return amount of trailing decimal zeros
596
597	_dec(obj) decrement object by one (input is >= 1)
598	_inc(obj) increment object by one
599
600	Input strings come in as unsigned but with prefix (aka as '123', '0xabc'
601	or '0b1101').
602
603	Testing of input parameter validity is done by the caller, so you need not to
604	worry about underflow (C<_sub()>, C<_dec()>) nor division by zero or similiar
605	cases.
606
607	=head1 LICENSE
608
609	This program is free software; you may redistribute it and/or modify it under
610	the same terms as Perl itself.
611
612	=head1 AUTHORS
613
614	Original math code by Mark Biggar, rewritten by Tels L<http://bloodgate.com/>
615	in late 2000, 2001.
616	Seperated from BigInt and shaped API with the help of John Peacock.
617
618	=head1 SEE ALSO
619
620	L<Math::BigInt>, L<Math::BigFloat>.
621
622	=cut