lib/nio/rtnlzr.rb in nio-0.2.3 vs lib/nio/rtnlzr.rb in nio-0.2.4
- old
+ new
@@ -1,13 +1,8 @@
# Rationalization of floating point numbers.
#--
-# Copyright (C) 2003-2005, Javier Goizueta <javier@goizueta.info>
-#
-# This program is free software; you can redistribute it and/or
-# modify it under the terms of the GNU General Public License
-# as published by the Free Software Foundation; either version 2
-# of the License, or (at your option) any later version.
+
#++
#
# Author:: Javier Goizueta (mailto:javier@goizueta.info)
# Copyright:: Copyright (c) 2002-2004 Javier Goizueta & Joe Horn
# License:: Distributes under the GPL license
@@ -25,34 +20,36 @@
# * BigDecimal#nio_r
require 'nio/tools'
-require 'nio/flttol'
+require 'flt/tolerance'
require 'rational'
require 'bigdecimal'
+require 'flt'
+
class Float
# Conversion to Rational preserving the exact value of the number.
def nio_xr
return Rational(self.to_i,1) if self.modulo(1)==0
if !self.finite?
return Rational(0,0) if self.nan?
return self<0 ? Rational(-1,0) : Rational(1,0)
end
-
+
f,e = Math.frexp(self)
-
+
if e < Float::MIN_EXP
bits = e+Float::MANT_DIG-Float::MIN_EXP
else
- bits = [Float::MANT_DIG,e].max
+ bits = [Float::MANT_DIG,e].max
#return Rational(self.to_i,1) if bits<e
- end
+ end
p = Math.ldexp(f,bits)
e = bits - e
if e<Float::MAX_EXP
q = Math.ldexp(1,e)
else
@@ -76,10 +73,17 @@
q = b**(e)
return Rational(p,q)
end
end
+class Flt::Num
+
+ def nio_xr
+ to_r
+ end
+end
+
class Integer
def nio_xr
return Rational(self,1)
end
@@ -88,68 +92,67 @@
class Rational
def nio_xr
return self
end
-
+
# helper method to return both the numerator and denominator
def nio_num_den
return [numerator,denominator]
end
end
class Float
- # Conversion to Rational. The optional argument must be one of:
- # - a Nio::Tolerance that defines the admisible tolerance;
- # in that case, the smallest denominator rational within the
- # tolerance will be found (which may take a long time for
- # small tolerances.)
- # - an integer that defines a maximum value for the denominator.
- # in which case, the best approximation with that maximum
- # denominator will be returned.
- def nio_r(tol = Nio::Tolerance.big_epsilon)
+
+ def nio_r(tol = Flt.Tolerance(:big_epsilon))
case tol
when Integer
Rational(*Nio::Rtnlzr.max_denominator(self,tol,Float))
else
- Rational(*Nio::Rtnlzr.new(Nio::Tol(tol)).rationalize(self))
+ Rational(*Nio::Rtnlzr.new(tol).rationalize(self))
end
end
end
class BigDecimal
- # Conversion to Rational. The optional argument must be one of:
- # - a Nio::BigTolerance that defines the admisible tolerance;
- # in that case, the smallest denominator rational within the
- # tolerance will be found (which may take a long time for
- # small tolerances.)
- # - an integer that defines a maximum value for the denominator.
- # in which case, the best approximation with that maximum
- # denominator will be returned.
+
def nio_r(tol = nil)
- tol ||= BigTolerance.decimals([precs[0],Float::DIG].max,:sig)
+ tol ||= Flt.Tolerance([precs[0],Float::DIG].max,:sig_decimals)
case tol
when Integer
Rational(*Nio::Rtnlzr.max_denominator(self,tol,BigDecimal))
else
- Rational(*Nio::Rtnlzr.new(Nio::BigTol(tol)).rationalize(self))
+ Rational(*Nio::Rtnlzr.new(tol).rationalize(self))
end
end
end
+class Flt::Num
+
+ def nio_r(tol = nil)
+ tol ||= Flt.Tolerance(Rational(1,2),:ulps)
+ case tol
+ when Integer
+ Rational(*Nio::Rtnlzr.max_denominator(self,tol,num_class))
+ else
+ Rational(*Nio::Rtnlzr.new(tol).rationalize(self))
+ end
+ end
+end
+
module Nio
- # This class provides conversion of fractions
+ # This class provides conversion of fractions
# (as approximate floating point numbers)
# to rational numbers.
class Rtnlzr
include StateEquivalent
# Create Rationalizator with given tolerance.
- def initialize(tol=Tolerance.new)
+ def initialize(tol=Flt.Tolerance(:epsilon))
@tol = tol
end
# Rationalization method that finds the fraction with
# smallest denominator fraction within the tolerance distance
@@ -157,12 +160,12 @@
#
# It uses the algorithm which has been found most efficient, rationalize_Knuth.
def rationalize(x)
rationalize_Knuth(x)
end
-
- # This algorithm is derived from exercise 39 of 4.5.3 in
+
+ # This algorithm is derived from exercise 39 of 4.5.3 in
# "The Art of Computer Programming", by Donald E. Knuth.
def rationalize_Knuth(x)
num_tol = @tol.kind_of?(Numeric)
@@ -173,18 +176,18 @@
negans=false
if x<0
negans = true
x = -x
end
- dx = num_tol ? @tol : @tol.get_value(x)
+ dx = num_tol ? @tol : @tol.value(x)
+
-
x = x.nio_xr
dx = dx.nio_xr
xp,xq = (x-dx).nio_num_den
yp,yq = (x+dx).nio_num_den
-
+
a = []
fin,odd = false,false
while !fin && xp!=0 && yp!=0
odd = !odd
xp,xq = xq,xp
@@ -203,17 +206,17 @@
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}
num,den = q,p
-
+
num = -num if negans
end
return num,den
-
+
end
# This is algorithm PDQ2 by Joe Horn.
def rationalize_Horn(x)
@@ -225,13 +228,13 @@
negans=false
if x<0
negans = true
x = -x
end
- dx = num_tol ? @tol : @tol.get_value(x)
+ dx = num_tol ? @tol : @tol.value(x)
+
-
z,t = x,dx # renaming
a,b = t.nio_xr.nio_num_den
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
@@ -255,23 +258,23 @@
if b*(n0*y-d0*x).abs <= a*d0*y
lo = mid
else
hi = mid
end
- end until hi-lo <= 1
+ end until hi-lo <= 1
x = cn - pn*lo
- y = cd - pd*lo
- end
-
+ y = cd - pd*lo
+ end
+
num,den = x,y # renaming
-
+
num = -num if negans
end
return num,den
-
+
end
# This is from a RPL program by Tony Hutchins (PDR6).
def rationalize_HornHutchins(x)
@@ -283,13 +286,13 @@
negans=false
if x<0
negans = true
x = -x
end
- dx = num_tol ? @tol : @tol.get_value(x)
+ dx = num_tol ? @tol : @tol.value(x)
+
-
z,t = x,dx # renaming
a,b = t.nio_xr.nio_num_den
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
@@ -313,23 +316,23 @@
if b*(n0*y-d0*x).abs <= a*d0*y
lo = mid
else
hi = mid
end
- end until hi-lo <= 1
+ end until hi-lo <= 1
x = cn - pn*lo
- y = cd - pd*lo
- end
-
+ y = cd - pd*lo
+ end
+
num,den = x,y # renaming
-
+
num = -num if negans
end
return num,den
-
+
end
end
# Best fraction given maximum denominator
# Algorithm Copyright (c) 1991 by Joseph K. Horn.
@@ -338,38 +341,41 @@
# arithmetic which limits the magnitude and precision of the results, specially
# using Float values.
def Rtnlzr.max_denominator(f, max_den=1000000000, num_class=nil)
return nil if max_den<1
num_class ||= f.class
+ context = num_class.context
return mth.ip(f),1 if mth.fp(f)==0
- one = Nio.numeric_cast(1, num_class)
-
- sign = f<0
- f = -f if sign
+ cast = lambda{|x| context.Num(x)}
- a,b,c = 0,1,f
- while b<max_den and c!=0
- cc = one/c
- a,b,c = b, mth.ip(cc)*b+a, mth.fp(cc)
- end
+ one = cast.call(1)
-
- if b>max_den
- b -= a*mth.ceil((b-max_den)/Float(a))
- end
-
+ sign = f<0
+ f = -f if sign
- f1,f2 = [a,b].collect{|x| mth.abs(mth.rnd(x*f)/Nio.numeric_cast(x, num_class)-f)}
+ a,b,c = 0,1,f
+ while b<max_den and c!=0
+ cc = one/c
+ a,b,c = b, mth.ip(cc)*b+a, mth.fp(cc)
+ end
- a = f1>f2 ? b : a
- num,den = mth.rnd(a*f).to_i,a
- den = 1 if mth.abs(den)<1
-
- num = -num if sign
+ if b>max_den
+ b -= a*mth.ceil(cast.call(b-max_den)/a)
+ end
+
+ f1,f2 = [a,b].collect{|x| mth.abs(cast.call(mth.rnd(x*f))/x-f)}
+
+ a = f1>f2 ? b : a
+
+ num,den = mth.rnd(a*f).to_i,a
+ den = 1 if mth.abs(den)<1
+
+ num = -num if sign
+
return num,den
end
class Rtnlzr
private
@@ -379,25 +385,26 @@
# y =x.modulo(1); return x<0 ? -y : y;
x-ip(x)
end
def self.ip(x)
- # x.to_i.to_f
+ # Note that ceil, floor return an Integer for Float and Flt::Num, but not for BigDecimal
(x<0 ? x.ceil : x.floor).to_i
end
def self.rnd(x)
- #x.round.to_i
- x.round
+ # Note that round returns an Integer for Float and Flt::Num, but not for BigDecimal
+ x.round.to_i
end
def self.abs(x)
x.abs
end
def self.ceil(x)
+ # Note that ceil returns an Integer for Float and Flt::Num, but not for BigDecimal
x.ceil.to_i
- end
+ end
end
def self.mth; Mth; end
end
module_function