calc.pm

source of perl for linux application,

    }
  my $rem = $src % $BASE_LEN;		# remainder to shift
  $src = int($src / $BASE_LEN);		# source
  if ($rem == 0)
    {
    splice (@$x,0,$src);		# even faster, 38.4 => 39.3
    }
  else
    {
    my $len = scalar @$x - $src;	# elems to go
    my $vd; my $z = '0'x $BASE_LEN;
    $x->[scalar @$x] = 0;		# avoid || 0 test inside loop
    while ($dst < $len)
      {
      $vd = $z.$x->[$src];
      $vd = substr($vd,-$BASE_LEN,$BASE_LEN-$rem);
      $src++;
      $vd = substr($z.$x->[$src],-$rem,$rem) . $vd;
      $vd = substr($vd,-$BASE_LEN,$BASE_LEN) if length($vd) > $BASE_LEN;
      $x->[$dst] = int($vd);
      $dst++;
      }
    splice (@$x,$dst) if $dst > 0;		# kill left-over array elems
    pop @$x if $x->[-1] == 0 && @$x > 1;	# kill last element if 0
    } # else rem == 0
  $x;
  }

sub _lsft
  {
  my ($c,$x,$y,$n) = @_;

  if ($n != 10)
    {
    $n = _new($c,$n); return _mul($c,$x, _pow($c,$n,$y));
    }

  # shortcut (faster) for shifting by 10) since we are in base 10eX
  # multiples of $BASE_LEN:
  my $src = scalar @$x;			# source
  my $len = _num($c,$y);		# shift-len as normal int
  my $rem = $len % $BASE_LEN;		# remainder to shift
  my $dst = $src + int($len/$BASE_LEN);	# destination
  my $vd;				# further speedup
  $x->[$src] = 0;			# avoid first ||0 for speed
  my $z = '0' x $BASE_LEN;
  while ($src >= 0)
    {
    $vd = $x->[$src]; $vd = $z.$vd;
    $vd = substr($vd,-$BASE_LEN+$rem,$BASE_LEN-$rem);
    $vd .= $src > 0 ? substr($z.$x->[$src-1],-$BASE_LEN,$rem) : '0' x $rem;
    $vd = substr($vd,-$BASE_LEN,$BASE_LEN) if length($vd) > $BASE_LEN;
    $x->[$dst] = int($vd);
    $dst--; $src--;
    }
  # set lowest parts to 0
  while ($dst >= 0) { $x->[$dst--] = 0; }
  # fix spurios last zero element
  splice @$x,-1 if $x->[-1] == 0;
  $x;
  }

sub _pow
  {
  # power of $x to $y
  # ref to array, ref to array, return ref to array
  my ($c,$cx,$cy) = @_;

  if (scalar @$cy == 1 && $cy->[0] == 0)
    {
    splice (@$cx,1); $cx->[0] = 1;		# y == 0 => x => 1
    return $cx;
    }
  if ((scalar @$cx == 1 && $cx->[0] == 1) ||	#    x == 1
      (scalar @$cy == 1 && $cy->[0] == 1))	# or y == 1
    {
    return $cx;
    }
  if (scalar @$cx == 1 && $cx->[0] == 0)
    {
    splice (@$cx,1); $cx->[0] = 0;		# 0 ** y => 0 (if not y <= 0)
    return $cx;
    }

  my $pow2 = _one();

  my $y_bin = _as_bin($c,$cy); $y_bin =~ s/^0b//;
  my $len = length($y_bin);
  while (--$len > 0)
    {
    _mul($c,$pow2,$cx) if substr($y_bin,$len,1) eq '1';		# is odd?
    _mul($c,$cx,$cx);
    }

  _mul($c,$cx,$pow2);
  $cx;
  }

