Parent

Class Index [+]

Quicksearch

Flt::Support::Formatter

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: @v - The number to be formatted = @f*@b**@e @b - The numberic base of the input floating-point representation of @v @f - The significand or characteristic (fraction) @e - The exponent

Quotients of integers will be used to hold the magnitudes: @s is the denominator of all fractions @r numerator of @v: @v = @r/@s @m_m numerator of the distance from the rounding-range lower limit, l, to @v: @m_m/@s = (@v - l) @m_p numerator of the distance from @v to the rounding-range upper limit, u: @m_p/@s = (u - @v) All numbers in the randound-range are rounded to @v (with the given precision p) @k scale factor that is applied to the quotients @r/@s, @m_m/@s and @m_p/@s to put the first 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+] @roundl is true if the lower limit of the rounding range is closed (i.e., if l rounds to @v) @roundh is true if the upper limit of the rounding range is closed (i.e., if u rounds to @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))

@output_b is the output base @output_min_e is the output minimum exponent p is the input floating point precision

Constants

ESTIMATE_FLOAT_LOG_B

Attributes

round_up[R]

Public Class Methods

new(input_b, input_min_e, output_b) click to toggle source

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.

     # File lib/flt/support.rb, line 832
832:       def initialize(input_b, input_min_e, output_b)
833:         @b = input_b
834:         @min_e = input_min_e
835:         @output_b = output_b
836:         # result of last operation
837:         @adjusted_digits = @digits = nil
838:         # for "all-digits" mode results (which are truncated, rather than rounded),
839:         # round_up contains information to round the result:
840:         # * it is nil if the rest of digits are zero (the result is exact)
841:         # * it is :lo if there exist non-zero digits beyond the significant ones (those returned), but
842:         #   the value is below the tie (the value must be rounded up only for :up rounding mode)
843:         # * it is :tie if there exists exactly one nonzero digit after the significant and it is radix/2,
844:         #   for round-to-nearest it is atie.
845:         # * it is :hi otherwise (the value should be rounded-up except for the :down mode)
846:         @round_up = nil
847:       end

Public Instance Methods

adjusted_digits(round_mode) click to toggle source

Access rounded result of format operation: scaling (position of radix point) and digits

      # File lib/flt/support.rb, line 983
 983:       def adjusted_digits(round_mode)
 984:         round_mode = Support.simplified_round_mode(round_mode, @minus)
 985:         if @adjusted_digits.nil? && !@digits.nil?
 986:           increment = (@round_up && (round_mode != :down)) &&
 987:                       ((round_mode == :up) ||
 988:                       (@round_up == :hi) ||
 989:                       ((@round_up == :tie) &&
 990:                        ((round_mode==:half_up) || ((round_mode==:half_even) && ((@digits.last % 2)==1)))))
 991:           # increment = (@round_up == :tie) || (@round_up == :hi) # old behaviour (:half_up)
 992:           if increment
 993:             base = @output_b
 994:             dec_pos = @k
 995:             digits = @digits.dup
 996:             # carry = increment ? 1 : 0
 997:             # digits = digits.reverse.map{|d| d += carry; d>=base ? 0 : (carry=0;d)}.reverse
 998:             # if carry != 0
 999:             #   digits.unshift carry
1000:             #   dec_pos += 1
1001:             # end
1002:             i = digits.size - 1
1003:             while i>=0
1004:               digits[i] += 1
1005:               if digits[i] == base
1006:                 digits[i] = 0
1007:               else
1008:                 break
1009:               end
1010:               i -= 1
1011:             end
1012:             if i<0
1013:               dec_pos += 1
1014:               digits.unshift 1
1015:             end
1016:             @adjusted_k = dec_pos
1017:             @adjusted_digits = digits
1018:           else
1019:             @adjusted_k = @k
1020:             @adjusted_digits = @digits
1021:           end
1022:         end
1023:         return @adjusted_k, @adjusted_digits
1024:       end
b_power(n) click to toggle source
      # File lib/flt/support.rb, line 1073
1073:       def b_power(n)
1074:         @b**n
1075:       end
digits() click to toggle source

Access result of format operation: scaling (position of radix point) and digits

     # File lib/flt/support.rb, line 975
975:       def digits
976:         return @k, @digits
977:       end
format(v, f, e, round_mode, p=nil, all=false) click to toggle source

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.

     # File lib/flt/support.rb, line 886
