Commit 725ece9a4a2e4e50b43f7da16b7bd04f1779766d

Authored by Georg Hopp
1 parent 41d4d98b

Generics for Vector and Transform

  1 +//
  2 +// A «continued fraction» is a representation of a fraction as a vector
  3 +// of integrals… Irrational fractions will result in infinite most of the
  4 +// time repetitive vectors. They can be used to get a resonable approximation
  5 +// for sqrt on fractionals.
  6 +//
  7 +// Georg Hopp <georg@steffers.org>
  8 +//
  9 +// Copyright © 2019 Georg Hopp
  10 +//
  11 +// This program is free software: you can redistribute it and/or modify
  12 +// it under the terms of the GNU General Public License as published by
  13 +// the Free Software Foundation, either version 3 of the License, or
  14 +// (at your option) any later version.
  15 +//
  16 +// This program is distributed in the hope that it will be useful,
  17 +// but WITHOUT ANY WARRANTY; without even the implied warranty of
  18 +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  19 +// GNU General Public License for more details.
  20 +//
  21 +// You should have received a copy of the GNU General Public License
  22 +// along with this program. If not, see <http://www.gnu.org/licenses/>.
  23 +//
  24 +use crate::Fractional;
  25 +
  26 +#[derive(Debug)]
  27 +pub struct Continuous (Vec<i64>);
  28 +
  29 +impl Continuous {
  30 + // calculate a sqrt as continued fraction sequence. Taken from:
  31 + // https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#
  32 + // Continued_fraction_expansion
  33 + pub fn sqrt(x :i64, a0 :i64) -> Self {
  34 + fn inner(mut v :Vec<i64>,
  35 + x :i64,
  36 + a0 :i64,
  37 + mn :i64,
  38 + dn :i64,
  39 + an :i64) -> Vec<i64> {
  40 + let mn_1 = dn * an - mn;
  41 + let dn_1 = (x - mn_1 * mn_1) / dn;
  42 + let an_1 = (a0 + mn_1) / dn_1;
  43 +
  44 + v.push(an);
  45 +
  46 + // The convergence criteria „an_1 == 2 * a0“ is not good for
  47 + // very small x thus I decided to break the iteration at constant
  48 + // time. Which is the 10 below.
  49 + match v.len() {
  50 + 10 => v,
  51 + _ => inner(v, x, a0, mn_1, dn_1, an_1),
  52 + }
  53 + }
  54 +
  55 + Continuous(inner(Vec::new(), x, a0, 0, 1, a0))
  56 + }
  57 +}
  58 +
  59 +impl From<&Fractional> for Continuous {
  60 + // general continous fraction form of a fractional...
  61 + fn from(x :&Fractional) -> Self {
  62 + fn inner(mut v :Vec<i64>, f :Fractional) -> Vec<i64> {
  63 + let Fractional(n, d) = f;
  64 + let a = n / d;
  65 + let Fractional(_n, _d) = f - a.into();
  66 +
  67 + v.push(a);
  68 + match _n {
  69 + 1 => { v.push(_d); v },
  70 + _ => inner(v, Fractional(_d, _n)),
  71 + }
  72 + }
  73 +
  74 + Continuous(inner(Vec::new(), *x))
  75 + }
  76 +}
  77 +
  78 +impl Into<Fractional> for &Continuous {
  79 + fn into(self) -> Fractional {
  80 + let Continuous(c) = self;
  81 + let Fractional(n, d) = c.iter().rev().fold( Fractional(0, 1)
  82 + , |acc, x| {
  83 + let Fractional(an, ad) = acc + (*x).into();
  84 + Fractional(ad, an)
  85 + });
  86 + Fractional(d, n)
  87 + }
  88 +}
