Module: Flt::Num::AuxiliarFunctions

Included in:: Flt::Num, Flt::Num
Defined in:: lib/flt/num.rb

Constant Summary

EXP_INC =

LOG_PREC_INC =

LOG_RADIX_INC =

LOG_RADIX_EXTRA =

LOG2_MULT = TODO: K=100? K=64? ...

LOG2_LB_CORRECTION = (1..15).map{|i| (LOG2_MULT*Math.log(16.0/i)/Math.log(2)).ceil}

[ # (1..15).map{|i| (LOG2_MULT*Math.log(16.0/i)/Math.log(2)).ceil}
  400, 300, 242, 200, 168, 142, 120, 100, 84, 68, 55, 42, 30, 20, 10
  # for LOG2_MULT=64: 256, 192, 155, 128, 108, 91, 77, 64, 54, 44, 35, 27, 20, 13, 6
]

LOG10_MULT =

LOG10_LB_CORRECTION = (1..9).map_hash{|i| LOG10_MULT - (LOG10_MULT*Math.log10(i)).floor}

{ # (1..9).map_hash{|i| LOG10_MULT - (LOG10_MULT*Math.log10(i)).floor}
  '1'=> 100, '2'=> 70, '3'=> 53, '4'=> 40, '5'=> 31,
  '6'=> 23, '7'=> 16, '8'=> 10, '9'=> 5
}

Class Attribute Summary (collapse)

+ (Object) log_radix_digits readonly

Returns the value of attribute log_radix_digits.

Instance Method Summary (collapse)

- (Object) _convert(x, error = true)

Convert a numeric value to decimal (internal use).
- (Object) _div_nearest(a, b)

Closest integer to a/b, a and b positive integers; rounds to even in the case of a tie.
- (Object) _exp(c, e, p)

Compute an approximation to exp(c*radix**e), with p decimal places of precision.
- (Object) _iexp(x, m, l = 8)

Given integers x and M, M > 0, such that x/M is small in absolute value, compute an integer approximation to M*exp(x/M).
- (Object) _ilog(x, m, l = 8)

Integer approximation to M*log(x/M), with absolute error boundable in terms only of x/M.
- (Object) _log(c, e, p)

Given integers c, e and p with c > 0, compute an integer approximation to radix**p * log(c*radix**e), with an absolute error of at most 1.
- (Object) _log_radix_digits(p)

Given an integer p >= 0, return floor(radix**p)*log(radix).
- (Object) _log_radix_lb(c)
- (Object) _log_radix_mult
- (Object) _normalize(op1, op2, prec = 0)

Normalizes op1, op2 to have the same exp and length of coefficient.
- (Object) _number_of_digits(v)
- (Object) _parser(txt)

Parse numeric text literals (internal use).
- (Object) _power(xc, xe, yc, ye, p)

Given integers xc, xe, yc and ye representing Num x = xc*radix**xe and y = yc*radix**ye, compute x**y.
- (Object) _rshift_nearest(x, shift)

Given an integer x and a nonnegative integer shift, return closest integer to x / 2**shift; use round-to-even in case of a tie.
- (Object) _sqrt_nearest(n, a)

Closest integer to the square root of the positive integer n.
- (Object) log10_lb(c)

Compute a lower bound for LOG10_MULT*log10(c) for a positive integer c.
- (Object) log2_lb(c)

Compute a lower bound for LOG2_MULT*log10(c) for a positive integer c.

Class Attribute Details

+ (`Object`) log_radix_digits (readonly)

Returns the value of attribute log_radix_digits



4298
4299
4300

# File 'lib/flt/num.rb', line 4298

def log_radix_digits
  @log_radix_digits
end

Instance Method Details

- (`Object`) _convert(x, error = true)

Convert a numeric value to decimal (internal use)

# File 'lib/flt/num.rb', line 3993

def _convert(x, error=true)
  case x
  when num_class
    x
  when *num_class.context.coercible_types
    num_class.new(x)
  else
    raise TypeError, "Unable to convert #{x.class} to #{num_class}" if error
    nil
  end
end

- (`Object`) _div_nearest(a, b)

Closest integer to a/b, a and b positive integers; rounds to even in the case of a tie.

# File 'lib/flt/num.rb', line 4287

def _div_nearest(a, b)
  q, r = a.divmod(b)
  q + (((2*r + (q&1)) > b) ? 1 : 0)
end

- (`Object`) _exp(c, e, p)

