trigonometry.rs 4.1 KB
//
// Some trigonometic functions with Fractions results.
// 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 <georg@steffers.org>
//
// Copyright © 2019 Georg Hopp
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program.  If not, see <http://www.gnu.org/licenses/>.
//
use std::cmp::Ordering;
use crate::{Fractional};

pub const PI :Fractional = Fractional(355, 113); // This is a really close
                                                 // fractional approximation
                                                 // for pi

// 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;

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),
    }
}

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 {
        0  ..=179 => TANTAB[d as usize],
        180..=359 => TANTAB[d as usize - 180],
        _         => tan(d % 360),
    }
}

// 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<Fractional> =
        (0..=90).map(|x| _sin(x)).collect();
}

// 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<Fractional> =
        (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 {
        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),
    }
}