bigint.pm

1. 记录每个帖子的访问人情况

  return
   wantarray ? ($x->round(@r),$self->bzero(@r)):$x->round(@r) if $x->is_zero();
 
  # Is $x in the interval [0, $y) (aka $x <= $y) ?
  my $cmp = $CALC->_acmp($x->{value},$y->{value});
  if (($cmp < 0) and (($x->{sign} eq $y->{sign}) or !wantarray))
    {
    return $upgrade->bdiv($upgrade->new($x),$upgrade->new($y),@r)
     if defined $upgrade;

    return $x->bzero()->round(@r) unless wantarray;
    my $t = $x->copy();      # make copy first, because $x->bzero() clobbers $x
    return ($x->bzero()->round(@r),$t);
    }
  elsif ($cmp == 0)
    {
    # shortcut, both are the same, so set to +/- 1
    $x->__one( ($x->{sign} ne $y->{sign} ? '-' : '+') ); 
    return $x unless wantarray;
    return ($x->round(@r),$self->bzero(@r));
    }
  return $upgrade->bdiv($upgrade->new($x),$upgrade->new($y),@r)
   if defined $upgrade;
   
  # calc new sign and in case $y == +/- 1, return $x
  my $xsign = $x->{sign};				# keep
  $x->{sign} = ($x->{sign} ne $y->{sign} ? '-' : '+'); 
  # check for / +-1 (cant use $y->is_one due to '-'
  if ($CALC->_is_one($y->{value}))
    {
    return wantarray ? ($x->round(@r),$self->bzero(@r)) : $x->round(@r); 
    }

  if (wantarray)
    {
    my $rem = $self->bzero(); 
    ($x->{value},$rem->{value}) = $CALC->_div($x->{value},$y->{value});
    $x->{sign} = '+' if $CALC->_is_zero($x->{value});
    $rem->{_a} = $x->{_a};
    $rem->{_p} = $x->{_p};
    $x->round(@r); 
    if (! $CALC->_is_zero($rem->{value}))
      {
      $rem->{sign} = $y->{sign};
      $rem = $y-$rem if $xsign ne $y->{sign};	# one of them '-'
      }
    else
      {
      $rem->{sign} = '+';			# dont leave -0
      }
    return ($x,$rem->round(@r));
    }

  $x->{value} = $CALC->_div($x->{value},$y->{value});
  $x->{sign} = '+' if $CALC->_is_zero($x->{value});

  $x->round(@r) if !exists $x->{_f} || $x->{_f} & MB_NEVER_ROUND == 0;
  $x;
  }

###############################################################################
# modulus functions

sub bmod 
  {
  # modulus (or remainder)
  # (BINT or num_str, BINT or num_str) return BINT
  
  # set up parameters
  my ($self,$x,$y,@r) = (ref($_[0]),@_);
  # objectify is costly, so avoid it
  if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
    {
    ($self,$x,$y,@r) = objectify(2,@_);
    }

  return $x if $x->modify('bmod');
  $r[3] = $y;					# no push!
  if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/) || $y->is_zero())
    {
    my ($d,$r) = $self->_div_inf($x,$y);
    $x->{sign} = $r->{sign};
    $x->{value} = $r->{value};
    return $x->round(@r);
    }

  if ($CALC->can('_mod'))
    {
    # calc new sign and in case $y == +/- 1, return $x
    $x->{value} = $CALC->_mod($x->{value},$y->{value});
    if (!$CALC->_is_zero($x->{value}))
      {
      my $xsign = $x->{sign};
      $x->{sign} = $y->{sign};
      if ($xsign ne $y->{sign})
        {
        my $t = $CALC->_copy($x->{value});		# copy $x
        $x->{value} = $CALC->_copy($y->{value});	# copy $y to $x
        $x->{value} = $CALC->_sub($y->{value},$t,1); 	# $y-$x
        }
      }
    else
      {
      $x->{sign} = '+';				# dont leave -0
      }
    $x->round(@r) if !exists $x->{_f} || $x->{_f} & MB_NEVER_ROUND == 0;
    return $x;
    }
  my ($t,$rem) = $self->bdiv($x->copy(),$y,@r);	# slow way (also rounds)
  # modify in place
  foreach (qw/value sign _a _p/)
    {
    $x->{$_} = $rem->{$_};
    }
  $x;
  }

