calc.pm

source of perl for linux application,

  print "l =  $l " if DEBUG;

  splice @$x,$l;		# keep ref($x), but modify it

  # we make the first part of the guess not '1000...0' but int(sqrt($lastelem))
  # that gives us:
  # 14400 00000 => sqrt(14400) => guess first digits to be 120
  # 144000 000000 => sqrt(144000) => guess 379

  print "$lastelem (elems $elems) => " if DEBUG;
  $lastelem = $lastelem / 10 if ($elems & 1 == 1);		# odd or even?
  my $g = sqrt($lastelem); $g =~ s/\.//;			# 2.345 => 2345
  $r -= 1 if $elems & 1 == 0;					# 70 => 7

  # padd with zeros if result is too short
  $x->[$l--] = int(substr($g . '0' x $r,0,$r+1));
  print "now ",$x->[-1] if DEBUG;
  print " would have been ", int('1' . '0' x $r),"\n" if DEBUG;

  # If @$x > 1, we could compute the second elem of the guess, too, to create
  # an even better guess. Not implemented yet. Does it improve performance?
  $x->[$l--] = 0 while ($l >= 0);	# all other digits of guess are zero

  print "start x= ",_str($c,$x),"\n" if DEBUG;
  my $two = _two();
  my $last = _zero();
  my $lastlast = _zero();
  $steps = 0 if DEBUG;
  while (_acmp($c,$last,$x) != 0 && _acmp($c,$lastlast,$x) != 0)
    {
    $steps++ if DEBUG;
    $lastlast = _copy($c,$last);
    $last = _copy($c,$x);
    _add($c,$x, _div($c,_copy($c,$y),$x));
    _div($c,$x, $two );
    print " x= ",_str($c,$x),"\n" if DEBUG;
    }
  print "\nsteps in sqrt: $steps, " if DEBUG;
  _dec($c,$x) if _acmp($c,$y,_mul($c,_copy($c,$x),$x)) < 0;	# overshot? 
  print " final ",$x->[-1],"\n" if DEBUG;
  $x;
  }

sub _root
  {
  # take n'th root of $x in place (n >= 3)
  my ($c,$x,$n) = @_;
 
  if (scalar @$x == 1)
    {
    if (scalar @$n > 1)
      {
      # result will always be smaller than 2 so trunc to 1 at once
      $x->[0] = 1;
      }
    else
      {
      # fits into one Perl scalar, so result can be computed directly
      # cannot use int() here, because it rounds wrongly (try 
      # (81 ** 3) ** (1/3) to see what I mean)
      #$x->[0] = int( $x->[0] ** (1 / $n->[0]) );
      # round to 8 digits, then truncate result to integer
      $x->[0] = int ( sprintf ("%.8f", $x->[0] ** (1 / $n->[0]) ) );
      }
    return $x;
    } 

  # we know now that X is more than one element long

  # if $n is a power of two, we can repeatedly take sqrt($X) and find the
  # proper result, because sqrt(sqrt($x)) == root($x,4)
  my $b = _as_bin($c,$n);
  if ($b =~ /0b1(0+)$/)
    {
    my $count = CORE::length($1);	# 0b100 => len('00') => 2
    my $cnt = $count;			# counter for loop
    unshift (@$x, 0);			# add one element, together with one
					# more below in the loop this makes 2
    while ($cnt-- > 0)
      {
      # 'inflate' $X by adding one element, basically computing
      # $x * $BASE * $BASE. This gives us more $BASE_LEN digits for result
      # since len(sqrt($X)) approx == len($x) / 2.
      unshift (@$x, 0);
      # calculate sqrt($x), $x is now one element to big, again. In the next
      # round we make that two, again.
      _sqrt($c,$x);
      }
    # $x is now one element to big, so truncate result by removing it
    splice (@$x,0,1);
    } 
  else
    {
    # trial computation by starting with 2,4,8,16 etc until we overstep
    my $step;
    my $trial = _two();

    # while still to do more than X steps
    do
      {
      $step = _two();
      while (_acmp($c, _pow($c, _copy($c, $trial), $n), $x) < 0)
        {
        _mul ($c, $step, [2]);
        _add ($c, $trial, $step);
        }

      # hit exactly?
      if (_acmp($c, _pow($c, _copy($c, $trial), $n), $x) == 0)
        {
        @$x = @$trial;			# make copy while preserving ref to $x
        return $x;
        }
      # overstepped, so go back on step
      _sub($c, $trial, $step);
      } while (scalar @$step > 1 || $step->[0] > 128);

    # reset step to 2
    $step = _two();
    # add two, because $trial cannot be exactly the result (otherwise we would
    # alrady have found it)
    _add($c, $trial, $step);
 
    # and now add more and more (2,4,6,8,10 etc)
    while (_acmp($c, _pow($c, _copy($c, $trial), $n), $x) < 0)
      {
      _add ($c, $trial, $step);
      }

    # hit not exactly? (overstepped)
    if (_acmp($c, _pow($c, _copy($c, $trial), $n), $x) > 0)
      {
      _dec($c,$trial);
      }

    # hit not exactly? (overstepped)
    # 80 too small, 81 slightly too big, 82 too big
    if (_acmp($c, _pow($c, _copy($c, $trial), $n), $x) > 0)
      {
      _dec ($c, $trial); 
      }

    @$x = @$trial;			# make copy while preserving ref to $x
    return $x;
    }
  $x; 
  }

