Class: Flt::Support::Rationalizer

Inherits:
Object
  • Object
show all
Extended by:
AuxiliarFunctions
Includes:
AuxiliarFunctions
Defined in:
lib/flt/support/rationalizer.rb,
lib/flt/support/rationalizer_extra.rb

Overview

This class provides efficient conversion of fraction (as approximate floating point numbers) to rational numbers.

Defined Under Namespace

Modules: AuxiliarFunctions

Class Method Summary collapse

Instance Method Summary collapse

Constructor Details

#initialize(tol = Tolerance(:epsilon)) ⇒ Rationalizer

Create Rationalizator with given tolerance.



113
114
115
 
# File 'lib/flt/support/rationalizer.rb', line 113

def initialize(tol=Tolerance(:epsilon))
  @tol = tol
end

Class Method Details

.[](*args) ⇒ Object 



117
118
119
 
# File 'lib/flt/support/rationalizer.rb', line 117

def self.[](*args)
  new *args
end

.max_denominator(f, max_den = 1000000000, num_class = nil) ⇒ Object

Best fraction given maximum denominator Algorithm Copyright © 1991 by Joseph K. Horn.

The implementation of this method uses floating point arithmetic which limits the magnitude and precision of the results, specially using Float values.



242
243
244
245
246
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
 
# File 'lib/flt/support/rationalizer.rb', line 242

def self.max_denominator(f, max_den=1000000000, num_class=nil)
  return rationalize_special(f) if special?(f)
  return nil if max_den < 1
  num_class ||= f.class
  context = num_class.context
  return ip(f),1 if fp(f) == 0

  cast = lambda{|x| context.Num(x)}

  one = cast[1]

   sign = f < 0
   f = -f if sign

   a,b,c = 0,1,f
   while b < max_den && c != 0
     cc = one/c
     a,b,c = b, ip(cc)*b+a, fp(cc)
   end

   if b>max_den
     b -= a*ceil(cast[b-max_den]/a)
   end

   f1,f2 = [a,b].collect{|x| abs(cast[rnd(x*f)]/x-f)}

   a = f1 > f2 ? b : a

   num,den = rnd(a*f).to_i,a
   den = 1 if abs(den) < 1

   num = -num if sign

  return num,den
end

.rationalize_special(x) ⇒ Object 



278
279
280
281
282
283
284
 
# File 'lib/flt/support/rationalizer.rb', line 278

def self.rationalize_special(x)
  if x.nan?
    [0, 0]
  else
    [sign(x), 0]
  end
end

.to_r(x) ⇒ Object

Exact conversion to rational. Ruby provides this method for all numeric types since version 1.9.1, but before that it wasn’t available for Float or BigDecimal. This methods supports old Ruby versions.



10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
 
# File 'lib/flt/support/rationalizer.rb', line 10

def self.to_r(x)
  if x.respond_to?(:to_r)
    x.to_r
  else
    case x
    when Float
      # Float did not had a #to_r method until Ruby 1.9.1
      return Rational(x.to_i, 1) if x.modulo(1) == 0
      if !x.finite?
        return Rational(0, 0) if x.nan?
        return x < 0 ? Rational(-1, 0) : Rational(1, 0)
      end

      f, e = Math.frexp(x)

      if e < Float::MIN_EXP
         bits = e + Float::MANT_DIG - Float::MIN_EXP
      else
         bits = [Float::MANT_DIG,e].max
         # return Rational(x.to_i, 1) if bits < e
      end
        p = Math.ldexp(f, bits)
        e = bits - e
        if e < Float::MAX_EXP
          q = Math.ldexp(1, e)
        else
          q = Float::RADIX**e
        end
      Rational(p.to_i, q.to_i)

    when BigDecimal
      # BigDecimal probably didn't have #to_r at some point
      s, f, b, e = x.split
      p = f.to_i
      p = -p if s < 0
      e = f.size - e
      if e < 0
        p *= b**(-e)
        e = 0
      end
      q = b**(e)
      Rational(p,q)

    else
      x.to_r
    end
  end
end

Instance Method Details

#exact_binary_rationalization(x) ⇒ Object

An exact rationalization method for binary floating point that yields smallest fractions when possible and is not too slow



147
148
149
150
151
152
153
154
 
# File 'lib/flt/support/rationalizer_extra.rb', line 147

def exact_binary_rationalization(x)
  p, q = x, 1
  while p.modulo(1) != 0
    p *= 2.0
    q <<= 1 # q *= 2
  end
  Rational(p.to_i, q)
end

#exact_float_rationalization(x) ⇒ Object