sub bmodinv
  {
  # Modular inverse.  given a number which is (hopefully) relatively
  # prime to the modulus, calculate its inverse using Euclid's
  # alogrithm.  If the number is not relatively prime to the modulus
  # (i.e. their gcd is not one) then NaN is returned.

  # set up parameters
  my ($self,$x,$y,@r) = (ref($_[0]),@_);
  # objectify is costly, so avoid it
  if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
    {
    ($self,$x,$y,@r) = objectify(2,@_);
    }

  return $x if $x->modify('bmodinv');

  return $x->bnan()
        if ($y->{sign} ne '+'                           # -, NaN, +inf, -inf
         || $x->is_zero()                               # or num == 0
         || $x->{sign} !~ /^[+-]$/                      # or num NaN, inf, -inf
        );

  # put least residue into $x if $x was negative, and thus make it positive
  $x->bmod($y) if $x->{sign} eq '-';

  if ($CALC->can('_modinv'))
    {
    my $sign;
    ($x->{value},$sign) = $CALC->_modinv($x->{value},$y->{value});
    $x->bnan() if !defined $x->{value};                 # in case no GCD found
    return $x if !defined $sign;                        # already real result
    $x->{sign} = $sign;                                 # flip/flop see below
    $x->bmod($y);                                       # calc real result
    return $x;
    }
  my ($u, $u1) = ($self->bzero(), $self->bone());
  my ($a, $b) = ($y->copy(), $x->copy());

  # first step need always be done since $num (and thus $b) is never 0
  # Note that the loop is aligned so that the check occurs between #2 and #1
  # thus saving us one step #2 at the loop end. Typical loop count is 1. Even
  # a case with 28 loops still gains about 3% with this layout.
  my $q;
  ($a, $q, $b) = ($b, $a->bdiv($b));                    # step #1
  # Euclid's Algorithm (calculate GCD of ($a,$b) in $a and also calculate
  # two values in $u and $u1, we use only $u1 afterwards)
  my $sign = 1;                                         # flip-flop
  while (!$b->is_zero())                                # found GCD if $b == 0
    {
    # the original algorithm had:
    # ($u, $u1) = ($u1, $u->bsub($u1->copy()->bmul($q))); # step #2
    # The following creates exact the same sequence of numbers in $u1,
    # except for the sign ($u1 is now always positive). Since formerly
    # the sign of $u1 was alternating between '-' and '+', the $sign
    # flip-flop will take care of that, so that at the end of the loop
    # we have the real sign of $u1. Keeping numbers positive gains us
    # speed since badd() is faster than bsub() and makes it possible
    # to have the algorithmn in Calc for even more speed.

    ($u, $u1) = ($u1, $u->badd($u1->copy()->bmul($q))); # step #2
    $sign = - $sign;                                    # flip sign

    ($a, $q, $b) = ($b, $a->bdiv($b));                  # step #1 again
    }

  # If the gcd is not 1, then return NaN! It would be pointless to
  # have called bgcd to check this first, because we would then be
  # performing the same Euclidean Algorithm *twice*.
  return $x->bnan() unless $a->is_one();

  $u1->bneg() if $sign != 1;                            # need to flip?

  $u1->bmod($y);                                        # calc result
  $x->{value} = $u1->{value};                           # and copy over to $x
  $x->{sign} = $u1->{sign};                             # to modify in place
  $x;
  }

sub bmodpow
  {
  # takes a very large number to a very large exponent in a given very
  # large modulus, quickly, thanks to binary exponentation.  supports
  # negative exponents.
  my ($self,$num,$exp,$mod,@r) = objectify(3,@_);

  return $num if $num->modify('bmodpow');

  # check modulus for valid values
  return $num->bnan() if ($mod->{sign} ne '+'		# NaN, - , -inf, +inf
                       || $mod->is_zero());

  # check exponent for valid values
  if ($exp->{sign} =~ /\w/) 
    {
    # i.e., if it's NaN, +inf, or -inf...
    return $num->bnan();
    }

  $num->bmodinv ($mod) if ($exp->{sign} eq '-');

  # check num for valid values (also NaN if there was no inverse but $exp < 0)
  return $num->bnan() if $num->{sign} !~ /^[+-]$/;

  if ($CALC->can('_modpow'))
    {
    # $mod is positive, sign on $exp is ignored, result also positive
    $num->{value} = $CALC->_modpow($num->{value},$exp->{value},$mod->{value});
    return $num;
    }

  # in the trivial case,
  return $num->bzero(@r) if $mod->is_one();
  return $num->bone('+',@r) if $num->is_zero() or $num->is_one();

  # $num->bmod($mod);           # if $x is large, make it smaller first
  my $acc = $num->copy();	# but this is not really faster...

  $num->bone(); # keep ref to $num

  my $expbin = $exp->as_bin(); $expbin =~ s/^[-]?0b//; # ignore sign and prefix
  my $len = length($expbin);
  while (--$len >= 0)
    {
    if( substr($expbin,$len,1) eq '1')
      {
      $num->bmul($acc)->bmod($mod);
      }
    $acc->bmul($acc)->bmod($mod);
    }

  $num;
  }

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