##############################################################################
# binary stuff

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

  # the shortcut makes equal, large numbers _really_ fast, and makes only a
  # very small performance drop for small numbers (e.g. something with less
  # than 32 bit) Since we optimize for large numbers, this is enabled.
  return $x if _acmp($c,$x,$y) == 0;		# shortcut
  
  my $m = _one(); my ($xr,$yr);
  my $mask = $AND_MASK;

  my $x1 = $x;
  my $y1 = _copy($c,$y);			# make copy
  $x = _zero();
  my ($b,$xrr,$yrr);
  use integer;
  while (!_is_zero($c,$x1) && !_is_zero($c,$y1))
    {
    ($x1, $xr) = _div($c,$x1,$mask);
    ($y1, $yr) = _div($c,$y1,$mask);

    # make ints() from $xr, $yr
    # this is when the AND_BITS are greater than $BASE and is slower for
    # small (<256 bits) numbers, but faster for large numbers. Disabled
    # due to KISS principle

#    $b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
#    $b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
#    _add($c,$x, _mul($c, _new( $c, ($xrr & $yrr) ), $m) );
    
    # 0+ due to '&' doesn't work in strings
    _add($c,$x, _mul($c, [ 0+$xr->[0] & 0+$yr->[0] ], $m) );
    _mul($c,$m,$mask);
    }
  $x;
  }

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

  return _zero() if _acmp($c,$x,$y) == 0;	# shortcut (see -and)

  my $m = _one(); my ($xr,$yr);
  my $mask = $XOR_MASK;

  my $x1 = $x;
  my $y1 = _copy($c,$y);			# make copy
  $x = _zero();
  my ($b,$xrr,$yrr);
  use integer;
  while (!_is_zero($c,$x1) && !_is_zero($c,$y1))
    {
    ($x1, $xr) = _div($c,$x1,$mask);
    ($y1, $yr) = _div($c,$y1,$mask);
    # make ints() from $xr, $yr (see _and())
    #$b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
    #$b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
    #_add($c,$x, _mul($c, _new( $c, ($xrr ^ $yrr) ), $m) );

    # 0+ due to '^' doesn't work in strings
    _add($c,$x, _mul($c, [ 0+$xr->[0] ^ 0+$yr->[0] ], $m) );
    _mul($c,$m,$mask);
    }
  # the loop stops when the shorter of the two numbers is exhausted
  # the remainder of the longer one will survive bit-by-bit, so we simple
  # multiply-add it in
  _add($c,$x, _mul($c, $x1, $m) ) if !_is_zero($c,$x1);
  _add($c,$x, _mul($c, $y1, $m) ) if !_is_zero($c,$y1);
  
  $x;
  }

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

  return $x if _acmp($c,$x,$y) == 0;		# shortcut (see _and)

  my $m = _one(); my ($xr,$yr);
  my $mask = $OR_MASK;

  my $x1 = $x;
  my $y1 = _copy($c,$y);			# make copy
  $x = _zero();
  my ($b,$xrr,$yrr);
  use integer;
  while (!_is_zero($c,$x1) && !_is_zero($c,$y1))
    {
    ($x1, $xr) = _div($c,$x1,$mask);
    ($y1, $yr) = _div($c,$y1,$mask);
    # make ints() from $xr, $yr (see _and())
#    $b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
#    $b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
#    _add($c,$x, _mul($c, _new( $c, ($xrr | $yrr) ), $m) );
    
    # 0+ due to '|' doesn't work in strings
    _add($c,$x, _mul($c, [ 0+$xr->[0] | 0+$yr->[0] ], $m) );
    _mul($c,$m,$mask);
    }
  # the loop stops when the shorter of the two numbers is exhausted
  # the remainder of the longer one will survive bit-by-bit, so we simple
  # multiply-add it in
  _add($c,$x, _mul($c, $x1, $m) ) if !_is_zero($c,$x1);
  _add($c,$x, _mul($c, $y1, $m) ) if !_is_zero($c,$y1);
  
  $x;
  }

