feo_math/
imaginary.rs

1//! Light implementations for imaginary numbers 
2//! 
3//! Quaternions are under rotations.
4//! 
5//! TODO
6//! 
7use std::ops::{Add, Div, Mul, Neg, Sub};
8
9use crate::{F32Fmt, One, SignOps, Three, Two, Zero};
10
11use super::Construct;
12
13pub trait ImaginaryConstruct<T>: Construct<T> { /* todo */ }
14
15// /// Imaginary construct type .0 n .1
16// pub struct Imaginary<T>(char, T);
17
18// pub struct Complex<T>(T, Vec<Imaginary<T>>);
19
20/// For now
21#[derive(Clone, Copy)]
22pub struct Simple2DComplexNumber<T>(pub T, pub T);
23
24impl<T> Simple2DComplexNumber<T> {
25    pub fn is_real (self) -> bool 
26    where T: Zero + PartialEq {
27        self.1 == T::ZERO
28    }
29
30    pub fn to_real(self) -> Result<T, ()> 
31    where T: Zero + PartialEq {
32        if self.1 == T::ZERO {
33            Ok(self.0)            
34        } else {
35            Err(())
36        }
37    }
38
39    pub fn cbrt (self) -> [Self; 3]
40    where T: F32Fmt + PartialEq + Add<T, Output = T> + Sub<T, Output = T> + Mul<T, Output = T> + One + Two + Three + SignOps + Div<T, Output = T> + Copy {
41        let Simple2DComplexNumber(real, imag) = self;
42        let r = (real * real + imag * imag).sqrt();
43        let angle_pfi = T::atan2_mul(imag, real, T::ONE);
44        
45        let cbrt_r = r.cbrt();
46        // Answers: 
47        // 1)
48        // r.cbrt() * (e^(iTheta/3)) = r.cbrt() * (e^(i(Theta/3))
49        // Eulers formula: (e^(i(Theta/3)) = cos(Theta/3) + sin(Theta/3)i
50        // r.cbrt() * cos(Theta/3) + sin(Theta/3)i
51        // r.cbrt() * cos(Theta/3) + (r.cbrt() * sin(Theta/3))i
52        // 2)
53        // r.cbrt() * (e^(iTheta/3 + 2iPI/3)) = r.cbrt() * (e^(i((Theta+2PI)/3))
54        // ... Eulers formula
55        // r.cbrt() * cos((Theta+2PI)/3) + sin((Theta+2PI)/3)i
56        // (r.cbrt() * cos((Theta+2PI)/3)) + (r.cbrt() * sin((Theta+2PI)/3))i
57        // 3)
58        // r.cbrt() * (e^(iTheta/3 - 2iPI/3)) = r.cbrt() * (e^(i((Theta-2PI)/3))
59        // ... Eulers formula
60        // r.cbrt() * cos((Theta-2PI)/3) + sin((Theta-2PI)/3)i
61        // (r.cbrt() * cos((Theta-2PI)/3)) + (r.cbrt() * sin((Theta-2PI)/3))i
62
63        let block0 = angle_pfi / T::THREE;
64        let answer0 = Simple2DComplexNumber(T::cos_mul(block0, cbrt_r), T::sin_mul(block0, cbrt_r));
65        let block1 = (angle_pfi + T::f32_const_mul(T::TWO, std::f32::consts::PI)) / T::THREE;
66        let answer1 = Simple2DComplexNumber(T::cos_mul(block1, cbrt_r), T::sin_mul(block1, cbrt_r));
67        let block2 = (angle_pfi - T::f32_const_mul(T::TWO, std::f32::consts::PI)) / T::THREE;
68        let answer2 = Simple2DComplexNumber(T::cos_mul(block2, cbrt_r), T::sin_mul(block2, cbrt_r));
69
70        [answer0, answer1, answer2]
71    }
72}
73
74impl<T> Add<Simple2DComplexNumber<T>> for Simple2DComplexNumber<T> 
75where T: Add<T, Output = T> + Copy {
76    type Output = Self;
77
78    fn add(self, rhs: Simple2DComplexNumber<T>) -> Self::Output {
79        Simple2DComplexNumber(self.0 + rhs.0, self.1 + rhs.1)
80    }
81}
82
83impl<T> Add<T> for Simple2DComplexNumber<T> 
84where T: Add<T, Output = T> + Copy {
85    type Output = Self;
86
87    fn add(self, rhs: T) -> Self::Output {
88        Simple2DComplexNumber(self.0 + rhs, self.1)
89    }
90}
91
92impl<T> Sub<Simple2DComplexNumber<T>> for Simple2DComplexNumber<T> 
93where T: Sub<T, Output = T> + Copy {
94    type Output = Self;
95
96    fn sub(self, rhs: Simple2DComplexNumber<T>) -> Self::Output {
97        Simple2DComplexNumber(self.0 - rhs.0, self.1 - rhs.1)
98    }
99}
100
101impl<T> Sub<T> for Simple2DComplexNumber<T> 
102where T: Sub<T, Output = T> + Copy {
103    type Output = Self;
104
105    fn sub(self, rhs: T) -> Self::Output {
106        Simple2DComplexNumber(self.0 - rhs, self.1)
107    }
108}
109
110impl<T> Mul<Simple2DComplexNumber<T>> for Simple2DComplexNumber<T> 
111where T: Add<T, Output = T> + Sub<T, Output = T> + Mul<T, Output = T> + Sub<T, Output = T> + Copy {
112    type Output = Self;
113
114    fn mul(self, rhs: Simple2DComplexNumber<T>) -> Self::Output {
115        Simple2DComplexNumber(self.0 * rhs.0 - self.1 * rhs.1, self.0 * rhs.1 + self.1 * rhs.0)
116    }
117}
118
119impl<T> Mul<T> for Simple2DComplexNumber<T> 
120where T: Mul<T, Output = T> + Copy {
121    type Output = Self;
122
123    fn mul(self, rhs: T) -> Self::Output {
124        Simple2DComplexNumber(self.0 * rhs, self.1 * rhs)
125    }
126}
127
128impl<T> Div<Simple2DComplexNumber<T>> for Simple2DComplexNumber<T> 
129where T: Add<T, Output = T> + Sub<T, Output = T> + Mul<T, Output = T> + Div<T, Output = T> + Copy {
130    type Output = Self;
131
132    fn div(self, rhs: Simple2DComplexNumber<T>) -> Self::Output {
133        let numerator = (self.0 * rhs.0 +/*-*/ self.1 * /*-*/rhs.1, self.1 * rhs.0 -/*+*/ self.0 * /*-*/rhs.1);
134        let real_denominator = rhs.0 * rhs.0 +/*-*/ rhs.1 * /*-*/rhs.1;
135        Simple2DComplexNumber(numerator.0 / real_denominator, numerator.1 / real_denominator)
136    }
137}
138
139impl<T> Div<T> for Simple2DComplexNumber<T> 
140where T: Div<T, Output = T> + Copy {
141    type Output = Self;
142
143    fn div(self, rhs: T) -> Self::Output {
144        Simple2DComplexNumber(self.0 / rhs, self.1 / rhs)
145    }
146}
147
148impl<T> Neg for Simple2DComplexNumber<T> 
149where T: Neg<Output = T> + Copy {
150    type Output = Self;
151
152    fn neg(self) -> Self::Output {
153        Simple2DComplexNumber(-self.0, -self.1)
154    }
155}
156
157#[derive(Debug, PartialEq, Clone, Copy)]
158pub enum PossibleImag<T> {
159    Complex(T, T),
160    Imag(T),
161    Real(T),
162}
163impl<T> PossibleImag<T> {
164    pub fn new(real: T, imag: T) -> Self 
165    where T: Zero + PartialEq {
166        if imag == T::ZERO {
167            PossibleImag::Real(real)
168        } else if real == T::ZERO {
169            PossibleImag::Imag(imag)
170        } else {
171            PossibleImag::Complex(real, imag)
172        }
173    }
174
175    pub fn do_sqrt(other: T) -> Self 
176    where T: F32Fmt + SignOps + Clone {
177        let n = other.clone().abs().sqrt();
178        if other.ptsignum() == -1 {
179            PossibleImag::Imag(n)
180        } else {
181            PossibleImag::Real(n)
182        }
183    }
184
185    pub fn cbrt(self) -> [Self; 3]
186    where T: Zero + F32Fmt + PartialEq + Add<T, Output = T> + Sub<T, Output = T> + Mul<T, Output = T> + One + Two + Three + SignOps + Div<T, Output = T> + Copy {
187        let (real, imag) = match self {
188            PossibleImag::Complex(r, i) => (r, i),
189            PossibleImag::Imag(i) => (T::ZERO, i),
190            PossibleImag::Real(r) => (r, T::ZERO)
191        };
192        let r = (real * real + imag * imag).sqrt();
193        let angle_pfi = T::atan2_mul(imag, real, T::ONE);
194        
195        let cbrt_r = r.cbrt();
196        // Answers: 
197        // 1)
198        // r.cbrt() * (e^(iTheta/3)) = r.cbrt() * (e^(i(Theta/3))
199        // Eulers formula: (e^(i(Theta/3)) = cos(Theta/3) + sin(Theta/3)i
200        // r.cbrt() * cos(Theta/3) + sin(Theta/3)i
201        // r.cbrt() * cos(Theta/3) + (r.cbrt() * sin(Theta/3))i
202        // 2)
203        // r.cbrt() * (e^(iTheta/3 + 2iPI/3)) = r.cbrt() * (e^(i((Theta+2PI)/3))
204        // ... Eulers formula
205        // r.cbrt() * cos((Theta+2PI)/3) + sin((Theta+2PI)/3)i
206        // (r.cbrt() * cos((Theta+2PI)/3)) + (r.cbrt() * sin((Theta+2PI)/3))i
207        // 3)
208        // r.cbrt() * (e^(iTheta/3 - 2iPI/3)) = r.cbrt() * (e^(i((Theta-2PI)/3))
209        // ... Eulers formula
210        // r.cbrt() * cos((Theta-2PI)/3) + sin((Theta-2PI)/3)i
211        // (r.cbrt() * cos((Theta-2PI)/3)) + (r.cbrt() * sin((Theta-2PI)/3))i
212
213        let block0 = angle_pfi/T::THREE;
214        let answer0 = PossibleImag::new(T::cos_mul(block0, cbrt_r), T::sin_mul(block0, cbrt_r));
215        let block1 = (angle_pfi + T::f32_const_mul(T::TWO, std::f32::consts::PI)) / T::THREE;
216        let answer1 = PossibleImag::new(T::cos_mul(block1, cbrt_r), T::sin_mul(block1, cbrt_r));
217        let block2 = (angle_pfi - T::f32_const_mul(T::TWO, std::f32::consts::PI)) / T::THREE;
218        let answer2 = PossibleImag::new(T::cos_mul(block2, cbrt_r), T::sin_mul(block2, cbrt_r));
219
220        [answer0, answer1, answer2]
221    }
222}
223
224impl<T: Add<T, Output = T> + PartialEq + Zero> Add<PossibleImag<T>> for PossibleImag<T> {
225    type Output = PossibleImag<T>;
226
227    fn add(self, rhs: PossibleImag<T>) -> Self::Output {
228        let (mut real, mut imag) = match self {
229            PossibleImag::Complex(r, i) => (r, i),
230            PossibleImag::Imag(i) => (T::ZERO, i),
231            PossibleImag::Real(r) => (r, T::ZERO)
232        };
233
234        match rhs {
235            PossibleImag::Complex(r, i) => {
236                real = real + r;
237                imag = imag + i;
238            },
239            PossibleImag::Imag(i) => imag = imag + i,
240            PossibleImag::Real(r) => real = real + r
241        }
242
243        match (real, imag) {
244            (real, zero) if zero == T::ZERO => PossibleImag::Real(real),
245            (zero, imag) if zero == T::ZERO => PossibleImag::Imag(imag),
246            (real, imag) => PossibleImag::Complex(real, imag)
247        }
248        
249    }
250}
251
252impl<T: Add<T, Output = T>> Add<T> for PossibleImag<T> {
253    type Output = PossibleImag<T>;
254
255    fn add(self, rhs: T) -> Self::Output {
256        match self {
257            PossibleImag::Complex(r, i) => PossibleImag::Complex(r + rhs, i),
258            PossibleImag::Imag(i) => PossibleImag::Complex(rhs, i),
259            PossibleImag::Real(r) => PossibleImag::Real(r + rhs)
260        }
261    }
262}
263
264impl<T: Sub<T, Output = T> + PartialEq + Zero> Sub<PossibleImag<T>> for PossibleImag<T> {
265    type Output = PossibleImag<T>;
266
267    fn sub(self, rhs: PossibleImag<T>) -> Self::Output {
268        let (mut real, mut imag) = match self {
269            PossibleImag::Complex(r, i) => (r, i),
270            PossibleImag::Imag(i) => (T::ZERO, i),
271            PossibleImag::Real(r) => (r, T::ZERO)
272        };
273
274        match rhs {
275            PossibleImag::Complex(r, i) => {
276                real = real - r;
277                imag = imag - i;
278            },
279            PossibleImag::Imag(i) => imag = imag - i,
280            PossibleImag::Real(r) => real = real - r
281        }
282
283        match (real, imag) {
284            (real, zero) if zero == T::ZERO  => PossibleImag::Real(real),
285            (zero, imag) if zero == T::ZERO  => PossibleImag::Imag(imag),
286            (real, imag) => PossibleImag::Complex(real, imag)
287        }
288    }
289}
290
291impl<T: Sub<T, Output = T> + Neg<Output = T>> Sub<T> for PossibleImag<T> {
292    type Output = PossibleImag<T>;
293
294    fn sub(self, rhs: T) -> Self::Output {
295        match self {
296            PossibleImag::Complex(r, i) => PossibleImag::Complex(r - rhs, i),
297            PossibleImag::Imag(i) => PossibleImag::Complex(-rhs, i),
298            PossibleImag::Real(r) => PossibleImag::Real(r - rhs)
299        }
300    }
301}
302
303impl<T> Mul<PossibleImag<T>> for PossibleImag<T> 
304where T: Add<T, Output = T> + Sub<T, Output = T> + Mul<T, Output = T> + Neg<Output = T> + Zero + PartialEq + Copy {
305    type Output = PossibleImag<T>;
306
307    fn mul(self, rhs: PossibleImag<T>) -> Self::Output {
308        let (mut real, mut imag) = match self {
309            PossibleImag::Complex(r, i) => (r, i),
310            PossibleImag::Imag(i) => (T::ZERO, i),
311            PossibleImag::Real(r) => (r, T::ZERO)
312        };
313
314        match rhs {
315            PossibleImag::Complex(r, i) => {
316                let tmp_real = real * r - imag * i;
317                imag = real * i + imag * r;
318                real = tmp_real;
319            },
320            PossibleImag::Imag(i) => {
321                let tmp_real = - imag * i;
322                imag = real * i;
323                real = tmp_real;
324            },
325            PossibleImag::Real(r) => {
326                real = real * r;
327                imag = imag * r;
328            }
329        }
330
331        match (real, imag) {
332            (real, zero) if zero == T::ZERO => PossibleImag::Real(real),
333            (zero, imag) if zero == T::ZERO => PossibleImag::Imag(imag),
334            _ => PossibleImag::Complex(real, imag)
335        }
336    }
337}
338
339impl<T: Mul<T, Output = T>> Mul<T> for PossibleImag<T> 
340where T: Copy {
341    type Output = PossibleImag<T>;
342
343    fn mul(self, rhs: T) -> Self::Output {
344        match self {
345            PossibleImag::Complex(r, i) => PossibleImag::Complex(r * rhs, i * rhs),
346            PossibleImag::Imag(i) => PossibleImag::Imag(i * rhs),
347            PossibleImag::Real(r) => PossibleImag::Real(r * rhs)
348        }
349    }
350}
351
352impl<T> Div<PossibleImag<T>> for PossibleImag<T> 
353where T: Add<T, Output = T> + Sub<T, Output = T> + Mul<T, Output = T> + Div<T, Output = T> + Neg<Output = T> + Zero + PartialEq + Copy {
354    type Output = PossibleImag<T>;
355
356    fn div(self, rhs: PossibleImag<T>) -> Self::Output {
357        let (mut real, mut imag) = match self {
358            PossibleImag::Complex(r, i) => (r, i),
359            PossibleImag::Imag(i) => (T::ZERO, i),
360            PossibleImag::Real(r) => (r, T::ZERO)
361        };
362
363        match rhs {
364            PossibleImag::Complex(r, i) => {
365                { // new denominator
366                    let tmp_real = real * r +/*-*/ imag * /*-*/i;
367                    imag = imag * r -/*+*/ real * /*-*/i;
368                    real = tmp_real;
369                }
370                let real_denominator = r * r +/*-*/ i * /*-*/i;
371                imag = imag / real_denominator;
372                real = real / real_denominator;
373            },
374            PossibleImag::Imag(i) => {
375                { // new denominator
376                    let tmp_real = imag * i;
377                    imag = real * i;
378                    real = tmp_real;
379                }
380                let real_denominator = - i * i;
381                imag = imag / real_denominator;
382                real = real / real_denominator;
383            },
384            PossibleImag::Real(r) => {
385                real = real / r;
386                imag = imag / r;
387            }
388        }
389
390        match (real, imag) {
391            (real, zero) if zero == T::ZERO => PossibleImag::Real(real),
392            (zero, imag) if zero == T::ZERO => PossibleImag::Imag(imag),
393            _ => PossibleImag::Complex(real, imag)
394        }
395    }
396}
397
398impl<T: Div<T, Output = T>> Div<T> for PossibleImag<T>  
399where T: Copy {
400    type Output = PossibleImag<T>;
401
402    fn div(self, rhs: T) -> Self::Output {
403        match self {
404            PossibleImag::Complex(r, i) => PossibleImag::Complex(r / rhs, i / rhs),
405            PossibleImag::Imag(i) => PossibleImag::Imag(i / rhs),
406            PossibleImag::Real(r) => PossibleImag::Real(r / rhs)
407        }
408    }
409}
410
411impl<T> Neg for PossibleImag<T> 
412where T: Neg<Output = T> {
413    type Output = PossibleImag<T>;
414
415    fn neg(self) -> Self::Output {
416        match self {
417            PossibleImag::Real(r) => PossibleImag::Real(-r),
418            PossibleImag::Imag(i) => PossibleImag::Imag(-i),
419            PossibleImag::Complex(r, i) => PossibleImag::Complex(-r, -i),
420        }
421    }
422}