diff options
Diffstat (limited to 'lib/pure/rationals.nim')
| -rw-r--r-- | lib/pure/rationals.nim | 82 |
1 files changed, 33 insertions, 49 deletions
diff --git a/lib/pure/rationals.nim b/lib/pure/rationals.nim index c2ba2b1f3..7fb24c26f 100644 --- a/lib/pure/rationals.nim +++ b/lib/pure/rationals.nim @@ -39,47 +39,13 @@ proc toRational*[T:SomeInteger](x: T): Rational[T] = result.num = x result.den = 1 -proc toRationalSub(x: float, n: int): Rational[int] = - var - a = 0'i64 - b, c, d = 1'i64 - result = 0 // 1 # rational 0 - while b <= n and d <= n: - let ac = (a+c) - let bd = (b+d) - # scale by 1000 so not overflow for high precision - let mediant = (ac.float/1000) / (bd.float/1000) - if x == mediant: - if bd <= n: - result.num = ac.int - result.den = bd.int - return result - elif d > b: - result.num = c.int - result.den = d.int - return result - else: - result.num = a.int - result.den = b.int - return result - elif x > mediant: - a = ac - b = bd - else: - c = ac - d = bd - if (b > n): - return initRational(c.int, d.int) - return initRational(a.int, b.int) - -proc toRational*(x: float, n: int = high(int)): Rational[int] = - ## Calculate the best rational numerator and denominator +proc toRational*(x: float, n: int = high(int32)): Rational[int] = + ## Calculates the best rational numerator and denominator ## that approximates to `x`, where the denominator is ## smaller than `n` (default is the largest possible - ## int to give maximum resolution) + ## int to give maximum resolution). ## - ## The algorithm is based on the Farey sequence named - ## after John Farey + ## The algorithm is based on the theory of continued fractions. ## ## .. code-block:: Nim ## import math, rationals @@ -88,13 +54,24 @@ proc toRational*(x: float, n: int = high(int)): Rational[int] = ## let x = toRational(PI, t) ## let newPI = x.num / x.den ## echo x, " ", newPI, " error: ", PI - newPI, " ", t - if x > 1: - result = toRationalSub(1.0/x, n) - swap(result.num, result.den) - elif x == 1.0: - result = 1 // 1 - else: - result = toRationalSub(x, n) + + # David Eppstein / UC Irvine / 8 Aug 1993 + # With corrections from Arno Formella, May 2008 + var + m11, m22 = 1 + m12, m21 = 0 + ai = int(x) + x = x + while m21 * ai + m22 <= n: + swap m12, m11 + swap m22, m21 + m11 = m12 * ai + m11 + m21 = m22 * ai + m21 + if x == float(ai): break # division by zero + x = 1/(x - float(ai)) + if x > float(high(int32)): break # representation failure + ai = int(x) + result = m11 // m21 proc toFloat*[T](x: Rational[T]): float = ## Convert a rational number `x` to a float. @@ -346,7 +323,14 @@ when isMainModule: assert abs(toFloat(y) - 0.4814814814814815) < 1.0e-7 assert toInt(z) == 0 - assert toRational(0.98765432) == 12345679 // 12500000 - assert toRational(0.1, 1000000) == 1 // 10 - assert toRational(0.9, 1000000) == 9 // 10 - #assert toRational(PI) == 80143857 // 25510582 + assert toRational(0.98765432) == 2111111029 // 2137499919 + assert toRational(PI) == 817696623 // 260280919 + assert toRational(0.1) == 1 // 10 + assert toRational(0.9) == 9 // 10 + + assert toRational(0.0) == 0 // 1 + assert toRational(-0.25) == 1 // -4 + assert toRational(3.2) == 16 // 5 + assert toRational(0.33) == 33 // 100 + assert toRational(0.22) == 11 // 50 + assert toRational(10.0) == 10 // 1 |