886:       def format(v, f, e, round_mode, p=nil, all=false)
887:         context = v.class.context
888:         # TODO: consider removing parameters f,e and using v.split instead
889:         @minus = (context.sign(v)==1)
890:         @v = context.copy_sign(v, 1) # don't use context.abs(v) because it rounds (and may overflow also)
891:         @f = f.abs
892:         @e = e
893:         @round_mode = round_mode
894:         @all_digits = all
895:         p ||= context.precision
896: 
897:         # adjust the rounding mode to work only with positive numbers
898:         @round_mode = Support.simplified_round_mode(@round_mode, @minus)
899: 
900:         # determine the high,low inclusion flags of the rounding limits
901:         case @round_mode
902:           when :half_even
903:             # rounding rage is (v-m-,v+m+) if v is odd and [v+m-,v+m+] if even
904:             @round_l = @round_h = ((@f%2)==0)
905:           when :up
906:             # rounding rage is (v-,v]
907:             # ceiling is treated here assuming f>0
908:             @round_l, @round_h = false, true
909:           when :down
910:             # rounding rage is [v,v+)
911:             # floor is treated here assuming f>0
912:             @round_l, @round_h = true, false
913:           when :half_up
914:             # rounding rage is [v+m-,v+m+)
915:             @round_l, @round_h = true, false
916:           when :half_down
917:             # rounding rage is (v+m-,v+m+]
918:             @round_l, @round_h = false, true
919:           else # :nearest
920:             # Here assume only that round-to-nearest will be used, but not which variant of it
921:             # The result is valid for any rounding (to nearest) but may produce more digits
922:             # than stricly necessary for specific rounding modes.
923:             # That is, enough digits are generated so that when the result is
924:             # converted to floating point with the specified precision and
925:             # correct rounding (to nearest), the result is the original number.
926:             # rounding range is (v+m-,v+m+)
927:             @round_l = @round_h = false
928:         end
929: 
930:         # TODO: use context.next_minus, next_plus instead of direct computing, don't require min_e & ps
931:         # Now compute the working quotients @r/@s, @m_p/@s = (v+ - @v), @m_m/@s = (@v - v-) and scale them.
932:         if @e >= 0
933:           if @f != b_power(p-1)
934:             be = b_power(@e)
935:             @r, @s, @m_p, @m_m = @f*be*2, 2, be, be
936:           else
937:             be = b_power(@e)
938:             be1 = be*@b
939:             @r, @s, @m_p, @m_m = @f*be1*2, @b*2, be1, be
940:           end
941:         else
942:           if @e==@min_e or @f != b_power(p-1)
943:             @r, @s, @m_p, @m_m = @f*2, b_power(-@e)*2, 1, 1
944:           else
945:             @r, @s, @m_p, @m_m = @f*@b*2, b_power(1-@e)*2, @b, 1
946:           end
947:         end
948:         @k = 0
949:         @context = context
950:         scale_optimized!
951: 
952: 
953:         # The value to be formatted is @v=@r/@s; m- = @m_m/@s = (@v - v-)/@s; m+ = @m_p/@s = (v+ - @v)/@s
954:         # Now adjust @m_m, @m_p so that they define the rounding range
955:         case @round_mode
956:         when :up
957:           # ceiling is treated here assuming @f>0
958:           # rounding range is -v,@v
959:           @m_m, @m_p = @m_m*2, 0
960:         when :down
961:           # floor is treated here assuming #f>0
962:           # rounding range is @v,v+
963:           @m_m, @m_p = 0, @m_p*2
964:         else
965:           # rounding range is v-,v+
966:           # @m_m, @m_p = @m_m, @m_p
967:         end
968: 
969:         # Now m_m, m_p define the rounding range
970:         all ? generate_max : generate
971: 
972:       end
generate() click to toggle source
      # File lib/flt/support.rb, line 1121
1121:       def generate
1122:         list = []
1123:         r, s, m_p, m_m, = @r, @s, @m_p, @m_m
1124:         loop do
1125:           d,r = (r*@output_b).divmod(s)
1126:           m_p *= @output_b
1127:           m_m *= @output_b
1128:           tc1 = @round_l ? (r<=m_m) : (r<m_m)
1129:           tc2 = @round_h ? (r+m_p >= s) : (r+m_p > s)
1130: 
1131:           if not tc1
1132:             if not tc2
1133:               list << d
1134:             else
1135:               list << d+1
1136:               break
1137:             end
1138:           else
1139:             if not tc2
1140:               list << d
1141:               break
1142:             else
1143:               if r*2 < s
1144:                 list << d
1145:                 break
1146:               else
1147:                 list << d+1
1148:                 break
1149:               end
1150:             end
1151:           end
1152: 
1153:         end
1154:         @digits = list
1155:       end
generate_max() click to toggle source
      # File lib/flt/support.rb, line 1081
1081:       def generate_max
1082:         @round_up = false
1083:         list = []
1084:         r, s, m_p, m_m, = @r, @s, @m_p, @m_m
1085:         n_iters, rs = 0, []
1086:         loop do
1087:           if (n_iters > 10000)
1088:             raise "Infinite digit sequence." if rs.include?(r)
1089:             rs << r
1090:           else
1091:             n_iters += 1
1092:           end
1093: 
1094:           d,r = (r*@output_b).divmod(s)
1095: 
1096:           m_p *= @output_b
1097:           m_m *= @output_b
1098: 
1099:           list << d
1100: 
1101:           tc1 = @round_l ? (r<=m_m) : (r<m_m)
1102:           tc2 = @round_h ? (r+m_p >= s) : (r+m_p > s)
1103: 
1104:           if tc1 && tc2
1105:             if r != 0
1106:               r *= 2
1107:               if r > s
1108:                 @round_up = :hi
1109:               elsif r == s
1110:                 @round_up = :tie
1111:               else
1112:                 @rund_up = :lo
1113:               end
1114:             end
1115:             break
1116:           end
1117:         end
1118:         @digits = list
1119:       end
output_b_power(n) click to toggle source
      # File lib/flt/support.rb, line 1077
