From 725ece9a4a2e4e50b43f7da16b7bd04f1779766d Mon Sep 17 00:00:00 2001 From: Georg Hopp Date: Fri, 13 Dec 2019 01:24:51 +0100 Subject: [PATCH] Generics for Vector and Transform --- fractional/src/continuous.rs | 88 +++++++++++ fractional/src/fractional.rs | 147 ++----------------- fractional/src/lib.rs | 10 +- fractional/src/main.rs | 44 +++--- fractional/src/transform.rs | 102 ++++++++----- fractional/src/trigonometry.rs | 258 +++++++++++++++++++++++++-------- fractional/src/vector.rs | 50 ++++--- 7 files changed, 424 insertions(+), 275 deletions(-) create mode 100644 fractional/src/continuous.rs diff --git a/fractional/src/continuous.rs b/fractional/src/continuous.rs new file mode 100644 index 0000000..361c84d --- /dev/null +++ b/fractional/src/continuous.rs @@ -0,0 +1,88 @@ +// +// A «continued fraction» is a representation of a fraction as a vector +// of integrals… Irrational fractions will result in infinite most of the +// time repetitive vectors. They can be used to get a resonable approximation +// for sqrt on fractionals. +// +// Georg Hopp +// +// 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 . +// +use crate::Fractional; + +#[derive(Debug)] +pub struct Continuous (Vec); + +impl Continuous { + // calculate a sqrt as continued fraction sequence. Taken from: + // https://en.wikipedia.org/wiki/Methods_of_computing_square_roots# + // Continued_fraction_expansion + pub fn sqrt(x :i64, a0 :i64) -> Self { + fn inner(mut v :Vec, + x :i64, + a0 :i64, + mn :i64, + dn :i64, + an :i64) -> Vec { + let mn_1 = dn * an - mn; + let dn_1 = (x - mn_1 * mn_1) / dn; + let an_1 = (a0 + mn_1) / dn_1; + + v.push(an); + + // The convergence criteria „an_1 == 2 * a0“ is not good for + // very small x thus I decided to break the iteration at constant + // time. Which is the 10 below. + match v.len() { + 10 => v, + _ => inner(v, x, a0, mn_1, dn_1, an_1), + } + } + + Continuous(inner(Vec::new(), x, a0, 0, 1, a0)) + } +} + +impl From<&Fractional> for Continuous { + // general continous fraction form of a fractional... + fn from(x :&Fractional) -> Self { + fn inner(mut v :Vec, f :Fractional) -> Vec { + let Fractional(n, d) = f; + let a = n / d; + let Fractional(_n, _d) = f - a.into(); + + v.push(a); + match _n { + 1 => { v.push(_d); v }, + _ => inner(v, Fractional(_d, _n)), + } + } + + Continuous(inner(Vec::new(), *x)) + } +} + +impl Into for &Continuous { + fn into(self) -> Fractional { + let Continuous(c) = self; + let Fractional(n, d) = c.iter().rev().fold( Fractional(0, 1) + , |acc, x| { + let Fractional(an, ad) = acc + (*x).into(); + Fractional(ad, an) + }); + Fractional(d, n) + } +} diff --git a/fractional/src/fractional.rs b/fractional/src/fractional.rs index 007da4a..92c8922 100644 --- a/fractional/src/fractional.rs +++ b/fractional/src/fractional.rs @@ -31,10 +31,6 @@ use std::num::TryFromIntError; #[derive(Debug, Eq, Clone, Copy)] pub struct Fractional (pub i64, pub i64); -pub type Continuous = Vec; - -pub type Error = &'static str; - #[inline] fn hcf(x :i64, y :i64) -> i64 { match y { @@ -43,36 +39,8 @@ fn hcf(x :i64, y :i64) -> i64 { } } -// calculate a sqrt as continued fraction sequence. Taken from: -// https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Continued_fraction_expansion -fn sqrt_cfrac(x :i64, a0 :i64) -> Continuous { - let v :Continuous = Vec::new(); - - fn inner(mut v :Continuous, - x :i64, - a0 :i64, - mn :i64, - dn :i64, - an :i64) -> Continuous { - let mn_1 = dn * an - mn; - let dn_1 = (x - mn_1 * mn_1) / dn; - let an_1 = (a0 + mn_1) / dn_1; - - v.push(an); - match v.len() { - 10 => v, - _ => inner(v, x, a0, mn_1, dn_1, an_1), - } - - // This convergence criterium is not good for very small x - // thus I decided to break the iteration at constant time. -// match an_1 == 2 * a0 { -// true => v, -// _ => inner(v, x, a0, mn_1, dn_1, an_1), -// } - } - - inner(v, x, a0, 0, 1, a0) +pub fn from_vector(xs: &Vec) -> Vec { + xs.iter().map(|x| Fractional(*x, 1)).collect() } impl Fractional { @@ -98,76 +66,6 @@ impl Fractional { Self(n / hcf(n, d), d / hcf(n, d)) } } - - #[inline] - pub fn numerator(self) -> i64 { - self.0 - } - - #[inline] - pub fn denominator(self) -> i64 { - self.1 - } - - // 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/ - pub fn sqrt(self) -> Result { - // 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 { - (&sqrt_cfrac(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) - } -} - -impl fmt::Display for Fractional { - fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { - write!(f, "({}/{})", self.0, self.1) - } } impl From for Fractional { @@ -176,6 +74,12 @@ impl From for Fractional { } } +impl From for Fractional { + fn from(x: i32) -> Self { + Self(x as i64, 1) + } +} + impl TryFrom for Fractional { type Error = &'static str; @@ -188,10 +92,6 @@ impl TryFrom for Fractional { } } -pub fn from_vector(xs: &Vec) -> Vec { - xs.iter().map(|x| Fractional(*x, 1)).collect() -} - impl TryInto for Fractional { type Error = TryFromIntError; @@ -212,34 +112,9 @@ impl TryInto<(i32, i32)> for Fractional { } } -impl Into for Fractional { - // general continous fraction form of a fractional... - fn into(self) -> Continuous { - let v :Continuous = Vec::new(); - - fn inner(mut v :Continuous, f :Fractional) -> Continuous { - let Fractional(n, d) = f; - let a = n / d; - let Fractional(_n, _d) = f - a.into(); - - v.push(a); - match _n { - 1 => { v.push(_d); v }, - _ => inner(v, Fractional(_d, _n)), - } - } - - inner(v, self) - } -} - -impl From<&Continuous> for Fractional { - fn from(x :&Continuous) -> Self { - let Self(n, d) = x.iter().rev().fold(Fractional(0, 1), |acc, x| { - let Self(an, ad) = acc + (*x).into(); - Self(ad, an) - }); - Self(d, n) +impl fmt::Display for Fractional { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + write!(f, "({}/{})", self.0, self.1) } } diff --git a/fractional/src/lib.rs b/fractional/src/lib.rs index b18c369..49ccdd6 100644 --- a/fractional/src/lib.rs +++ b/fractional/src/lib.rs @@ -20,11 +20,13 @@ // extern crate lazy_static; +pub type Error = &'static str; + +pub mod continuous; pub mod fractional; +pub mod transform; pub mod trigonometry; pub mod vector; -pub mod transform; -use fractional::{Fractional}; -use trigonometry::{sin, cos}; -use vector::{Vector}; +use fractional::Fractional; +use vector::Vector; diff --git a/fractional/src/main.rs b/fractional/src/main.rs index 604734e..93bb2d3 100644 --- a/fractional/src/main.rs +++ b/fractional/src/main.rs @@ -22,8 +22,9 @@ use std::convert::{TryFrom, TryInto, Into}; use std::num::TryFromIntError; use std::f64::consts::PI as FPI; -use fractional::fractional::{Fractional, from_vector, Continuous}; -use fractional::trigonometry::{sin, cos, tan, PI}; +use fractional::continuous::Continuous; +use fractional::fractional::{Fractional, from_vector}; +use fractional::trigonometry::Trig; use fractional::vector::{Vector}; use fractional::transform::{translate, rotate_x, rotate_y, rotate_z, rotate_v}; @@ -51,8 +52,8 @@ fn continuous() { let d = Fractional(45, 16); let e = Fractional(16, 45); - let dc :Continuous = d.into(); - let ec :Continuous = e.into(); + let dc :Continuous = (&d).into(); + let ec :Continuous = (&e).into(); println!("cont frac of d : {} => {:?}", d, dc); println!("cont frac of e : {} => {:?}", e, ec); @@ -82,21 +83,22 @@ fn sqrt() { } fn pi() { - let pir :f64 = PI.try_into().unwrap(); - let pit :(i32, i32) = PI.try_into().unwrap(); - let pi2r :f64 = (PI * PI).try_into().unwrap(); + let pi = Fractional::pi(); + let pir :f64 = pi.try_into().unwrap(); + let pit :(i32, i32) = pi.try_into().unwrap(); + let pi2r :f64 = (pi * pi).try_into().unwrap(); println!(" Rust π : {}" , FPI); - println!(" π : {} {}" , PI, pir); + println!(" π : {} {}" , pi, pir); println!(" π as tuple : {:?}" , pit); println!(" Rust π² : {}" , FPI * FPI); - println!(" π² : {} {}" , PI * PI, pi2r); + println!(" π² : {} {}" , pi * pi, pi2r); } fn _sin() { for d in [ 0, 45, 90, 135, 180, 225, 270, 315 , 9, 17, 31, 73, 89, 123, 213, 312, 876 ].iter() { - let s = sin(*d as i32); + let s = Fractional::sin(*d as i32); let sr :f64 = s.try_into().unwrap(); let _s = f64::sin(*d as f64 * FPI / 180.0); @@ -107,22 +109,22 @@ fn _sin() { fn _tan() { for d in [ 0, 45, 90, 135, 180, 225, 270, 315 , 9, 17, 31, 73, 89, 123, 213, 312, 876 ].iter() { - let s = tan(*d as i32); - let sr :f64 = s.try_into().unwrap(); - let _s = f64::tan(*d as f64 * FPI / 180.0); + let t = Fractional::tan(*d as i32); + let tr :f64 = t.try_into().unwrap(); + let _t = f64::tan(*d as f64 * FPI / 180.0); - println!("{:>14} : {} {} / {}", format!("tan {}", d), s, sr, _s); + println!("{:>14} : {} {} / {}", format!("tan {}", d), t, tr, _t); } } fn _cos() { for d in [ 0, 45, 90, 135, 180, 225, 270, 315 , 9, 17, 31, 73, 89, 123, 213, 312, 876 ].iter() { - let s = cos(*d as i32); - let sr :f64 = s.try_into().unwrap(); - let _s = f64::cos(*d as f64 * FPI / 180.0); + let c = Fractional::cos(*d as i32); + let cr :f64 = c.try_into().unwrap(); + let _c = f64::cos(*d as f64 * FPI / 180.0); - println!("{:>14} : {} {} / {}", format!("cos {}", d), s, sr, _s); + println!("{:>14} : {} {} / {}", format!("cos {}", d), c, cr, _c); } } @@ -152,8 +154,8 @@ fn _vector() { } fn _transform() { - let v = Vector(1.into(), 1.into(), 1.into()); - let v1 = Vector(1.into(), 2.into(), 3.into()); + let v = Vector(Fractional(1,1), Fractional(1,1), Fractional(1,1)); + let v1 = Vector(Fractional(1,1), Fractional(2,1), Fractional(3,1)); let mt = translate(v); println!("{:>14} : {}", "Vector v1", v1); @@ -184,7 +186,7 @@ fn _transform() { } println!(); - let v3 = Vector(1.into(), 0.into(), 1.into()); + let v3 = Vector(Fractional(1,1), Fractional(0,1), Fractional(1,1)); println!("{:>14} : {}", "Vector v3", v3); for d in [ 30, 45, 60, 90, 120, 135, 150, 180 , 210, 225, 240, 270, 300, 315, 330 ].iter() { diff --git a/fractional/src/transform.rs b/fractional/src/transform.rs index 2469f2f..591c383 100644 --- a/fractional/src/transform.rs +++ b/fractional/src/transform.rs @@ -18,15 +18,19 @@ // You should have received a copy of the GNU General Public License // along with this program. If not, see . // -use std::ops::{Mul}; -use crate::{Fractional, cos, sin, Vector}; - -pub struct TMatrix( (Fractional, Fractional, Fractional, Fractional) - , (Fractional, Fractional, Fractional, Fractional) - , (Fractional, Fractional, Fractional, Fractional) - , (Fractional, Fractional, Fractional, Fractional) ); - -pub fn translate(v :Vector) -> TMatrix { +use std::ops::{Add, Sub, Neg, Mul, Div}; +use crate::Vector; +use crate::trigonometry::Trig; + +pub struct TMatrix( (T, T, T, T) + , (T, T, T, T) + , (T, T, T, T) + , (T, T, T, T) ) + where T: Add + Sub + Neg + Mul + Div + Trig + From + Copy; + +pub fn translate(v :Vector) -> TMatrix +where T: Add + Sub + Neg + + Mul + Div + Trig + From + Copy { let Vector(x, y, z) = v; TMatrix( (1.into(), 0.into(), 0.into(), 0.into()) @@ -35,49 +39,71 @@ pub fn translate(v :Vector) -> TMatrix { , ( x, y, z, 1.into()) ) } -pub fn rotate_x(a :i32) -> TMatrix { +pub fn rotate_x(a :i32) -> TMatrix +where T: Add + Sub + Neg + + Mul + Div + Trig + From + Copy { + let sin :T = Trig::sin(a); + let cos :T = Trig::cos(a); + TMatrix( (1.into(), 0.into(), 0.into(), 0.into()) - , (0.into(), cos(a), -sin(a), 0.into()) - , (0.into(), sin(a), cos(a), 0.into()) + , (0.into(), cos , -sin , 0.into()) + , (0.into(), sin , cos , 0.into()) , (0.into(), 0.into(), 0.into(), 1.into()) ) } -pub fn rotate_y(a :i32) -> TMatrix { - TMatrix( ( cos(a), 0.into(), sin(a), 0.into()) +pub fn rotate_y(a :i32) -> TMatrix +where T: Add + Sub + Neg + + Mul + Div + Trig + From + Copy { + let sin :T = Trig::sin(a); + let cos :T = Trig::cos(a); + + TMatrix( (cos , 0.into(), sin , 0.into()) , (0.into(), 1.into(), 0.into(), 0.into()) - , ( -sin(a), 0.into(), cos(a), 0.into()) + , (-sin , 0.into(), cos , 0.into()) , (0.into(), 0.into(), 0.into(), 1.into()) ) } -pub fn rotate_z(a :i32) -> TMatrix { - TMatrix( ( cos(a), -sin(a), 0.into(), 0.into()) - , ( sin(a), cos(a), 0.into(), 0.into()) +pub fn rotate_z(a :i32) -> TMatrix +where T: Add + Sub + Neg + + Mul + Div + Trig + From + Copy { + let sin :T = Trig::sin(a); + let cos :T = Trig::cos(a); + + TMatrix( (cos , -sin , 0.into(), 0.into()) + , (sin , cos , 0.into(), 0.into()) , (0.into(), 0.into(), 1.into(), 0.into()) , (0.into(), 0.into(), 0.into(), 1.into()) ) } -pub fn rotate_v(v :&Vector, a :i32) -> TMatrix { +pub fn rotate_v(v :&Vector, a :i32) -> TMatrix +where T: Add + Sub + Neg + + Mul + Div + Trig + From + Copy { let Vector(x, y, z) = *v; - let zero :Fractional = 0.into(); - let one :Fractional = 1.into(); + let sin :T = Trig::sin(a); + let cos :T = Trig::cos(a); + + let zero :T = 0.into(); + let one :T = 1.into(); - TMatrix( ( (one - cos(a)) * x * x + cos(a) - , (one - cos(a)) * x * y - sin(a) * z - , (one - cos(a)) * x * z + sin(a) * y + TMatrix( ( (one - cos) * x * x + cos + , (one - cos) * x * y - sin * z + , (one - cos) * x * z + sin * y , zero ) - , ( (one - cos(a)) * x * y + sin(a) * z - , (one - cos(a)) * y * y + cos(a) - , (one - cos(a)) * y * z - sin(a) * x + , ( (one - cos) * x * y + sin * z + , (one - cos) * y * y + cos + , (one - cos) * y * z - sin * x , zero ) - , ( (one - cos(a)) * x * z - sin(a) * y - , (one - cos(a)) * y * z + sin(a) * x - , (one - cos(a)) * z * z + cos(a) + , ( (one - cos) * x * z - sin * y + , (one - cos) * y * z + sin * x + , (one - cos) * z * z + cos , zero ) , (0.into(), 0.into(), 0.into(), 1.into()) ) } -pub fn scale(v :Vector) -> TMatrix { +pub fn scale(v :Vector) -> TMatrix +where T: Add + Sub + Neg + + Mul + Div + Trig + From + Copy { let Vector(x, y, z) = v; TMatrix( ( x, 0.into(), 0.into(), 0.into()) @@ -86,8 +112,10 @@ pub fn scale(v :Vector) -> TMatrix { , (0.into(), 0.into(), 0.into(), 1.into()) ) } -impl TMatrix { - pub fn apply(&self, v :&Vector) -> Vector { +impl TMatrix +where T: Add + Sub + Neg + + Mul + Div + Trig + From + Copy { + pub fn apply(&self, v :&Vector) -> Vector { let TMatrix( (a11, a12, a13, a14) , (a21, a22, a23, a24) , (a31, a32, a33, a34) @@ -97,13 +125,15 @@ impl TMatrix { let v = Vector( a11 * x + a21 * y + a31 * z + a41 , a12 * x + a22 * y + a32 * z + a42 , a13 * x + a23 * y + a33 * z + a43 ); - let Fractional(wn, wd) = a14 * x + a24 * y + a34 * z + a44; + let w = a14 * x + a24 * y + a34 * z + a44; - v.mul(&Fractional(wd, wn)) + v.mul(&w.recip()) } } -impl Mul for TMatrix { +impl Mul for TMatrix +where T: Add + Sub + Neg + + Mul + Div + Trig + From + Copy { type Output = Self; // ATTENTION: This is not commutative, nor assoziative. diff --git a/fractional/src/trigonometry.rs b/fractional/src/trigonometry.rs index 65166b8..580bcfd 100644 --- a/fractional/src/trigonometry.rs +++ b/fractional/src/trigonometry.rs @@ -26,90 +26,230 @@ // along with this program. If not, see . // use std::cmp::Ordering; -use crate::{Fractional}; +use std::ops::Neg; +use std::marker::Sized; +use crate::{Fractional, Error}; +use crate::continuous::Continuous; -pub const PI :Fractional = Fractional(355, 113); // This is a really close - // fractional approximation - // for pi +pub trait Trig { + fn pi() -> Self; + fn recip(self) -> Self; + fn sqrt(self) -> Result where Self: Sized; + fn sintab() -> Vec where Self: Sized; + fn tantab() -> Vec where Self: Sized; + + fn sin(d :i32) -> Self + where Self: Sized + Neg + 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(d % 360), + } + } + + fn cos(d :i32) -> Self + where Self: Sized + Neg + 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(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(d % 360), + } + } +} // Try to keep precision as high as possible while having a denominator -// as small as possible. -// The values are taken by triel and error. +// as small as possible. The values are taken by try 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), +// This is a really close fractional approximation for pi. +impl Trig for Fractional { + fn pi() -> Self { + Fractional(355, 113) } -} -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), + fn recip(self) -> Self { + let Fractional(n, d) = self; + Fractional(d, n) } -} -pub fn tan(d :i32) -> Fractional { - match d { - 0 ..=179 => TANTAB[d as usize], - 180..=359 => TANTAB[d as usize - 180], - _ => tan(d % 360), + // 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 { + // 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) } -} -// 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 sintab() -> Vec { + // 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(); + } -// 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 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 { + // 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 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() } } -// 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), +impl Trig for f64 { + fn pi() -> Self { + std::f64::consts::PI + } + + fn recip(self) -> Self { + self.recip() + } + + fn sqrt(self) -> Result { + let x = self.sqrt(); + match x.is_nan() { + true => Err("sqrt on negative undefined"), + false => Ok(x), + } + } + + fn sintab() -> Vec { + lazy_static::lazy_static! { + static ref SINTAB :Vec = + (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 { + // 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 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 10000 that -// provides the minimal precision criteria. +// 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 10000000.0 which should be fine. - let s = (f(f64::to_radians(d as f64)) * p as f64).round() as i64; + // 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), diff --git a/fractional/src/vector.rs b/fractional/src/vector.rs index 1c307be..43bfb54 100644 --- a/fractional/src/vector.rs +++ b/fractional/src/vector.rs @@ -18,29 +18,32 @@ // You should have received a copy of the GNU General Public License // along with this program. If not, see . // -use std::ops::{Add, Sub, Neg, Mul}; +use std::ops::{Add, Sub, Neg, Mul, Div}; use std::fmt; -use crate::{Fractional}; +use crate::trigonometry::Trig; #[derive(Debug, Eq, Clone, Copy)] -pub struct Vector(pub Fractional, pub Fractional, pub Fractional); +pub struct Vector(pub T, pub T, pub T) + where T: Add + Sub + Neg + Mul + Div + Trig + Copy; -impl Vector { - pub fn x(self) -> Fractional { self.0 } - pub fn y(self) -> Fractional { self.1 } - pub fn z(self) -> Fractional { self.2 } +impl Vector +where T: Add + Sub + Neg + + Mul + Div + Trig + Copy { + pub fn x(self) -> T { self.0 } + pub fn y(self) -> T { self.1 } + pub fn z(self) -> T { self.2 } - pub fn abs(self) -> Fractional { + pub fn abs(self) -> T { let Vector(x, y, z) = self; (x * x + y * y + z * z).sqrt().unwrap() } - pub fn mul(self, s :&Fractional) -> Self { + pub fn mul(self, s :&T) -> Self { let Vector(x, y, z) = self; Vector(x * *s, y * *s, z * *s) } - pub fn dot(self, other :Self) -> Fractional { + pub fn dot(self, other :Self) -> T { let Vector(x1, y1, z1) = self; let Vector(x2, y2, z2) = other; @@ -48,24 +51,25 @@ impl Vector { } pub fn norm(self) -> Self { - let Fractional(n, d) = self.abs(); // TODO This can result in 0 or inf Vectors… // Maybe we need to handle zero and inf abs here… - self.mul(&Fractional(d, n)) + self.mul(&self.abs()) } - pub fn distance(self, other :Self) -> Fractional { + pub fn distance(self, other :Self) -> T { (self - other).abs() } } -impl fmt::Display for Vector { +impl fmt::Display for Vector +where T: Add + Sub + Neg + Mul + Div + Trig + fmt::Display + Copy { fn fmt(&self, f :&mut fmt::Formatter<'_>) -> fmt::Result { write!(f, "({}, {}, {})", self.0, self.1, self.2) } } -impl PartialEq for Vector { +impl PartialEq for Vector +where T: Add + Sub + Neg + Mul + Div + Trig + PartialEq + Copy { fn eq(&self, other :&Self) -> bool { let Vector(x1, y1, z1) = self; let Vector(x2, y2, z2) = other; @@ -73,7 +77,9 @@ impl PartialEq for Vector { } } -impl Add for Vector { +impl Add for Vector +where T: Add + Sub + Neg + + Mul + Div + Trig + Copy { type Output = Self; fn add(self, other :Self) -> Self { @@ -83,7 +89,9 @@ impl Add for Vector { } } -impl Sub for Vector { +impl Sub for Vector +where T: Add + Sub + Neg + + Mul + Div + Trig + Copy { type Output = Self; fn sub(self, other :Self) -> Self { @@ -91,7 +99,9 @@ impl Sub for Vector { } } -impl Neg for Vector { +impl Neg for Vector +where T: Add + Sub + Neg + + Mul + Div + Trig + Copy { type Output = Self; fn neg(self) -> Self { @@ -100,7 +110,9 @@ impl Neg for Vector { } } -impl Mul for Vector { +impl Mul for Vector +where T: Add + Sub + Neg + + Mul + Div + Trig + Copy { type Output = Self; fn mul(self, other :Self) -> Self {