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

package Math::BigInt;
$VERSION='0.01';

use overload
'+'	=>	sub {new Math::BigInt &badd},
'-'	=>	sub {new Math::BigInt
		       $_[2]? bsub($_[1],${$_[0]}) : bsub(${$_[0]},$_[1])},
'<=>'	=>	sub {$_[2]? bcmp($_[1],${$_[0]}) : bcmp(${$_[0]},$_[1])},
'cmp'	=>	sub {$_[2]? ($_[1] cmp ${$_[0]}) : (${$_[0]} cmp $_[1])},
'*'	=>	sub {new Math::BigInt &bmul},
'/'	=>	sub {new Math::BigInt 
		       $_[2]? scalar bdiv($_[1],${$_[0]}) :
			 scalar bdiv(${$_[0]},$_[1])},
'%'	=>	sub {new Math::BigInt
		       $_[2]? bmod($_[1],${$_[0]}) : bmod(${$_[0]},$_[1])},
'**'	=>	sub {new Math::BigInt
		       $_[2]? bpow($_[1],${$_[0]}) : bpow(${$_[0]},$_[1])},
'neg'	=>	sub {new Math::BigInt &bneg},
'abs'	=>	sub {new Math::BigInt &babs},
'<<'	=>	sub {new Math::BigInt
		       $_[2]? blsft($_[1],${$_[0]}) : blsft(${$_[0]},$_[1])},
'>>'	=>	sub {new Math::BigInt
		       $_[2]? brsft($_[1],${$_[0]}) : brsft(${$_[0]},$_[1])},
'&'	=>	sub {new Math::BigInt &band},
'|'	=>	sub {new Math::BigInt &bior},
'^'	=>	sub {new Math::BigInt &bxor},
'~'	=>	sub {new Math::BigInt &bnot},

qw(
""	stringify
0+	numify)			# Order of arguments unsignificant
;

$NaNOK=1;

sub new {
  my($class) = shift;
  my($foo) = bnorm(shift);
  die "Not a number initialized to Math::BigInt" if !$NaNOK && $foo eq "NaN";
  bless \$foo, $class;
}
sub stringify { "${$_[0]}" }
sub numify { 0 + "${$_[0]}" }	# Not needed, additional overhead
				# comparing to direct compilation based on
				# stringify
sub import {
  shift;
  return unless @_;
  die "unknown import: @_" unless @_ == 1 and $_[0] eq ':constant';
  overload::constant integer => sub {Math::BigInt->new(shift)};
}

$zero = 0;

# overcome a floating point problem on certain osnames (posix-bc, os390)
BEGIN {
    my $x = 100000.0;
    my $use_mult = int($x*1e-5)*1e5 == $x ? 1 : 0;
}

# normalize string form of number.   Strip leading zeros.  Strip any
#   white space and add a sign, if missing.
# Strings that are not numbers result the value 'NaN'.

