In Files

Class Index [+]

Quicksearch

Flt::Num::AuxiliarFunctions

Constants

EXP_INC
LOG_PREC_INC
LOG_RADIX_INC
LOG_RADIX_EXTRA
LOG2_MULT
LOG2_LB_CORRECTION
LOG10_MULT
LOG10_LB_CORRECTION

Attributes

log_radix_digits[R]

Public Instance Methods

_convert(x, error=true) click to toggle source

Convert a numeric value to decimal (internal use)

      # File lib/flt/num.rb, line 3945
3945:     def _convert(x, error=true)
3946:       case x
3947:       when num_class
3948:         x
3949:       when *num_class.context.coercible_types
3950:         num_class.new(x)
3951:       else
3952:         raise TypeError, "Unable to convert #{x.class} to #{num_class}" if error
3953:         nil
3954:       end
3955:     end
_div_nearest(a, b) click to toggle source

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 4239
4239:     def _div_nearest(a, b)
4240:       q, r = a.divmod(b)
4241:       q + (((2*r + (q&1)) > b) ? 1 : 0)
4242:     end
_exp(c, e, p) click to toggle source

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 4042
4042:     def _exp(c, e, p)
4043:         # we'll call iexp with M = radix**(p+2), giving p+3 digits of precision
4044:         p += EXP_INC
4045: 
4046:         # compute log(radix) with extra precision = adjusted exponent of c*radix**e
4047:         # TODO: without the .abs tests fail because c is negative: c should not be negative!!
4048:         extra = [0, e + _number_of_digits(c.abs) - 1].max
4049:         q = p + extra
4050: 
4051:         # compute quotient c*radix**e/(log(radix)) = c*radix**(e+q)/(log(radix)*radix**q),
4052:         # rounding down
4053:         shift = e+q
4054:         if shift >= 0
4055:             cshift = c*num_class.int_radix_power(shift)
4056:         else
4057:             cshift = c/num_class.int_radix_power(-shift)
4058:         end
4059:         quot, rem = cshift.divmod(_log_radix_digits(q))
4060: 
4061:         # reduce remainder back to original precision
4062:         rem = _div_nearest(rem, num_class.int_radix_power(extra))
4063: 
4064:         # for radix=10: error in result of _iexp < 120;  error after division < 0.62
4065:         r = _div_nearest(_iexp(rem, num_class.int_radix_power(p)), num_class.int_radix_power(EXP_INC+1)), quot - p + (EXP_INC+1)
4066:         return r
4067:     end
_iexp(x, m, l=8) click to toggle source

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 4123
4123:     def _iexp(x, m, l=8)
4124: 
4125:         # Algorithm: to compute exp(z) for a real number z, first divide z
4126:         # by a suitable power R of 2 so that |z/2**R| < 2**-L.  Then
4127:         # compute expm1(z/2**R) = exp(z/2**R) - 1 using the usual Taylor
4128:         # series
4129:         #
4130:         #     expm1(x) = x + x**2/2! + x**3/3! + ...
4131:         #
4132:         # Now use the identity
4133:         #
4134:         #     expm1(2x) = expm1(x)*(expm1(x)+2)
4135:         #
4136:         # R times to compute the sequence expm1(z/2**R),
4137:         # expm1(z/2**(R-1)), ... , exp(z/2), exp(z).
4138: 
4139:         # Find R such that x/2**R/M <= 2**-L
4140:         r = _nbits((x<<l)/m)
4141: 
4142:         # Taylor series.  (2**L)**T > M
4143:         t = -(-num_class.radix*_number_of_digits(m)/(3*l)).to_i
4144:         y = _div_nearest(x, t)
4145:         mshift = m<<r
4146:         (1...t).to_a.reverse.each do |i|
4147:             y = _div_nearest(x*(mshift + y), mshift * i)
4148:         end
4149: 
4150:         # Expansion
4151:         (0...r).to_a.reverse.each do |k|
4152:             mshift = m<<(k+2)
4153:             y = _div_nearest(y*(y+mshift), mshift)
4154:         end
4155: 
4156:         return m+y
4157:     end
_ilog(x, m, l = 8) click to toggle source

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 4168
4168:     def _ilog(x, m, l = 8)
4169:       # The basic algorithm is the following: let log1p be the function
4170:       # log1p(x) = log(1+x).  Then log(x/M) = log1p((x-M)/M).  We use
4171:       # the reduction
4172:       #
4173:       #    log1p(y) = 2*log1p(y/(1+sqrt(1+y)))
4174:       #
4175:       # repeatedly until the argument to log1p is small (< 2**-L in
4176:       # absolute value).  For small y we can use the Taylor series
4177:       # expansion
4178:       #
4179:       #    log1p(y) ~ y - y**2/2 + y**3/3 - ... - (-y)**T/T
4180:       #
4181:       # truncating at T such that y**T is small enough.  The whole
4182:       # computation is carried out in a form of fixed-point arithmetic,
4183:       # with a real number z being represented by an integer
4184:       # approximation to z*M.  To avoid loss of precision, the y below
4185:       # is actually an integer approximation to 2**R*y*M, where R is the
4186:       # number of reductions performed so far.
4187: 
4188:       y = x-m
4189:       # argument reduction; R = number of reductions performed
4190:       r = 0
4191:       # while (r <= l && y.abs << l-r >= m ||
4192:       #        r > l and y.abs>> r-l >= m)
4193:       while (((r <= l) && ((y.abs << (l-r)) >= m)) ||
4194:              ((r > l) && ((y.abs>>(r-l)) >= m)))
4195:           y = _div_nearest((m*y) << 1,
4196:                            m + _sqrt_nearest(m*(m+_rshift_nearest(y, r)), m))
4197:           r += 1
4198:       end
4199: 
4200:       # Taylor series with T terms
4201:       t = -(10*_number_of_digits(m)/(3*l)).to_i
4202:       yshift = _rshift_nearest(y, r)
4203:       w = _div_nearest(m, t)
4204:       # (1...t).reverse_each do |k| # Ruby 1.9
4205:       (1...t).to_a.reverse.each do |k|
4206:          w = _div_nearest(m, k) - _div_nearest(yshift*w, m)
4207:       end
4208: 
4209:       return _div_nearest(w*y, m)
4210:     end
_log(c, e, p) click to toggle source

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 4073
4073:     def _log(c, e, p)
4074:         # Increase precision by 2. The precision increase is compensated
4075:         # for at the end with a division
4076:         p += LOG_PREC_INC
4077: 
4078:         # rewrite c*radix**e as d*radix**f with either f >= 0 and 1 <= d <= radix,
4079:         # or f <= 0 and 1/radix <= d <= 1.  Then we can compute radix**p * log(c*radix**e)
4080:         # as radix**p * log(d) + radix**p*f * log(radix).
4081:         l = _number_of_digits(c)
4082:         f = e+l - ((e+l >= 1) ? 1 : 0)
4083: 
4084:         # compute approximation to radix**p*log(d), with error < 27 for radix=10
4085:         if p > 0
4086:             k = e+p-f
4087:             if k >= 0
4088:                 c *= num_class.int_radix_power(k)
4089:             else
4090:                 c = _div_nearest(c, num_class.int_radix_power(-k))  # error of <= 0.5 in c for radix=10
4091:             end
4092: 
4093:             # _ilog magnifies existing error in c by a factor of at most radix
4094:             log_d = _ilog(c, num_class.int_radix_power(p)) # error < 5 + 22 = 27 for radix=10
4095:         else
4096:             # p <= 0: just approximate the whole thing by 0; error < 2.31 for radix=10
4097:             log_d = 0
4098:         end
4099: 
4100:         # compute approximation to f*radix**p*log(radix), with error < 11 for radix=10.
4101:         if f
4102:             extra = _number_of_digits(f.abs) - 1
4103:             if p + extra >= 0
4104:                 # for radix=10:
4105:                 # error in f * _log10_digits(p+extra) < |f| * 1 = |f|
4106:                 # after division, error < |f|/10**extra + 0.5 < 10 + 0.5 < 11
4107:                 f_log_r = _div_nearest(f*_log_radix_digits(p+extra), num_class.int_radix_power(extra))
4108:             else
4109:                 f_log_r = 0
4110:             end
4111:         else
4112:             f_log_r = 0
4113:         end
4114: 
4115:         # error in sum < 11+27 = 38; error after division < 0.38 + 0.5 < 1 for radix=10
4116:         return _div_nearest(f_log_r + log_d, num_class.int_radix_power(LOG_PREC_INC)) # extra radix factor for base 2 ???
4117:     end
_log_radix_digits(p) click to toggle source

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

      # File lib/flt/num.rb, line 4256
