calc.pm

网页留言本,比一般的留言簿管用

package Math::BigInt::Calc;

use 5.005;
use strict;
# use warnings;	# dont use warnings for older Perls

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

$VERSION = '0.32';

# 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 1eX where X is determined
# automatically at loading time to be the maximum possible value

# todo:
# - fully remove funky $# stuff (maybe)

# 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). Other platforms
# BS2000, some Crays need USE_DIV instead.
# The BEGIN block is used to determine which of the two variants gives the
# correct result.

##############################################################################
# global constants, flags and accessory
 
# constants for easier life
my $nan = 'NaN';
my ($MBASE,$BASE,$RBASE,$BASE_LEN,$MAX_VAL,$BASE_LEN2,$BASE_LEN_SMALL);
my ($AND_BITS,$XOR_BITS,$OR_BITS);
my ($AND_MASK,$XOR_MASK,$OR_MASK);
my ($LEN_CONVERT);

sub _base_len 
  {
  # set/get the BASE_LEN and assorted other, connected values
  # used only be the testsuite, set is used only by the BEGIN block below
  shift;

  my $b = shift;
  if (defined $b)
    {
    # find whether we can use mul or div or none in mul()/div()
    # (in last case reduce BASE_LEN_SMALL)
    $BASE_LEN_SMALL = $b+1;
    my $caught = 0;
    while (--$BASE_LEN_SMALL > 5)
      {
      $MBASE = int("1e".$BASE_LEN_SMALL);
      $RBASE = abs('1e-'.$BASE_LEN_SMALL);		# see USE_MUL
      $caught = 0;
      $caught += 1 if (int($MBASE * $RBASE) != 1);	# should be 1
      $caught += 2 if (int($MBASE / $MBASE) != 1);	# should be 1
      last if $caught != 3;
      }
    # BASE_LEN is used for anything else than mul()/div()
    $BASE_LEN = $BASE_LEN_SMALL;
    $BASE_LEN = shift if (defined $_[0]);		# one more arg?
    $BASE = int("1e".$BASE_LEN);

    $BASE_LEN2 = int($BASE_LEN_SMALL / 2);		# for mul shortcut
    $MBASE = int("1e".$BASE_LEN_SMALL);
    $RBASE = abs('1e-'.$BASE_LEN_SMALL);		# see USE_MUL
    $MAX_VAL = $MBASE-1;
    $LEN_CONVERT = 0;
    $LEN_CONVERT = 1 if $BASE_LEN_SMALL != $BASE_LEN;

    #print "BASE_LEN: $BASE_LEN MAX_VAL: $MAX_VAL BASE: $BASE RBASE: $RBASE ";
    #print "BASE_LEN_SMALL: $BASE_LEN_SMALL MBASE: $MBASE\n";

    undef &_mul;
    undef &_div;

    if ($caught & 1 != 0)
      {
      # must USE_MUL
      *{_mul} = \&_mul_use_mul;
      *{_div} = \&_div_use_mul;
      }
    else		# $caught must be 2, since it can't be 1 nor 3
      {
      # can USE_DIV instead
      *{_mul} = \&_mul_use_div;
      *{_div} = \&_div_use_div;
      }
    }
  return $BASE_LEN unless wantarray;
  return ($BASE_LEN, $AND_BITS, $XOR_BITS, $OR_BITS, $BASE_LEN_SMALL, $MAX_VAL);
  }