sub _as_hex
  {
  # convert a decimal number to hex (ref to array, return ref to string)
  my ($c,$x) = @_;

  # fits into one element (handle also 0x0 case)
  return sprintf("0x%x",$x->[0]) if @$x == 1;

  my $x1 = _copy($c,$x);

  my $es = '';
  my ($xr, $h, $x10000);
  if ($] >= 5.006)
    {
    $x10000 = [ 0x10000 ]; $h = 'h4';
    }
  else
    {
    $x10000 = [ 0x1000 ]; $h = 'h3';
    }
  while (@$x1 != 1 || $x1->[0] != 0)		# _is_zero()
    {
    ($x1, $xr) = _div($c,$x1,$x10000);
    $es .= unpack($h,pack('V',$xr->[0]));
    }
  $es = reverse $es;
  $es =~ s/^[0]+//;   # strip leading zeros
  '0x' . $es;					# return result prepended with 0x
  }

sub _as_bin
  {
  # convert a decimal number to bin (ref to array, return ref to string)
  my ($c,$x) = @_;

  # fits into one element (and Perl recent enough), handle also 0b0 case
  # handle zero case for older Perls
  if ($] <= 5.005 && @$x == 1 && $x->[0] == 0)
    {
    my $t = '0b0'; return $t;
    }
  if (@$x == 1 && $] >= 5.006)
    {
    my $t = sprintf("0b%b",$x->[0]);
    return $t;
    }
  my $x1 = _copy($c,$x);

  my $es = '';
  my ($xr, $b, $x10000);
  if ($] >= 5.006)
    {
    $x10000 = [ 0x10000 ]; $b = 'b16';
    }
  else
    {
    $x10000 = [ 0x1000 ]; $b = 'b12';
    }
  while (!(@$x1 == 1 && $x1->[0] == 0))		# _is_zero()
    {
    ($x1, $xr) = _div($c,$x1,$x10000);
    $es .= unpack($b,pack('v',$xr->[0]));
    }
  $es = reverse $es;
  $es =~ s/^[0]+//;   # strip leading zeros
  '0b' . $es;					# return result prepended with 0b
  }

sub _as_oct
  {
  # convert a decimal number to octal (ref to array, return ref to string)
  my ($c,$x) = @_;

  # fits into one element (handle also 0 case)
  return sprintf("0%o",$x->[0]) if @$x == 1;

  my $x1 = _copy($c,$x);

  my $es = '';
  my $xr;
  my $x1000 = [ 0100000 ];
  while (@$x1 != 1 || $x1->[0] != 0)		# _is_zero()
    {
    ($x1, $xr) = _div($c,$x1,$x1000);
    $es .= reverse sprintf("%05o", $xr->[0]);
    }
  $es = reverse $es;
  $es =~ s/^[0]+//;   # strip leading zeros
  '0' . $es;					# return result prepended with 0
  }

