From 7f9b67443ed79dac2fc20b86fab675821298bfc8 Mon Sep 17 00:00:00 2001 From: Georg Hopp Date: Wed, 4 Dec 2019 17:05:27 +0100 Subject: [PATCH] Add tangens and improve fractional rep --- fractional/src/trigonometry.rs | 133 +++++++++++++++++++-------------- 1 file changed, 77 insertions(+), 56 deletions(-) diff --git a/fractional/src/trigonometry.rs b/fractional/src/trigonometry.rs index cc809a1..65166b8 100644 --- a/fractional/src/trigonometry.rs +++ b/fractional/src/trigonometry.rs @@ -1,10 +1,12 @@ // // Some trigonometic functions with Fractions results. -// Currently only sin and cos are implemented. As I was unable -// to find a really good integral approximation for them I -// implement them as tables which are predefined using the -// floating point function f64::sin and then transformed into -// a fraction of a given PRECISION. +// Currently only sin, cos and tan are implemented. +// As I was unable to find a really good integral approximation for them I +// implement them as a table which is predefined using the floating point +// function f64::sin and then transformed into a fraction of a given +// PRECISION. +// These approximations are quite good and for a few edge cases +// even better than the floating point implementations. // // Georg Hopp // @@ -23,75 +25,94 @@ // You should have received a copy of the GNU General Public License // along with this program. If not, see . // +use std::cmp::Ordering; use crate::{Fractional}; pub const PI :Fractional = Fractional(355, 113); // This is a really close // fractional approximation // for pi -const PRECISION :i64 = 100000; +// Try to keep precision as high as possible while having a denominator +// as small as possible. +// The values are taken by triel and error. +const PRECISION :i64 = 1000000; +const MAX_DENOMINATOR :i64 = 7000; -#[inline] -pub fn rad(d: u32) -> f64 { - use std::f64::consts::PI; - d as f64 * PI / 180.0 -} - -pub fn sin(d: i32) -> Fractional { - // hold sin Fractionals from 0 to 89 ... - lazy_static::lazy_static! { - static ref SINTAB :Vec = - (0..90).map(|x| _sin(x)).collect(); +pub fn sin(d :i32) -> Fractional { + match d { + 0 ..=90 => SINTAB[d as usize], + 91 ..=180 => SINTAB[180 - d as usize], + 181..=270 => -SINTAB[d as usize - 180], + 271..=359 => -SINTAB[360 - d as usize], + _ => sin(d % 360), } +} - // fractional sin from f64 sin. (From 1° to 89°) - fn _sin(d: u32) -> Fractional { - match d { - 0 => Fractional(0, 1), - _ => { - // This is undefined behaviour for very large f64, but our f64 - // is always between 0.0 and 10000.0 which should be fine. - let s = (f64::sin(rad(d)) * PRECISION as f64).round() as i64; - Fractional(s, PRECISION).reduce() - } - } +pub fn cos(d :i32) -> Fractional { + match d { + 0 ..=90 => SINTAB[90 - d as usize], + 91 ..=180 => -SINTAB[90 - (180 - d as usize)], + 181..=270 => -SINTAB[90 - (d as usize - 180)], + 271..=359 => SINTAB[90 - (360 - d as usize)], + _ => cos(d % 360), } +} +pub fn tan(d :i32) -> Fractional { match d { - 90 => Fractional(1, 1), - 180 => SINTAB[0], - 270 => -Fractional(1, 1), - 1..=89 => SINTAB[d as usize], - 91..=179 => SINTAB[180 - d as usize], - 181..=269 => -SINTAB[d as usize - 180], - 271..=359 => -SINTAB[360 - d as usize], - _ => sin(d % 360), + 0 ..=179 => TANTAB[d as usize], + 180..=359 => TANTAB[d as usize - 180], + _ => tan(d % 360), } } -pub fn cos(d: i32) -> Fractional { - lazy_static::lazy_static! { - static ref COSTAB :Vec = - (0..90).map(|x| _cos(x)).collect(); - } +// hold sin Fractionals from 0 to 89 ... +// luckily with a bit of index tweeking this can also be used for +// cosine values. +lazy_static::lazy_static! { + static ref SINTAB :Vec = + (0..=90).map(|x| _sin(x)).collect(); +} - fn _cos(d: u32) -> Fractional { - match d { - 0 => Fractional(1, 1), - _ => { - let s = (f64::cos(rad(d)) * PRECISION as f64).round() as i64; - Fractional(s, PRECISION).reduce() - } - } +// This table exists only because the sin(α) / cos(α) method +// yields very large unreducable denominators in a lot of cases. +lazy_static::lazy_static! { + static ref TANTAB :Vec = + (0..180).map(|x| _tan(x)).collect(); +} + +// fractional sin from f64 sin. (From 0° to 90°) +fn _sin(d: u32) -> Fractional { + match d { + 0 => Fractional(0, 1), + 90 => Fractional(1, 1), + _ => reduce(d, PRECISION, &f64::sin), } +} +// fractional tan from f64 tan. (From 0° to 179°) +fn _tan(d: u32) -> Fractional { match d { - 90 | 270 => Fractional(0, 1), - 180 => -COSTAB[0], - 1..=89 => COSTAB[d as usize], - 91..=179 => -COSTAB[180 - d as usize], - 181..=269 => -COSTAB[d as usize - 180], - 271..=359 => COSTAB[360 - d as usize], - _ => cos(d % 360), + 0 => Fractional(0, 1), + 45 => Fractional(1, 1), + 90 => Fractional(1, 0), // although they are both inf and -inf. + 135 => -Fractional(1, 1), + _ => reduce(d, PRECISION, &f64::tan), + } +} + +// search for a fraction with a denominator less than 10000 that +// provides the minimal precision criteria. +// !! With f = &f64::tan and d close to the inf boundarys of tan +// we get very large numerators because the numerator becomes a +// multiple of the denominator. +fn reduce(d :u32, p :i64, f :&dyn Fn(f64) -> f64) -> Fractional { + // This is undefined behaviour for very large f64, but our f64 + // is always between 0.0 and 10000000.0 which should be fine. + let s = (f(f64::to_radians(d as f64)) * p as f64).round() as i64; + let Fractional(n, dn) = Fractional(s, p).reduce(); + match dn.abs().cmp(&MAX_DENOMINATOR) { + Ordering::Less => Fractional(n, dn), + _ => reduce(d, p + 1, f), } }