Compute an approximation to exp(c*radix**e), with p decimal places of precision. Returns integers d, f such that:

radix**(p-1) <= d <= radix**p, and
(d-1)*radix**f < exp(c*radix**e) < (d+1)*radix**f

In other words, d*radix**f is an approximation to exp(c*radix**e) with p digits of precision, and with an error in d of at most 1. This is almost, but not quite, the same as the error being < 1ulp: when d

radix**(p-1) the error could be up to radix ulp.

# File 'lib/flt/num.rb', line 4090

def _exp(c, e, p)
    # we'll call iexp with M = radix**(p+2), giving p+3 digits of precision
    p += EXP_INC

    # compute log(radix) with extra precision = adjusted exponent of c*radix**e
    # TODO: without the .abs tests fail because c is negative: c should not be negative!!
    extra = [0, e + _number_of_digits(c.abs) - 1].max
    q = p + extra

    # compute quotient c*radix**e/(log(radix)) = c*radix**(e+q)/(log(radix)*radix**q),
    # rounding down
    shift = e+q
    if shift >= 0
        cshift = c*num_class.int_radix_power(shift)
    else
        cshift = c/num_class.int_radix_power(-shift)
    end
    quot, rem = cshift.divmod(_log_radix_digits(q))

    # reduce remainder back to original precision
    rem = _div_nearest(rem, num_class.int_radix_power(extra))

    # for radix=10: error in result of _iexp < 120;  error after division < 0.62
    r = _div_nearest(_iexp(rem, num_class.int_radix_power(p)), num_class.int_radix_power(EXP_INC+1)), quot - p + (EXP_INC+1)
    return r
end

- (`Object`) _iexp(x, m, l = 8)

Given integers x and M, M > 0, such that x/M is small in absolute value, compute an integer approximation to M*exp(x/M). For redix=10, and 0 <= x/M <= 2.4, the absolute error in the result is bounded by 60 (and is usually much smaller).

# File 'lib/flt/num.rb', line 4171

def _iexp(x, m, l=8)

    # Algorithm: to compute exp(z) for a real number z, first divide z
    # by a suitable power R of 2 so that |z/2**R| < 2**-L.  Then
    # compute expm1(z/2**R) = exp(z/2**R) - 1 using the usual Taylor
    # series
    #
    #     expm1(x) = x + x**2/2! + x**3/3! + ...
    #
    # Now use the identity
    #
    #     expm1(2x) = expm1(x)*(expm1(x)+2)
    #
    # R times to compute the sequence expm1(z/2**R),
    # expm1(z/2**(R-1)), ... , exp(z/2), exp(z).

    # Find R such that x/2**R/M <= 2**-L
    r = _nbits((x<<l)/m)

    # Taylor series.  (2**L)**T > M
    t = -(-num_class.radix*_number_of_digits(m)/(3*l)).to_i
    y = _div_nearest(x, t)
    mshift = m<<r
    (1...t).to_a.reverse.each do |i|
        y = _div_nearest(x*(mshift + y), mshift * i)
    end

    # Expansion
    (0...r).to_a.reverse.each do |k|
        mshift = m<<(k+2)
        y = _div_nearest(y*(y+mshift), mshift)
    end

    return m+y
end

- (`Object`) _ilog(x, m, l = 8)

Integer approximation to M*log(x/M), with absolute error boundable in terms only of x/M.

Given positive integers x and M, return an integer approximation to M * log(x/M). For radix=10, L = 8 and 0.1 <= x/M <= 10 the difference between the approximation and the exact result is at most 22. For L = 8 and 1.0 <= x/M <= 10.0 the difference is at most 15. In both cases these are upper bounds on the error; it will usually be much smaller.

# File 'lib/flt/num.rb', line 4216