... ...
... ... @@ -31,10 +31,6 @@ use std::num::TryFromIntError;
31 31 #[derive(Debug, Eq, Clone, Copy)]
32 32 pub struct Fractional (pub i64, pub i64);
33 33
34   -pub type Continuous = Vec<i64>;
35   -
36   -pub type Error = &'static str;
37   -
38 34 #[inline]
39 35 fn hcf(x :i64, y :i64) -> i64 {
40 36 match y {
... ... @@ -43,36 +39,8 @@ fn hcf(x :i64, y :i64) -> i64 {
43 39 }
44 40 }
45 41
46   -// calculate a sqrt as continued fraction sequence. Taken from:
47   -// https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Continued_fraction_expansion
48   -fn sqrt_cfrac(x :i64, a0 :i64) -> Continuous {
49   - let v :Continuous = Vec::new();
50   -
51   - fn inner(mut v :Continuous,
52   - x :i64,
53   - a0 :i64,
54   - mn :i64,
55   - dn :i64,
56   - an :i64) -> Continuous {
57   - let mn_1 = dn * an - mn;
58   - let dn_1 = (x - mn_1 * mn_1) / dn;
59   - let an_1 = (a0 + mn_1) / dn_1;
60   -
61   - v.push(an);
62   - match v.len() {
63   - 10 => v,
64   - _ => inner(v, x, a0, mn_1, dn_1, an_1),
65   - }
66   -
67   - // This convergence criterium is not good for very small x
68   - // thus I decided to break the iteration at constant time.
69   -// match an_1 == 2 * a0 {
70   -// true => v,
71   -// _ => inner(v, x, a0, mn_1, dn_1, an_1),
72   -// }
73   - }
74   -
75   - inner(v, x, a0, 0, 1, a0)
  42 +pub fn from_vector(xs: &Vec<i64>) -> Vec<Fractional> {
  43 + xs.iter().map(|x| Fractional(*x, 1)).collect()
76 44 }
77 45
78 46 impl Fractional {
... ... @@ -98,76 +66,6 @@ impl Fractional {
98 66 Self(n / hcf(n, d), d / hcf(n, d))
99 67 }
100 68 }
101   -
102   - #[inline]
103   - pub fn numerator(self) -> i64 {
104   - self.0
105   - }
106   -
107   - #[inline]
108   - pub fn denominator(self) -> i64 {
109   - self.1
110   - }
111   -
112   - // This is a really bad approximation of sqrt for a fractional...
113   - // for (9/3) it will result 3 which if way to far from the truth,
114   - // which is ~1.7320508075
115   - // BUT we can use this value as starting guess for creating a
116   - // continous fraction for the sqrt... and create a much better
117   - // fractional representation of the sqrt.
118   - // So, if inner converges, but is not a perfect square (does not
119   - // end up in an Ordering::Equal - which is the l > h case)
120   - // we use the l - 1 as starting guess for sqrt_cfrac.
121   - // taken from:
122   - // https://www.geeksforgeeks.org/square-root-of-an-integer/
123   - pub fn sqrt(self) -> Result<Self, Error> {
124   - // find the sqrt of x in O(log x/2).
125   - // This stops if a perfect sqare was found. Else it passes
126   - // the found value as starting guess to the continous fraction
127   - // sqrt function.
128   - fn floor_sqrt(x :i64) -> Fractional {
129   - fn inner(l :i64, h :i64, x :i64) -> Fractional {
130   - if l > h {
131   - (&sqrt_cfrac(x, l - 1)).into()
132   - } else {
133   - let m = (l + h) / 2;
134   - match x.cmp(&(m * m)) {
135   - Ordering::Equal => m.into(),
136   - Ordering::Less => inner(l, m - 1, x),
137   - Ordering::Greater => inner(m + 1, h, x),
138   - }
139   - }
140   - }
141   -
142   - match x {
143   - 0 => 0.into(),
144   - 1 => 1.into(),
145   - _ => inner(1, x / 2, x),
146   - }
147   - }
148   -
149   - let Fractional(n, d) = self;
150   -
151   - let n = match n.cmp(&0) {
152   - Ordering::Equal => 0.into(),
153   - Ordering::Less => return Err("sqrt on negative undefined"),
154   - Ordering::Greater => floor_sqrt(n),
155   - };
156   -
157   - let d = match d.cmp(&0) {
158   - Ordering::Equal => 0.into(),
159   - Ordering::Less => return Err("sqrt on negative undefined"),
160   - Ordering::Greater => floor_sqrt(d),
161   - };
162   -
163   - Ok(n / d)
164   - }
165   -}
166   -
167   -impl fmt::Display for Fractional {
168   - fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
169   - write!(f, "({}/{})", self.0, self.1)
170   - }
171 69 }
172 70
173 71 impl From<i64> for Fractional {
... ... @@ -176,6 +74,12 @@ impl From<i64> for Fractional {
176 74 }
177 75 }
178 76
  77 +impl From<i32> for Fractional {
  78 + fn from(x: i32) -> Self {
  79 + Self(x as i64, 1)
  80 + }
  81 +}
  82 +