1077:       def output_b_power(n)
1078:         @output_b**n
1079:       end
scale!() click to toggle source

using local vars instead of instance vars: it makes a difference in performance

      # File lib/flt/support.rb, line 1051
1051:       def scale!
1052:         r, s, m_p, m_m, k,output_b = @r, @s, @m_p, @m_m, @k,@output_b
1053:         loop do
1054:           if (@round_h ? (r+m_p >= s) : (r+m_p > s)) # k is too low
1055:             s *= output_b
1056:             k += 1
1057:           elsif (@round_h ? ((r+m_p)*output_b<s) : ((r+m_p)*output_b<=s)) # k is too high
1058:             r *= output_b
1059:             m_p *= output_b
1060:             m_m *= output_b
1061:             k -= 1
1062:           else
1063:             @s = s
1064:             @r = r
1065:             @m_p = m_p
1066:             @m_m = m_m
1067:             @k = k
1068:             break
1069:           end
1070:         end
1071:       end
scale_fixup!() click to toggle source

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

      # File lib/flt/support.rb, line 1231
1231:       def scale_fixup!
1232:         if (@round_h ? (@r+@m_p >= @s) : (@r+@m_p > @s)) # too low?
1233:           @s *= @output_b
1234:           @k += 1
1235:         end
1236:       end
scale_optimized!() click to toggle source

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.)
      # File lib/flt/support.rb, line 1165
1165:       def scale_optimized!
1166:         context = @context # @v.class.context
1167:         return scale! if context.zero?(@v)
1168: 
1169:         # 1. compute estimated_scale
1170: 
1171:         # 1.1. try to use Float logarithms (Math.log)
1172:         v = @v
1173:         v_abs = context.copy_sign(v, 1) # don't use v.abs because it rounds (and may overflow also)
1174:         v_flt = v_abs.to_f
1175:         b = @output_b
1176:         log_b = ESTIMATE_FLOAT_LOG_B[b]
1177:         log_b = ESTIMATE_FLOAT_LOG_B[b] = 1.0/Math.log(b) if log_b.nil?
1178:         estimated_scale = nil
1179:         fixup = false
1180:         begin
1181:           l = ((b==10) ? Math.log10(v_flt) : Math.log(v_flt)*log_b)
1182:           estimated_scale =(l - 1E-10).ceil
1183:           fixup = true
1184:         rescue
1185:           # rescuing errors is more efficient than checking (v_abs < Float::MAX.to_i) && (v_flt > Float::MIN) when v is a Flt
1186:         else
1187:           # estimated_scale = nil
1188:         end
1189: 
1190:         # 1.2. Use Flt::DecNum logarithm
1191:         if estimated_scale.nil?
1192:           v.to_decimal_exact(:precision=>12) if v.is_a?(BinNum)
1193:           if v.is_a?(DecNum)
1194:             l = nil
1195:             DecNum.context(:precision=>12) do
1196:               case b
1197:               when 10
1198:                 l = v_abs.log10
1199:               else
1200:                 l = v_abs.ln/Flt.DecNum(b).ln
1201:               end
1202:             end
1203:             l -= Flt.DecNum(1,1,10)
1204:             estimated_scale = l.ceil
1205:             fixup = true
1206:           end
1207:         end
1208: 
1209:         # 1.3 more rough Float aproximation
1210:           # TODO: optimize denominator, correct numerator for more precision with first digit or part
1211:           # of the coefficient (like _log_10_lb)
1212:         estimated_scale ||= (v.adjusted_exponent.to_f * Math.log(v.class.context.radix) * log_b).ceil
1213: 
1214:         if estimated_scale >= 0
1215:           @k = estimated_scale
1216:           @s *= output_b_power(estimated_scale)
1217:         else
1218:           sc = output_b_power(-estimated_scale)
1219:           @k = estimated_scale
1220:           @r *= sc
1221:           @m_p *= sc
1222:           @m_m *= sc
1223:         end
1224:         fixup ? scale_fixup! : scale!
1225: 
1226:       end
scale_original!(really=false) click to toggle source

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.

      # File lib/flt/support.rb, line 1035
1035:       def scale_original!(really=false)
1036:         loop do
1037:           if (@round_h ? (@r+@m_p >= @s) : (@r+@m_p > @s)) # k is too low
1038:             @s *= @output_b
1039:             @k += 1
1040:           elsif (@round_h ? ((@r+@m_p)*@output_b<@s) : ((@r+@m_p)*@output_b<=@s)) # k is too high
1041:             @r *= @output_b
1042:             @m_p *= @output_b
1043:             @m_m *= @output_b
1044:             @k -= 1
1045:           else
1046:             break
1047:           end
1048:         end
1049:       end

Disabled; run with --debug to generate this.

[Validate]

Generated with the Darkfish Rdoc Generator 1.1.6.