feo_math/
polynomial.rs

1use std::{fmt, ops::{Add, Div, Mul, Sub, Neg}};
2
3use crate::{F32Fmt, One, SignOps, Three, Two, imaginary::PossibleImag};
4
5// use super::Equation; todo
6
7use crate::Zero;
8
9#[derive(Clone, Copy)]
10pub struct CubicEquation<T>(pub T, pub T, pub T, pub T); // enough for now got a lot to do still
11
12impl<T> fmt::Debug for CubicEquation<T> 
13where T: fmt::Debug {
14    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
15        write!(f, "{:?}x^3 + {:?}x^2 + {:?}x + {:?}", self.0, self.1, self.2, self.3)
16    }
17}
18
19impl<T> CubicEquation<T>{ // for now 
20    pub fn new(pow_3: T, pow_2: T, var: T, scl: T) -> Self {
21        CubicEquation(pow_3, pow_2, var, scl)
22    }
23    /// returns a tuple. The first element of which is a Vec containing the real solutions
24    /// and the second element of which is a vec containing the imaginary solutions.
25    pub fn solve(&self) -> [PossibleImag<T>; 3]
26    where T: Add<T, Output = T> + Sub<T, Output = T> + Mul<T, Output = T> + Div<T, Output = T> + Neg<Output = T> + SignOps + Zero + One + Two + Three + Copy + F32Fmt + PartialEq + fmt::Debug {
27        // 0 = ax^3 + bx^2 + cx + d
28        // Cubic formula
29
30        let CubicEquation(a, b, c, d) = *self;
31
32        let block0: T = -(b * b * b) / (a * T::THREE * a * T::THREE * a * T::THREE) + (b * c) / (a * a * T::THREE * T::TWO) - d / (a * T::TWO);
33        let block1: T = c / (a * T::THREE) - (b * b) / (a * a * T::THREE * T::THREE);
34        let block2 = PossibleImag::do_sqrt(block0 * block0 + block1 * block1 * block1);
35
36        let part1 = (block2 + block0).cbrt();
37        let part2 = (-block2 + block0).cbrt();
38        
39        let block3 = b / (a * T::THREE);
40        [
41            // part1[0] + part2[0] - block3, 
42            // part1[0] + part2[1] - block3, 
43            part1[0] + part2[2] - block3, // yes // I DONT KNOW WHY THESE ARE THE ANSWERS but they seem to be correct ): I'll be back.
44            // part1[1] + part2[0] - block3,  
45            part1[1] + part2[1] - block3, // yes
46            // part1[1] + part2[2] - block3, 
47            part1[2] + part2[0] - block3, // yes
48            // part1[2] + part2[1] - block3, 
49            // part1[2] + part2[2] - block3,
50        ] // part1 + part2
51    }
52}