Module: Flt::Num::AuxiliarFunctions
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
-
.log_radix_digits ⇒ Object
readonly
Returns the value of attribute log_radix_digits.
Class Method Summary collapse
-
._convert(x, error = true) ⇒ Object
Convert a numeric value to decimal (internal use).
-
._div_nearest(a, b) ⇒ Object
Closest integer to a/b, a and b positive integers; rounds to even in the case of a tie.
-
._exp(c, e, p) ⇒ Object
Compute an approximation to exp(c*radix**e), with p decimal places of precision.
-
._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).
-
._ilog(x, m, l = 8) ⇒ Object
Integer approximation to M*log(x/M), with absolute error boundable in terms only of x/M.
-
._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.
-
._log_radix_digits(p) ⇒ Object
Given an integer p >= 0, return floor(radix**p)*log(radix).
- ._log_radix_lb(c) ⇒ Object
- ._log_radix_mult ⇒ Object
-
._normalize(op1, op2, prec = 0) ⇒ Object
Normalizes op1, op2 to have the same exp and length of coefficient.
- ._number_of_digits(v) ⇒ Object
-
._parser(txt) ⇒ Object
Parse numeric text literals (internal use).
-
._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.
-
._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.
-
._sqrt_nearest(n, a) ⇒ Object
Closest integer to the square root of the positive integer n.
-
.log10_lb(c) ⇒ Object
Compute a lower bound for LOG10_MULT*log10© for a positive integer c.
-
.log2_lb(c) ⇒ Object
Compute a lower bound for LOG2_MULT*log10© for a positive integer c.
Class Attribute Details
.log_radix_digits ⇒ Object (readonly)
Returns the value of attribute log_radix_digits.
4311 4312 4313 |
# File 'lib/flt/num.rb', line 4311 def log_radix_digits @log_radix_digits end |
Class Method Details
._convert(x, error = true) ⇒ Object
Convert a numeric value to decimal (internal use)
4007 4008 4009 4010 4011 4012 4013 4014 4015 4016 4017 |
# File 'lib/flt/num.rb', line 4007 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.
4300 4301 4302 4303 |
# File 'lib/flt/num.rb', line 4300 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.
4104 4105 4106 4107 4108 4109 4110 4111 4112 4113 4114 4115 4116 4117 4118 4119 4120 4121 4122 4123 4124 4125 4126 4127 4128 4129 |
# File 'lib/flt/num.rb', line 4104 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).
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 4211 4212 4213 4214 4215 4216 4217 4218 4219 |
# File 'lib/flt/num.rb', line 4185 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.
4230 4231 4232 4233 4234 4235 4236 4237 4238 4239 4240 4241 4242 4243 4244 4245 4246 4247 4248 4249 4250 4251 4252 4253 4254 4255 4256 4257 4258 4259 4260 4261 4262 4263 4264 4265 4266 4267 4268 4269 4270 4271 |
# File 'lib/flt/num.rb', line 4230 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.
4135 4136 4137 4138 4139 4140 4141 4142 4143 4144 4145 4146 4147 4148 4149 4150 4151 4152 4153 4154 4155 4156 4157 4158 4159 4160 4161 4162 4163 4164 4165 4166 4167 4168 4169 4170 4171 4172 4173 4174 4175 4176 4177 4178 4179 |
# File 'lib/flt/num.rb', line 4135 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).
4317 4318 4319 4320 4321 4322 4323 4324 4325 4326 4327 4328 4329 4330 4331 4332 4333 4334 4335 4336 4337 4338 4339 4340 4341 4342 4343 4344 |
# File 'lib/flt/num.rb', line 4317 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
4381 4382 4383 4384 4385 4386 4387 4388 4389 4390 |
# File 'lib/flt/num.rb', line 4381 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
4370 4371 4372 4373 4374 4375 4376 4377 4378 4379 |
# File 'lib/flt/num.rb', line 4370 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.
4029 4030 4031 4032 4033 4034 4035 4036 4037 4038 4039 4040 4041 4042 4043 4044 4045 4046 4047 |
# File 'lib/flt/num.rb', line 4029 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
4392 4393 4394 |
# File 'lib/flt/num.rb', line 4392 def _number_of_digits(v) _ndigits(v, num_class.radix) end |
._parser(txt) ⇒ Object
Parse numeric text literals (internal use)
4020 4021 4022 4023 4024 4025 4026 |
# File 'lib/flt/num.rb', line 4020 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.
4061 4062 4063 4064 4065 4066 4067 4068 4069 4070 4071 4072 4073 4074 4075 4076 4077 4078 4079 4080 4081 4082 4083 4084 4085 4086 4087 4088 4089 4090 4091 |
# File 'lib/flt/num.rb', line 4061 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.
4292 4293 4294 4295 4296 |
# File 'lib/flt/num.rb', line 4292 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.
4277 4278 4279 4280 4281 4282 4283 4284 4285 4286 4287 4288 |
# File 'lib/flt/num.rb', line 4277 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.
4364 4365 4366 4367 4368 |
# File 'lib/flt/num.rb', line 4364 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.
4352 4353 4354 4355 4356 |
# File 'lib/flt/num.rb', line 4352 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 |