4256:     def _log_radix_digits(p)
4257:       # digits are stored as a string, for quick conversion to
4258:       # integer in the case that we've already computed enough
4259:       # digits; the stored digits should always bge correct
4260:       # (truncated, not rounded to nearest).
4261:       raise ArgumentError, "p should be nonnegative" if p<0
4262:       stored_digits = (AuxiliarFunctions.log_radix_digits[num_class.radix] || "")
4263:       if p >= stored_digits.length
4264:           digits = nil
4265:           # compute p+3, p+6, p+9, ... digits; continue until at
4266:           # least one of the extra digits is nonzero
4267:           extra = LOG_RADIX_EXTRA
4268:           loop do
4269:             # compute p+extra digits, correct to within 1ulp
4270:             m = num_class.int_radix_power(p+extra+LOG_RADIX_INC)
4271:             digits = _div_nearest(_ilog(num_class.radix*m, m), num_class.int_radix_power(LOG_RADIX_INC)).to_s(num_class.radix)
4272:             break if digits[-extra..1] != '0'*extra
4273:             extra += LOG_RADIX_EXTRA
4274:           end
4275:           # if the radix < e (i.e. only for radix==2), we must prefix with a 0 because log(radix)<1
4276:           # BUT THIS REDUCES PRECISION BY ONE? : may be avoid prefix and adjust scaling in the caller
4277:           prefix = num_class.radix==2 ? '0' : ''
4278:           # keep all reliable digits so far; remove trailing zeros
4279:           # and next nonzero digit
4280:           AuxiliarFunctions.log_radix_digits[num_class.radix] = prefix + digits.sub(/0*$/,'')[0...1]
4281:       end
4282:       return (AuxiliarFunctions.log_radix_digits[num_class.radix][0..p]).to_i(num_class.radix)
4283:     end
_log_radix_lb(c) click to toggle source
      # File lib/flt/num.rb, line 4320
