From b138cadbcb21f99f0f7c48c891fa41d3bf6b309c Mon Sep 17 00:00:00 2001 From: Hungry Blue Dev <46458389+hungrybluedev@users.noreply.github.com> Date: Sun, 17 May 2020 14:30:29 +0530 Subject: [PATCH] math.fractions: add approximation.v and tests --- vlib/math/fractions/approximations.v | 119 ++++++++++++++ vlib/math/fractions/approximations_test.v | 180 ++++++++++++++++++++++ vlib/math/fractions/fraction_test.v | 3 + 3 files changed, 302 insertions(+) create mode 100644 vlib/math/fractions/approximations.v create mode 100644 vlib/math/fractions/approximations_test.v diff --git a/vlib/math/fractions/approximations.v b/vlib/math/fractions/approximations.v new file mode 100644 index 0000000000..bd325774f0 --- /dev/null +++ b/vlib/math/fractions/approximations.v @@ -0,0 +1,119 @@ +// Copyright (c) 2019-2020 Alexander Medvednikov. All rights reserved. +// Use of this source code is governed by an MIT license +// that can be found in the LICENSE file. +module fractions + +import math + +const ( + default_eps = 1.0e-4 + max_iterations = 50 + zero = fraction(0, 1) +) + +// ------------------------------------------------------------------------ +// Unwrapped evaluation methods for fast evaluation of continued fractions. +// ------------------------------------------------------------------------ +// We need these functions because the evaluation of continued fractions +// always has to be done from the end. Also, the numerator-denominator pairs +// are generated from front to end. This means building a result from a +// previous one isn't possible. So we need unrolled versions to ensure that +// we don't take too much of a performance penalty by calling eval_cf +// several times. +// ------------------------------------------------------------------------ +// eval_1 returns the result of evaluating a continued fraction series of length 1 +fn eval_1(whole i64, d []i64) Fraction { + return fraction(whole * d[0] + 1, d[0]) +} + +// eval_2 returns the result of evaluating a continued fraction series of length 2 +fn eval_2(whole i64, d []i64) Fraction { + den := d[0] * d[1] + 1 + return fraction(whole * den + d[1], den) +} + +// eval_3 returns the result of evaluating a continued fraction series of length 3 +fn eval_3(whole i64, d []i64) Fraction { + d1d2_plus_n2 := d[1] * d[2] + 1 + den := d[0] * d1d2_plus_n2 + d[2] + return fraction(whole * den + d1d2_plus_n2, den) +} + +// eval_cf evaluates a continued fraction series and returns a Fraction. +fn eval_cf(whole i64, den []i64) Fraction { + count := den.len + // Offload some small-scale calculations + // to dedicated functions + match count { + 1 { + return eval_1(whole, den) + } + 2 { + return eval_2(whole, den) + } + 3 { + return eval_3(whole, den) + } + else { + last := count - 1 + mut n := 1 + mut d := den[last] + // The calculations are done from back to front + for index := count - 2; index >= 0; index-- { + t := d + d = den[index] * d + n + n = t + } + return fraction(d * whole + n, d) + } + } +} + +// approximate returns a Fraction that approcimates the given value to +// within the default epsilon value (1.0e-4). This means the result will +// be accurate to 3 places after the decimal. +pub fn approximate(val f64) Fraction { + return approximate_with_eps(val, default_eps) +} + +// approximate_with_eps returns a Fraction +pub fn approximate_with_eps(val, eps f64) Fraction { + if val == 0.0 { + return zero + } + if eps < 0.0 { + panic('Epsilon value cannot be negative.') + } + if math.fabs(val) > math.max_i64 { + panic('Value out of range.') + } + // The integer part is separated first. Then we process the fractional + // part to generate numerators and denominators in tandem. + whole := i64(val) + mut frac := val - whole + // Quick exit for integers + if frac == 0.0 { + return fraction(whole, 1) + } + mut d := []i64{} + mut partial := zero + // We must complete the approximation within the maximum number of + // itertations allowed. If we can't panic. + // Empirically tested: the hardest constant to approximate is the + // golden ratio (math.phi) and for f64s, it only needs 38 iterations. + for _ in 0 .. max_iterations { + // We calculate the reciprocal. That's why the numerator is + // always 1. + frac = 1.0 / frac + den := i64(frac) + d << den + // eval_cf is called often so it needs to be performant + partial = eval_cf(whole, d) + // Check if we're done + if math.fabs(val - partial.f64()) < eps { + return partial + } + frac -= den + } + panic("Couldn\'t converge. Please create an issue on https://github.com/vlang/v") +} diff --git a/vlib/math/fractions/approximations_test.v b/vlib/math/fractions/approximations_test.v new file mode 100644 index 0000000000..af1406071b --- /dev/null +++ b/vlib/math/fractions/approximations_test.v @@ -0,0 +1,180 @@ +// Copyright (c) 2019-2020 Alexander Medvednikov. All rights reserved. +// Use of this source code is governed by an MIT license +// that can be found in the LICENSE file. +import fractions +import math + +fn test_half() { + float_val := 0.5 + fract_val := fractions.approximate(float_val) + assert fract_val.equals(fractions.fraction(1, 2)) +} + +fn test_third() { + float_val := 1.0 / 3.0 + fract_val := fractions.approximate(float_val) + assert fract_val.equals(fractions.fraction(1, 3)) +} + +fn test_minus_one_twelfth() { + float_val := -1.0 / 12.0 + fract_val := fractions.approximate(float_val) + assert fract_val.equals(fractions.fraction(-1, 12)) +} + +fn test_zero() { + float_val := 0.0 + println('Pre') + fract_val := fractions.approximate(float_val) + println('Post') + assert fract_val.equals(fractions.fraction(0, 1)) +} + +fn test_minus_one() { + float_val := -1.0 + fract_val := fractions.approximate(float_val) + assert fract_val.equals(fractions.fraction(-1, 1)) +} + +fn test_thirty_three() { + float_val := 33.0 + fract_val := fractions.approximate(float_val) + assert fract_val.equals(fractions.fraction(33, 1)) +} + +fn test_millionth() { + float_val := 1.0 / 1000000.0 + fract_val := fractions.approximate(float_val) + assert fract_val.equals(fractions.fraction(1, 1000000)) +} + +fn test_minus_27_by_57() { + float_val := -27.0 / 57.0 + fract_val := fractions.approximate(float_val) + assert fract_val.equals(fractions.fraction(-27, 57)) +} + +fn test_29_by_104() { + float_val := 29.0 / 104.0 + fract_val := fractions.approximate(float_val) + assert fract_val.equals(fractions.fraction(29, 104)) +} + +fn test_140710_232() { + float_val := 140710.232 + fract_val := fractions.approximate(float_val) + // Approximation will match perfectly for upto 3 places after the decimal + // The result will be within default_eps of original value + assert fract_val.f64() == float_val +} + +fn test_pi_1_digit() { + assert fractions.approximate_with_eps(math.pi, 5.0e-2).equals(fractions.fraction(22, 7)) +} + +fn test_pi_2_digits() { + assert fractions.approximate_with_eps(math.pi, 5.0e-3).equals(fractions.fraction(22, 7)) +} + +fn test_pi_3_digits() { + assert fractions.approximate_with_eps(math.pi, 5.0e-4).equals(fractions.fraction(333, 106)) +} + +fn test_pi_4_digits() { + assert fractions.approximate_with_eps(math.pi, 5.0e-5).equals(fractions.fraction(355, 113)) +} + +fn test_pi_5_digits() { + assert fractions.approximate_with_eps(math.pi, 5.0e-6).equals(fractions.fraction(355, 113)) +} + +fn test_pi_6_digits() { + assert fractions.approximate_with_eps(math.pi, 5.0e-7).equals(fractions.fraction(355, 113)) +} + +fn test_pi_7_digits() { + assert fractions.approximate_with_eps(math.pi, 5.0e-8).equals(fractions.fraction(103993, + 33102)) +} + +fn test_pi_8_digits() { + assert fractions.approximate_with_eps(math.pi, 5.0e-9).equals(fractions.fraction(103993, + 33102)) +} + +fn test_pi_9_digits() { + assert fractions.approximate_with_eps(math.pi, 5.0e-10).equals(fractions.fraction(104348, + 33215)) +} + +fn test_pi_10_digits() { + assert fractions.approximate_with_eps(math.pi, 5.0e-11).equals(fractions.fraction(312689, + 99532)) +} + +fn test_pi_11_digits() { + assert fractions.approximate_with_eps(math.pi, 5.0e-12).equals(fractions.fraction(1146408, + 364913)) +} + +fn test_pi_12_digits() { + assert fractions.approximate_with_eps(math.pi, 5.0e-13).equals(fractions.fraction(4272943, + 1360120)) +} + +fn test_phi_1_digit() { + assert fractions.approximate_with_eps(math.phi, 5.0e-2).equals(fractions.fraction(5, 3)) +} + +fn test_phi_2_digits() { + assert fractions.approximate_with_eps(math.phi, 5.0e-3).equals(fractions.fraction(21, 13)) +} + +fn test_phi_3_digits() { + assert fractions.approximate_with_eps(math.phi, 5.0e-4).equals(fractions.fraction(55, 34)) +} + +fn test_phi_4_digits() { + assert fractions.approximate_with_eps(math.phi, 5.0e-5).equals(fractions.fraction(233, + 144)) +} + +fn test_phi_5_digits() { + assert fractions.approximate_with_eps(math.phi, 5.0e-6).equals(fractions.fraction(610, + 377)) +} + +fn test_phi_6_digits() { + assert fractions.approximate_with_eps(math.phi, 5.0e-7).equals(fractions.fraction(1597, + 987)) +} + +fn test_phi_7_digits() { + assert fractions.approximate_with_eps(math.phi, 5.0e-8).equals(fractions.fraction(6765, + 4181)) +} + +fn test_phi_8_digits() { + assert fractions.approximate_with_eps(math.phi, 5.0e-9).equals(fractions.fraction(17711, + 10946)) +} + +fn test_phi_9_digits() { + assert fractions.approximate_with_eps(math.phi, 5.0e-10).equals(fractions.fraction(75025, + 46368)) +} + +fn test_phi_10_digits() { + assert fractions.approximate_with_eps(math.phi, 5.0e-11).equals(fractions.fraction(196418, + 121393)) +} + +fn test_phi_11_digits() { + assert fractions.approximate_with_eps(math.phi, 5.0e-12).equals(fractions.fraction(514229, + 317811)) +} + +fn test_phi_12_digits() { + assert fractions.approximate_with_eps(math.phi, 5.0e-13).equals(fractions.fraction(2178309, + 1346269)) +} \ No newline at end of file diff --git a/vlib/math/fractions/fraction_test.v b/vlib/math/fractions/fraction_test.v index 2ff34b62bf..afaeb2c7ae 100644 --- a/vlib/math/fractions/fraction_test.v +++ b/vlib/math/fractions/fraction_test.v @@ -1,3 +1,6 @@ +// Copyright (c) 2019-2020 Alexander Medvednikov. All rights reserved. +// Use of this source code is governed by an MIT license +// that can be found in the LICENSE file. import math.fractions // (Old) results are verified using https://www.calculatorsoup.com/calculators/math/fractions.php