def _ilog(x, m, l = 8)
  # The basic algorithm is the following: let log1p be the function
  # log1p(x) = log(1+x).  Then log(x/M) = log1p((x-M)/M).  We use
  # the reduction
  #
  #    log1p(y) = 2*log1p(y/(1+sqrt(1+y)))
  #
  # repeatedly until the argument to log1p is small (< 2**-L in
  # absolute value).  For small y we can use the Taylor series
  # expansion
  #
  #    log1p(y) ~ y - y**2/2 + y**3/3 - ... - (-y)**T/T
  #
  # truncating at T such that y**T is small enough.  The whole
  # computation is carried out in a form of fixed-point arithmetic,
  # with a real number z being represented by an integer
  # approximation to z*M.  To avoid loss of precision, the y below
  # is actually an integer approximation to 2**R*y*M, where R is the
  # number of reductions performed so far.

  y = x-m
  # argument reduction; R = number of reductions performed
  r = 0
  # while (r <= l && y.abs << l-r >= m ||
  #        r > l and y.abs>> r-l >= m)
  while (((r <= l) && ((y.abs << (l-r)) >= m)) ||
         ((r > l) && ((y.abs>>(r-l)) >= m)))
      y = _div_nearest((m*y) << 1,
                       m + _sqrt_nearest(m*(m+_rshift_nearest(y, r)), m))
      r += 1
  end

  # Taylor series with T terms
  t = -(-10*_number_of_digits(m)/(3*l)).to_i
  yshift = _rshift_nearest(y, r)
  w = _div_nearest(m, t)
  # (1...t).reverse_each do |k| # Ruby 1.9
  (1...t).to_a.reverse.each do |k|
     w = _div_nearest(m, k) - _div_nearest(yshift*w, m)
  end

  return _div_nearest(w*y, m)
end

- (`Object`) _log(c, e, p)

Given integers c, e and p with c > 0, compute an integer approximation to radix**p * log(c*radix**e), with an absolute error of at most 1. Assumes that c*radix**e is not exactly 1.

# File 'lib/flt/num.rb', line 4121

def _log(c, e, p)
    # Increase precision by 2. The precision increase is compensated
    # for at the end with a division
    p += LOG_PREC_INC

    # rewrite c*radix**e as d*radix**f with either f >= 0 and 1 <= d <= radix,
    # or f <= 0 and 1/radix <= d <= 1.  Then we can compute radix**p * log(c*radix**e)
    # as radix**p * log(d) + radix**p*f * log(radix).
    l = _number_of_digits(c)
    f = e+l - ((e+l >= 1) ? 1 : 0)

    # compute approximation to radix**p*log(d), with error < 27 for radix=10
    if p > 0
        k = e+p-f
        if k >= 0
            c *= num_class.int_radix_power(k)
        else
            c = _div_nearest(c, num_class.int_radix_power(-k))  # error of <= 0.5 in c for radix=10
        end

        # _ilog magnifies existing error in c by a factor of at most radix
        log_d = _ilog(c, num_class.int_radix_power(p)) # error < 5 + 22 = 27 for radix=10
    else
        # p <= 0: just approximate the whole thing by 0; error < 2.31 for radix=10
        log_d = 0
    end

    # compute approximation to f*radix**p*log(radix), with error < 11 for radix=10.
    if f
        extra = _number_of_digits(f.abs) - 1
        if p + extra >= 0
            # for radix=10:
            # error in f * _log10_digits(p+extra) < |f| * 1 = |f|
            # after division, error < |f|/10**extra + 0.5 < 10 + 0.5 < 11
            f_log_r = _div_nearest(f*_log_radix_digits(p+extra), num_class.int_radix_power(extra))
        else
            f_log_r = 0
        end
    else
        f_log_r = 0
    end

    # error in sum < 11+27 = 38; error after division < 0.38 + 0.5 < 1 for radix=10
    return _div_nearest(f_log_r + log_d, num_class.int_radix_power(LOG_PREC_INC)) # extra radix factor for base 2 ???
end

- (`Object`) _log_radix_digits(p)

Given an integer p >= 0, return floor(radix**p)*log(radix).

Raises:

(ArgumentError)

# File 'lib/flt/num.rb', line 4304

def _log_radix_digits(p)
  # digits are stored as a string, for quick conversion to
  # integer in the case that we've already computed enough
  # digits; the stored digits should always bge correct
  # (truncated, not rounded to nearest).
  raise ArgumentError, "p should be nonnegative" if p<0
  stored_digits = (AuxiliarFunctions.log_radix_digits[num_class.radix] || "")
  if p >= stored_digits.length
      digits = nil
      # compute p+3, p+6, p+9, ... digits; continue until at
      # least one of the extra digits is nonzero
      extra = LOG_RADIX_EXTRA
      loop do
        # compute p+extra digits, correct to within 1ulp
        m = num_class.int_radix_power(p+extra+LOG_RADIX_INC)
        digits = _div_nearest(_ilog(num_class.radix*m, m), num_class.int_radix_power(LOG_RADIX_INC)).to_s(num_class.radix)
        break if digits[-extra..-1] != '0'*extra
        extra += LOG_RADIX_EXTRA
      end
      # if the radix < e (i.e. only for radix==2), we must prefix with a 0 because log(radix)<1
      # BUT THIS REDUCES PRECISION BY ONE? : may be avoid prefix and adjust scaling in the caller
      prefix = num_class.radix==2 ? '0' : ''
      # keep all reliable digits so far; remove trailing zeros
      # and next nonzero digit
      AuxiliarFunctions.log_radix_digits[num_class.radix] = prefix + digits.sub(/0*$/,'')[0...-1]
  end
  return (AuxiliarFunctions.log_radix_digits[num_class.radix][0..p]).to_i(num_class.radix)