4320:     def _log_radix_lb(c)
4321:       case num_class.radix
4322:       when 10
4323:         log10_lb(c)
4324:       when 2
4325:         log2_lb(c)
4326:       else
4327:         raise ArgumentError, "_log_radix_lb not implemented for base #{num_class.radix}"
4328:       end
4329:     end
_log_radix_mult() click to toggle source
      # File lib/flt/num.rb, line 4309
4309:     def _log_radix_mult
4310:       case num_class.radix
4311:       when 10
4312:         LOG10_MULT
4313:       when 2
4314:         LOG2_MULT
4315:       else
4316:         raise ArgumentError, "_log_radix_mult not implemented for base #{num_class.radix}"
4317:       end
4318:     end
_normalize(op1, op2, prec=0) click to toggle source

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

      # File lib/flt/num.rb, line 3967
3967:     def _normalize(op1, op2, prec=0)
3968:       if op1.exponent < op2.exponent
3969:         swap = true
3970:         tmp,other = op2,op1
3971:       else
3972:         swap = false
3973:         tmp,other = op1,op2
3974:       end
3975:       tmp_len = tmp.number_of_digits
3976:       other_len = other.number_of_digits
3977:       exp = tmp.exponent + [1, tmp_len - prec - 2].min
3978:       if (other_len+other.exponent-1 < exp) && prec>0
3979:         other = num_class.new([other.sign, 1, exp])
3980:       end
3981:       tmp = Num(tmp.sign,
3982:                         num_class.int_mult_radix_power(tmp.coefficient, tmp.exponent-other.exponent),
3983:                         other.exponent)
3984:       return swap ? [other, tmp] : [tmp, other]
3985:     end
_number_of_digits(v) click to toggle source
      # File lib/flt/num.rb, line 4331
