Module: Flt::Num::AuxiliarFunctions

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

Constant Summary collapse

EXP_INC = 
4
LOG_PREC_INC = 
4
LOG_RADIX_INC = 
2
LOG_RADIX_EXTRA = 
3
LOG2_MULT =

TODO: K=100? K=64? …

100
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 = 
100
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

Class Method Summary collapse

Class Attribute Details

.log_radix_digits ⇒ Object (readonly)

Returns the value of attribute log_radix_digits.



4250
4251
4252
 
# File 'lib/flt/num.rb', line 4250

def log_radix_digits
  @log_radix_digits
end

Class Method Details

._convert(x, error = true) ⇒ Object

Convert a numeric value to decimal (internal use)



3945
3946
3947
3948
3949
3950
3951
3952
3953
3954
3955
 
# File 'lib/flt/num.rb', line 3945

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

._div_nearest(a, b) ⇒ Object

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



4239
4240
4241
4242
 
# File 'lib/flt/num.rb', line 4239

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

._exp(c, e, p) ⇒ Object

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.



4042
4043
4044
4045
4046
4047
4048
4049
4050
4051
4052
4053
4054
4055
4056
4057
4058
4059
4060
4061
4062
4063
4064
4065
4066
4067
 
# File 'lib/flt/num.rb', line 4042

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

._iexp(x, m, l = 8) ⇒ Object

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).



4123
4124
4125
4126
4127
4128
4129
4130
4131
4132
4133
4134
4135
4136
4137
4138
4139
4140
4141
4142
4143
4144
4145
4146
4147
4148
4149
4150
4151
4152
4153
4154
4155
4156
4157
 
# File 'lib/flt/num.rb', line 4123

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

._ilog(x, m, l = 8) ⇒ Object

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.



4168
4169
4170
4171
4172
4173
4174
4175
4176
4177
4178
4179
4180
4181
4182
4183
4184
4185
4186
4187
4188
4189
4190
4191
4192
4193
4194
4195
4196
4197
4198
4199
4200
4201
4202
4203
4204
4205
4206
4207
4208
4209
4210
 
# File 'lib/flt/num.rb', line 4168

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

._log(c, e, p) ⇒ Object

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.



4073
4074
4075
4076
4077
4078
4079
4080
4081
4082
4083
4084
4085
4086
4087
4088
4089
4090
4091
4092
4093
4094
4095
4096
4097
4098
4099
4100
4101
4102
4103
4104
4105
4106
4107
4108
4109
4110
4111
4112
4113
4114
4115
4116
4117
 
# File 'lib/flt/num.rb', line 4073

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

._log_radix_digits(p) ⇒ Object

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

Raises:

  • (ArgumentError) 


4256
4257
4258
4259
4260
4261
4262
4263
4264
4265
4266
4267
4268
4269
4270
4271
4272
4273
4274
4275
4276
4277
4278
4279
4280
4281
4282
4283
 
# File 'lib/flt/num.rb', line 4256

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

._log_radix_lb(c) ⇒ Object 



4320
4321
4322
4323
4324
4325
4326
4327
4328
4329
 
# File 'lib/flt/num.rb', line 4320

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

._log_radix_mult ⇒ Object 



4309
4310
4311
4312
4313
4314
4315
4316
4317
4318
 
# File 'lib/flt/num.rb', line 4309

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

._normalize(op1, op2, prec = 0) ⇒ Object

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



3967
3968
3969
3970
3971
3972
3973
3974
3975
3976
3977
3978
3979
3980
3981
3982
3983
3984
3985
 
# File 'lib/flt/num.rb', line 3967

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

._number_of_digits(v) ⇒ Object 



4331
4332
4333
 
# File 'lib/flt/num.rb', line 4331

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

._parser(txt) ⇒ Object

Parse numeric text literals (internal use)



3958
3959
3960
3961
3962
3963
3964
 
# File 'lib/flt/num.rb', line 3958

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

._power(xc, xe, yc, ye, p) ⇒ Object

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.



3999
4000
4001
4002
4003
4004
4005
4006
4007
4008
4009
4010
4011
4012
4013
4014
4015
4016
4017
4018
4019
4020
4021
4022
4023
4024
4025
4026
4027
4028
4029
 
# File 'lib/flt/num.rb', line 3999

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

._rshift_nearest(x, shift) ⇒ Object

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.



4231
4232
4233
4234
4235
 
# File 'lib/flt/num.rb', line 4231

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

._sqrt_nearest(n, a) ⇒ Object

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.



4216
4217
4218
4219
4220
4221
4222
4223
4224
4225
4226
4227
 
# File 'lib/flt/num.rb', line 4216

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

.log10_lb(c) ⇒ Object

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

Raises:

  • (ArgumentError) 


4303
4304
4305
4306
4307
 
# File 'lib/flt/num.rb', line 4303

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

.log2_lb(c) ⇒ Object

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

Raises:

  • (ArgumentError) 


4291
4292
4293
4294
4295
 
# File 'lib/flt/num.rb', line 4291

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