An a here’s a shorter implementation relying on the semantics of the power operator, but which is somewhat slow:



158
159
160
161
162
163
 
# File 'lib/flt/support/rationalizer_extra.rb', line 158

def exact_float_rationalization(x)
  f,e = Math.frexp(x)
  f = Math.ldexp(f, Float::MANT_DIG)
  e -= Float::MANT_DIG
  return Rational(f.to_i*(Float::RADIX**e.to_i), 1)
end

#rationalize(x) ⇒ Object

Rationalization method that finds the fraction with smallest denominator fraction within the tolerance distance of an approximate (floating point) number.



125
126
127
128
 
# File 'lib/flt/support/rationalizer.rb', line 125

def rationalize(x)
  # Use the algorithm which has been found most efficient, rationalize_Knuth.
  rationalize_Knuth(x)
end

#rationalize_Horn(x) ⇒ Object

This is algorithm PDQ2 by Joe Horn.



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
198
 
# File 'lib/flt/support/rationalizer.rb', line 165

def rationalize_Horn(x)
  rationalization(x) do |z, t|
    a,b = num_den(t)
    n0,d0 = (n,d = num_den(z))
    cn,x,pn,cd,y,pd,lo,hi,mid,q,r = 1,1,0,0,0,1,0,1,1,0,0
    begin
      q,r = n.divmod(d)
      x = q*cn+pn
      y = q*cd+pd
      pn = cn
      cn = x
      pd = cd
      cd = y
      n,d = d,r
    end until b*(n0*y-d0*x).abs <= a*d0*y

    if q > 1
      hi = q
      begin
        mid = (lo + hi).div(2)
        x = cn - pn*mid
        y = cd - pd*mid
        if b*(n0*y - d0*x).abs <= a*d0*y
          lo = mid
        else
          hi = mid
        end
      end until hi - lo <= 1
      x = cn - pn*lo
      y = cd - pd*lo
    end
    [x, y]
  end
end

#rationalize_Horn_Hutchins(x) ⇒ Object

Smallest denominator rationalization procedure by Joe Horn and Tony Hutchins; this is the most efficient method as implemented in RPL. Tony Hutchins has come up with PDR6, an improvement over PDQ2; though benchmarking does not show any speed improvement under Ruby.



50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
 
# File 'lib/flt/support/rationalizer_extra.rb', line 50

def rationalize_Horn_Hutchins(x)
  rationalization(x) do |x, dx|
    a,b = num_den(dx)
    n,d = num_den(x)
    pc,ce = n,-d
    pc,cd = 1,0
    t = a*b
    begin
      tt = (-pe).div(ce)
      pd,cd = cd,pd+tt*cd
      pe,ce = ce,pe+tt*ce
    end until b*ce.abs <= t*cd
    tt = t * (pe<0 ? -1 : (pe>0 ? +1 : 0))
    tt = (tt*d+b*ce).div(tt*pd+b*pe)
    [(n*cd-ce-(n*pd-pe)*tt)/d, tt/(cd-tt*pd)]
  end
end

#rationalize_Horn_simple(x, smallest_denominator = false) ⇒ Object

Simple Rationalization by Joe Horn



9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
 
# File 'lib/flt/support/rationalizer_extra.rb', line 9

def rationalize_Horn_simple(x, smallest_denominator = false)
  rationalization(x) do |z, t|
    a,b = num_den(t)
    n0,d0 = (n,d = z.nio_xr.nio_num_den)
    cn,x,pn,cd,y,pd,lo,hi,mid,q,r = 1,1,0,0,0,1,0,1,1,0,0
    begin
      q,r = n.divmod(d)
      x = q*cn+pn
      y = q*cd+pd
      pn = cn
      cn = x
      pd = cd
      cd = y
      n,d = d,r
    end until b*(n0*y-d0*x).abs <= a*d0*y
    if smallest_denominator
      if q>1
        hi = q
        begin
          mid = (lo+hi).div(2)
          x = cn-pn*mid
          y = cd-pd*mid
          if b*(n0*y-d0*x).abs <= a*d0*y
            lo = mid
          else
            hi = mid
          end
        end until hi-lo <= 1
        x = cn - pn*lo
        y = cd - pd*lo
      end
    end
    [x, y]
  end
end

#rationalize_HornHutchins(x) ⇒ Object

This is from a RPL program by Tony Hutchins (PDR6).



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
226
227
228
229
230
231
232
233
234
 
# File 'lib/flt/support/rationalizer.rb', line 201

