//
// 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::Div;
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 round(&self) -> i32;
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 })
},
}
}
fn cot(d :i32) -> Self
where Self: Sized + Copy + From<i32> + Div<Output = Self> {
Into::<Self>::into(1) / Self::tan(d)
}
}
// 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)
}
fn round(&self) -> i32 {
let Fractional(n, d) = self;
(n / d) as i32
}
// 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),
_ => generate(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),
_ => generate(d, PRECISION, &f64::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 generate(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),
_ => generate(d, p + 1, f),
}
}
impl Trig for f64 {
fn pi() -> Self {
std::f64::consts::PI
}
fn recip(self) -> Self {
self.recip()
}
fn round(&self) -> i32 {
f64::round(*self) as i32
}
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()
}
}