Class: Flt::Support::Reader

Inherits:
Object show all
Defined in:
lib/flt/support/reader.rb

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

Instance Method Summary collapse

Constructor Details

#initialize(options = {}) ⇒ Reader

There are two 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 r<v_r
    # z = z
  elsif r>v_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<rp_n) || k==min_e || k==max_e
        z, exact = Reader.ratio_float(context,u,v,k,round_mode)
        @exact = exact
        if context.respond_to?(:exception)
          if k==min_e
            context.exception(Num::Subnormal) if z.subnormal?
            context.exception(Num::Underflow,"Input literal out of range") if z.zero? && f!=0
          elsif k==max_e
            if !context.exact? && z.coefficient > 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_1
       u *= r
       k -= 1
     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 exact_mode = 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