BEGIN
  {
  # from Daniel Pfeiffer: determine largest group of digits that is precisely
  # multipliable with itself plus carry
  # Test now changed to expect the proper pattern, not a result off by 1 or 2
  my ($e, $num) = 3;	# lowest value we will use is 3+1-1 = 3
  do 
    {
    $num = ('9' x ++$e) + 0;
    $num *= $num + 1.0;
    } while ("$num" =~ /9{$e}0{$e}/);	# must be a certain pattern
  $e--; 				# last test failed, so retract one step
  # the limits below brush the problems with the test above under the rug:
  # the test should be able to find the proper $e automatically
  $e = 5 if $^O =~ /^uts/;	# UTS get's some special treatment
  $e = 5 if $^O =~ /^unicos/;	# unicos is also problematic (6 seems to work
				# there, but we play safe)
  $e = 5 if $] < 5.006;		# cap, for older Perls
  $e = 7 if $e > 7;		# cap, for VMS, OS/390 and other 64 bit systems
				# 8 fails inside random testsuite, so take 7

  # determine how many digits fit into an integer and can be safely added 
  # together plus carry w/o causing an overflow

  # this below detects 15 on a 64 bit system, because after that it becomes
  # 1e16  and not 1000000 :/ I can make it detect 18, but then I get a lot of
  # test failures. Ugh! (Tomake detect 18: uncomment lines marked with *)
  use integer;
  my $bi = 5;			# approx. 16 bit
  $num = int('9' x $bi);
  # $num = 99999; # *
  # while ( ($num+$num+1) eq '1' . '9' x $bi)	# *
  while ( int($num+$num+1) eq '1' . '9' x $bi)
    {
    $bi++; $num = int('9' x $bi);
    # $bi++; $num *= 10; $num += 9;	# *
    }
  $bi--;				# back off one step
  # by setting them equal, we ignore the findings and use the default
  # one-size-fits-all approach from former versions
  $bi = $e;				# XXX, this should work always

  __PACKAGE__->_base_len($e,$bi);	# set and store

  # find out how many bits _and, _or and _xor can take (old default = 16)
  # I don't think anybody has yet 128 bit scalars, so let's play safe.
  local $^W = 0;	# don't warn about 'nonportable number'
  $AND_BITS = 15; $XOR_BITS = 15; $OR_BITS  = 15;

  # find max bits, we will not go higher than numberofbits that fit into $BASE
  # to make _and etc simpler (and faster for smaller, slower for large numbers)
  my $max = 16;
  while (2 ** $max < $BASE) { $max++; }
  {
    no integer;
    $max = 16 if $] < 5.006;	# older Perls might not take >16 too well
  }
  my ($x,$y,$z);
  do {
    $AND_BITS++;
    $x = oct('0b' . '1' x $AND_BITS); $y = $x & $x;
    $z = (2 ** $AND_BITS) - 1;
    } while ($AND_BITS < $max && $x == $z && $y == $x);
  $AND_BITS --;						# retreat one step
  do {
    $XOR_BITS++;
    $x = oct('0b' . '1' x $XOR_BITS); $y = $x ^ 0;
    $z = (2 ** $XOR_BITS) - 1;
    } while ($XOR_BITS < $max && $x == $z && $y == $x);
  $XOR_BITS --;						# retreat one step
  do {
    $OR_BITS++;
    $x = oct('0b' . '1' x $OR_BITS); $y = $x | $x;
    $z = (2 ** $OR_BITS) - 1;
    } while ($OR_BITS < $max && $x == $z && $y == $x);
  $OR_BITS --;						# retreat one step
  
  }

##############################################################################
# convert between the "small" and the "large" representation

sub _to_large
  {
  # take an array in base $BASE_LEN_SMALL and convert it in-place to $BASE_LEN
  my ($c,$x) = @_;

#  print "_to_large $BASE_LEN_SMALL => $BASE_LEN\n";

  return $x if $LEN_CONVERT == 0 ||		# nothing to converconvertor
	 @$x == 1;				# only one element => early out
  
  #     12345    67890    12345    67890   contents
  # to      3        2        1        0   index 
  #             123456  7890123  4567890   contents 

#  # faster variant
#  my @d; my $str = '';
#  my $z = '0' x $BASE_LEN_SMALL;
#  foreach (@$x)
#    {
#    # ... . 04321 . 000321
#    $str = substr($z.$_,-$BASE_LEN_SMALL,$BASE_LEN_SMALL) . $str;
#    if (length($str) > $BASE_LEN)
#      {
#      push @d, substr($str,-$BASE_LEN,$BASE_LEN);	# extract one piece
#      substr($str,-$BASE_LEN,$BASE_LEN) = '';		# remove it
#      }
#    }
#  push @d, $str if $str !~ /^0*$/;			# extract last piece
#  @$x = @d;
#  $x->[-1] = int($x->[-1]);			# strip leading zero
#  $x;

  my $ret = "";
  my $l = scalar @$x;		# number of parts
  $l --; $ret .= int($x->[$l]); $l--;
  my $z = '0' x ($BASE_LEN_SMALL-1);                            
  while ($l >= 0)
    {
    $ret .= substr($z.$x->[$l],-$BASE_LEN_SMALL); 
    $l--;
    }
  my $str = _new($c,\$ret);			# make array
  @$x = @$str;					# clobber contents of $x
  $x->[-1] = int($x->[-1]);			# strip leading zero
  }

sub _to_small
  {
  # take an array in base $BASE_LEN and convert it in-place to $BASE_LEN_SMALL
  my ($c,$x) = @_;

  return $x if $LEN_CONVERT == 0;		# nothing to do
  return $x if @$x == 1 && length(int($x->[0])) <= $BASE_LEN_SMALL;

  my $d = _str($c,$x);
  my $il = length($$d)-1;
  ## this leaves '00000' instead of int 0 and will be corrected after any op
  # clobber contents of $x
  @$x = reverse(unpack("a" . ($il % $BASE_LEN_SMALL+1) 
      . ("a$BASE_LEN_SMALL" x ($il / $BASE_LEN_SMALL)), $$d));	

  $x->[-1] = int($x->[-1]);			# strip leading zero
  }

