easel3d/fractional/src/trigonometry.rs

//
// 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 std::ops::Neg;
use std::marker::Sized;
use crate::{Fractional, Error};
use crate::continuous::Continuous;

pub trait Trig {
    fn pi()        -> Self;
    fn recip(self) -> Self;
    fn sqrt(self)  -> Result<Self, Error> where Self: Sized;
    fn sintab()    -> Vec<Self> where Self: Sized;
    fn tantab()    -> Vec<Self> where Self: Sized;

    fn sin(d :i32) -> Self
        where Self: Sized + Neg<Output = Self> + Copy {
        match d {
            0  ..=90  => Self::sintab()[d as usize],
            91 ..=180 => Self::sintab()[180 - d as usize],
            181..=270 => -Self::sintab()[d as usize - 180],
            271..=359 => -Self::sintab()[360 - d as usize],
            _         => {
                Self::sin(if d < 0 { d % 360 + 360 } else { d % 360 })
            },
        }
    }

    fn cos(d :i32) -> Self
        where Self: Sized + Neg<Output = Self> + Copy {
        match d {
            0  ..=90  => Self::sintab()[90 - d as usize],
            91 ..=180 => -Self::sintab()[90 - (180 - d as usize)],
            181..=270 => -Self::sintab()[90 - (d as usize - 180)],
            271..=359 => Self::sintab()[90 - (360 - d as usize)],
            _         => {
                Self::cos(if d < 0 { d % 360 + 360 } else { d % 360 })
            },
        }
    }

    fn tan(d :i32) -> Self where Self: Sized + Copy {
        match d {
            0  ..=179 => Self::tantab()[d as usize],
            180..=359 => Self::tantab()[d as usize - 180],
            _         => {
                Self::tan(if d < 0 { d % 360 + 360 } else { d % 360 })
            },
        }
    }
}

// Try to keep precision as high as possible while having a denominator
// as small as possible. The values are taken by try and error.
const PRECISION       :i64 = 1000000;
const MAX_DENOMINATOR :i64 = 7000;

// This is a really close fractional approximation for pi.
impl Trig for Fractional {
    fn pi() -> Self {
        Fractional(355, 113)
    }

    fn recip(self) -> Self {
        let Fractional(n, d) = self;
        Fractional(d, n)
    }

    // This is a really bad approximation of sqrt for a fractional...
    // for (9/3) it will result 3 which if way to far from the truth,
    // which is ~1.7320508075
    // BUT we can use this value as starting guess for creating a 
    // continous fraction for the sqrt... and create a much better
    // fractional representation of the sqrt.
    // So, if inner converges, but is not a perfect square (does not
    // end up in an Ordering::Equal - which is the l > h case)
    // we use the l - 1 as starting guess for sqrt_cfrac.
    // taken from:
    // https://www.geeksforgeeks.org/square-root-of-an-integer/
    fn sqrt(self) -> Result<Self, Error> {
        // find the sqrt of x in O(log x/2).
        // This stops if a perfect sqare was found. Else it passes
        // the found value as starting guess to the continous fraction
        // sqrt function.
        fn floor_sqrt(x :i64) -> Fractional {
            fn inner(l :i64, h :i64, x :i64) -> Fractional {
                if l > h {
                    (&Continuous::sqrt(x, l - 1)).into()
                } else {
                    let m = (l + h) / 2;
                    match x.cmp(&(m * m)) {
                        Ordering::Equal   => m.into(),
                        Ordering::Less    => inner(l, m - 1, x),
                        Ordering::Greater => inner(m + 1, h, x),
                    }
                }
            }

            match x {
                0 => 0.into(),
                1 => 1.into(),
                _ => inner(1, x / 2, x),
            }
        }

        let Fractional(n, d) = self;

        let n = match n.cmp(&0) {
            Ordering::Equal   => 0.into(),
            Ordering::Less    => return Err("sqrt on negative undefined"),
            Ordering::Greater => floor_sqrt(n),
        };

        let d = match d.cmp(&0) {
            Ordering::Equal   => 0.into(),
            Ordering::Less    => return Err("sqrt on negative undefined"),
            Ordering::Greater => floor_sqrt(d),
        };

        Ok(n / d)
    }

    fn sintab() -> Vec<Self> {
        // 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();
        }

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

        SINTAB.to_vec()
    }

    fn tantab() -> Vec<Self> {
        // 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 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),
            }
        }

        TANTAB.to_vec()
    }
}

impl Trig for f64 {
    fn pi() -> Self {
        std::f64::consts::PI
    }

    fn recip(self) -> Self {
        self.recip()
    }

    fn sqrt(self) -> Result<Self, Error> {
        let x = self.sqrt();
        match x.is_nan() {
            true  => Err("sqrt on negative undefined"),
            false => Ok(x),
        }
    }

    fn sintab() -> Vec<Self> {
        lazy_static::lazy_static! {
            static ref SINTAB :Vec<f64> =
                (0..=90).map(|x| _sin(x)).collect();
        }

        // f64 sin. (From 0° to 90°)
        fn _sin(d: u32) -> f64 {
            match d {
                0  => 0.0,
                90 => 1.0,
                _  => (d as f64).to_radians().sin(),
            }
        }

        SINTAB.to_vec()
    }

    fn tantab() -> Vec<Self> {
        // 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<f64> =
                (0..180).map(|x| _tan(x)).collect();
        }

        // fractional tan from f64 tan. (From 0° to 179°)
        fn _tan(d: u32) -> f64 {
            match d {
                0   => 0.0,
                45  => 1.0,
                90  => std::f64::INFINITY,
                135 => -1.0,
                _   => (d as f64).to_radians().tan(),
            }
        }

        TANTAB.to_vec()
    }
}

// search for a fraction with a denominator less than MAX_DENOMINATOR 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 1000000.0 which should be fine.
    let s = (f((d as f64).to_radians()) * 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),
    }
}