end

- (`Object`) _log_radix_lb(c)

# File 'lib/flt/num.rb', line 4368

def _log_radix_lb(c)
  case num_class.radix
  when 10
    log10_lb(c)
  when 2
    log2_lb(c)
  else
    raise ArgumentError, "_log_radix_lb not implemented for base #{num_class.radix}"
  end
end

- (`Object`) _log_radix_mult

# File 'lib/flt/num.rb', line 4357

def _log_radix_mult
  case num_class.radix
  when 10
    LOG10_MULT
  when 2
    LOG2_MULT
  else
    raise ArgumentError, "_log_radix_mult not implemented for base #{num_class.radix}"
  end
end

- (`Object`) _normalize(op1, op2, prec = 0)

Normalizes op1, op2 to have the same exp and length of coefficient. Used for addition.

# File 'lib/flt/num.rb', line 4015

def _normalize(op1, op2, prec=0)
  if op1.exponent < op2.exponent
    swap = true
    tmp,other = op2,op1
  else
    swap = false
    tmp,other = op1,op2
  end
  tmp_len = tmp.number_of_digits
  other_len = other.number_of_digits
  exp = tmp.exponent + [-1, tmp_len - prec - 2].min
  if (other_len+other.exponent-1 < exp) && prec>0
    other = num_class.new([other.sign, 1, exp])
  end
  tmp = Num(tmp.sign,
                    num_class.int_mult_radix_power(tmp.coefficient, tmp.exponent-other.exponent),
                    other.exponent)
  return swap ? [other, tmp] : [tmp, other]
end

- (`Object`) _number_of_digits(v)



4379
4380
4381

# File 'lib/flt/num.rb', line 4379

def _number_of_digits(v)
  _ndigits(v, num_class.radix)
end

- (`Object`) _parser(txt)

Parse numeric text literals (internal use)

# File 'lib/flt/num.rb', line 4006

def _parser(txt)
  md = /^\s*([-+])?(?:(?:(\d+)(?:\.(\d*))?|\.(\d+))(?:E([-+]?\d+))?|Inf(?:inity)?|(s)?NaN(\d*))\s*$/i.match(txt)
  if md
    OpenStruct.new :sign=>md[1], :int=>md[2], :frac=>md[3], :onlyfrac=>md[4], :exp=>md[5],
                   :signal=>md[6], :diag=>md[7]
  end
end

- (`Object`) _power(xc, xe, yc, ye, p)

Given integers xc, xe, yc and ye representing Num x = xc*radix**xe and y = yc*radix**ye, compute x**y. Returns a pair of integers (c, e) such that:

radix**(p-1) <= c <= radix**p, and
(c-1)*radix**e < x**y < (c+1)*radix**e

in other words, c*radix**e is an approximation to x**y with p digits of precision, and with an error in c of at most 1. (This is almost, but not quite, the same as the error being < 1ulp: when c

radix**(p-1) we can only guarantee error < radix ulp.)

We assume that: x is positive and not equal to 1, and y is nonzero.

# File 'lib/flt/num.rb', line 4047

def _power(xc, xe, yc, ye, p)
  # Find b such that radix**(b-1) <= |y| <= radix**b
  b = _number_of_digits(yc.abs) + ye

  # log(x) = lxc*radix**(-p-b-1), to p+b+1 places after the decimal point
  lxc = _log(xc, xe, p+b+1)

  # compute product y*log(x) = yc*lxc*radix**(-p-b-1+ye) = pc*radix**(-p-1)
  shift = ye-b
  if shift >= 0
      pc = lxc*yc*num_class.int_radix_power(shift)
  else
      pc = _div_nearest(lxc*yc, num_class.int_radix_power(-shift))
  end

  if pc == 0
      # we prefer a result that isn't exactly 1; this makes it
      # easier to compute a correctly rounded result in __pow__
      if (_number_of_digits(xc) + xe >= 1) == (yc > 0) # if x**y > 1:
          coeff, exp = num_class.int_radix_power(p-1)+1, 1-p
      else
          coeff, exp = num_class.int_radix_power(p)-1, -p
      end
  else
      coeff, exp = _exp(pc, -(p+1), p+1)
      coeff = _div_nearest(coeff, num_class.radix)
      exp += 1
  end

  return coeff, exp