sub _nok
  {
  # n over k
  # ref to array, return ref to array
  my ($c,$n,$k) = @_;

  # ( 7 )    7!          7*6*5 * 4*3*2*1   7 * 6 * 5
  # ( - ) = --------- =  --------------- = ---------
  # ( 3 )   3! (7-3)!    3*2*1 * 4*3*2*1   3 * 2 * 1 

  # compute n - k + 2 (so we start with 5 in the example above)
  my $x = _copy($c,$n);

  _sub($c,$n,$k);
  if (!_is_one($c,$n))
    {
    _inc($c,$n);
    my $f = _copy($c,$n); _inc($c,$f);		# n = 5, f = 6, d = 2
    my $d = _two($c);
    while (_acmp($c,$f,$x) <= 0)		# f < n ?
      {
      # n = (n * f / d) == 5 * 6 / 2 => n == 3
      $n = _mul($c,$n,$f); $n = _div($c,$n,$d);
      # f = 7, d = 3
      _inc($c,$f); _inc($c,$d);
      }
    }
  else 
    {
    # keep ref to $n and set it to 1
    splice (@$n,1); $n->[0] = 1;
    }
  $n;
  }

my @factorials = (
  1,
  1,
  2,
  2*3,
  2*3*4,
  2*3*4*5,
  2*3*4*5*6,
  2*3*4*5*6*7,
);