4331:     def _number_of_digits(v)
4332:       _ndigits(v, num_class.radix)
4333:     end
_parser(txt) click to toggle source

Parse numeric text literals (internal use)

      # File lib/flt/num.rb, line 3958
3958:     def _parser(txt)
3959:       md = /^\s*([-+])?(?:(?:(\d+)(?:\.(\d*))?|\.(\d+))(?:E([-+]?\d+))?|Inf(?:inity)?|(s)?NaN(\d*))\s*$/.match(txt)
3960:       if md
3961:         OpenStruct.new :sign=>md[1], :int=>md[2], :frac=>md[3], :onlyfrac=>md[4], :exp=>md[5],
3962:                        :signal=>md[6], :diag=>md[7]
3963:       end
3964:     end
_power(xc, xe, yc, ye, p) click to toggle source

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 3999
3999:     def _power(xc, xe, yc, ye, p)
4000:       # Find b such that radix**(b-1) <= |y| <= radix**b
4001:       b = _number_of_digits(yc.abs) + ye
4002: 
4003:       # log(x) = lxc*radix**(-p-b-1), to p+b+1 places after the decimal point
4004:       lxc = _log(xc, xe, p+b+1)
4005: 
4006:       # compute product y*log(x) = yc*lxc*radix**(-p-b-1+ye) = pc*radix**(-p-1)
4007:       shift = ye-b
4008:       if shift >= 0
4009:           pc = lxc*yc*num_class.int_radix_power(shift)
4010:       else
4011:           pc = _div_nearest(lxc*yc, num_class.int_radix_power(-shift))
4012:       end
4013: 
4014:       if pc == 0
4015:           # we prefer a result that isn't exactly 1; this makes it
4016:           # easier to compute a correctly rounded result in __pow__
4017:           if (_number_of_digits(xc) + xe >= 1) == (yc > 0) # if x**y > 1:
4018:               coeff, exp = num_class.int_radix_power(p-1)+1, 1-p
4019:           else
4020:               coeff, exp = num_class.int_radix_power(p)-1, -p
4021:           end
4022:       else
4023:           coeff, exp = _exp(pc, -(p+1), p+1)
4024:           coeff = _div_nearest(coeff, num_class.radix)
4025:           exp += 1
4026:       end
4027: 
4028:       return coeff, exp
4029:     end
_rshift_nearest(x, shift) click to toggle source

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 4231
4231:     def _rshift_nearest(x, shift)
4232:         b, q = (1 << shift), (x >> shift)
4233:         return q + (((2*(x & (b-1)) + (q&1)) > b) ? 1 : 0)
4234:         #return q + (2*(x & (b-1)) + (((q&1) > b) ? 1 : 0))
4235:     end
_sqrt_nearest(n, a) click to toggle source

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 4216
4216:     def _sqrt_nearest(n, a)
4217: 
4218:         if n <= 0 or a <= 0
4219:             raise ArgumentError, "Both arguments to _sqrt_nearest should be positive."
4220:         end
4221: 
4222:         b=0
4223:         while a != b
4224:             b, a = a, a--n/a>>1 # ??
4225:         end
4226:         return a
4227:     end
log10_lb(c) click to toggle source

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

      # File lib/flt/num.rb, line 4303
4303:     def log10_lb(c)
4304:         raise ArgumentError, "The argument to _log10_lb should be nonnegative." if c <= 0
4305:         str_c = c.to_s
4306:         return LOG10_MULT*str_c.length - LOG10_LB_CORRECTION[str_c[0,1]]
4307:     end
log2_lb(c) click to toggle source

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

      # File lib/flt/num.rb, line 4291
4291:     def log2_lb(c)
4292:         raise ArgumentError, "The argument to _log2_lb should be nonnegative." if c <= 0
4293:         str_c = c.to_s(16)
4294:         return LOG2_MULT*4*str_c.length - LOG2_LB_CORRECTION[str_c[0,1].to_i(16)-1]
4295:     end

Disabled; run with --debug to generate this.

[Validate]

Generated with the Darkfish Rdoc Generator 1.1.6.