Inherits:
Object
show all
Defined in:

## Overview

Clinger algorithms to read floating point numbers from text literals with correct rounding. from his paper: “How to Read Floating Point Numbers Accurately” (William D. Clinger)

## Class Method Summary collapse

• Minimum & maximum adjusted exponent for numbers in base to be in the range of Floats.

• Given exact positive integers u and v with beta**(n-1) <= u/v < beta**n and exact integer k, returns the floating point number closest to u/v * beta**n (beta is the floating-point radix).

## Instance Method Summary collapse

• Algorithm M to read floating point numbers from text literals with correct rounding from his paper: “How to Read Floating Point Numbers Accurately” (William D. Clinger).

• Fast for Float.

• constructor

There are three different reading approaches, selected by the :mode parameter: * :fixed (the destination context defines the resulting precision) input is rounded as specified by the context; if the context precision is ‘exact’, the exact input value will be represented in the destination base, which can lead to a Inexact exception (or a NaN result and an Inexact flag) * :free The input precision is preserved, and the destination context precision is ignored; in this case the result can be converted back to the original number (with the same precision) a rounding mode for the back conversion may be passed; otherwise any round-to-nearest is assumed.

• Given exact integers f and e, with f nonnegative, returns the floating-point number closest to f * eb**e (eb is the input radix).

## Constructor Details

### #initialize(options = {}) ⇒ Reader

There are three different reading approaches, selected by the :mode parameter:

• :fixed (the destination context defines the resulting precision) input is rounded as specified by the context; if the context precision is ‘exact’, the exact input value will be represented in the destination base, which can lead to a Inexact exception (or a NaN result and an Inexact flag)

• :free The input precision is preserved, and the destination context precision is ignored; in this case the result can be converted back to the original number (with the same precision) a rounding mode for the back conversion may be passed; otherwise any round-to-nearest is assumed. (to increase the precision of the result the input precision must be increased –adding trailing zeros)

• :short is like :free, but the minumum number of digits that preserve the original value are generated (with :free, all significant digits are generated)

For the fixed mode there are three conversion algorithms available that can be selected with the :algorithm parameter:

• :A Arithmetic algorithm, using correctly rounded Flt::Num arithmetic.

• :M The Clinger Algorithm M is the slowest method, but it was the first implemented and testes and is kept as a reference for testing.

• :R The Clinger Algorithm R, which requires an initial approximation is currently only implemented for Float and is the fastest by far.

 ``` 47 48 49 50 51``` ```# File 'lib/flt/support/reader.rb', line 47 def initialize(options={}) @exact = nil @algorithm = options[:algorithm] @mode = options[:mode] || :fixed end```

## Class Method Details

### .float_min_max_adj_exp(base, normalized = false) ⇒ Object

Minimum & maximum adjusted exponent for numbers in base to be in the range of Floats

 ``` 234 235 236 237 238 239 240 241 242 243 244``` ```# File 'lib/flt/support/reader.rb', line 234 def float_min_max_adj_exp(base, normalized=false) k = normalized ? base : -base unless min_max = @float_min_max_exp_values[k] max_exp = (Math.log(Float::MAX)/Math.log(base)).floor e = Float::MIN_EXP e -= Float::MANT_DIG unless normalized min_exp = (e*Math.log(Float::RADIX)/Math.log(base)).ceil @float_min_max_exp_values[k] = min_max = [min_exp, max_exp] end min_max.map{|exp| exp - 1} # adjust end```

### .ratio_float(context, u, v, k, round_mode) ⇒ Object

Given exact positive integers u and v with beta**(n-1) <= u/v < beta**n and exact integer k, returns the floating point number closest to u/v * beta**n (beta is the floating-point radix)

 ``` 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422``` ```# File 'lib/flt/support/reader.rb', line 397 def self.ratio_float(context, u, v, k, round_mode) # since this handles only positive numbers and ceiling and floor # are not symmetrical, they should have been swapped before calling this. q = u.div v r = u-q*v v_r = v-r z = context.Num(+1,q,k) exact = (r==0) if round_mode == :down # z = z elsif (round_mode == :up) && r>0 z = context.next_plus(z) elsif rv_r z = context.next_plus(z) else # tie if (round_mode == :half_down) || (round_mode == :half_even && ((q%2)==0)) || (round_mode == :down) # z = z else z = context.next_plus(z) end end return z, exact end```

## Instance Method Details

### #_alg_m(context, round_mode, sign, f, e, eb, n) ⇒ Object