sub _fac
  {
  # factorial of $x
  # ref to array, return ref to array
  my ($c,$cx) = @_;

  if ((@$cx == 1) && ($cx->[0] <= 7))
    {
    $cx->[0] = $factorials[$cx->[0]];		# 0 => 1, 1 => 1, 2 => 2 etc.
    return $cx;
    }

  if ((@$cx == 1) && 		# we do this only if $x >= 12 and $x <= 7000
      ($cx->[0] >= 12 && $cx->[0] < 7000))
    {

  # Calculate (k-j) * (k-j+1) ... k .. (k+j-1) * (k + j)
  # See http://blogten.blogspot.com/2007/01/calculating-n.html
  # The above series can be expressed as factors:
  #   k * k - (j - i) * 2
  # We cache k*k, and calculate (j * j) as the sum of the first j odd integers

  # This will not work when N exceeds the storage of a Perl scalar, however,
  # in this case the algorithm would be way to slow to terminate, anyway.

  # As soon as the last element of $cx is 0, we split it up and remember
  # how many zeors we got so far. The reason is that n! will accumulate
  # zeros at the end rather fast.
  my $zero_elements = 0;

  # If n is even, set n = n -1
  my $k = _num($c,$cx); my $even = 1;
  if (($k & 1) == 0)
    {
    $even = $k; $k --;
    }
  # set k to the center point
  $k = ($k + 1) / 2;
#  print "k $k even: $even\n";
  # now calculate k * k
  my $k2 = $k * $k;
  my $odd = 1; my $sum = 1;
  my $i = $k - 1;
  # keep reference to x
  my $new_x = _new($c, $k * $even);
  @$cx = @$new_x;
  if ($cx->[0] == 0)
    {
    $zero_elements ++; shift @$cx;
    }
#  print STDERR "x = ", _str($c,$cx),"\n";
  my $BASE2 = int(sqrt($BASE))-1;
  my $j = 1; 
  while ($j <= $i)
    {
    my $m = ($k2 - $sum); $odd += 2; $sum += $odd; $j++;
    while ($j <= $i && ($m < $BASE2) && (($k2 - $sum) < $BASE2))
      {
      $m *= ($k2 - $sum);
      $odd += 2; $sum += $odd; $j++;
#      print STDERR "\n k2 $k2 m $m sum $sum odd $odd\n"; sleep(1);
      }
    if ($m < $BASE)
      {
      _mul($c,$cx,[$m]);
      }
    else
      {
      _mul($c,$cx,$c->_new($m));
      }
    if ($cx->[0] == 0)
      {
      $zero_elements ++; shift @$cx;
      }
#    print STDERR "Calculate $k2 - $sum = $m (x = ", _str($c,$cx),")\n";
    }
  # multiply in the zeros again
  unshift @$cx, (0) x $zero_elements; 
  return $cx;
  }

  # go forward until $base is exceeded
  # limit is either $x steps (steps == 100 means a result always too high) or
  # $base.
  my $steps = 100; $steps = $cx->[0] if @$cx == 1;
  my $r = 2; my $cf = 3; my $step = 2; my $last = $r;
  while ($r*$cf < $BASE && $step < $steps)
    {
    $last = $r; $r *= $cf++; $step++;
    }
  if ((@$cx == 1) && $step == $cx->[0])
    {
    # completely done, so keep reference to $x and return
    $cx->[0] = $r;
    return $cx;
    }
  
  # now we must do the left over steps
  my $n;					# steps still to do
  if (scalar @$cx == 1)
    {
    $n = $cx->[0];
    }
  else
    {
    $n = _copy($c,$cx);
    }

  # Set $cx to the last result below $BASE (but keep ref to $x)
  $cx->[0] = $last; splice (@$cx,1);
  # As soon as the last element of $cx is 0, we split it up and remember
  # how many zeors we got so far. The reason is that n! will accumulate
  # zeros at the end rather fast.
  my $zero_elements = 0;

  # do left-over steps fit into a scalar?
  if (ref $n eq 'ARRAY')
    {
    # No, so use slower inc() & cmp()
    # ($n is at least $BASE here)
    my $base_2 = int(sqrt($BASE)) - 1;
    #print STDERR "base_2: $base_2\n"; 
    while ($step < $base_2)
      {
      if ($cx->[0] == 0)
        {
        $zero_elements ++; shift @$cx;
        }
      my $b = $step * ($step + 1); $step += 2;
      _mul($c,$cx,[$b]);
      }
    $step = [$step];
    while (_acmp($c,$step,$n) <= 0)
      {
      if ($cx->[0] == 0)
        {
        $zero_elements ++; shift @$cx;
        }
      _mul($c,$cx,$step); _inc($c,$step);
      }
    }
  else
    {
    # Yes, so we can speed it up slightly
  
#    print "# left over steps $n\n";

    my $base_4 = int(sqrt(sqrt($BASE))) - 2;
    #print STDERR "base_4: $base_4\n";
    my $n4 = $n - 4; 
    while ($step < $n4 && $step < $base_4)
      {
      if ($cx->[0] == 0)
        {
        $zero_elements ++; shift @$cx;
        }
      my $b = $step * ($step + 1); $step += 2; $b *= $step * ($step + 1); $step += 2;
      _mul($c,$cx,[$b]);
      }
    my $base_2 = int(sqrt($BASE)) - 1;
    my $n2 = $n - 2; 
    #print STDERR "base_2: $base_2\n"; 
    while ($step < $n2 && $step < $base_2)
      {
      if ($cx->[0] == 0)
        {
        $zero_elements ++; shift @$cx;
        }
      my $b = $step * ($step + 1); $step += 2;
      _mul($c,$cx,[$b]);
      }
    # do what's left over
    while ($step <= $n)
      {
      _mul($c,$cx,[$step]); $step++;
      if ($cx->[0] == 0)
        {
        $zero_elements ++; shift @$cx;
        }
      }
    }
  # multiply in the zeros again
  unshift @$cx, (0) x $zero_elements;
  $cx;			# return result
  }

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