sub _from_oct
  {
  # convert a octal number to decimal (string, return ref to array)
  my ($c,$os) = @_;

  # for older Perls, play safe
  my $m = [ 0100000 ];
  my $d = 5;					# 5 digits at a time

  my $mul = _one();
  my $x = _zero();

  my $len = int( (length($os)-1)/$d );		# $d digit parts, w/o the '0'
  my $val; my $i = -$d;
  while ($len >= 0)
    {
    $val = substr($os,$i,$d);			# get oct digits
    $val = CORE::oct($val);
    $i -= $d; $len --;
    my $adder = [ $val ];
    _add ($c, $x, _mul ($c, $adder, $mul ) ) if $val != 0;
    _mul ($c, $mul, $m ) if $len >= 0; 		# skip last mul
    }
  $x;
  }

sub _from_hex
  {
  # convert a hex number to decimal (string, return ref to array)
  my ($c,$hs) = @_;

  my $m = _new($c, 0x10000000);			# 28 bit at a time (<32 bit!)
  my $d = 7;					# 7 digits at a time
  if ($] <= 5.006)
    {
    # for older Perls, play safe
    $m = [ 0x10000 ];				# 16 bit at a time (<32 bit!)
    $d = 4;					# 4 digits at a time
    }

  my $mul = _one();
  my $x = _zero();

  my $len = int( (length($hs)-2)/$d );		# $d digit parts, w/o the '0x'
  my $val; my $i = -$d;
  while ($len >= 0)
    {
    $val = substr($hs,$i,$d);			# get hex digits
    $val =~ s/^0x// if $len == 0;		# for last part only because
    $val = CORE::hex($val);			# hex does not like wrong chars
    $i -= $d; $len --;
    my $adder = [ $val ];
    # if the resulting number was to big to fit into one element, create a
    # two-element version (bug found by Mark Lakata - Thanx!)
    if (CORE::length($val) > $BASE_LEN)
      {
      $adder = _new($c,$val);
      }
    _add ($c, $x, _mul ($c, $adder, $mul ) ) if $val != 0;
    _mul ($c, $mul, $m ) if $len >= 0; 		# skip last mul
    }
  $x;
  }

sub _from_bin
  {
  # convert a hex number to decimal (string, return ref to array)
  my ($c,$bs) = @_;

  # instead of converting X (8) bit at a time, it is faster to "convert" the
  # number to hex, and then call _from_hex.

  my $hs = $bs;
  $hs =~ s/^[+-]?0b//;					# remove sign and 0b
  my $l = length($hs);					# bits
  $hs = '0' x (8-($l % 8)) . $hs if ($l % 8) != 0;	# padd left side w/ 0
  my $h = '0x' . unpack('H*', pack ('B*', $hs));	# repack as hex
  
  $c->_from_hex($h);
  }

##############################################################################
# special modulus functions

sub _modinv
  {
  # modular inverse
  my ($c,$x,$y) = @_;

  my $u = _zero($c); my $u1 = _one($c);
  my $a = _copy($c,$y); my $b = _copy($c,$x);

  # Euclid's Algorithm for bgcd(), only that we calc bgcd() ($a) and the
  # result ($u) at the same time. See comments in BigInt for why this works.
  my $q;
  ($a, $q, $b) = ($b, _div($c,$a,$b));		# step 1
  my $sign = 1;
  while (!_is_zero($c,$b))
    {
    my $t = _add($c, 				# step 2:
       _mul($c,_copy($c,$u1), $q) ,		#  t =  u1 * q
       $u );					#     + u
    $u = $u1;					#  u = u1, u1 = t
    $u1 = $t;
    $sign = -$sign;
    ($a, $q, $b) = ($b, _div($c,$a,$b));	# step 1
    }

  # if the gcd is not 1, then return NaN
  return (undef,undef) unless _is_one($c,$a);
 
  ($u1, $sign == 1 ? '+' : '-');
  }

sub _modpow
  {
  # modulus of power ($x ** $y) % $z
  my ($c,$num,$exp,$mod) = @_;

  # in the trivial case,
  if (_is_one($c,$mod))
    {
    splice @$num,0,1; $num->[0] = 0;
    return $num;
    }
  if ((scalar @$num == 1) && (($nu