sub bfac
  {
  # (BINT or num_str, BINT or num_str) return BINT
  # compute factorial numbers
  # modifies first argument
  my ($self,$x,@r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);

  return $x if $x->modify('bfac');
 
  return $x->bnan() if $x->{sign} ne '+';	# inf, NnN, <0 etc => NaN
  return $x->bone('+',@r) if $x->is_zero() || $x->is_one();	# 0 or 1 => 1

  if ($CALC->can('_fac'))
    {
    $x->{value} = $CALC->_fac($x->{value});
    return $x->round(@r);
    }

  my $n = $x->copy();
  $x->bone();
  # seems we need not to temp. clear A/P of $x since the result is the same
  my $f = $self->new(2);
  while ($f->bacmp($n) < 0)
    {
    $x->bmul($f); $f->binc();
    }
  $x->bmul($f,@r);			# last step and also round
  }
 
sub bpow 
  {
  # (BINT or num_str, BINT or num_str) return BINT
  # compute power of two numbers -- stolen from Knuth Vol 2 pg 233
  # modifies first argument
  
  # set up parameters
  my ($self,$x,$y,@r) = (ref($_[0]),@_);
  # objectify is costly, so avoid it
  if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
    {
    ($self,$x,$y,@r) = objectify(2,@_);
    }

  return $x if $x->modify('bpow');

  return $upgrade->bpow($upgrade->new($x),$y,@r)
   if defined $upgrade && !$y->isa($self);

  $r[3] = $y;					# no push!
  return $x if $x->{sign} =~ /^[+-]inf$/;	# -inf/+inf ** x
  return $x->bnan() if $x->{sign} eq $nan || $y->{sign} eq $nan;
  return $x->bone('+',@r) if $y->is_zero();
  return $x->round(@r) if $x->is_one() || $y->is_one();
  if ($x->{sign} eq '-' && $CALC->_is_one($x->{value}))
    {
    # if $x == -1 and odd/even y => +1/-1
    return $y->is_odd() ? $x->round(@r) : $x->babs()->round(@r);
    # my Casio FX-5500L has a bug here: -1 ** 2 is -1, but -1 * -1 is 1;
    }
  # 1 ** -y => 1 / (1 ** |y|)
  # so do test for negative $y after above's clause
  return $x->bnan() if $y->{sign} eq '-';
  return $x->round(@r) if $x->is_zero();  # 0**y => 0 (if not y <= 0)

  if ($CALC->can('_pow'))
    {
    $x->{value} = $CALC->_pow($x->{value},$y->{value});
    $x->round(@r) if !exists $x->{_f} || $x->{_f} & MB_NEVER_ROUND == 0;
    return $x;
    }

# based on the assumption that shifting in base 10 is fast, and that mul
# works faster if numbers are small: we count trailing zeros (this step is
# O(1)..O(N), but in case of O(N) we save much more time due to this),
# stripping them out of the multiplication, and add $count * $y zeros
# afterwards like this:
# 300 ** 3 == 300*300*300 == 3*3*3 . '0' x 2 * 3 == 27 . '0' x 6
# creates deep recursion since brsft/blsft use bpow sometimes.
#  my $zeros = $x->_trailing_zeros();
#  if ($zeros > 0)
#    {
#    $x->brsft($zeros,10);	# remove zeros
#    $x->bpow($y);		# recursion (will not branch into here again)
#    $zeros = $y * $zeros; 	# real number of zeros to add
#    $x->blsft($zeros,10);
#    return $x->round(@r);
#    }

  my $pow2 = $self->__one();
  my $y_bin = $y->as_bin(); $y_bin =~ s/^0b//;
  my $len = length($y_bin);
  while (--$len > 0)
    {
    $pow2->bmul($x) if substr($y_bin,$len,1) eq '1';	# is odd?
    $x->bmul($x);
    }
  $x->bmul($pow2);
  $x->round(@r) if !exists $x->{_f} || $x->{_f} & MB_NEVER_ROUND == 0;
  $x;
  }

