Class: Flt::Support::Formatter
Overview
Burger and Dybvig free formatting algorithm, from their paper: “Printing Floating-Point Numbers Quickly and Accurately” (Robert G. Burger, R. Kent Dybvig)
This algorithm formats arbitrary base floating point numbers as decimal text literals. The floating-point (with fixed precision) is interpreted as an approximated value, representing any value in its ‘rounding-range’ (the interval where all values round to the floating-point value, with the given precision and rounding mode). An alternative approach which is not taken here would be to represent the exact floating-point value with some given precision and rounding mode requirements; that can be achieved with Clinger algorithm (which may fail for exact precision).
The variables used by the algorithm are stored in instance variables: Quotients of integers will be used to hold the magnitudes: All numbers in the randound-range are rounded to @v (with the given precision p) significant digit right after the radix point. @b**@k is the first power of @b >= u
The rounding range of @v is the interval of values that round to @v under the runding-mode. If the rounding mode is one of the round-to-nearest variants (even, up, down), then it is ((v+v-)/2 = (@v-@m_m)/@s, (v+v+)/2 = (@v+@m_)/2) whith the boundaries open or closed as explained below. In this case:
@m_m/@s = (@v - (v + v-)/2) where v- = @v.next_minus is the lower adjacent to v floating point value
@m_p/@s = ((v + v+)/2 - @v) where v+ = @v.next_plus is the upper adjacent to v floating point value
If the rounding is directed, then the rounding interval is either (v-, @v] or [@v, v+] if @roundh, then @k is the minimum @k with (@r+@m_p)/@s <= @output_b**@k
@k = ceil(logB((@r+@m_p)/2)) with lobB the @output_b base logarithm
if @roundh, then @k is the minimum @k with (@r+@m_p)/@s < @output_b**@k
@k = 1+floor(logB((@r+@m_p)/2))
p is the input floating point precision
Constant Summary collapse
- ESTIMATE_FLOAT_LOG_B =
{2=>1/Math.log(2), 10=>1/Math.log(10), 16=>1/Math.log(16)}
Instance Attribute Summary collapse
-
#round_up ⇒ Object
readonly
Returns the value of attribute round_up.
Instance Method Summary collapse
-
#adjusted_digits(round_mode) ⇒ Object
Access rounded result of format operation: scaling (position of radix point) and digits.
- #b_power(n) ⇒ Object
-
#digits ⇒ Object
Access result of format operation: scaling (position of radix point) and digits.
-
#format(v, f, e, round_mode, p = nil, all = false) ⇒ Object
This method converts v = f*b**e into a sequence of output_b-base digits, so that if the digits are converted back to a floating-point value of precision p (correctly rounded), the result is v.
- #generate ⇒ Object
- #generate_max ⇒ Object
-
#initialize(input_b, input_min_e, output_b) ⇒ Formatter
constructor
This Object-oriented implementation is slower than the functional one for two reasons: * The overhead of object creation * The use of instance variables instead of local variables But if scale is optimized or local variables are used in the inner loops, then this implementation is on par with the functional one for Float and it is more efficient for Flt types, where the variables passed as parameters hold larger objects.
- #output_b_power(n) ⇒ Object
-
#scale! ⇒ Object
using local vars instead of instance vars: it makes a difference in performance.
-
#scale_fixup! ⇒ Object
fix up scaling (final step): specialized version of scale! This performs a single up scaling step, i.e.
-
#scale_optimized! ⇒ Object
scale_o1! is an optimized version of scale!; it requires an additional parameters with the floating-point number v=r/s.
-
#scale_original!(really = false) ⇒ Object
Given r/s = v (number to convert to text), m_m/s = (v - v-)/s, m_p/s = (v+ - v)/s Scale the fractions so that the first significant digit is right after the radix point, i.e.
Constructor Details
#initialize(input_b, input_min_e, output_b) ⇒ Formatter
This Object-oriented implementation is slower than the functional one for two reasons:
-
The overhead of object creation
-
The use of instance variables instead of local variables
But if scale is optimized or local variables are used in the inner loops, then this implementation is on par with the functional one for Float and it is more efficient for Flt types, where the variables passed as parameters hold larger objects.
832 833 834 835 836 837 838 839 840 841 842 843 844 845 846 847 |
# File 'lib/flt/support.rb', line 832 def initialize(input_b, input_min_e, output_b) @b = input_b @min_e = input_min_e @output_b = output_b # result of last operation @adjusted_digits = @digits = nil # for "all-digits" mode results (which are truncated, rather than rounded), # round_up contains information to round the result: # * it is nil if the rest of digits are zero (the result is exact) # * it is :lo if there exist non-zero digits beyond the significant ones (those returned), but # the value is below the tie (the value must be rounded up only for :up rounding mode) # * it is :tie if there exists exactly one nonzero digit after the significant and it is radix/2, # for round-to-nearest it is atie. # * it is :hi otherwise (the value should be rounded-up except for the :down mode) @round_up = nil end |
Instance Attribute Details
#round_up ⇒ Object (readonly)
Returns the value of attribute round_up.
979 980 981 |
# File 'lib/flt/support.rb', line 979 def round_up @round_up end |
Instance Method Details
#adjusted_digits(round_mode) ⇒ Object
Access rounded result of format operation: scaling (position of radix point) and digits
983 984 985 986 987 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 1019 1020 1021 1022 1023 1024 |
# File 'lib/flt/support.rb', line 983 def adjusted_digits(round_mode) round_mode = Support.simplified_round_mode(round_mode, @minus) if @adjusted_digits.nil? && !@digits.nil? increment = (@round_up && (round_mode != :down)) && ((round_mode == :up) || (@round_up == :hi) || ((@round_up == :tie) && ((round_mode==:half_up) || ((round_mode==:half_even) && ((@digits.last % 2)==1))))) # increment = (@round_up == :tie) || (@round_up == :hi) # old behaviour (:half_up) if increment base = @output_b dec_pos = @k digits = @digits.dup # carry = increment ? 1 : 0 # digits = digits.reverse.map{|d| d += carry; d>=base ? 0 : (carry=0;d)}.reverse # if carry != 0 # digits.unshift carry # dec_pos += 1 # end i = digits.size - 1 while i>=0 digits[i] += 1 if digits[i] == base digits[i] = 0 else break end i -= 1 end if i<0 dec_pos += 1 digits.unshift 1 end @adjusted_k = dec_pos @adjusted_digits = digits else @adjusted_k = @k @adjusted_digits = @digits end end return @adjusted_k, @adjusted_digits end |
#b_power(n) ⇒ Object
1073 1074 1075 |
# File 'lib/flt/support.rb', line 1073 def b_power(n) @b**n end |
#digits ⇒ Object
Access result of format operation: scaling (position of radix point) and digits
975 976 977 |
# File 'lib/flt/support.rb', line 975 def digits return @k, @digits end |
#format(v, f, e, round_mode, p = nil, all = false) ⇒ Object
This method converts v = f*b**e into a sequence of output_b-base digits, so that if the digits are converted back to a floating-point value of precision p (correctly rounded), the result is v. If round_mode is not nil, just enough digits to produce v using that rounding is used; otherwise enough digits to produce v with any rounding are delivered.
If the all
parameter is true, all significant digits are generated without rounding, i.e. all digits that, if used on input, cannot arbitrarily change while preserving the parsed value of the floating point number. Since the digits are not rounded more digits may be needed to assure round-trip value preservation. This is useful to reflect the precision of the floating point value in the output; in particular trailing significant zeros are shown. But note that, for directed rounding and base conversion this may need to produce an infinite number of digits, in which case an exception will be raised. This is specially frequent for the :up rounding mode, in which any number with a finite number of nonzero digits equal to or less than the precision will haver and infinite sequence of zero significant digits.
With :down rounding (truncation) this could be used to show the exact value of the floating point but beware: when used with directed rounding, if the value has not an exact representation in the output base this will lead to an infinite loop. formatting ‘0.1’ (as a decimal floating-point number) in base 2 with :down rounding
When the all
parameters is used the result is not rounded (is truncated), and the round_up flag is set to indicate that nonzero digits exists beyond the returned digits; the possible values of the round_up flag are:
-
nil : the rest of digits are zero (the result is exact)
-
:lo : there exist non-zero digits beyond the significant ones (those returned), but the value is below the tie (the value must be rounded up only for :up rounding mode)
-
:tie : there exists exactly one nonzero digit after the significant and it is radix/2, for round-to-nearest it is atie.
-
:hi : the value is closer to the rounded-up value (incrementing the last significative digit.)
Note that the round_mode here is not the rounding mode applied to the output; it is the rounding mode that applied to input preserves the original floating-point value (with the same precision as input). should be rounded-up.
886 887 888 889 890 891 892 893 894 895 896 897 898 899 900 901 902 903 904 905 906 907 908 909 910 911 912 913 914 915 916 917 918 919 920 921 922 923 924 925 926 927 928 929 930 931 932 933 934 935 936 937 938 939 940 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 |
# File 'lib/flt/support.rb', line 886 def format(v, f, e, round_mode, p=nil, all=false) context = v.class.context # TODO: consider removing parameters f,e and using v.split instead @minus = (context.sign(v)==-1) @v = context.copy_sign(v, +1) # don't use context.abs(v) because it rounds (and may overflow also) @f = f.abs @e = e @round_mode = round_mode @all_digits = all p ||= context.precision # adjust the rounding mode to work only with positive numbers @round_mode = Support.simplified_round_mode(@round_mode, @minus) # determine the high,low inclusion flags of the rounding limits case @round_mode when :half_even # rounding rage is (v-m-,v+m+) if v is odd and [v+m-,v+m+] if even @round_l = @round_h = ((@f%2)==0) when :up # rounding rage is (v-,v] # ceiling is treated here assuming f>0 @round_l, @round_h = false, true when :down # rounding rage is [v,v+) # floor is treated here assuming f>0 @round_l, @round_h = true, false when :half_up # rounding rage is [v+m-,v+m+) @round_l, @round_h = true, false when :half_down # rounding rage is (v+m-,v+m+] @round_l, @round_h = false, true else # :nearest # Here assume only that round-to-nearest will be used, but not which variant of it # The result is valid for any rounding (to nearest) but may produce more digits # than stricly necessary for specific rounding modes. # That is, enough digits are generated so that when the result is # converted to floating point with the specified precision and # correct rounding (to nearest), the result is the original number. # rounding range is (v+m-,v+m+) @round_l = @round_h = false end # TODO: use context.next_minus, next_plus instead of direct computing, don't require min_e & ps # Now compute the working quotients @r/@s, @m_p/@s = (v+ - @v), @m_m/@s = (@v - v-) and scale them. if @e >= 0 if @f != b_power(p-1) be = b_power(@e) @r, @s, @m_p, @m_m = @f*be*2, 2, be, be else be = b_power(@e) be1 = be*@b @r, @s, @m_p, @m_m = @f*be1*2, @b*2, be1, be end else if @e==@min_e or @f != b_power(p-1) @r, @s, @m_p, @m_m = @f*2, b_power(-@e)*2, 1, 1 else @r, @s, @m_p, @m_m = @f*@b*2, b_power(1-@e)*2, @b, 1 end end @k = 0 @context = context scale_optimized! # The value to be formatted is @v=@r/@s; m- = @m_m/@s = (@v - v-)/@s; m+ = @m_p/@s = (v+ - @v)/@s # Now adjust @m_m, @m_p so that they define the rounding range case @round_mode when :up # ceiling is treated here assuming @f>0 # rounding range is -v,@v @m_m, @m_p = @m_m*2, 0 when :down # floor is treated here assuming #f>0 # rounding range is @v,v+ @m_m, @m_p = 0, @m_p*2 else # rounding range is v-,v+ # @m_m, @m_p = @m_m, @m_p end # Now m_m, m_p define the rounding range all ? generate_max : generate end |
#generate ⇒ Object
1121 1122 1123 1124 1125 1126 1127 1128 1129 1130 1131 1132 1133 1134 1135 1136 1137 1138 1139 1140 1141 1142 1143 1144 1145 1146 1147 1148 1149 1150 1151 1152 1153 1154 1155 |
# File 'lib/flt/support.rb', line 1121 def generate list = [] r, s, m_p, m_m, = @r, @s, @m_p, @m_m loop do d,r = (r*@output_b).divmod(s) m_p *= @output_b m_m *= @output_b tc1 = @round_l ? (r<=m_m) : (r<m_m) tc2 = @round_h ? (r+m_p >= s) : (r+m_p > s) if not tc1 if not tc2 list << d else list << d+1 break end else if not tc2 list << d break else if r*2 < s list << d break else list << d+1 break end end end end @digits = list end |
#generate_max ⇒ Object
1081 1082 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 1118 1119 |
# File 'lib/flt/support.rb', line 1081 def generate_max @round_up = false list = [] r, s, m_p, m_m, = @r, @s, @m_p, @m_m n_iters, rs = 0, [] loop do if (n_iters > 10000) raise "Infinite digit sequence." if rs.include?(r) rs << r else n_iters += 1 end d,r = (r*@output_b).divmod(s) m_p *= @output_b m_m *= @output_b list << d tc1 = @round_l ? (r<=m_m) : (r<m_m) tc2 = @round_h ? (r+m_p >= s) : (r+m_p > s) if tc1 && tc2 if r != 0 r *= 2 if r > s @round_up = :hi elsif r == s @round_up = :tie else @rund_up = :lo end end break end end @digits = list end |
#output_b_power(n) ⇒ Object
1077 1078 1079 |
# File 'lib/flt/support.rb', line 1077 def output_b_power(n) @output_b**n end |
#scale! ⇒ Object
using local vars instead of instance vars: it makes a difference in performance
1051 1052 1053 1054 1055 1056 1057 1058 1059 1060 1061 1062 1063 1064 1065 1066 1067 1068 1069 1070 1071 |
# File 'lib/flt/support.rb', line 1051 def scale! r, s, m_p, m_m, k,output_b = @r, @s, @m_p, @m_m, @k,@output_b loop do if (@round_h ? (r+m_p >= s) : (r+m_p > s)) # k is too low s *= output_b k += 1 elsif (@round_h ? ((r+m_p)*output_b<s) : ((r+m_p)*output_b<=s)) # k is too high r *= output_b m_p *= output_b m_m *= output_b k -= 1 else @s = s @r = r @m_p = m_p @m_m = m_m @k = k break end end end |
#scale_fixup! ⇒ Object
fix up scaling (final step): specialized version of scale! This performs a single up scaling step, i.e. behaves like scale2, but the input must be at most one step down from the final result
1231 1232 1233 1234 1235 1236 |
# File 'lib/flt/support.rb', line 1231 def scale_fixup! if (@round_h ? (@r+@m_p >= @s) : (@r+@m_p > @s)) # too low? @s *= @output_b @k += 1 end end |
#scale_optimized! ⇒ Object
scale_o1! is an optimized version of scale!; it requires an additional parameters with the floating-point number v=r/s
It uses a Float estimate of ceil(logB(v)) that may need to adjusted one unit up TODO: find easy to use estimate; determine max distance to correct value and use it for fixing,
or use the general scale! for fixing (but remembar to multiply by exptt(...))
(determine when Math.log is aplicable, etc.)
1165 1166 1167 1168 1169 1170 1171 1172 1173 1174 1175 1176 1177 1178 1179 1180 1181 1182 1183 1184 1185 1186 1187 1188 1189 1190 1191 1192 1193 1194 1195 1196 1197 1198 1199 1200 1201 1202 1203 1204 1205 1206 1207 1208 1209 1210 1211 1212 1213 1214 1215 1216 1217 1218 1219 1220 1221 1222 1223 1224 1225 1226 |
# File 'lib/flt/support.rb', line 1165 def scale_optimized! context = @context # @v.class.context return scale! if context.zero?(@v) # 1. compute estimated_scale # 1.1. try to use Float logarithms (Math.log) v = @v v_abs = context.copy_sign(v, +1) # don't use v.abs because it rounds (and may overflow also) v_flt = v_abs.to_f b = @output_b log_b = ESTIMATE_FLOAT_LOG_B[b] log_b = ESTIMATE_FLOAT_LOG_B[b] = 1.0/Math.log(b) if log_b.nil? estimated_scale = nil fixup = false begin l = ((b==10) ? Math.log10(v_flt) : Math.log(v_flt)*log_b) estimated_scale =(l - 1E-10).ceil fixup = true rescue # rescuing errors is more efficient than checking (v_abs < Float::MAX.to_i) && (v_flt > Float::MIN) when v is a Flt else # estimated_scale = nil end # 1.2. Use Flt::DecNum logarithm if estimated_scale.nil? v.to_decimal_exact(:precision=>12) if v.is_a?(BinNum) if v.is_a?(DecNum) l = nil DecNum.context(:precision=>12) do case b when 10 l = v_abs.log10 else l = v_abs.ln/Flt.DecNum(b).ln end end l -= Flt.DecNum(+1,1,-10) estimated_scale = l.ceil fixup = true end end # 1.3 more rough Float aproximation # TODO: optimize denominator, correct numerator for more precision with first digit or part # of the coefficient (like _log_10_lb) estimated_scale ||= (v.adjusted_exponent.to_f * Math.log(v.class.context.radix) * log_b).ceil if estimated_scale >= 0 @k = estimated_scale @s *= output_b_power(estimated_scale) else sc = output_b_power(-estimated_scale) @k = estimated_scale @r *= sc @m_p *= sc @m_m *= sc end fixup ? scale_fixup! : scale! end |
#scale_original!(really = false) ⇒ Object
Given r/s = v (number to convert to text), m_m/s = (v - v-)/s, m_p/s = (v+ - v)/s Scale the fractions so that the first significant digit is right after the radix point, i.e. find k = ceil(logB((r+m_p)/s)), the smallest integer such that (r+m_p)/s <= B^k if k>=0 return:
r=r, s=s*B^k, m_p=m_p, m_m=m_m
if k<0 return:
r=r*B^k, s=s, m_p=m_p*B^k, m_m=m_m*B^k
scale! is a general iterative method using only (multiprecision) integer arithmetic.
1035 1036 1037 1038 1039 1040 1041 1042 1043 1044 1045 1046 1047 1048 1049 |
# File 'lib/flt/support.rb', line 1035 def scale_original!(really=false) loop do if (@round_h ? (@r+@m_p >= @s) : (@r+@m_p > @s)) # k is too low @s *= @output_b @k += 1 elsif (@round_h ? ((@r+@m_p)*@output_b<@s) : ((@r+@m_p)*@output_b<=@s)) # k is too high @r *= @output_b @m_p *= @output_b @m_m *= @output_b @k -= 1 else break end end end |