[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;


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