Algorithm M to read floating point numbers from text literals with correct rounding from his paper: “How to Read Floating Point Numbers Accurately” (William D. Clinger)

 ``` 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392``` ```# File 'lib/flt/support/reader.rb', line 356 def _alg_m(context, round_mode, sign, f, e, eb, n) if e<0 u,v,k = f,eb**(-e),0 else u,v,k = f*(eb**e),1,0 end min_e = context.etiny max_e = context.etop rp_n = context.int_radix_power(n) rp_n_1 = context.int_radix_power(n-1) r = context.radix loop do x = u.div(v) # bottleneck if (x>=rp_n_1 && x context.maximum_coefficient context.exception(Num::Overflow,"Input literal out of range") end end context.exception Num::Inexact if !exact end return z.copy_sign(sign) elsif x=rp_n v *= r k += 1 end end end```

### #_alg_r(z0, context, round_mode, sign, f, e, eb, n) ⇒ Object

Fast for Float

 ``` 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225``` ```# File 'lib/flt/support/reader.rb', line 199 def _alg_r(z0, context, round_mode, sign, f, e, eb, n) # Fast for Float #raise InvalidArgument, "Reader Algorithm R only supports base 2" if context.radix != 2 @z = z0 @r = context.radix @rp_n_1 = context.int_radix_power(n-1) @round_mode = round_mode ret = nil loop do m, k = context.to_int_scale(@z) # TODO: replace call to compare by setting the parameters in local variables, # then insert the body of compare here; # then eliminate innecesary instance variables if e >= 0 && k >= 0 ret = compare m, f*eb**e, m*@r**k, context elsif e >= 0 && k < 0 ret = compare m, f*eb**e*@r**(-k), m, context elsif e < 0 && k >= 0 ret = compare m, f, m*@r**k*eb**(-e), context else # e < 0 && k < 0 ret = compare m, f*@r**(-k), m*eb**(-e), context end break if ret end ret && context.copy_sign(ret, sign) # TODO: normalize? end```

### #_alg_r_approx(context, round_mode, sign, f, e, eb, n) ⇒ Object

 ``` 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197``` ```# File 'lib/flt/support/reader.rb', line 148 def _alg_r_approx(context, round_mode, sign, f, e, eb, n) return nil if context.radix != Float::RADIX || context.exact? || context.precision > Float::MANT_DIG # Compute initial approximation; if Float uses IEEE-754 binary arithmetic, the approximation # is good enough to be adjusted in just one step. @good_approx = true ndigits = Support::AuxiliarFunctions._ndigits(f, eb) adj_exp = e + ndigits - 1 min_exp, max_exp = Reader.float_min_max_adj_exp(eb) if adj_exp >= min_exp && adj_exp <= max_exp if eb==2 z0 = Math.ldexp(f,e) elsif eb==10 unless Flt.float_correctly_rounded? min_exp_norm, max_exp_norm = Reader.float_min_max_adj_exp(eb, true) @good_approx = false return nil if e <= min_exp_norm end z0 = Float("#{f}E#{e}") else ff = f ee = e min_exp_norm, max_exp_norm = Reader.float_min_max_adj_exp(eb, true) if e <= min_exp_norm # avoid loss of precision due to gradual underflow return nil if e <= min_exp @good_approx = false ff = Float(f)*Float(eb)**(e-min_exp_norm-1) ee = min_exp_norm + 1 end # if ee < 0 # z0 = Float(ff)/Float(eb**(-ee)) # else # z0 = Float(ff)*Float(eb**ee) # end z0 = Float(ff)*Float(eb)**ee end if z0 && context.num_class != Float @good_approx = false z0 = context.Num(z0).plus(context) # context.plus(z0) ? else z0 = context.Num(z0) end end end```