sub _log_int
  {
  # calculate integer log of $x to base $base
  # ref to array, ref to array - return ref to array
  my ($c,$x,$base) = @_;

  # X == 0 => NaN
  return if (scalar @$x == 1 && $x->[0] == 0);
  # BASE 0 or 1 => NaN
  return if (scalar @$base == 1 && $base->[0] < 2);
  my $cmp = _acmp($c,$x,$base); # X == BASE => 1
  if ($cmp == 0)
    {
    splice (@$x,1); $x->[0] = 1;
    return ($x,1)
    }
  # X < BASE
  if ($cmp < 0)
    {
    splice (@$x,1); $x->[0] = 0;
    return ($x,undef);
    }

  my $x_org = _copy($c,$x);		# preserve x
  splice(@$x,1); $x->[0] = 1;		# keep ref to $x

  # Compute a guess for the result based on:
  # $guess = int ( length_in_base_10(X) / ( log(base) / log(10) ) )
  my $len = _len($c,$x_org);
  my $log = log($base->[-1]) / log(10);

  # for each additional element in $base, we add $BASE_LEN to the result,
  # based on the observation that log($BASE,10) is BASE_LEN and
  # log(x*y) == log(x) + log(y):
  $log += ((scalar @$base)-1) * $BASE_LEN;

  # calculate now a guess based on the values obtained above:
  my $res = int($len / $log);

  $x->[0] = $res;
  my $trial = _pow ($c, _copy($c, $base), $x);
  my $a = _acmp($c,$trial,$x_org);

#  print STDERR "# trial ", _str($c,$x)," was: $a (0 = exact, -1 too small, +1 too big)\n";

  # found an exact result?
  return ($x,1) if $a == 0;

  if ($a > 0)
    {
    # or too big
    _div($c,$trial,$base); _dec($c, $x);
    while (($a = _acmp($c,$trial,$x_org)) > 0)
      {
#      print STDERR "# big _log_int at ", _str($c,$x), "\n"; 
      _div($c,$trial,$base); _dec($c, $x);
      }
    # result is now exact (a == 0), or too small (a < 0)
    return ($x, $a == 0 ? 1 : 0);
    }

  # else: result was to small
  _mul($c,$trial,$base);

  # did we now get the right result?
  $a = _acmp($c,$trial,$x_org);

  if ($a == 0)				# yes, exactly
    {
    _inc($c, $x);
    return ($x,1); 
    }
  return ($x,0) if $a > 0;  

  # Result still too small (we should come here only if the estimate above
  # was very off base):
 
  # Now let the normal trial run obtain the real result
  # Simple loop that increments $x by 2 in each step, possible overstepping
  # the real result

  my $base_mul = _mul($c, _copy($c,$base), $base);	# $base * $base

  while (($a = _acmp($c,$trial,$x_org)) < 0)
    {
#    print STDERR "# small _log_int at ", _str($c,$x), "\n"; 
    _mul($c,$trial,$base_mul); _add($c, $x, [2]);
    }

  my $exact = 1;
  if ($a > 0)
    {
    # overstepped the result
    _dec($c, $x);
    _div($c,$trial,$base);
    $a = _acmp($c,$trial,$x_org);
    if ($a > 0)
      {
      _dec($c, $x);
      }
    $exact = 0 if $a != 0;		# a = -1 => not exact result, a = 0 => exact
    }
  
  ($x,$exact);				# return result
  }

# for debugging:
  use constant DEBUG => 0;
  my $steps = 0;
  sub steps { $steps };

sub _sqrt
  {
  # square-root of $x in place
  # Compute a guess of the result (by rule of thumb), then improve it via
  # Newton's method.
  my ($c,$x) = @_;

  if (scalar @$x == 1)
    {
    # fits into one Perl scalar, so result can be computed directly
    $x->[0] = int(sqrt($x->[0]));
    return $x;
    } 
  my $y = _copy($c,$x);
  # hopefully _len/2 is < $BASE, the -1 is to always undershot the guess
  # since our guess will "grow"
  my $l = int((_len($c,$x)-1) / 2);	

  my $lastelem = $x->[-1];					# for guess
  my $elems = scalar @$x - 1;
  # not enough digits, but could have more?
  if ((length($lastelem) <= 3) && ($elems > 1))
    {
    # right-align with zero pad
    my $len = length($lastelem) & 1;
    print "$lastelem => " if DEBUG;
    $lastelem .= substr($x->[-2] . '0' x $BASE_LEN,0,$BASE_LEN);
    # former odd => make odd again, or former even to even again
    $lastelem = $lastelem / 10 if (length($lastelem) & 1) != $len;
    print "$lastelem\n" if DEBUG;
    }

  # construct $x (instead of _lsft($c,$x,$l,10)
  my $r = $l % $BASE_LEN;	# 10000 00000 00000 00000 ($BASE_LEN=5)
  $l = int($l / $BASE_LEN);