sub bnorm { #(num_str) return num_str
    local($_) = @_;
    s/\s+//g;                           # strip white space
    if (s/^([+-]?)0*(\d+)$/$1$2/) {     # test if number
	substr($_,$[,0) = '+' unless $1; # Add missing sign
	s/^-0/+0/;
	$_;
    } else {
	'NaN';
    }
}

# Convert a number from string format to internal base 100000 format.
#   Assumes normalized value as input.
sub internal { #(num_str) return int_num_array
    local($d) = @_;
    ($is,$il) = (substr($d,$[,1),length($d)-2);
    substr($d,$[,1) = '';
    ($is, reverse(unpack("a" . ($il%5+1) . ("a5" x ($il/5)), $d)));
}

# Convert a number from internal base 100000 format to string format.
#   This routine scribbles all over input array.
sub external { #(int_num_array) return num_str
    $es = shift;
    grep($_ > 9999 || ($_ = substr('0000'.$_,-5)), @_);   # zero pad
    &bnorm(join('', $es, reverse(@_)));    # reverse concat and normalize
}

# Negate input value.
sub bneg { #(num_str) return num_str
    local($_) = &bnorm(@_);
    return $_ if $_ eq '+0' or $_ eq 'NaN';
    vec($_,0,8) ^= ord('+') ^ ord('-');
    $_;
}

# Returns the absolute value of the input.
sub babs { #(num_str) return num_str
    &abs(&bnorm(@_));
}

sub abs { # post-normalized abs for internal use
    local($_) = @_;
    s/^-/+/;
    $_;
}

# Compares 2 values.  Returns one of undef, <0, =0, >0. (suitable for sort)
sub bcmp { #(num_str, num_str) return cond_code
    local($x,$y) = (&bnorm($_[$[]),&bnorm($_[$[+1]));
    if ($x eq 'NaN') {
	undef;
    } elsif ($y eq 'NaN') {
	undef;
    } else {
	&cmp($x,$y) <=> 0;
    }
}

sub cmp { # post-normalized compare for internal use
    local($cx, $cy) = @_;
    
    return 0 if ($cx eq $cy);

    local($sx, $sy) = (substr($cx, 0, 1), substr($cy, 0, 1));
    local($ld);

    if ($sx eq '+') {
      return  1 if ($sy eq '-' || $cy eq '+0');
      $ld = length($cx) - length($cy);
      return $ld if ($ld);
      return $cx cmp $cy;
    } else { # $sx eq '-'
      return -1 if ($sy eq '+');
      $ld = length($cy) - length($cx);
      return $ld if ($ld);
      return $cy cmp $cx;
    }
}

sub badd { #(num_str, num_str) return num_str
    local(*x, *y); ($x, $y) = (&bnorm($_[$[]),&bnorm($_[$[+1]));
    if ($x eq 'NaN') {
	'NaN';
    } elsif ($y eq 'NaN') {
	'NaN';
    } else {
	@x = &internal($x);             # convert to internal form
	@y = &internal($y);
	local($sx, $sy) = (shift @x, shift @y); # get signs
	if ($sx eq $sy) {
	    &external($sx, &add(*x, *y)); # if same sign add
	} else {
	    ($x, $y) = (&abs($x),&abs($y)); # make abs
	    if (&cmp($y,$x) > 0) {
		&external($sy, &sub(*y, *x));
	    } else {
		&external($sx, &sub(*x, *y));
	    }
	}
    }
}

sub bsub { #(num_str, num_str) return num_str
    &badd($_[$[],&bneg($_[$[+1]));    
}

# GCD -- Euclids algorithm Knuth Vol 2 pg 296
sub bgcd { #(num_str, num_str) return num_str
    local($x,$y) = (&bnorm($_[$[]),&bnorm($_[$[+1]));
    if ($x eq 'NaN' || $y eq 'NaN') {
	'NaN';
    } else {
	($x, $y) = ($y,&bmod($x,$y)) while $y ne '+0';
	$x;
    }
}

# 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
sub add { #(int_num_array, int_num_array) return int_num_array
    local(*x, *y) = @_;
    $car = 0;
    for $x (@x) {
	last unless @y || $car;
	$x -= 1e5 if $car = (($x += (@y ? shift(@y) : 0) + $car) >= 1e5) ? 1 : 0;
    }
    for $y (@y) {
	last unless $car;
	$y -= 1e5 if $car = (($y += $car) >= 1e5) ? 1 : 0;
    }
    (@x, @y, $car);
}

# subtract base 1e5 numbers -- stolen from Knuth Vol 2 pg 232, $x > $y
sub sub { #(int_num_array, int_num_array) return int_num_array
    local(*sx, *sy) = @_;
    $bar = 0;
    for $sx (@sx) {
	last unless @sy || $bar;
	$sx += 1e5 if $bar = (($sx -= (@sy ? shift(@sy) : 0) + $bar) < 0);
    }
    @sx;
}

# multiply two numbers -- stolen from Knuth Vol 2 pg 233
sub bmul { #(num_str, num_str) return num_str
    local(*x, *y); ($x, $y) = (&bnorm($_[$[]), &bnorm($_[$[+1]));
    if ($x eq 'NaN') {
	'NaN';
    } elsif ($y eq 'NaN') {
	'NaN';
    } else {
	@x = &internal($x);
	@y = &internal($y);
	&external(&mul(*x,*y));
    }
}

# multiply two numbers in internal representation
# destroys the arguments, supposes that two arguments are different
sub mul { #(*int_num_array, *int_num_array) return int_num_array
    local(*x, *y) = (shift, shift);
    local($signr) = (shift @x ne shift @y) ? '-' : '+';
    @prod = ();
    for $x (@x) {
      ($car, $cty) = (0, $[);
      for $y (@y) {
	$prod = $x * $y + ($prod[$cty] || 0) + $car;
        if ($use_mult) {
	$prod[$cty++] =
	  $prod - ($car = int($prod * 1e-5)) * 1e5;
        }
        else {
	$prod[$cty++] =
	  $prod - ($car = int($prod / 1e5)) * 1e5;
        }
      }
      $prod[$cty] += $car if $car;
      $x = shift @prod;
    }
    ($signr, @x, @prod);
}

# modulus
sub bmod { #(num_str, num_str) return num_str
    (&bdiv(@_))[$[+1];
}

sub bdiv { #(dividend: num_str, divisor: num_str) return num_str
    local (*x, *y); ($x, $y) = (&bnorm($_[$[]), &bnorm($_[$[+1]));
    return wantarray ? ('NaN','NaN') : 'NaN'
	if ($x eq 'NaN' || $y eq 'NaN' || $y eq '+0');
    return wantarray ? ('+0',$x) : '+0' if (&cmp(&abs($x),&abs($y)) < 0);
    @x = &internal($x); @y = &internal($y);
    $srem = $y[$[];
    $sr = (shift @x ne shift @y) ? '-' : '+';
    $car = $bar = $prd = 0;
    if (($dd = int(1e5/($y[$#y]+1))) != 1) {
	for $x (@x) {
	    $x = $x * $dd + $car;
            if ($use_mult) {
	    $x -= ($car = int($x * 1e-5)) * 1e5;
            }
            else {
	    $x -= ($car = int($x / 1e5)) * 1e5;
            }
	}
	push(@x, $car); $car = 0;
	for $y (@y) {
	    $y = $y * $dd + $car;
            if ($use_mult) {
	    $y -= ($car = int($y * 1e-5)) * 1e5;
            }
            else {
	    $y -= ($car = int($y / 1e5)) * 1e5;
            }
	}
    }
    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;
	$q = (($u0 == $v1) ? 99999 : int(($u0*1e5+$u1)/$v1));
	--$q while ($v2*$q > ($u0*1e5+$u1-$q*$v1)*1e5+$u2);
	if ($q) {
	    ($car, $bar) = (0,0);
	    for ($y = $[, $x = $#x-$#y+$[-1; $y <= $#y; ++$y,++$x) {
		$prd = $q * $y[$y] + $car;
                if ($use_mult) {
		$prd -= ($car = int($prd * 1e-5)) * 1e5;
                }
                else {
		$prd -= ($car = int($prd / 1e5)) * 1e5;
                }
		$x[$x] += 1e5 if ($bar = (($x[$x] -= $prd + $bar) < 0));
	    }
	    if ($x[$#x] < $car + $bar) {
		$car = 0; --$q;
		for ($y = $[, $x = $#x-$#y+$[-1; $y <= $#y; ++$y,++$x) {
		    $x[$x] -= 1e5
			if ($car = (($x[$x] += $y[$y] + $car) > 1e5));
		}
	    }   
	}
	pop(@x); unshift(@q, $q);
    }
    if (wantarray) {
	@d = ();
	if ($dd != 1) {
	    $car = 0;
	    for $x (reverse @x) {
		$prd = $car * 1e5 + $x;
		$car = $prd - ($tmp = int($prd / $dd)) * $dd;
		unshift(@d, $tmp);
	    }
	}
	else {
	    @d = @x;
	}
	(&external($sr, @q), &external($srem, @d, $zero));
    } else {
	&external($sr, @q);
    }
}

# compute power of two numbers -- stolen from Knuth Vol 2 pg 233
sub bpow { #(num_str, num_str) return num_str
    local(*x, *y); ($x, $y) = (&bnorm($_[$[]), &bnorm($_[$[+1]));
    if ($x eq 'NaN') {
	'NaN';
    } elsif ($y eq 'NaN') {
	'NaN';
    } elsif ($x eq '+1') {
	'+1';
    } elsif ($x eq '-1') {
	&bmod($x,2) ? '-1': '+1';
    } elsif ($y =~ /^-/) {
	'NaN';
    } elsif ($x eq '+0' && $y eq '+0') {
	'NaN';
    } else {
	@x = &internal($x);
	local(@pow2)=@x;
	local(@pow)=&internal("+1");
	local($y1,$res,@tmp1,@tmp2)=(1); # need tmp to send to mul
	while ($y ne '+0') {
	  ($y,$res)=&bdiv($y,2);
	  if ($res ne '+0') {@tmp=@pow2; @pow=&mul(*pow,*tmp);}
	  if ($y ne '+0') {@tmp=@pow2;@pow2=&mul(*pow2,*tmp);}
	}
	&external(@pow);
    }
}

# compute x << y, y >= 0
sub blsft { #(num_str, num_str) return num_str
    &bmul($_[$[], &bpow(2, $_[$[+1]));
}

# compute x >> y, y >= 0
sub brsft { #(num_str, num_str) return num_str
    &bdiv($_[$[], &bpow(2, $_[$[+1]));
}

# compute x & y
sub band { #(num_str, num_str) return num_str
    local($x,$y,$r,$m,$xr,$yr) = (&bnorm($_[$[]),&bnorm($_[$[+1]),0,1);
    if ($x eq 'NaN' || $y eq 'NaN') {
	'NaN';
    } else {
	while ($x ne '+0' && $y ne '+0') {
	    ($x, $xr) = &bdiv($x, 0x10000);
	    ($y, $yr) = &bdiv($y, 0x10000);
	    $r = &badd(&bmul(int $xr & $yr, $m), $r);
	    $m = &bmul($m, 0x10000);
	}
	$r;
    }
}

# compute x | y
sub bior { #(num_str, num_str) return num_str
    local($x,$y,$r,$m,$xr,$yr) = (&bnorm($_[$[]),&bnorm($_[$[+1]),0,1);
    if ($x eq 'NaN' || $y eq 'NaN') {
	'NaN';
    } else {
	while ($x ne '+0' || $y ne '+0') {
	    ($x, $xr) = &bdiv($x, 0x10000);
	    ($y, $yr) = &bdiv($y, 0x10000);
	    $r = &badd(&bmul(int $xr | $yr, $m), $r);
	    $m = &bmul($m, 0x10000);
	}
	$r;
    }
}

# compute x ^ y
sub bxor { #(num_str, num_str) return num_str
    local($x,$y,$r,$m,$xr,$yr) = (&bnorm($_[$[]),&bnorm($_[$[+1]),0,1);
    if ($x eq 'NaN' || $y eq 'NaN') {
	'NaN';
    } else {
	while ($x ne '+0' || $y ne '+0') {
	    ($x, $xr) = &bdiv($x, 0x10000);
	    ($y, $yr) = &bdiv($y, 0x10000);
	    $r = &badd(&bmul(int $xr ^ $yr, $m), $r);
	    $m = &bmul($m, 0x10000);
	}
	$r;
    }
}

# represent ~x as twos-complement number
sub bnot { #(num_str) return num_str
    &bsub(-1,$_[$[]);
}

1;
__END__

=head1 NAME

Math::BigInt - Arbitrary size integer math package

=head1 SYNOPSIS

  use Math::BigInt;
  $i = Math::BigInt->new($string);

  $i->bneg return BINT               negation
  $i->babs return BINT               absolute value
  $i->bcmp(BINT) return CODE         compare numbers (undef,<0,=0,>0)
  $i->badd(BINT) return BINT         addition
  $i->bsub(BINT) return BINT         subtraction
  $i->bmul(BINT) return BINT         multiplication
  $i->bdiv(BINT) return (BINT,BINT)  division (quo,rem) just quo if scalar
  $i->bmod(BINT) return BINT         modulus
  $i->bgcd(BINT) return BINT         greatest common divisor
  $i->bnorm return BINT              normalization
  $i->blsft(BINT) return BINT        left shift
  $i->brsft(BINT) return (BINT,BINT) right shift (quo,rem) just quo if scalar
  $i->band(BINT) return BINT         bit-wise and
  $i->bior(BINT) return BINT         bit-wise inclusive or
  $i->bxor(BINT) return BINT         bit-wise exclusive or
  $i->bnot return BINT               bit-wise not

=head1 DESCRIPTION

All basic math operations are overloaded if you declare your big
integers as

  $i = new Math::BigInt '123 456 789 123 456 789';


=over 2

=item Canonical notation

Big integer value are strings of the form C</^[+-]\d+$/> with leading
zeros suppressed.

=item Input

Input values to these routines may be strings of the form
C</^\s*[+-]?[\d\s]+$/>.

=item Output

Output values always always in canonical form

=back

Actual math is done in an internal format consisting of an array
whose first element is the sign (/^[+-]$/) and whose remaining 
elements are base 100000 digits with the least significant digit first.
The string 'NaN' is used to represent the result when input arguments 
are not numbers, as well as the result of dividing by zero.

=head1 EXAMPLES

   '+0'                            canonical zero value
   '   -123 123 123'               canonical value '-123123123'
   '1 23 456 7890'                 canonical value '+1234567890'


=head1 Autocreating constants

After C<use Math::BigInt ':constant'> all the integer decimal constants
in the given scope are converted to C<Math::BigInt>.  This conversion
happens at compile time.

In particular

  perl -MMath::BigInt=:constant -e 'print 2**100'

print the integer value of C<2**100>.  Note that without conversion of 
constants the expression 2**100 will be calculated as floating point number.

=head1 BUGS

The current version of this module is a preliminary version of the
real thing that is currently (as of perl5.002) under development.

=head1 AUTHOR

Mark Biggar, overloaded interface by Ilya Zakharevich.

=cut
Commit	Line	Data
a0d0e21e	1	package Math::BigInt;
3deb277d	2	$VERSION='0.01';
a0d0e21e	3
a5f75d66	4	use overload
748a9306	5	'+' => sub {new Math::BigInt &badd},
748a9306	6	'-' => sub {new Math::BigInt
a0d0e21e	7	$_[2]? bsub($_[1],${$_[0]}) : bsub(${$_[0]},$_[1])},
5d7098d5	8	'<=>' => sub {$_[2]? bcmp($_[1],${$_[0]}) : bcmp(${$_[0]},$_[1])},
5d7098d5	9	'cmp' => sub {$_[2]? ($_[1] cmp ${$_[0]}) : (${$_[0]} cmp $_[1])},
748a9306	10	'*' => sub {new Math::BigInt &bmul},
748a9306	11	'/' => sub {new Math::BigInt
a0d0e21e	12	$_[2]? scalar bdiv($_[1],${$_[0]}) :
a0d0e21e	13	scalar bdiv(${$_[0]},$_[1])},
748a9306	14	'%' => sub {new Math::BigInt
a0d0e21e	15	$_[2]? bmod($_[1],${$_[0]}) : bmod(${$_[0]},$_[1])},
748a9306	16	'**' => sub {new Math::BigInt
a0d0e21e	17	$_[2]? bpow($_[1],${$_[0]}) : bpow(${$_[0]},$_[1])},
748a9306	18	'neg' => sub {new Math::BigInt &bneg},
748a9306	19	'abs' => sub {new Math::BigInt &babs},
e16b8f49	20	'<<' => sub {new Math::BigInt
	21	$_[2]? blsft($_[1],${$_[0]}) : blsft(${$_[0]},$_[1])},
	22	'>>' => sub {new Math::BigInt
	23	$_[2]? brsft($_[1],${$_[0]}) : brsft(${$_[0]},$_[1])},
	24	'&' => sub {new Math::BigInt &band},
	25	'\|' => sub {new Math::BigInt &bior},
	26	'^' => sub {new Math::BigInt &bxor},
	27	'~' => sub {new Math::BigInt &bnot},
a0d0e21e	28
	29	qw(
	30	"" stringify
	31	0+ numify) # Order of arguments unsignificant
a5f75d66	32	;
a0d0e21e	33
748a9306	34	$NaNOK=1;
748a9306	35
a0d0e21e	36	sub new {
a5f75d66	37	my($class) = shift;
a5f75d66	38	my($foo) = bnorm(shift);
748a9306	39	die "Not a number initialized to Math::BigInt" if !$NaNOK && $foo eq "NaN";
a5f75d66	40	bless \$foo, $class;
a0d0e21e	41	}
	42	sub stringify { "${$_[0]}" }
	43	sub numify { 0 + "${$_[0]}" } # Not needed, additional overhead
	44	# comparing to direct compilation based on
	45	# stringify
b3ac6de7	46	sub import {
	47	shift;
	48	return unless @_;
	49	die "unknown import: @_" unless @_ == 1 and $_[0] eq ':constant';
	50	overload::constant integer => sub {Math::BigInt->new(shift)};
	51	}
a0d0e21e	52
a0d0e21e	53	$zero = 0;
a0d0e21e	54
1f45ae4a	55	# overcome a floating point problem on certain osnames (posix-bc, os390)
	56	BEGIN {
	57	my $x = 100000.0;
	58	my $use_mult = int($x1e-5)1e5 == $x ? 1 : 0;
	59	}
a5f75d66	60
a0d0e21e	61	# normalize string form of number. Strip leading zeros. Strip any
	62	# white space and add a sign, if missing.
	63	# Strings that are not numbers result the value 'NaN'.
	64
	65	sub bnorm { #(num_str) return num_str
	66	local($_) = @_;
	67	s/\s+//g; # strip white space
	68	if (s/^([+-]?)0*(\d+)$/$1$2/) { # test if number
	69	substr($_,$[,0) = '+' unless $1; # Add missing sign
	70	s/^-0/+0/;
	71	$_;
	72	} else {
	73	'NaN';
	74	}
	75	}
	76
	77	# Convert a number from string format to internal base 100000 format.
	78	# Assumes normalized value as input.
	79	sub internal { #(num_str) return int_num_array
	80	local($d) = @_;
	81	($is,$il) = (substr($d,$[,1),length($d)-2);
	82	substr($d,$[,1) = '';
	83	($is, reverse(unpack("a" . ($il%5+1) . ("a5" x ($il/5)), $d)));
	84	}
	85
	86	# Convert a number from internal base 100000 format to string format.
	87	# This routine scribbles all over input array.
	88	sub external { #(int_num_array) return num_str
	89	$es = shift;
	90	grep($_ > 9999 \|\| ($_ = substr('0000'.$_,-5)), @_); # zero pad
	91	&bnorm(join('', $es, reverse(@_))); # reverse concat and normalize
	92	}
	93
	94	# Negate input value.
	95	sub bneg { #(num_str) return num_str
	96	local($_) = &bnorm(@_);
e3c7ef20	97	return $_ if $_ eq '+0' or $_ eq 'NaN';
e3c7ef20	98	vec($_,0,8) ^= ord('+') ^ ord('-');
a0d0e21e	99	$_;
	100	}
	101
	102	# Returns the absolute value of the input.
	103	sub babs { #(num_str) return num_str
	104	&abs(&bnorm(@_));
	105	}
	106
	107	sub abs { # post-normalized abs for internal use
	108	local($_) = @_;
	109	s/^-/+/;
	110	$_;
	111	}
a5f75d66	112
a0d0e21e	113	# Compares 2 values. Returns one of undef, <0, =0, >0. (suitable for sort)
	114	sub bcmp { #(num_str, num_str) return cond_code
	115	local($x,$y) = (&bnorm($_[$[]),&bnorm($_[$[+1]));
	116	if ($x eq 'NaN') {
	117	undef;
	118	} elsif ($y eq 'NaN') {
	119	undef;
	120	} else {
e3c7ef20	121	&cmp($x,$y) <=> 0;
a0d0e21e	122	}
	123	}
	124
	125	sub cmp { # post-normalized compare for internal use
	126	local($cx, $cy) = @_;
1e2e1ae8	127
	128	return 0 if ($cx eq $cy);
	129
	130	local($sx, $sy) = (substr($cx, 0, 1), substr($cy, 0, 1));
	131	local($ld);
	132
	133	if ($sx eq '+') {
	134	return 1 if ($sy eq '-' \|\| $cy eq '+0');
	135	$ld = length($cx) - length($cy);
	136	return $ld if ($ld);
	137	return $cx cmp $cy;
	138	} else { # $sx eq '-'
	139	return -1 if ($sy eq '+');
	140	$ld = length($cy) - length($cx);
	141	return $ld if ($ld);
	142	return $cy cmp $cx;
	143	}
a0d0e21e	144	}
	145
	146	sub badd { #(num_str, num_str) return num_str
	147	local(x, y); ($x, $y) = (&bnorm($_[$[]),&bnorm($_[$[+1]));
	148	if ($x eq 'NaN') {
	149	'NaN';
	150	} elsif ($y eq 'NaN') {
	151	'NaN';
	152	} else {
	153	@x = &internal($x); # convert to internal form
	154	@y = &internal($y);
	155	local($sx, $sy) = (shift @x, shift @y); # get signs
	156	if ($sx eq $sy) {
	157	&external($sx, &add(x, y)); # if same sign add
	158	} else {
	159	($x, $y) = (&abs($x),&abs($y)); # make abs
	160	if (&cmp($y,$x) > 0) {
	161	&external($sy, &sub(y, x));
	162	} else {
	163	&external($sx, &sub(x, y));
	164	}
	165	}
	166	}
	167	}
	168
	169	sub bsub { #(num_str, num_str) return num_str
	170	&badd($_[$[],&bneg($_[$[+1]));
	171	}
	172
	173	# GCD -- Euclids algorithm Knuth Vol 2 pg 296
	174	sub bgcd { #(num_str, num_str) return num_str
	175	local($x,$y) = (&bnorm($_[$[]),&bnorm($_[$[+1]));
	176	if ($x eq 'NaN' \|\| $y eq 'NaN') {
	177	'NaN';
	178	} else {
	179	($x, $y) = ($y,&bmod($x,$y)) while $y ne '+0';
	180	$x;
	181	}
	182	}
a5f75d66	183
a0d0e21e	184	# routine to add two base 1e5 numbers
	185	# stolen from Knuth Vol 2 Algorithm A pg 231
	186	# there are separate routines to add and sub as per Kunth pg 233
	187	sub add { #(int_num_array, int_num_array) return int_num_array
	188	local(x, y) = @_;
	189	$car = 0;
	190	for $x (@x) {
	191	last unless @y \|\| $car;
20408e3c	192	$x -= 1e5 if $car = (($x += (@y ? shift(@y) : 0) + $car) >= 1e5) ? 1 : 0;
a0d0e21e	193	}
	194	for $y (@y) {
	195	last unless $car;
55497cff	196	$y -= 1e5 if $car = (($y += $car) >= 1e5) ? 1 : 0;
a0d0e21e	197	}
	198	(@x, @y, $car);
	199	}
	200
	201	# subtract base 1e5 numbers -- stolen from Knuth Vol 2 pg 232, $x > $y
	202	sub sub { #(int_num_array, int_num_array) return int_num_array
	203	local(sx, sy) = @_;
	204	$bar = 0;
	205	for $sx (@sx) {
20408e3c	206	last unless @sy \|\| $bar;
20408e3c	207	$sx += 1e5 if $bar = (($sx -= (@sy ? shift(@sy) : 0) + $bar) < 0);
a0d0e21e	208	}
	209	@sx;
	210	}
	211
	212	# multiply two numbers -- stolen from Knuth Vol 2 pg 233
	213	sub bmul { #(num_str, num_str) return num_str
	214	local(x, y); ($x, $y) = (&bnorm($_[$[]), &bnorm($_[$[+1]));
	215	if ($x eq 'NaN') {
	216	'NaN';
	217	} elsif ($y eq 'NaN') {
	218	'NaN';
	219	} else {
	220	@x = &internal($x);
	221	@y = &internal($y);
	222	&external(&mul(x,y));
	223	}
	224	}
	225
	226	# multiply two numbers in internal representation
	227	# destroys the arguments, supposes that two arguments are different
	228	sub mul { #(int_num_array, int_num_array) return int_num_array
	229	local(x, y) = (shift, shift);
	230	local($signr) = (shift @x ne shift @y) ? '-' : '+';
	231	@prod = ();
	232	for $x (@x) {
	233	($car, $cty) = (0, $[);
	234	for $y (@y) {
8b6a6e55	235	$prod = $x * $y + ($prod[$cty] \|\| 0) + $car;
1f45ae4a	236	if ($use_mult) {
a0d0e21e	237	$prod[$cty++] =
a0d0e21e	238	$prod - ($car = int($prod * 1e-5)) * 1e5;
1f45ae4a	239	}
	240	else {
	241	$prod[$cty++] =
	242	$prod - ($car = int($prod / 1e5)) * 1e5;
	243	}
a0d0e21e	244	}
	245	$prod[$cty] += $car if $car;
	246	$x = shift @prod;
	247	}
	248	($signr, @x, @prod);
	249	}
	250
	251	# modulus
	252	sub bmod { #(num_str, num_str) return num_str
	253	(&bdiv(@_))[$[+1];
	254	}
a5f75d66	255
a0d0e21e	256	sub bdiv { #(dividend: num_str, divisor: num_str) return num_str
	257	local (x, y); ($x, $y) = (&bnorm($_[$[]), &bnorm($_[$[+1]));
	258	return wantarray ? ('NaN','NaN') : 'NaN'
	259	if ($x eq 'NaN' \|\| $y eq 'NaN' \|\| $y eq '+0');
	260	return wantarray ? ('+0',$x) : '+0' if (&cmp(&abs($x),&abs($y)) < 0);
	261	@x = &internal($x); @y = &internal($y);
	262	$srem = $y[$[];
	263	$sr = (shift @x ne shift @y) ? '-' : '+';
	264	$car = $bar = $prd = 0;
	265	if (($dd = int(1e5/($y[$#y]+1))) != 1) {
	266	for $x (@x) {
	267	$x = $x * $dd + $car;
1f45ae4a	268	if ($use_mult) {
a0d0e21e	269	$x -= ($car = int($x * 1e-5)) * 1e5;
1f45ae4a	270	}
	271	else {
	272	$x -= ($car = int($x / 1e5)) * 1e5;
	273	}
a0d0e21e	274	}
	275	push(@x, $car); $car = 0;
	276	for $y (@y) {
	277	$y = $y * $dd + $car;
1f45ae4a	278	if ($use_mult) {
a0d0e21e	279	$y -= ($car = int($y * 1e-5)) * 1e5;
1f45ae4a	280	}
	281	else {
	282	$y -= ($car = int($y / 1e5)) * 1e5;
	283	}
a0d0e21e	284	}
	285	}
	286	else {
	287	push(@x, 0);
	288	}
5d7098d5	289	@q = (); ($v2,$v1) = @y[-2,-1];
0a6a0d52	290	$v2 = 0 unless $v2;
a0d0e21e	291	while ($#x > $#y) {
5d7098d5	292	($u2,$u1,$u0) = @x[-3..-1];
0a6a0d52	293	$u2 = 0 unless $u2;
a0d0e21e	294	$q = (($u0 == $v1) ? 99999 : int(($u0*1e5+$u1)/$v1));
	295	--$q while ($v2$q > ($u01e5+$u1-$q$v1)1e5+$u2);
	296	if ($q) {
	297	($car, $bar) = (0,0);
	298	for ($y = $[, $x = $#x-$#y+$[-1; $y <= $#y; ++$y,++$x) {
	299	$prd = $q * $y[$y] + $car;
1f45ae4a	300	if ($use_mult) {
a0d0e21e	301	$prd -= ($car = int($prd * 1e-5)) * 1e5;
1f45ae4a	302	}
	303	else {
	304	$prd -= ($car = int($prd / 1e5)) * 1e5;
	305	}
a0d0e21e	306	$x[$x] += 1e5 if ($bar = (($x[$x] -= $prd + $bar) < 0));
	307	}
	308	if ($x[$#x] < $car + $bar) {
	309	$car = 0; --$q;
	310	for ($y = $[, $x = $#x-$#y+$[-1; $y <= $#y; ++$y,++$x) {
	311	$x[$x] -= 1e5
	312	if ($car = (($x[$x] += $y[$y] + $car) > 1e5));
	313	}
	314	}
	315	}
	316	pop(@x); unshift(@q, $q);
	317	}
	318	if (wantarray) {
	319	@d = ();
	320	if ($dd != 1) {
	321	$car = 0;
	322	for $x (reverse @x) {
	323	$prd = $car * 1e5 + $x;
	324	$car = $prd - ($tmp = int($prd / $dd)) * $dd;
	325	unshift(@d, $tmp);
	326	}
	327	}
	328	else {
	329	@d = @x;
	330	}
	331	(&external($sr, @q), &external($srem, @d, $zero));
	332	} else {
	333	&external($sr, @q);
	334	}
	335	}
	336
	337	# compute power of two numbers -- stolen from Knuth Vol 2 pg 233
	338	sub bpow { #(num_str, num_str) return num_str
	339	local(x, y); ($x, $y) = (&bnorm($_[$[]), &bnorm($_[$[+1]));
	340	if ($x eq 'NaN') {
	341	'NaN';
	342	} elsif ($y eq 'NaN') {
	343	'NaN';
	344	} elsif ($x eq '+1') {
	345	'+1';
	346	} elsif ($x eq '-1') {
	347	&bmod($x,2) ? '-1': '+1';
	348	} elsif ($y =~ /^-/) {
	349	'NaN';
	350	} elsif ($x eq '+0' && $y eq '+0') {
	351	'NaN';
	352	} else {
	353	@x = &internal($x);
	354	local(@pow2)=@x;
	355	local(@pow)=&internal("+1");
	356	local($y1,$res,@tmp1,@tmp2)=(1); # need tmp to send to mul
	357	while ($y ne '+0') {
	358	($y,$res)=&bdiv($y,2);
	359	if ($res ne '+0') {@tmp=@pow2; @pow=&mul(pow,tmp);}
	360	if ($y ne '+0') {@tmp=@pow2;@pow2=&mul(pow2,tmp);}
	361	}
	362	&external(@pow);
	363	}
	364	}
	365
e16b8f49	366	# compute x << y, y >= 0
	367	sub blsft { #(num_str, num_str) return num_str
	368	&bmul($_[$[], &bpow(2, $_[$[+1]));
	369	}
	370
	371	# compute x >> y, y >= 0
	372	sub brsft { #(num_str, num_str) return num_str
	373	&bdiv($_[$[], &bpow(2, $_[$[+1]));
	374	}
	375
	376	# compute x & y
	377	sub band { #(num_str, num_str) return num_str
	378	local($x,$y,$r,$m,$xr,$yr) = (&bnorm($_[$[]),&bnorm($_[$[+1]),0,1);
	379	if ($x eq 'NaN' \|\| $y eq 'NaN') {
	380	'NaN';
	381	} else {
	382	while ($x ne '+0' && $y ne '+0') {
	383	($x, $xr) = &bdiv($x, 0x10000);
	384	($y, $yr) = &bdiv($y, 0x10000);
	385	$r = &badd(&bmul(int $xr & $yr, $m), $r);
	386	$m = &bmul($m, 0x10000);
	387	}
	388	$r;
	389	}
	390	}
	391
	392	# compute x \| y
	393	sub bior { #(num_str, num_str) return num_str
	394	local($x,$y,$r,$m,$xr,$yr) = (&bnorm($_[$[]),&bnorm($_[$[+1]),0,1);
	395	if ($x eq 'NaN' \|\| $y eq 'NaN') {
	396	'NaN';
	397	} else {
	398	while ($x ne '+0' \|\| $y ne '+0') {
	399	($x, $xr) = &bdiv($x, 0x10000);
	400	($y, $yr) = &bdiv($y, 0x10000);
	401	$r = &badd(&bmul(int $xr \| $yr, $m), $r);
	402	$m = &bmul($m, 0x10000);
	403	}
	404	$r;
	405	}
	406	}
	407
	408	# compute x ^ y
	409	sub bxor { #(num_str, num_str) return num_str
	410	local($x,$y,$r,$m,$xr,$yr) = (&bnorm($_[$[]),&bnorm($_[$[+1]),0,1);
	411	if ($x eq 'NaN' \|\| $y eq 'NaN') {
	412	'NaN';
	413	} else {
	414	while ($x ne '+0' \|\| $y ne '+0') {
	415	($x, $xr) = &bdiv($x, 0x10000);
	416	($y, $yr) = &bdiv($y, 0x10000);
	417	$r = &badd(&bmul(int $xr ^ $yr, $m), $r);
	418	$m = &bmul($m, 0x10000);
	419	}
	420	$r;
	421	}
	422	}
	423
	424	# represent ~x as twos-complement number
	425	sub bnot { #(num_str) return num_str
	426	&bsub(-1,$_[$[]);
	427	}
	428
a0d0e21e	429	1;
a5f75d66	430	__END__
	431
	432	=head1 NAME
	433
	434	Math::BigInt - Arbitrary size integer math package
	435
	436	=head1 SYNOPSIS
	437
	438	use Math::BigInt;
	439	$i = Math::BigInt->new($string);
	440
	441	$i->bneg return BINT negation
	442	$i->babs return BINT absolute value
	443	$i->bcmp(BINT) return CODE compare numbers (undef,<0,=0,>0)
	444	$i->badd(BINT) return BINT addition
	445	$i->bsub(BINT) return BINT subtraction
	446	$i->bmul(BINT) return BINT multiplication
	447	$i->bdiv(BINT) return (BINT,BINT) division (quo,rem) just quo if scalar
	448	$i->bmod(BINT) return BINT modulus
	449	$i->bgcd(BINT) return BINT greatest common divisor
	450	$i->bnorm return BINT normalization
e16b8f49	451	$i->blsft(BINT) return BINT left shift
	452	$i->brsft(BINT) return (BINT,BINT) right shift (quo,rem) just quo if scalar
	453	$i->band(BINT) return BINT bit-wise and
	454	$i->bior(BINT) return BINT bit-wise inclusive or
	455	$i->bxor(BINT) return BINT bit-wise exclusive or
	456	$i->bnot return BINT bit-wise not
a5f75d66	457
	458	=head1 DESCRIPTION
	459
	460	All basic math operations are overloaded if you declare your big
	461	integers as
	462
	463	$i = new Math::BigInt '123 456 789 123 456 789';
	464
	465
	466	=over 2
	467
	468	=item Canonical notation
	469
	470	Big integer value are strings of the form C</^[+-]\d+$/> with leading
	471	zeros suppressed.
	472
	473	=item Input
	474
	475	Input values to these routines may be strings of the form
	476	C</^\s*[+-]?[\d\s]+$/>.
	477
	478	=item Output
	479
	480	Output values always always in canonical form
	481
	482	=back
	483
	484	Actual math is done in an internal format consisting of an array
	485	whose first element is the sign (/^[+-]$/) and whose remaining
	486	elements are base 100000 digits with the least significant digit first.
	487	The string 'NaN' is used to represent the result when input arguments
	488	are not numbers, as well as the result of dividing by zero.
	489
	490	=head1 EXAMPLES
	491
	492	'+0' canonical zero value
	493	' -123 123 123' canonical value '-123123123'
	494	'1 23 456 7890' canonical value '+1234567890'
	495
	496
b3ac6de7	497	=head1 Autocreating constants
	498
	499	After C<use Math::BigInt ':constant'> all the integer decimal constants
e3c7ef20	500	in the given scope are converted to C<Math::BigInt>. This conversion
b3ac6de7	501	happens at compile time.
	502
	503	In particular
	504
	505	perl -MMath::BigInt=:constant -e 'print 2**100'
	506
8dcee03e	507	print the integer value of C<2**100>. Note that without conversion of
8dcee03e	508	constants the expression 2**100 will be calculated as floating point number.
b3ac6de7	509
a5f75d66	510	=head1 BUGS
	511
	512	The current version of this module is a preliminary version of the
	513	real thing that is currently (as of perl5.002) under development.
	514
	515	=head1 AUTHOR
	516
	517	Mark Biggar, overloaded interface by Ilya Zakharevich.
	518
	519	=cut