sub blsft 
  {
  # (BINT or num_str, BINT or num_str) return BINT
  # compute x << y, base n, y >= 0
 
  # set up parameters
  my ($self,$x,$y,$n,@r) = (ref($_[0]),@_);
  # objectify is costly, so avoid it
  if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
    {
    ($self,$x,$y,$n,@r) = objectify(2,@_);
    }

  return $x if $x->modify('blsft');
  return $x->bnan() if ($x->{sign} !~ /^[+-]$/ || $y->{sign} !~ /^[+-]$/);
  return $x->round(@r) if $y->is_zero();

  $n = 2 if !defined $n; return $x->bnan() if $n <= 0 || $y->{sign} eq '-';

  my $t; $t = $CALC->_lsft($x->{value},$y->{value},$n) if $CALC->can('_lsft');
  if (defined $t)
    {
    $x->{value} = $t; return $x->round(@r);
    }
  # fallback
  return $x->bmul( $self->bpow($n, $y, @r), @r );
  }

sub brsft 
  {
  # (BINT or num_str, BINT or num_str) return BINT
  # compute x >> y, base n, y >= 0
  
  # set up parameters
  my ($self,$x,$y,$n,@r) = (ref($_[0]),@_);
  # objectify is costly, so avoid it
  if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
    {
    ($self,$x,$y,$n,@r) = objectify(2,@_);
    }

  return $x if $x->modify('brsft');
  return $x->bnan() if ($x->{sign} !~ /^[+-]$/ || $y->{sign} !~ /^[+-]$/);
  return $x->round(@r) if $y->is_zero();
  return $x->bzero(@r) if $x->is_zero();		# 0 => 0

  $n = 2 if !defined $n; return $x->bnan() if $n <= 0 || $y->{sign} eq '-';

   # this only works for negative numbers when shifting in base 2
  if (($x->{sign} eq '-') && ($n == 2))
    {
    return $x->round(@r) if $x->is_one('-');	# -1 => -1
    if (!$y->is_one())
      {
      # although this is O(N*N) in calc (as_bin!) it is O(N) in Pari et al
      # but perhaps there is a better emulation for two's complement shift...
      # if $y != 1, we must simulate it by doing:
      # convert to bin, flip all bits, shift, and be done
      $x->binc();			# -3 => -2
      my $bin = $x->as_bin();
      $bin =~ s/^-0b//;			# strip '-0b' prefix
      $bin =~ tr/10/01/;		# flip bits
      # now shift
      if (CORE::length($bin) <= $y)
        {
	$bin = '0'; 			# shifting to far right creates -1
					# 0, because later increment makes 
					# that 1, attached '-' makes it '-1'
					# because -1 >> x == -1 !
        } 
      else
	{
	$bin =~ s/.{$y}$//;		# cut off at the right side
        $bin = '1' . $bin;		# extend left side by one dummy '1'
        $bin =~ tr/10/01/;		# flip bits back
	}
      my $res = $self->new('0b'.$bin);	# add prefix and convert back
      $res->binc();			# remember to increment
      $x->{value} = $res->{value};	# take over value
      return $x->round(@r);		# we are done now, magic, isn't?
      }
    $x->bdec();				# n == 2, but $y == 1: this fixes it
    }

  my $t; $t = $CALC->_rsft($x->{value},$y->{value},$n) if $CALC->can('_rsft');
  if (defined $t)
    {
    $x->{value} = $t;
    return $x->round(@r);
    }
  # fallback
  $x->bdiv($self->bpow($n,$y, @r), @r);
  $x;
  }

sub band 
  {
  #(BINT or num_str, BINT or num_str) return BINT
  # compute x & y
 
  # set up parameters
  my ($self,$x,$y,@r) = (ref($_[0]),@_);