math.fractions: add approximation.v and tests

2020-05-17 14:30:29 +05:30 · 2020-05-17 14:30:29 +05:30 · b138cadbcb
parent 02fb393747
commit b138cadbcb
3 changed files with 302 additions and 0 deletions
--- a/vlib/math/fractions/approximations.v
+++ 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")
+}
--- a/vlib/math/fractions/approximations_test.v
+++ 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))
+}
--- 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