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

package Math::BigInt;

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;


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