179 83 impl TryFrom<usize> for Fractional {
180 84 type Error = &'static str;
181 85
... ... @@ -188,10 +92,6 @@ impl TryFrom<usize> for Fractional {
188 92 }
189 93 }
190 94
191   -pub fn from_vector(xs: &Vec<i64>) -> Vec<Fractional> {
192   - xs.iter().map(|x| Fractional(*x, 1)).collect()
193   -}
194   -
195 95 impl TryInto<f64> for Fractional {
196 96 type Error = TryFromIntError;
197 97
... ... @@ -212,34 +112,9 @@ impl TryInto<(i32, i32)> for Fractional {
212 112 }
213 113 }
214 114
215   -impl Into<Continuous> for Fractional {
216   - // general continous fraction form of a fractional...
217   - fn into(self) -> Continuous {
218   - let v :Continuous = Vec::new();
219   -
220   - fn inner(mut v :Continuous, f :Fractional) -> Continuous {
221   - let Fractional(n, d) = f;
222   - let a = n / d;
223   - let Fractional(_n, _d) = f - a.into();
224   -
225   - v.push(a);
226   - match _n {
227   - 1 => { v.push(_d); v },
228   - _ => inner(v, Fractional(_d, _n)),
229   - }
230   - }
231   -
232   - inner(v, self)
233   - }
234   -}
235   -
236   -impl From<&Continuous> for Fractional {
237   - fn from(x :&Continuous) -> Self {
238   - let Self(n, d) = x.iter().rev().fold(Fractional(0, 1), |acc, x| {
239   - let Self(an, ad) = acc + (*x).into();
240   - Self(ad, an)
241   - });
242   - Self(d, n)
  115 +impl fmt::Display for Fractional {
  116 + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
  117 + write!(f, "({}/{})", self.0, self.1)
243 118 }
244 119 }
245 120
... ...
... ... @@ -20,11 +20,13 @@
20 20 //
21 21 extern crate lazy_static;
22 22
  23 +pub type Error = &'static str;
  24 +
  25 +pub mod continuous;
23 26 pub mod fractional;
  27 +pub mod transform;
24 28 pub mod trigonometry;
25 29 pub mod vector;
26   -pub mod transform;
27 30
28   -use fractional::{Fractional};
29   -use trigonometry::{sin, cos};
30   -use vector::{Vector};
  31 +use fractional::Fractional;
  32 +use vector::Vector;
... ...
... ... @@ -22,8 +22,9 @@ use std::convert::{TryFrom, TryInto, Into};
22 22 use std::num::TryFromIntError;
23 23 use std::f64::consts::PI as FPI;
24 24
25   -use fractional::fractional::{Fractional, from_vector, Continuous};
26   -use fractional::trigonometry::{sin, cos, tan, PI};
  25 +use fractional::continuous::Continuous;
  26 +use fractional::fractional::{Fractional, from_vector};
  27 +use fractional::trigonometry::Trig;
27 28 use fractional::vector::{Vector};
28 29 use fractional::transform::{translate, rotate_x, rotate_y, rotate_z, rotate_v};
29 30
... ... @@ -51,8 +52,8 @@ fn continuous() {
51 52 let d = Fractional(45, 16);
52 53 let e = Fractional(16, 45);
53 54
54   - let dc :Continuous = d.into();
55   - let ec :Continuous = e.into();
  55 + let dc :Continuous = (&d).into();
  56 + let ec :Continuous = (&e).into();
56 57
57 58 println!("cont frac of d : {} => {:?}", d, dc);
58 59 println!("cont frac of e : {} => {:?}", e, ec);
... ... @@ -82,21 +83,22 @@ fn sqrt() {
82 83 }
83 84
84 85 fn pi() {
85   - let pir :f64 = PI.try_into().unwrap();
86   - let pit :(i32, i32) = PI.try_into().unwrap();
87   - let pi2r :f64 = (PI * PI).try_into().unwrap();
  86 + let pi = Fractional::pi();
  87 + let pir :f64 = pi.try_into().unwrap();
  88 + let pit :(i32, i32) = pi.try_into().unwrap();
  89 + let pi2r :f64 = (pi * pi).try_into().unwrap();
88 90
89 91 println!(" Rust π : {}" , FPI);
90   - println!(" π : {} {}" , PI, pir);
  92 + println!(" π : {} {}" , pi, pir);
91 93 println!(" π as tuple : {:?}" , pit);
92 94 println!(" Rust π² : {}" , FPI * FPI);
93   - println!(" π² : {} {}" , PI * PI, pi2r);
  95 + println!(" π² : {} {}" , pi * pi, pi2r);
94 96 }
95 97
96 98 fn _sin() {
97 99 for d in [ 0, 45, 90, 135, 180, 225, 270, 315
98 100 , 9, 17, 31, 73, 89, 123, 213, 312, 876 ].iter() {
99   - let s = sin(*d as i32);
  101 + let s = Fractional::sin(*d as i32);
100 102 let sr :f64 = s.try_into().unwrap();
101 103 let _s = f64::sin(*d as f64 * FPI / 180.0);
102 104
... ... @@ -107,22 +109,22 @@ fn _sin() {
107 109 fn _tan() {
108 110 for d in [ 0, 45, 90, 135, 180, 225, 270, 315
109 111 , 9, 17, 31, 73, 89, 123, 213, 312, 876 ].iter() {
110   - let s = tan(*d as i32);
111   - let sr :f64 = s.try_into().unwrap();
112   - let _s = f64::tan(*d as f64 * FPI / 180.0);
  112 + let t = Fractional::tan(*d as i32);
  113 + let tr :f64 = t.try_into().unwrap();
  114 + let _t = f64::tan(*d as f64 * FPI / 180.0);
113 115
114   - println!("{:>14} : {} {} / {}", format!("tan {}", d), s, sr, _s);
  116 + println!("{:>14} : {} {} / {}", format!("tan {}", d), t, tr, _t);
115 117 }
116 118 }
117 119
118 120 fn _cos() {
119 121 for d in [ 0, 45, 90, 135, 180, 225, 270, 315
120 122 , 9, 17, 31, 73, 89, 123, 213, 312, 876 ].iter() {
121   - let s = cos(*d as i32);
122   - let sr :f64 = s.try_into().unwrap();
123   - let _s = f64::cos(*d as f64 * FPI / 180.0);
  123 + let c = Fractional::cos(*d as i32);
  124 + let cr :f64 = c.try_into().unwrap();
  125 + let _c = f64::cos(*d as f64 * FPI / 180.0);
124 126
125   - println!("{:>14} : {} {} / {}", format!("cos {}", d), s, sr, _s);
  127 + println!("{:>14} : {} {} / {}", format!("cos {}", d), c, cr, _c);
126 128 }
127 129 }
128 130
... ... @@ -152,8 +154,8 @@ fn _vector() {
152 154 }
153 155
154 156 fn _transform() {
155   - let v = Vector(1.into(), 1.into(), 1.into());
156   - let v1 = Vector(1.into(), 2.into(), 3.into());
  157 + let v = Vector(Fractional(1,1), Fractional(1,1), Fractional(1,1));
  158 + let v1 = Vector(Fractional(1,1), Fractional(2,1), Fractional(3,1));
157 159 let mt = translate(v);
158 160
159 161 println!("{:>14} : {}", "Vector v1", v1);
... ... @@ -184,7 +186,7 @@ fn _transform() {
184 186 }
185 187 println!();
186 188
187   - let v3 = Vector(1.into(), 0.into(), 1.into());
  189 + let v3 = Vector(Fractional(1,1), Fractional(0,1), Fractional(1,1));
188 190 println!("{:>14} : {}", "Vector v3", v3);
189 191 for d in [ 30, 45, 60, 90, 120, 135, 150, 180
190 192 , 210, 225, 240, 270, 300, 315, 330 ].iter() {
... ...
... ... @@ -18,15 +18,19 @@
18 18 // You should have received a copy of the GNU General Public License
19 19 // along with this program. If not, see <http://www.gnu.org/licenses/>.
20 20 //
21   -use std::ops::{Mul};
22   -use crate::{Fractional, cos, sin, Vector};
23   -
24   -pub struct TMatrix( (Fractional, Fractional, Fractional, Fractional)
25   - , (Fractional, Fractional, Fractional, Fractional)
26   - , (Fractional, Fractional, Fractional, Fractional)
27   - , (Fractional, Fractional, Fractional, Fractional) );
28   -
29   -pub fn translate(v :Vector) -> TMatrix {
  21 +use std::ops::{Add, Sub, Neg, Mul, Div};
  22 +use crate::Vector;
  23 +use crate::trigonometry::Trig;
  24 +
  25 +pub struct TMatrix<T>( (T, T, T, T)
  26 + , (T, T, T, T)
  27 + , (T, T, T, T)
  28 + , (T, T, T, T) )
  29 + where T: Add + Sub + Neg + Mul + Div + Trig + From<i32> + Copy;
  30 +
  31 +pub fn translate<T>(v :Vector<T>) -> TMatrix<T>
  32 +where T: Add<Output = T> + Sub<Output = T> + Neg<Output = T>
  33 + + Mul<Output = T> + Div<Output = T> + Trig + From<i32> + Copy {
30 34 let Vector(x, y, z) = v;
31 35
32 36 TMatrix( (1.into(), 0.into(), 0.into(), 0.into())
... ... @@ -35,49 +39,71 @@ pub fn translate(v :Vector) -> TMatrix {
35 39 , ( x, y, z, 1.into()) )
36 40 }
37 41
38   -pub fn rotate_x(a :i32) -> TMatrix {
  42 +pub fn rotate_x<T>(a :i32) -> TMatrix<T>
  43 +where T: Add<Output = T> + Sub<Output = T> + Neg<Output = T>
  44 + + Mul<Output = T> + Div<Output = T> + Trig + From<i32> + Copy {
  45 + let sin :T = Trig::sin(a);
  46 + let cos :T = Trig::cos(a);
  47 +
39 48 TMatrix( (1.into(), 0.into(), 0.into(), 0.into())
40   - , (0.into(), cos(a), -sin(a), 0.into())
41   - , (0.into(), sin(a), cos(a), 0.into())
  49 + , (0.into(), cos , -sin , 0.into())
  50 + , (0.into(), sin , cos , 0.into())
42 51 , (0.into(), 0.into(), 0.into(), 1.into()) )
43 52 }
44 53
45   -pub fn rotate_y(a :i32) -> TMatrix {
46   - TMatrix( ( cos(a), 0.into(), sin(a), 0.into())
  54 +pub fn rotate_y<T>(a :i32) -> TMatrix<T>
  55 +where T: Add<Output = T> + Sub<Output = T> + Neg<Output = T>
  56 + + Mul<Output = T> + Div<Output = T> + Trig + From<i32> + Copy {
  57 + let sin :T = Trig::sin(a);
  58 + let cos :T = Trig::cos(a);
  59 +
  60 + TMatrix( (cos , 0.into(), sin , 0.into())
47 61 , (0.into(), 1.into(), 0.into(), 0.into())
48   - , ( -sin(a), 0.into(), cos(a), 0.into())
  62 + , (-sin , 0.into(), cos , 0.into())
49 63 , (0.into(), 0.into(), 0.into(), 1.into()) )
50 64 }
51 65
52   -pub fn rotate_z(a :i32) -> TMatrix {
53   - TMatrix( ( cos(a), -sin(a), 0.into(), 0.into())
54   - , ( sin(a), cos(a), 0.into(), 0.into())
  66 +pub fn rotate_z<T>(a :i32) -> TMatrix<T>
  67 +where T: Add<Output = T> + Sub<Output = T> + Neg<Output = T>
  68 + + Mul<Output = T> + Div<Output = T> + Trig + From<i32> + Copy {
  69 + let sin :T = Trig::sin(a);
  70 + let cos :T = Trig::cos(a);
  71 +
  72 + TMatrix( (cos , -sin , 0.into(), 0.into())
  73 + , (sin , cos , 0.into(), 0.into())
55 74 , (0.into(), 0.into(), 1.into(), 0.into())
56 75 , (0.into(), 0.into(), 0.into(), 1.into()) )
57 76 }
58 77
59   -pub fn rotate_v(v :&Vector, a :i32) -> TMatrix {
  78 +pub fn rotate_v<T>(v :&Vector<T>, a :i32) -> TMatrix<T>
  79 +where T: Add<Output = T> + Sub<Output = T> + Neg<Output = T>
  80 + + Mul<Output = T> + Div<Output = T> + Trig + From<i32> + Copy {
60 81 let Vector(x, y, z) = *v;
61 82
62   - let zero :Fractional = 0.into();
63   - let one :Fractional = 1.into();
  83 + let sin :T = Trig::sin(a);
  84 + let cos :T = Trig::cos(a);
  85 +
  86 + let zero :T = 0.into();
  87 + let one :T = 1.into();
64 88
65   - TMatrix( ( (one - cos(a)) * x * x + cos(a)
66   - , (one - cos(a)) * x * y - sin(a) * z
67   - , (one - cos(a)) * x * z + sin(a) * y
  89 + TMatrix( ( (one - cos) * x * x + cos
  90 + , (one - cos) * x * y - sin * z
  91 + , (one - cos) * x * z + sin * y
68 92 , zero )
69   - , ( (one - cos(a)) * x * y + sin(a) * z
70   - , (one - cos(a)) * y * y + cos(a)
71   - , (one - cos(a)) * y * z - sin(a) * x
  93 + , ( (one - cos) * x * y + sin * z
  94 + , (one - cos) * y * y + cos
  95 + , (one - cos) * y * z - sin * x
72 96 , zero )
73   - , ( (one - cos(a)) * x * z - sin(a) * y
74   - , (one - cos(a)) * y * z + sin(a) * x
75   - , (one - cos(a)) * z * z + cos(a)
  97 + , ( (one - cos) * x * z - sin * y
  98 + , (one - cos) * y * z + sin * x
  99 + , (one - cos) * z * z + cos
76 100 , zero )
77 101 , (0.into(), 0.into(), 0.into(), 1.into()) )
78 102 }
79 103
80   -pub fn scale(v :Vector) -> TMatrix {
  104 +pub fn scale<T>(v :Vector<T>) -> TMatrix<T>
  105 +where T: Add<Output = T> + Sub<Output = T> + Neg<Output = T>
  106 + + Mul<Output = T> + Div<Output = T> + Trig + From<i32> + Copy {
81 107 let Vector(x, y, z) = v;
82 108
83 109 TMatrix( ( x, 0.into(), 0.into(), 0.into())
... ... @@ -86,8 +112,10 @@ pub fn scale(v :Vector) -> TMatrix {
86 112 , (0.into(), 0.into(), 0.into(), 1.into()) )
87 113 }
88 114
89   -impl TMatrix {
90   - pub fn apply(&self, v :&Vector) -> Vector {
  115 +impl<T> TMatrix<T>
  116 +where T: Add<Output = T> + Sub<Output = T> + Neg<Output = T>
  117 + + Mul<Output = T> + Div<Output = T> + Trig + From<i32> + Copy {
  118 + pub fn apply(&self, v :&Vector<T>) -> Vector<T> {
91 119 let TMatrix( (a11, a12, a13, a14)
92 120 , (a21, a22, a23, a24)
93 121 , (a31, a32, a33, a34)
... ... @@ -97,13 +125,15 @@ impl TMatrix {
97 125 let v = Vector( a11 * x + a21 * y + a31 * z + a41
98 126 , a12 * x + a22 * y + a32 * z + a42
99 127 , a13 * x + a23 * y + a33 * z + a43 );
100   - let Fractional(wn, wd) = a14 * x + a24 * y + a34 * z + a44;
  128 + let w = a14 * x + a24 * y + a34 * z + a44;
101 129
102   - v.mul(&Fractional(wd, wn))
  130 + v.mul(&w.recip())
103 131 }
104 132 }
105 133
106   -impl Mul for TMatrix {
  134 +impl<T> Mul for TMatrix<T>
  135 +where T: Add<Output = T> + Sub<Output = T> + Neg<Output = T>
  136 + + Mul<Output = T> + Div<Output = T> + Trig + From<i32> + Copy {
107 137 type Output = Self;
108 138
109 139 // ATTENTION: This is not commutative, nor assoziative.
... ...
... ... @@ -26,90 +26,230 @@
26 26 // along with this program. If not, see <http://www.gnu.org/licenses/>.
27 27 //
28 28 use std::cmp::Ordering;
29   -use crate::{Fractional};
  29 +use std::ops::Neg;
  30 +use std::marker::Sized;
  31 +use crate::{Fractional, Error};
  32 +use crate::continuous::Continuous;
30 33
31   -pub const PI :Fractional = Fractional(355, 113); // This is a really close
32   - // fractional approximation
33   - // for pi
  34 +pub trait Trig {
  35 + fn pi() -> Self;
  36 + fn recip(self) -> Self;
  37 + fn sqrt(self) -> Result<Self, Error> where Self: Sized;
  38 + fn sintab() -> Vec<Self> where Self: Sized;
  39 + fn tantab() -> Vec<Self> where Self: Sized;
  40 +
  41 + fn sin(d :i32) -> Self
  42 + where Self: Sized + Neg<Output = Self> + Copy {
  43 + match d {
  44 + 0 ..=90 => Self::sintab()[d as usize],
  45 + 91 ..=180 => Self::sintab()[180 - d as usize],
  46 + 181..=270 => -Self::sintab()[d as usize - 180],
  47 + 271..=359 => -Self::sintab()[360 - d as usize],
  48 + _ => Self::sin(d % 360),
  49 + }
  50 + }
  51 +
  52 + fn cos(d :i32) -> Self
  53 + where Self: Sized + Neg<Output = Self> + Copy {
  54 + match d {
  55 + 0 ..=90 => Self::sintab()[90 - d as usize],
  56 + 91 ..=180 => -Self::sintab()[90 - (180 - d as usize)],
  57 + 181..=270 => -Self::sintab()[90 - (d as usize - 180)],
  58 + 271..=359 => Self::sintab()[90 - (360 - d as usize)],
  59 + _ => Self::cos(d % 360),
  60 + }
  61 + }
  62 +
  63 + fn tan(d :i32) -> Self where Self: Sized + Copy {
  64 + match d {
  65 + 0 ..=179 => Self::tantab()[d as usize],
  66 + 180..=359 => Self::tantab()[d as usize - 180],
  67 + _ => Self::tan(d % 360),
  68 + }
  69 + }
  70 +}
34 71
35 72 // Try to keep precision as high as possible while having a denominator
36   -// as small as possible.
37   -// The values are taken by triel and error.
  73 +// as small as possible. The values are taken by try and error.
38 74 const PRECISION :i64 = 1000000;
39 75 const MAX_DENOMINATOR :i64 = 7000;
40 76
41   -pub fn sin(d :i32) -> Fractional {
42   - match d {
43   - 0 ..=90 => SINTAB[d as usize],
44   - 91 ..=180 => SINTAB[180 - d as usize],
45   - 181..=270 => -SINTAB[d as usize - 180],
46   - 271..=359 => -SINTAB[360 - d as usize],
47   - _ => sin(d % 360),
  77 +// This is a really close fractional approximation for pi.
  78 +impl Trig for Fractional {
  79 + fn pi() -> Self {
  80 + Fractional(355, 113)
48 81 }
49   -}
50 82
51   -pub fn cos(d :i32) -> Fractional {
52   - match d {
53   - 0 ..=90 => SINTAB[90 - d as usize],
54   - 91 ..=180 => -SINTAB[90 - (180 - d as usize)],
55   - 181..=270 => -SINTAB[90 - (d as usize - 180)],
56   - 271..=359 => SINTAB[90 - (360 - d as usize)],
57   - _ => cos(d % 360),
  83 + fn recip(self) -> Self {
  84 + let Fractional(n, d) = self;
  85 + Fractional(d, n)
58 86 }
59   -}
60 87
61   -pub fn tan(d :i32) -> Fractional {
62   - match d {
63   - 0 ..=179 => TANTAB[d as usize],
64   - 180..=359 => TANTAB[d as usize - 180],
65   - _ => tan(d % 360),
  88 + // This is a really bad approximation of sqrt for a fractional...
  89 + // for (9/3) it will result 3 which if way to far from the truth,
  90 + // which is ~1.7320508075
  91 + // BUT we can use this value as starting guess for creating a
  92 + // continous fraction for the sqrt... and create a much better
  93 + // fractional representation of the sqrt.
  94 + // So, if inner converges, but is not a perfect square (does not
  95 + // end up in an Ordering::Equal - which is the l > h case)
  96 + // we use the l - 1 as starting guess for sqrt_cfrac.
  97 + // taken from:
  98 + // https://www.geeksforgeeks.org/square-root-of-an-integer/
  99 + fn sqrt(self) -> Result<Self, Error> {
  100 + // find the sqrt of x in O(log x/2).
  101 + // This stops if a perfect sqare was found. Else it passes
  102 + // the found value as starting guess to the continous fraction
  103 + // sqrt function.
  104 + fn floor_sqrt(x :i64) -> Fractional {
  105 + fn inner(l :i64, h :i64, x :i64) -> Fractional {
  106 + if l > h {
  107 + (&Continuous::sqrt(x, l - 1)).into()
  108 + } else {
  109 + let m = (l + h) / 2;
  110 + match x.cmp(&(m * m)) {
  111 + Ordering::Equal => m.into(),
  112 + Ordering::Less => inner(l, m - 1, x),
  113 + Ordering::Greater => inner(m + 1, h, x),
  114 + }
  115 + }
  116 + }
  117 +
  118 + match x {
  119 + 0 => 0.into(),
  120 + 1 => 1.into(),
  121 + _ => inner(1, x / 2, x),
  122 + }
  123 + }
  124 +
  125 + let Fractional(n, d) = self;
  126 +
  127 + let n = match n.cmp(&0) {
  128 + Ordering::Equal => 0.into(),
  129 + Ordering::Less => return Err("sqrt on negative undefined"),
  130 + Ordering::Greater => floor_sqrt(n),
  131 + };
  132 +
  133 + let d = match d.cmp(&0) {
  134 + Ordering::Equal => 0.into(),
  135 + Ordering::Less => return Err("sqrt on negative undefined"),
  136 + Ordering::Greater => floor_sqrt(d),
  137 + };
  138 +
  139 + Ok(n / d)
66 140 }
67   -}
68 141
69   -// hold sin Fractionals from 0 to 89 ...
70   -// luckily with a bit of index tweeking this can also be used for
71   -// cosine values.
72   -lazy_static::lazy_static! {
73   - static ref SINTAB :Vec<Fractional> =
74   - (0..=90).map(|x| _sin(x)).collect();
75   -}
  142 + fn sintab() -> Vec<Self> {
  143 + // hold sin Fractionals from 0 to 89 ...
  144 + // luckily with a bit of index tweeking this can also be used for
  145 + // cosine values.
  146 + lazy_static::lazy_static! {
  147 + static ref SINTAB :Vec<Fractional> =
  148 + (0..=90).map(|x| _sin(x)).collect();
  149 + }
76 150
77   -// This table exists only because the sin(α) / cos(α) method
78   -// yields very large unreducable denominators in a lot of cases.
79   -lazy_static::lazy_static! {
80   - static ref TANTAB :Vec<Fractional> =
81   - (0..180).map(|x| _tan(x)).collect();
82   -}
  151 + // fractional sin from f64 sin. (From 0° to 90°)
  152 + fn _sin(d: u32) -> Fractional {
  153 + match d {
  154 + 0 => Fractional(0, 1),
  155 + 90 => Fractional(1, 1),
  156 + _ => reduce(d, PRECISION, &f64::sin),
  157 + }
  158 + }
83 159
84   -// fractional sin from f64 sin. (From 0° to 90°)
85   -fn _sin(d: u32) -> Fractional {
86   - match d {
87   - 0 => Fractional(0, 1),
88   - 90 => Fractional(1, 1),
89   - _ => reduce(d, PRECISION, &f64::sin),
  160 + SINTAB.to_vec()
  161 + }
  162 +
  163 + fn tantab() -> Vec<Self> {
  164 + // This table exists only because the sin(α) / cos(α) method
  165 + // yields very large unreducable denominators in a lot of cases.
  166 + lazy_static::lazy_static! {
  167 + static ref TANTAB :Vec<Fractional> =
  168 + (0..180).map(|x| _tan(x)).collect();
  169 + }
  170 +
  171 + // fractional tan from f64 tan. (From 0° to 179°)
  172 + fn _tan(d: u32) -> Fractional {
  173 + match d {
  174 + 0 => Fractional(0, 1),
  175 + 45 => Fractional(1, 1),
  176 + 90 => Fractional(1, 0), // although they are both inf and -inf.
  177 + 135 => -Fractional(1, 1),
  178 + _ => reduce(d, PRECISION, &f64::tan),
  179 + }
  180 + }
  181 +
  182 + TANTAB.to_vec()
90 183 }
91 184 }
92 185
93   -// fractional tan from f64 tan. (From 0° to 179°)
94   -fn _tan(d: u32) -> Fractional {
95   - match d {
96   - 0 => Fractional(0, 1),
97   - 45 => Fractional(1, 1),
98   - 90 => Fractional(1, 0), // although they are both inf and -inf.
99   - 135 => -Fractional(1, 1),
100   - _ => reduce(d, PRECISION, &f64::tan),
  186 +impl Trig for f64 {
  187 + fn pi() -> Self {
  188 + std::f64::consts::PI
  189 + }
  190 +
  191 + fn recip(self) -> Self {
  192 + self.recip()
  193 + }
  194 +
  195 + fn sqrt(self) -> Result<Self, Error> {
  196 + let x = self.sqrt();
  197 + match x.is_nan() {
  198 + true => Err("sqrt on negative undefined"),
  199 + false => Ok(x),
  200 + }
  201 + }
  202 +
  203 + fn sintab() -> Vec<Self> {
  204 + lazy_static::lazy_static! {
  205 + static ref SINTAB :Vec<f64> =
  206 + (0..=90).map(|x| _sin(x)).collect();
  207 + }
  208 +
  209 + // f64 sin. (From 0° to 90°)
  210 + fn _sin(d: u32) -> f64 {
  211 + match d {
  212 + 0 => 0.0,
  213 + 90 => 1.0,
  214 + _ => (d as f64).to_radians().sin(),
  215 + }
  216 + }
  217 +
  218 + SINTAB.to_vec()
  219 + }
  220 +
  221 + fn tantab() -> Vec<Self> {
  222 + // This table exists only because the sin(α) / cos(α) method
  223 + // yields very large unreducable denominators in a lot of cases.
  224 + lazy_static::lazy_static! {
  225 + static ref TANTAB :Vec<f64> =
  226 + (0..180).map(|x| _tan(x)).collect();
  227 + }
  228 +
  229 + // fractional tan from f64 tan. (From 0° to 179°)
  230 + fn _tan(d: u32) -> f64 {
  231 + match d {
  232 + 0 => 0.0,
  233 + 45 => 1.0,
  234 + 90 => std::f64::INFINITY,
  235 + 135 => -1.0,
  236 + _ => (d as f64).to_radians().tan(),
  237 + }
  238 + }
  239 +
  240 + TANTAB.to_vec()
101 241 }
102 242 }
103 243
104   -// search for a fraction with a denominator less than 10000 that
105   -// provides the minimal precision criteria.
  244 +// search for a fraction with a denominator less than MAX_DENOMINATOR that
  245 +// provides the minimal PRECISION criteria.
106 246 // !! With f = &f64::tan and d close to the inf boundarys of tan
107 247 // we get very large numerators because the numerator becomes a
108 248 // multiple of the denominator.
109 249 fn reduce(d :u32, p :i64, f :&dyn Fn(f64) -> f64) -> Fractional {
110 250 // This is undefined behaviour for very large f64, but our f64
111   - // is always between 0.0 and 10000000.0 which should be fine.
112   - let s = (f(f64::to_radians(d as f64)) * p as f64).round() as i64;
  251 + // is always between 0.0 and 1000000.0 which should be fine.
  252 + let s = (f((d as f64).to_radians()) * p as f64).round() as i64;
113 253 let Fractional(n, dn) = Fractional(s, p).reduce();
114 254 match dn.abs().cmp(&MAX_DENOMINATOR) {
115 255 Ordering::Less => Fractional(n, dn),
... ...
... ... @@ -18,29 +18,32 @@
18 18 // You should have received a copy of the GNU General Public License
19 19 // along with this program. If not, see <http://www.gnu.org/licenses/>.
20 20 //
21   -use std::ops::{Add, Sub, Neg, Mul};
  21 +use std::ops::{Add, Sub, Neg, Mul, Div};
22 22 use std::fmt;
23   -use crate::{Fractional};
  23 +use crate::trigonometry::Trig;
24 24
25 25 #[derive(Debug, Eq, Clone, Copy)]
26   -pub struct Vector(pub Fractional, pub Fractional, pub Fractional);
  26 +pub struct Vector<T>(pub T, pub T, pub T)
  27 + where T: Add + Sub + Neg + Mul + Div + Trig + Copy;
27 28
28   -impl Vector {
29   - pub fn x(self) -> Fractional { self.0 }
30   - pub fn y(self) -> Fractional { self.1 }
31   - pub fn z(self) -> Fractional { self.2 }
  29 +impl<T> Vector<T>
  30 +where T: Add<Output = T> + Sub<Output = T> + Neg<Output = T>
  31 + + Mul<Output = T> + Div<Output = T> + Trig + Copy {
  32 + pub fn x(self) -> T { self.0 }
  33 + pub fn y(self) -> T { self.1 }
  34 + pub fn z(self) -> T { self.2 }
32 35
33   - pub fn abs(self) -> Fractional {
  36 + pub fn abs(self) -> T {
34 37 let Vector(x, y, z) = self;
35 38 (x * x + y * y + z * z).sqrt().unwrap()
36 39 }
37 40
38   - pub fn mul(self, s :&Fractional) -> Self {
  41 + pub fn mul(self, s :&T) -> Self {
39 42 let Vector(x, y, z) = self;
40 43 Vector(x * *s, y * *s, z * *s)
41 44 }
42 45
43   - pub fn dot(self, other :Self) -> Fractional {
  46 + pub fn dot(self, other :Self) -> T {
44 47 let Vector(x1, y1, z1) = self;
45 48 let Vector(x2, y2, z2) = other;
46 49
... ... @@ -48,24 +51,25 @@ impl Vector {
48 51 }
49 52
50 53 pub fn norm(self) -> Self {
51   - let Fractional(n, d) = self.abs();
52 54 // TODO This can result in 0 or inf Vectors…
53 55 // Maybe we need to handle zero and inf abs here…
54   - self.mul(&Fractional(d, n))
  56 + self.mul(&self.abs())
55 57 }
56 58
57   - pub fn distance(self, other :Self) -> Fractional {
  59 + pub fn distance(self, other :Self) -> T {
58 60 (self - other).abs()
59 61 }
60 62 }
61 63
62   -impl fmt::Display for Vector {
  64 +impl<T> fmt::Display for Vector<T>
  65 +where T: Add + Sub + Neg + Mul + Div + Trig + fmt::Display + Copy {
63 66 fn fmt(&self, f :&mut fmt::Formatter<'_>) -> fmt::Result {
64 67 write!(f, "({}, {}, {})", self.0, self.1, self.2)
65 68 }
66 69 }
67 70
68   -impl PartialEq for Vector {
  71 +impl<T> PartialEq for Vector<T>
  72 +where T: Add + Sub + Neg + Mul + Div + Trig + PartialEq + Copy {
69 73 fn eq(&self, other :&Self) -> bool {
70 74 let Vector(x1, y1, z1) = self;
71 75 let Vector(x2, y2, z2) = other;
... ... @@ -73,7 +77,9 @@ impl PartialEq for Vector {
73 77 }
74 78 }
75 79
76   -impl Add for Vector {
  80 +impl<T> Add for Vector<T>
  81 +where T: Add<Output = T> + Sub<Output = T> + Neg<Output = T>
  82 + + Mul<Output = T> + Div<Output = T> + Trig + Copy {
77 83 type Output = Self;
78 84
79 85 fn add(self, other :Self) -> Self {
... ... @@ -83,7 +89,9 @@ impl Add for Vector {
83 89 }
84 90 }
85 91
86   -impl Sub for Vector {
  92 +impl<T> Sub for Vector<T>
  93 +where T: Add<Output = T> + Sub<Output = T> + Neg<Output = T>
  94 + + Mul<Output = T> + Div<Output = T> + Trig + Copy {
87 95 type Output = Self;
88 96
89 97 fn sub(self, other :Self) -> Self {
... ... @@ -91,7 +99,9 @@ impl Sub for Vector {
91 99 }
92 100 }
93 101
94   -impl Neg for Vector {
  102 +impl<T> Neg for Vector<T>
  103 +where T: Add<Output = T> + Sub<Output = T> + Neg<Output = T>
  104 + + Mul<Output = T> + Div<Output = T> + Trig + Copy {
95 105 type Output = Self;
96 106
97 107 fn neg(self) -> Self {
... ... @@ -100,7 +110,9 @@ impl Neg for Vector {
100 110 }
101 111 }
102 112
103   -impl Mul for Vector {
  113 +impl<T> Mul for Vector<T>
  114 +where T: Add<Output = T> + Sub<Output = T> + Neg<Output = T>
  115 + + Mul<Output = T> + Div<Output = T> + Trig + Copy {
104 116 type Output = Self;
105 117
106 118 fn mul(self, other :Self) -> Self {
... ...
Please register or login to post a comment