def rationalize_HornHutchins(x)
  rationalization(x) do |z, t|
    a,b = num_den(t)
    n0,d0 = (n,d = num_den(z))
    cn,x,pn,cd,y,pd,lo,hi,mid,q,r = 1,1,0,0,0,1,0,1,1,0,0
    begin
      q,r = n.divmod(d)
      x = q*cn+pn
      y = q*cd+pd
      pn = cn
      cn = x
      pd = cd
      cd = y
      n,d = d,r
    end until b*(n0*y-d0*x).abs <= a*d0*y

    if q > 1
      hi = q
      begin
        mid = (lo + hi).div(2)
        x = cn - pn*mid
        y = cd - pd*mid
        if b*(n0*y - d0*x).abs <= a*d0*y
          lo = mid
        else
          hi = mid
        end
      end until hi - lo <= 1
      x = cn - pn*lo
      y = cd - pd*lo
    end
    [x, y]
  end
end

#rationalize_Knuth(x) ⇒ Object

This algorithm is derived from exercise 39 of 4.5.3 in “The Art of Computer Programming”, by Donald E. Knuth.



132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
 
# File 'lib/flt/support/rationalizer.rb', line 132

def rationalize_Knuth(x)
  rationalization(x) do |x, dx|
    x     = to_r(x)
    dx    = to_r(dx)
    xp,xq = num_den(x-dx)
    yp,yq = num_den(x+dx)

    a = []
    fin, odd = false, false
    while !fin && xp != 0 && yp != 0
      odd = !odd
      xp,xq = xq,xp
      ax = xp.div(xq)
      xp -= ax*xq

      yp,yq = yq,yp
      ay = yp.div(yq)
      yp -= ay*yq

      if ax!=ay
        fin = true
        ax,xp,xq = ay,yp,yq if odd
      end
      a << ax # .to_i
    end
    a[-1] += 1 if xp != 0 && a.size > 0
    p,q = 1,0
    (1..a.size).each{|i| p, q = q+p*a[-i], p}
    [q, p]
  end
end

#rationalize_Knuth_Goizueta(x) ⇒ Object

Smallest denominator rationalization based on exercise 39 of cite[S 4.5.3]Knuth. This has been found the most efficient method (except for large tolerances) as implemented in Ruby. Here’s the rationalization procedure based on the exercise by Knuth. We need first to calculate the limits (x-dx, x+dx) of the range where we’ll look for the rational number. If we compute them using floating point and then convert then to fractions this method is always more efficient than the other procedures implemented here, but it may be less accurate. We can achieve perfect accuracy as the other methods by doing the substraction and addition with rationals, but then this method becomes less efficient than the others for a low number of iterations (low precision required).



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
 
# File 'lib/flt/support/rationalizer_extra.rb', line 79

def rationalize_Knuth_Goizueta(x)
  rationalization(x) do |x, dx|
    x = to_r(x)
    dx = to_r(dx)
    xp,xq = num_den(x-dx)
    yp,yq = num_den(x+dx)

    a = []
    fin,odd = false,false
    while !fin && xp!=0 && yp!=0
      odd = !odd
      xp,xq = xq,xp
      ax = xp.div(xq)
      xp -= ax*xq

      yp,yq = yq,yp
      ay = yp.div(yq)
      yp -= ay*yq

      if ax!=ay
        fin = true
        ax,xp,xq = ay,yp,yq if odd
      end
      a << ax # .to_i
    end
    a[-1] += 1 if xp!=0 && a.size>0
    p,q = 1,0
    (1..a.size).each{|i| p,q=q+p*a[-i],p}
    [q, p]
  end
end

#rationalize_Knuth_Goizueta_b(x) ⇒ Object

La siguiente variante realiza una iteración menos si xq<xp y una iteración más si xq>xp.



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
 
# File 'lib/flt/support/rationalizer_extra.rb', line 113

def rationalize_Knuth_Goizueta_b(x)
  rationalization(x) do |x, dx|
    x = to_r(x)
    dx = to_r(dx)
    xq,xp = num_den(x-dx)
    yq,yp = num_den(x+dx)

    a = []
    fin,odd = false,false
    while !fin && xp!=0 && yp!=0
      odd = !odd
      xp,xq = xq,xp
      ax = xp.div(xq)
      xp -= ax*xq

      yp,yq = yq,yp
      ay = yp.div(yq)
      yp -= ay*yq

      if ax!=ay
        fin = true
        ax,xp,xq = ay,yp,yq if odd
      end
      a << ax # .to_i
    end
    a[-1] += 1 if xp!=0 && a.size>0
    p,q = 1,0
    (1..a.size).each{|i| p,q=q+p*a[-i],p}
    [p, q]
  end
end