###############################################################################

sub _new
  {
  # (ref to string) return ref to num_array
  # Convert a number from string format (without sign) to internal base
  # 1ex format. Assumes normalized value as input.
  my $d = $_[1];
  my $il = length($$d)-1;
  # this leaves '00000' instead of int 0 and will be corrected after any op
  [ reverse(unpack("a" . ($il % $BASE_LEN+1) 
    . ("a$BASE_LEN" x ($il / $BASE_LEN)), $$d)) ];
  }                                                                             
  
BEGIN
  {
  $AND_MASK = __PACKAGE__->_new( \( 2 ** $AND_BITS ));
  $XOR_MASK = __PACKAGE__->_new( \( 2 ** $XOR_BITS ));
  $OR_MASK = __PACKAGE__->_new( \( 2 ** $OR_BITS ));
  }

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

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

sub _two
  {
  # create a two (used internally for shifting)
  [ 2 ];
  }

sub _copy
  {
  [ @{$_[1]} ];
  }

# catch and throw away
sub import { }

##############################################################################
# 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")
  my $ar = $_[1];
  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 .= int($ar->[$l]); $l--;
  # Interestingly, the pre-padd method uses more time
  # the old grep variant takes longer (14 to 10 sec)
  my $z = '0' x ($BASE_LEN-1);                            
  while ($l >= 0)
    {
    $ret .= substr($z.$ar->[$l],-$BASE_LEN); # fastest way I could think of
    $l--;
    }
  \$ret;
  }                                                                             

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

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

sub _add
  {
  # (ref to int_num_array, ref to int_num_array)
  # routine to add two base 1eX numbers
  # stolen from Knuth Vol 2 Algorithm A pg 231
  # there are separate routines to add and sub as per Knuth pg 233
  # This routine clobbers up array x, but not y.
 
  my ($c,$x,$y) = @_;

  return $x if (@$y == 1) && $y->[0] == 0;		# $x + 0 => $x
  if ((@$x == 1) && $x->[0] == 0)			# 0 + $y => $y->copy
    {
    # twice as slow as $x = [ @$y ], but necc. to retain $x as ref :(
    @$x = @$y; return $x;		
    }
 
  # 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 of 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++;
    }
  $x;
  }                                                                             

sub _inc
  {
  # (ref to int_num_array, ref to int_num_array)
  # routine to add 1 to a base 1eX numbers
  # This routine modifies array x
  my ($c,$x) = @_;

  for my $i (@$x)
    {
    return $x if (($i += 1) < $BASE);		# early out
    $i = 0;					# overflow, next
    }
  push @$x,1 if ($x->[-1] == 0);		# last overflowed, so extend
  $x;
  }                                                                             

sub _dec
  {
  # (ref to int_num_array, ref to int_num_array)
  # routine to add 1 to a base 1eX numbers
  # This routine modifies array x
  my ($c,$x) = @_;

  my $MAX = $BASE-1;				# since MAX_VAL based on MBASE
  for my $i (@$x)
    {
    last if (($i -= 1) >= 0);			# early out
    $i = $MAX;					# overflow, next
    }
  pop @$x if $x->[-1] == 0 && @$x > 1;		# last overflowed (but leave 0)
  $x;
  }                                                                             

sub _sub
  {
  # (ref to int_num_array, ref to int_num_array, swap)
  # subtract base 1eX numbers -- stolen from Knuth Vol 2 pg 232, $x > $y
  # subtract Y from X by modifying x in place
  my ($c,$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;
      $i += $BASE if $car = (($i -= ($sy->[$j] || 0) + $car) < 0); $j++;
      }
    # might leave leading zeros, so fix that
    return __strip_zeros($sx);
    }
  #print "case 1 (swap)\n";
  for $i (@$sx)
    {
    # we can't do an early out if $x is < than $y, since we
    # need to copy the high chunks from $y. Found by Bob Mathews.
    #last unless defined $sy->[$j] || $car;
    $sy->[$j] += $BASE
     if $car = (($sy->[$j] = $i-($sy->[$j]||0) - $car) < 0);
    $j++;
    }
  # might leave leading zeros, so fix that
  __strip_zeros($sy);
  }                                                                             

sub _mul_use_mul
  {
  # (ref to int_num_array, ref to int_num_array)
  # multiply two numbers in internal representation
  # modifies first arg, second need not be different from first
  my ($c,$xv,$yv) = @_;

  # shortcut for two very short numbers (improved by Nathan Zook)
  # works also if xv and yv are the same reference
  if ((@$xv == 1) && (@$yv == 1))
    {
    if (($xv->[0] *= $yv->[0]) >= $MBASE)