### #compare(m, x, y, context) ⇒ Object

 ``` 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352``` ```# File 'lib/flt/support/reader.rb', line 247 def compare(m, x, y, context) ret = nil d = x-y d2 = 2*m*d.abs # v = f*eb**e is the number to be approximated # z = m*@r**k is the current aproximation # the error of @z is eps = abs(v-z) = 1/2 * d2 / y # we have x, y integers such that x/y = v/z # so eps < 1/2 <=> d2 < y # d < 0 <=> x < y <=> v < z directed_rounding = [:up, :down].include?(@round_mode) if directed_rounding if @round_mode==:up ? (d <= 0) : (d < 0) # v <(=) z chk = (m == @rp_n_1) ? d2*@r : d2 if (@round_mode == :up) && (chk < 2*y) # eps < 1 ret = @z else @z = context.next_minus(@z) end else # @round_mode==:up ? (d > 0) : (d >= 0) # v >(=) z if (@round_mode == :down) && (d2 < 2*y) # eps < 1 ret = @z else @z = context.next_plus(@z) end end else if d2 < y # eps < 1/2 if (m == @rp_n_1) && (d < 0) && (y < @r*d2) # z has the minimum normalized significand, i.e. is a power of @r # and v < z # and @r*eps > 1/2 # On the left of z the ulp is 1/@r than the ulp on the right; if v < z we # must require an error @r times smaller. @z = context.next_minus(@z) else # unambiguous nearest ret = @z end elsif d2 == y # eps == 1/2 # round-to-nearest tie if @round_mode == :half_even if (m%2) == 0 # m is even if (m == @rp_n_1) && (d < 0) # z is power of @r and v < z; this wasn't really a tie because # there are closer values on the left @z = context.next_minus(@z) else # m is even => round tie to z ret = @z end elsif d < 0 # m is odd, v < z => round tie to prev ret = context.next_minus(@z) elsif d > 0 # m is odd, v > z => round tie to next ret = context.next_plus(@z) end elsif @round_mode == :half_up if d < 0 # v < z if (m == @rp_n_1) # this was not really a tie @z = context.next_minus(@z) else ret = @z end else # d > 0 # v >= z ret = context.next_plus(@z) end else # @round_mode == :half_down if d < 0 # v < z if (m == @rp_n_1) # this was not really a tie @z = context.next_minus(@z) else ret = context.next_minus(@z) end else # d < 0 # v > z ret = @z end end elsif d < 0 # eps > 1/2 and v < z @z = context.next_minus(@z) elsif d > 0 # eps > 1/2 and v > z @z = context.next_plus(@z) end end # Assume the initial approx is good enough (uses IEEE-754 arithmetic with round-to-nearest), # so we can avoid further iteration, except for directed rounding ret ||= @z unless directed_rounding || !@good_approx return ret end```

### #exact? ⇒ Boolean

Returns:

• (Boolean)
 ``` 53 54 55``` ```# File 'lib/flt/support/reader.rb', line 53 def exact? @exact end```

### #read(context, round_mode, sign, f, e, eb = 10) ⇒ Object

Given exact integers f and e, with f nonnegative, returns the floating-point number closest to f * eb**e (eb is the input radix)

If the context precision is exact an Inexact exception may occur (an NaN be returned) if an exact conversion is not possible.

round_mode: in :fixed mode it specifies how to round the result (to the context precision); it is passed separate from context for flexibility. in :free mode it specifies what rounding would be used to convert back the output to the input base eb (using the same precision that f has).

 ``` 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146``` ```# File 'lib/flt/support/reader.rb', line 68 def read(context, round_mode, sign, f, e, eb=10) @exact = true case @mode when :free, :short all_digits = (@mode == :free) # for free mode, (any) :nearest rounding is used by default Num.convert(Num[eb].Num(sign, f, e), context.num_class, :rounding=>round_mode||:nearest, :all_digits=>all_digits) when :fixed if context.exact? a,b = [eb, context.radix].sort m = (Math.log(b)/Math.log(a)).round if b == a**m # conmensurable bases if eb > context.radix n = AuxiliarFunctions._ndigits(f, eb)*m else n = (AuxiliarFunctions._ndigits(f, eb)+m-1)/m end else # inconmesurable bases; exact result may not be possible x = Num[eb].Num(sign, f, e) x = Num.convert_exact(x, context.num_class, context) @exact = !x.nan? return x end else n = context.precision end if round_mode == :nearest # :nearest is not meaningful here in :fixed mode; replace it if [:half_even, :half_up, :half_down].include?(context.rounding) round_mode = context.rounding else round_mode = :half_even end end # for fixed mode, use the context rounding by default round_mode ||= context.rounding alg = @algorithm if (context.radix == 2 && alg.nil?) || alg==:R z0 = _alg_r_approx(context, round_mode, sign, f, e, eb, n) alg = z0 && :R end alg ||= :A case alg when :M, :R round_mode = Support.simplified_round_mode(round_mode, sign == -1) case alg when :M _alg_m(context, round_mode, sign, f, e, eb, n) when :R _alg_r(z0, context, round_mode, sign, f, e, eb, n) end else # :A # direct arithmetic conversion if round_mode == context.rounding x = Num.convert_exact(Num[eb].Num(sign, f, e), context.num_class, context) x = context.normalize(x) unless !context.respond_to?(:normalize) || context.exact? x else if context.num_class == Float float = true context = BinNum::FloatContext end x = context.num_class.context(context) do |local_context| local_context.rounding = round_mode Num.convert_exact(Num[eb].Num(sign, f, e), local_context.num_class, local_context) end if float x = x.to_f else x = context.normalize(x) unless context.exact? end x end end end end```