diff options
Diffstat (limited to 'v_windows/v/old/vlib/math/fractions/approximations.v')
-rw-r--r-- | v_windows/v/old/vlib/math/fractions/approximations.v | 119 |
1 files changed, 119 insertions, 0 deletions
diff --git a/v_windows/v/old/vlib/math/fractions/approximations.v b/v_windows/v/old/vlib/math/fractions/approximations.v new file mode 100644 index 0000000..dd4f855 --- /dev/null +++ b/v_windows/v/old/vlib/math/fractions/approximations.v @@ -0,0 +1,119 @@ +// Copyright (c) 2019-2021 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 := i64(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, fractions.default_eps) +} + +// approximate_with_eps returns a Fraction +pub fn approximate_with_eps(val f64, eps f64) Fraction { + if val == 0.0 { + return fractions.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 - f64(whole) + // Quick exit for integers + if frac == 0.0 { + return fraction(whole, 1) + } + mut d := []i64{} + mut partial := fractions.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 .. fractions.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 -= f64(den) + } + panic("Couldn't converge. Please create an issue on https://github.com/vlang/v") +} |