end

- (`Object`) _rshift_nearest(x, shift)

Given an integer x and a nonnegative integer shift, return closest integer to x / 2**shift; use round-to-even in case of a tie.

# File 'lib/flt/num.rb', line 4279

def _rshift_nearest(x, shift)
    b, q = (1 << shift), (x >> shift)
    return q + (((2*(x & (b-1)) + (q&1)) > b) ? 1 : 0)
    #return q + (2*(x & (b-1)) + (((q&1) > b) ? 1 : 0))
end

- (`Object`) _sqrt_nearest(n, a)

Closest integer to the square root of the positive integer n. a is an initial approximation to the square root. Any positive integer will do for a, but the closer a is to the square root of n the faster convergence will be.

# File 'lib/flt/num.rb', line 4264

def _sqrt_nearest(n, a)

    if n <= 0 or a <= 0
        raise ArgumentError, "Both arguments to _sqrt_nearest should be positive."
    end

    b=0
    while a != b
        b, a = a, a--n/a>>1 # ??
    end
    return a
end

- (`Object`) log10_lb(c)

Compute a lower bound for LOG10_MULT*log10(c) for a positive integer c.

Raises:

(ArgumentError)

# File 'lib/flt/num.rb', line 4351

def log10_lb(c)
    raise ArgumentError, "The argument to _log10_lb should be nonnegative." if c <= 0
    str_c = c.to_s
    return LOG10_MULT*str_c.length - LOG10_LB_CORRECTION[str_c[0,1]]
end

- (`Object`) log2_lb(c)

Compute a lower bound for LOG2_MULT*log10(c) for a positive integer c.

Raises:

(ArgumentError)

# File 'lib/flt/num.rb', line 4339

def log2_lb(c)
    raise ArgumentError, "The argument to _log2_lb should be nonnegative." if c <= 0
    str_c = c.to_s(16)
    return LOG2_MULT*4*str_c.length - LOG2_LB_CORRECTION[str_c[0,1].to_i(16)-1]
end

Module: Flt::Num::AuxiliarFunctions

Constant Summary

Class Attribute Summary (collapse)

Instance Method Summary (collapse)

Class Attribute Details

+ (Object) log_radix_digits (readonly)

Instance Method Details

- (Object) _convert(x, error = true)

- (Object) _div_nearest(a, b)

- (Object) _exp(c, e, p)

radix**(p-1) the error could be up to radix ulp.

- (Object) _iexp(x, m, l = 8)

- (Object) _ilog(x, m, l = 8)

- (Object) _log(c, e, p)

- (Object) _log_radix_digits(p)

- (Object) _log_radix_lb(c)

- (Object) _log_radix_mult

- (Object) _normalize(op1, op2, prec = 0)

- (Object) _number_of_digits(v)

- (Object) _parser(txt)

- (Object) _power(xc, xe, yc, ye, p)

radix**(p-1) we can only guarantee error < radix ulp.)

- (Object) _rshift_nearest(x, shift)

- (Object) _sqrt_nearest(n, a)

- (Object) log10_lb(c)

- (Object) log2_lb(c)

+ (`Object`) log_radix_digits (readonly)

- (`Object`) _convert(x, error = true)

- (`Object`) _div_nearest(a, b)

- (`Object`) _exp(c, e, p)

- (`Object`) _iexp(x, m, l = 8)

- (`Object`) _ilog(x, m, l = 8)

- (`Object`) _log(c, e, p)

- (`Object`) _log_radix_digits(p)

- (`Object`) _log_radix_lb(c)

- (`Object`) _log_radix_mult

- (`Object`) _normalize(op1, op2, prec = 0)

- (`Object`) _number_of_digits(v)

- (`Object`) _parser(txt)

- (`Object`) _power(xc, xe, yc, ye, p)

- (`Object`) _rshift_nearest(x, shift)

- (`Object`) _sqrt_nearest(n, a)

- (`Object`) log10_lb(c)

- (`Object`) log2_lb(c)