Module: Flt::DecNum::AuxiliarFunctions
- Included in:
- Flt::DecNum
- Defined in:
- lib/flt/dec_num.rb
Overview
:nodoc:
Constant Summary collapse
- LOG10_LB_CORRECTION =
(1..9).map_hash{|i| 100 - (100*Math.log10(i)).floor}
{ # (1..9).map_hash{|i| 100 - (100*Math.log10(i)).floor} '1'=> 100, '2'=> 70, '3'=> 53, '4'=> 40, '5'=> 31, '6'=> 23, '7'=> 16, '8'=> 10, '9'=> 5}
- NUMBER_OF_DIGITS_MAX_VALID_LOG =
10**(Float::DIG-1)
Class Attribute Summary collapse
-
.log10_digits ⇒ Object
Returns the value of attribute log10_digits.
Class Method Summary collapse
-
._dexp(c, e, p) ⇒ Object
Compute an approximation to exp(c*10**e), with p decimal places of precision.
-
._div_nearest(a, b) ⇒ Object
Closest integer to a/b, a and b positive integers; rounds to even in the case of a tie.
-
._dlog(c, e, p) ⇒ Object
Given integers c, e and p with c > 0, compute an integer approximation to 10**p * log(c*10**e), with an absolute error of at most 1.
-
._dlog10(c, e, p) ⇒ Object
Given integers c, e and p with c > 0, p >= 0, compute an integer approximation to 10**p * log10(c*10**e), with an absolute error of at most 1.
-
._dpower(xc, xe, yc, ye, p) ⇒ Object
Given integers xc, xe, yc and ye representing Decimals x = xc*10**xe and y = yc*10**ye, compute x**y.
-
._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.
-
._log10_digits(p) ⇒ Object
Given an integer p >= 0, return floor(10**p)*log(10).
-
._log10_lb(c) ⇒ Object
Compute a lower bound for 100*log10© for a positive integer c.
-
._number_of_digits(i) ⇒ Object
number of bits in a nonnegative integer.
-
._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.
-
.dexp(c, e, p) ⇒ Object
Compute an approximation to exp(c*10**e), with p decimal places of precision.
Class Attribute Details
.log10_digits ⇒ Object
Returns the value of attribute log10_digits.
1122 1123 1124 |
# File 'lib/flt/dec_num.rb', line 1122 def log10_digits @log10_digits end |
Class Method Details
._dexp(c, e, p) ⇒ Object
Compute an approximation to exp(c*10**e), with p decimal places of precision. Returns integers d, f such that:
10**(p-1) <= d <= 10**p, and
(d-1)*10**f < exp(c*10**e) < (d+1)*10**f
In other words, d*10**f is an approximation to exp(c*10**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
10**(p-1) the error could be up to 10 ulp.
874 875 876 877 878 879 880 881 882 883 884 885 886 887 888 889 890 891 892 893 894 895 896 897 898 |
# File 'lib/flt/dec_num.rb', line 874 def _dexp(c, e, p) # we'll call iexp with M = 10**(p+2), giving p+3 digits of precision p += 2 # compute log(10) with extra precision = adjusted exponent of c*10**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*10**e/(log(10)) = c*10**(e+q)/(log(10)*10**q), # rounding down shift = e+q if shift >= 0 cshift = c*10**shift else cshift = c/10**-shift end quot, rem = cshift.divmod(_log10_digits(q)) # reduce remainder back to original precision rem = _div_nearest(rem, 10**extra) # error in result of _iexp < 120; error after division < 0.62 return _div_nearest(_iexp(rem, 10**p), 1000), quot - p + 3 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.
902 903 904 905 |
# File 'lib/flt/dec_num.rb', line 902 def _div_nearest(a, b) q, r = a.divmod(b) q + (((2*r + (q&1)) > b) ? 1 : 0) end |
._dlog(c, e, p) ⇒ Object
Given integers c, e and p with c > 0, compute an integer approximation to 10**p * log(c*10**e), with an absolute error of at most 1. Assumes that c*10**e is not exactly 1.
1033 1034 1035 1036 1037 1038 1039 1040 1041 1042 1043 1044 1045 1046 1047 1048 1049 1050 1051 1052 1053 1054 1055 1056 1057 1058 1059 1060 1061 1062 1063 1064 1065 1066 1067 1068 1069 1070 1071 1072 1073 1074 1075 1076 1077 |
# File 'lib/flt/dec_num.rb', line 1033 def _dlog(c, e, p) # Increase precision by 2. The precision increase is compensated # for at the end with a division by 100. p += 2 # rewrite c*10**e as d*10**f with either f >= 0 and 1 <= d <= 10, # or f <= 0 and 0.1 <= d <= 1. Then we can compute 10**p * log(c*10**e) # as 10**p * log(d) + 10**p*f * log(10). l = _number_of_digits(c) f = e+l - ((e+l >= 1) ? 1 : 0) # compute approximation to 10**p*log(d), with error < 27 if p > 0 k = e+p-f if k >= 0 c *= 10**k else c = _div_nearest(c, 10**-k) # error of <= 0.5 in c end # _ilog magnifies existing error in c by a factor of at most 10 log_d = _ilog(c, 10**p) # error < 5 + 22 = 27 else # p <= 0: just approximate the whole thing by 0; error < 2.31 log_d = 0 end # compute approximation to f*10**p*log(10), with error < 11. if f extra = _number_of_digits(f.abs) - 1 if p + extra >= 0 # error in f * _log10_digits(p+extra) < |f| * 1 = |f| # after division, error < |f|/10**extra + 0.5 < 10 + 0.5 < 11 f_log_ten = _div_nearest(f*_log10_digits(p+extra), 10**extra) else f_log_ten = 0 end else f_log_ten = 0 end # error in sum < 11+27 = 38; error after division < 0.38 + 0.5 < 1 return _div_nearest(f_log_ten + log_d, 100) end |
._dlog10(c, e, p) ⇒ Object
Given integers c, e and p with c > 0, p >= 0, compute an integer approximation to 10**p * log10(c*10**e), with an absolute error of at most 1. Assumes that c*10**e is not exactly 1.
988 989 990 991 992 993 994 995 996 997 998 999 1000 1001 1002 1003 1004 1005 1006 1007 1008 1009 1010 1011 1012 1013 1014 1015 1016 1017 1018 |
# File 'lib/flt/dec_num.rb', line 988 def _dlog10(c, e, p) # increase precision by 2; compensate for this by dividing # final result by 100 p += 2 # write c*10**e as d*10**f with either: # f >= 0 and 1 <= d <= 10, or # f <= 0 and 0.1 <= d <= 1. # Thus for c*10**e close to 1, f = 0 l = _number_of_digits(c) f = e+l - ((e+l >= 1) ? 1 : 0) if p > 0 m = 10**p k = e+p-f if k >= 0 c *= 10**k else c = _div_nearest(c, 10**-k) end log_d = _ilog(c, m) # error < 5 + 22 = 27 log_10 = _log10_digits(p) # error < 1 log_d = _div_nearest(log_d*m, log_10) log_tenpower = f*m # exact else log_d = 0 # error < 2.31 log_tenpower = _div_nearest(f, 10**-p) # error < 0.5 end return _div_nearest(log_tenpower+log_d, 100) end |
._dpower(xc, xe, yc, ye, p) ⇒ Object
Given integers xc, xe, yc and ye representing Decimals x = xc*10**xe and y = yc*10**ye, compute x**y. Returns a pair of integers (c, e) such that:
10**(p-1) <= c <= 10**p, and
(c-1)*10**e < x**y < (c+1)*10**e
in other words, c*10**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
10**(p-1) we can only guarantee error < 10ulp.)
We assume that: x is positive and not equal to 1, and y is nonzero.
832 833 834 835 836 837 838 839 840 841 842 843 844 845 846 847 848 849 850 851 852 853 854 855 856 857 858 859 860 861 862 |
# File 'lib/flt/dec_num.rb', line 832 def _dpower(xc, xe, yc, ye, p) # Find b such that 10**(b-1) <= |y| <= 10**b b = _number_of_digits(yc.abs) + ye # log(x) = lxc*10**(-p-b-1), to p+b+1 places after the decimal point lxc = _dlog(xc, xe, p+b+1) # compute product y*log(x) = yc*lxc*10**(-p-b-1+ye) = pc*10**(-p-1) shift = ye-b if shift >= 0 pc = lxc*yc*10**shift else pc = _div_nearest(lxc*yc, 10**-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 = 10**(p-1)+1, 1-p else coeff, exp = 10**p-1, -p end else coeff, exp = _dexp(pc, -(p+1), p+1) coeff = _div_nearest(coeff, 10) exp += 1 end return coeff, exp 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 0 <= x/M <= 2.4, the absolute error in the result is bounded by 60 (and is usually much smaller).
1083 1084 1085 1086 1087 1088 1089 1090 1091 1092 1093 1094 1095 1096 1097 1098 1099 1100 1101 1102 1103 1104 1105 1106 1107 1108 1109 1110 1111 1112 1113 1114 1115 1116 1117 |
# File 'lib/flt/dec_num.rb', line 1083 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 = -(-10*_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 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.
941 942 943 944 945 946 947 948 949 950 951 952 953 954 955 956 957 958 959 960 961 962 963 964 965 966 967 968 969 970 971 972 973 974 975 976 977 978 979 980 981 982 983 |
# File 'lib/flt/dec_num.rb', line 941 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 |
._log10_digits(p) ⇒ Object
Given an integer p >= 0, return floor(10**p)*log(10).
1126 1127 1128 1129 1130 1131 1132 1133 1134 1135 1136 1137 1138 1139 1140 1141 1142 1143 1144 1145 1146 1147 1148 1149 |
# File 'lib/flt/dec_num.rb', line 1126 def _log10_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 be correct # (truncated, not rounded to nearest). raise ArgumentError, "p should be nonnegative" if p<0 if p >= AuxiliarFunctions.log10_digits.length digits = nil # compute p+3, p+6, p+9, ... digits; continue until at # least one of the extra digits is nonzero extra = 3 loop do # compute p+extra digits, correct to within 1ulp m = 10**(p+extra+2) digits = _div_nearest(_ilog(10*m, m), 100).to_s break if digits[-extra..-1] != '0'*extra extra += 3 end # keep all reliable digits so far; remove trailing zeros # and next nonzero digit AuxiliarFunctions.log10_digits = digits.sub(/0*$/,'')[0...-1] end return (AuxiliarFunctions.log10_digits[0...p+1]).to_i end |
._log10_lb(c) ⇒ Object
Compute a lower bound for 100*log10© for a positive integer c.
1021 1022 1023 1024 1025 |
# File 'lib/flt/dec_num.rb', line 1021 def _log10_lb(c) raise ArgumentError, "The argument to _log10_lb should be nonnegative." if c <= 0 str_c = c.to_s return 100*str_c.length - LOG10_LB_CORRECTION[str_c[0,1]] end |
._number_of_digits(i) ⇒ Object
number of bits in a nonnegative integer
1189 1190 1191 1192 1193 1194 1195 1196 1197 1198 |
# File 'lib/flt/dec_num.rb', line 1189 def _number_of_digits(i) raise TypeError, "The argument to _number_of_digits should be nonnegative." if i < 0 if i.is_a?(Fixnum) || (i > NUMBER_OF_DIGITS_MAX_VALID_LOG) # for short integers this is faster # note that here we return 1 for 0 i.to_s.length else (::Math.log10(i)+1).floor end 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.
926 927 928 929 930 |
# File 'lib/flt/dec_num.rb', line 926 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.
911 912 913 914 915 916 917 918 919 920 921 922 |
# File 'lib/flt/dec_num.rb', line 911 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 |
.dexp(c, e, p) ⇒ Object
Compute an approximation to exp(c*10**e), with p decimal places of precision.
Returns integers d, f such that:
10**(p-1) <= d <= 10**p, and
(d-1)*10**f < exp(c*10**e) < (d+1)*10**f
In other words, d*10**f is an approximation to exp(c*10**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
10**(p-1) the error could be up to 10 ulp.
1163 1164 1165 1166 1167 1168 1169 1170 1171 1172 1173 1174 1175 1176 1177 1178 1179 1180 1181 1182 1183 1184 1185 1186 |
# File 'lib/flt/dec_num.rb', line 1163 def dexp(c, e, p) # we'll call iexp with M = 10**(p+2), giving p+3 digits of precision p += 2 # compute log(10) with extra precision = adjusted exponent of c*10**e extra = [0, e + _number_of_digits(c) - 1].max q = p + extra # compute quotient c*10**e/(log(10)) = c*10**(e+q)/(log(10)*10**q), # rounding down shift = e+q if shift >= 0 cshift = c*10**shift else cshift = c/10**-shift end quot, rem = cshift.divmod(_log10_digits(q)) # reduce remainder back to original precision rem = _div_nearest(rem, 10**extra) # error in result of _iexp < 1s20; error after division < 0.62 return _div_nearest(_iexp(rem, 10**p), 1000), quot - p + 3 end |