computorv1/src/maths/solver.rs

use super::*;

fn degree_zero(equation: Vec<GaussianRational>) {
    println!("Polynomial degree: 0");
    if equation[0] == zero() {
        println!("Every real number is a solution ! Personally, I'd rather say that this equation is valid.")
    } else {
        println!("This equation has no solution ! Personally, I'd rather say that this equation is invalid.")
    }
}

fn degree_one(equation: Vec<GaussianRational>) {
    println!("Polynomial degree: 1");
    let x = minus_one() * equation[0] / equation[1];
    println!("The solution is");
    let mut format = x.format(false);
    if format == "" {
        format = String::from("0");
    }
    println!("x = {}", format);
}

#[derive(Debug)]
struct MySqrt {
    numerator_natural: i128,
    numerator_irrational: i128,
    denominator: i128,
}

#[derive(Debug)]
struct MyCompSqrt {
    sqrt: MySqrt,
    rational: Rational,
    sign: i128,
}

fn simplify_sqrt(n: i128) -> (i128, i128) {
    if n == 0 {
        return (0, 0);
    }
    let mut prime_factors = get_prime_factors(n);
    let (mut natural, mut irrational) = (1, 1);
    prime_factors = prime_factors.into_iter().rev().collect();
    while prime_factors.len() > 1 {
        let pop = prime_factors.pop().unwrap();
        if *prime_factors.last().unwrap() == pop {
            natural *= pop;
            prime_factors.pop();
        } else {
            irrational *= pop;
        }
    }
    if let Some(n) = prime_factors.pop() {
        irrational *= n;
    }
    (natural, irrational)
}

// Eric Naslund for the math https://math.stackexchange.com/questions/44406/how-do-i-get-the-square-root-of-a-complex-number
fn comp_sqrt(c: GaussianRational) -> MyCompSqrt {
    let real = c.real();
    let imag = c.imaginary();
    let n = real * real + imag * imag;
    let numerator_irrational = n.numerator * n.denominator;
    let (mut numerator_natural, numerator_irrational) = simplify_sqrt(numerator_irrational);

    let mut denominator = n.denominator;
    let gcd_var = gcd(numerator_natural, denominator);
    denominator /= gcd_var;
    numerator_natural /= gcd_var;
    let mut sqrt = MySqrt {
        numerator_natural,
        numerator_irrational,
        denominator,
    };

    let mut sign = 1;
    if c.imaginary.numerator < 0 {
        sign = -1;
    }

    let mut rational = c.real;
    rational = rational / 2;
    if sqrt.numerator_natural % 2 == 0 {
        sqrt.numerator_natural /= 2;
    } else {
        sqrt.denominator *= 2;
    }

    // Here we multiply by 4 to avoid dividing by 2 later and ease ascii art
    let mem = sqrt.denominator;
    sqrt.denominator *= rational.denominator * 4;
    sqrt.numerator_natural *= rational.denominator;
    rational.numerator *= mem;
    rational.denominator *= mem * 4;

    let gcd = gcd(gcd(sqrt.denominator, sqrt.numerator_natural), rational.numerator);

    sqrt.denominator /= gcd;
    sqrt.numerator_natural /= gcd;
    rational.numerator /= gcd;
    rational.denominator /= gcd;

    MyCompSqrt {
        sqrt,
        rational,
        sign,
    }
}

fn get_strings(left_part: Rational, sqrt_delta: &MySqrt) -> (String, String, String, String, String, String) {
    let mut a = left_part.numerator;
    let mut b = sqrt_delta.numerator_natural;
    if b < 0 {
        b = -b;
    }
    let c = sqrt_delta.numerator_irrational;
    let mut d = left_part.denominator;
    if d < 0 {
        a = -a;
        d = -d;
    }

    let mut string = String::from("");
    if b != 0 || c != 0 {
        if b != 1 && c != 1 {
            string = format!("{b} * sqrt({c})");
        } else if b != 1 {
            string = format!("{b}");
        } else if c != 1 {
            string = format!("sqrt({c})");
        } else {
            string = format!("1");
        }
    }

    let mut string_a = format!("{string}");
    let mut string_b = format!("{string}");
    if a != 0 {
        string_a = format!("{a} + {string}");
        string_b = format!("{} + {string}", -a);
    }


    let len = std::cmp::max(string_a.len(), format!("{d}").len()) + 2;
    let len_b = std::cmp::max(string_b.len(), format!("{d}").len()) + 2;

    let len1 = len - string_a.len();
    let len2 = len - d.to_string().len();
    let len1b = len_b - string_b.len();
    let len2b = len_b - d.to_string().len();

    let string1;
    let string2;
    let string3;
    let string4;
    let string5;
    let string6;

    if d != 1 {
        string1 = format!("{}{string_a}{}", " ".repeat(len1 / 2), " ".repeat(len1 / 2 + len1 % 2));
        string2 = format!("{}", "-".repeat(len));
        string3 = format!("{}{d}{}", " ".repeat(len2 / 2), " ".repeat(len2 / 2 + len2 % 2));
        string4 = format!("{}{string_b}{}", " ".repeat(len1b / 2), " ".repeat(len1b / 2 + len1b % 2));
        string5 = format!("{}", "-".repeat(len_b));
        string6 = format!("{}{d}{}", " ".repeat(len2b / 2), " ".repeat(len2b / 2 + len2b % 2));
    } else if string.len() > 0 {
        string1 = format!("{}", string.len());
        string2 = format!("{string_a}");
        string3 = format!("{}", string.len());
        string4 = format!("{}", string.len());
        string5 = format!("{string_b}");
        string6 = format!("{}", string.len());
    } else {
        string1 = format!(" ");
        string2 = format!("0");
        string3 = format!(" ");
        string4 = format!(" ");
        string5 = format!("0");
        string6 = format!(" ");
    }

    (string1, string2, string3, string4, string5, string6)
}

fn degree_two_complex(a: GaussianRational, b: GaussianRational, c: GaussianRational) {
    // Here we divide by a to avoid dividing by complex later and ease ascii art
    let b = b / a;
    let c = c / a;
    let delta = b * b - 4 * c;
    let sqrt_delta = comp_sqrt(delta);

    let b = b / GaussianRational::new(Rational::new(-2, 1), Rational::new(0, 1));

    let mut sign = '±';
    if sqrt_delta.sign < 0 {
        sign = '∓';
    }

    let (s1, s2, s3, s4, s5, s6) = get_strings(sqrt_delta.rational, &sqrt_delta.sqrt);
    let mut string1 = String::from("");
    let mut string2 = String::from("");
    let mut string3 = String::from("");

    if b.real.numerator != 0 {
        if b.real.denominator != 1 {
            let b_real_len = std::cmp::max(format!("{}", b.real.denominator).len(), format!("{}", b.real.numerator).len()) + 2;
            let space_num = b_real_len - format!("{}", b.real.numerator).len();
            let space_dem = b_real_len - format!("{}", b.real.denominator).len();
            string1 = format!("{}{}{}", " ".repeat(space_num / 2), b.real.numerator, " ".repeat(space_num / 2 + space_num % 2));
            string2 = format!("{}", "-".repeat(b_real_len));
            string3 = format!("{}{}{}", " ".repeat(space_dem / 2), b.real.denominator, " ".repeat(space_dem / 2 + space_dem % 2));
        } else {
            let string = format!("{}", b.real.numerator);
            string1 = format!("{}", " ".repeat(string.len()));
            string2 = format!("{string}");
            string3 = format!("{}", " ".repeat(string.len()));
        }
    }

    string1 = format!("{string1}         /{s1}\\  ");
    string2 = format!("{string2} ± sqrt | {s2} | ");
    string3 = format!("{string3}         \\{s3}/  ");

    if b.imaginary.numerator != 0 {
        if b.imaginary.denominator != 1 {
            let b_imag_len = std::cmp::max(format!("{}", b.imaginary.denominator).len(), format!("{}", b.imaginary.numerator).len()) + 2;
            let space_num = b_imag_len - format!("{}", b.imaginary.numerator).len();
            let space_dem = b_imag_len - format!("{}", b.imaginary.denominator).len();
            string1 = format!("{string1}  {}{}{}   ", " ".repeat(space_num / 2), b.imaginary.numerator, " ".repeat(space_num / 2 + space_num % 2));
            string2 = format!("{string2}+ {} i ", "-".repeat(b_imag_len));
            string3 = format!("{string3}  {}{}{}   ", " ".repeat(space_dem / 2), b.imaginary.denominator, " ".repeat(space_dem / 2 + space_dem % 2));
        } else {
            let string = format!("{}", b.imaginary.numerator);
            string1 = format!("{string1}  {}  ", " ".repeat(string.len()));
            string2 = format!("{string2}+ {string}i ");
            string3 = format!("{string3}  {}  ", " ".repeat(string.len()));
        }
    }

    string1 = format!("{string1}        /{s4}\\   ");
    string2 = format!("{string2}{sign} sqrt | {s5} | i");
    string3 = format!("{string3}        \\{s6}/   ");

    println!("The two solutions are:");
    println!("{string1}");
    println!("{string2}");
    println!("{string3}");
    println!("Approximations:");
    let sqrt = sqrt_delta.sqrt.numerator_natural as f64 * approx_sqrt(sqrt_delta.sqrt.numerator_irrational);
    let left = b.real.numerator as f64 / b.real.denominator as f64;
    let left_sqrt = approx_sqrt_f64((sqrt_delta.rational.numerator as f64 + sqrt) / sqrt_delta.rational.denominator as f64);

    let right = b.imaginary.numerator as f64 / b.imaginary.denominator as f64;
    let right_sqrt = sqrt_delta.sign as f64 * approx_sqrt_f64((-sqrt_delta.rational.numerator as f64 + sqrt) / sqrt_delta.rational.denominator as f64);
    println!("{}", format!("{} + {}i", left + left_sqrt, right + right_sqrt).split("+ -").collect::<Vec<_>>().join("- "));
    println!("{}", format!("{} + {}i", left - left_sqrt, right - right_sqrt).split("+ -").collect::<Vec<_>>().join("- "));
}

fn sqrt(n: Rational, a: Rational) -> MySqrt {
    let numerator_irrational = n.numerator * n.denominator;
    let (mut numerator_natural, numerator_irrational) = simplify_sqrt(numerator_irrational);

    numerator_natural *= a.denominator;
    let mut denominator = n.denominator * 2 * a.numerator;
    let gcd = gcd(numerator_natural, denominator);
    denominator /= gcd;
    numerator_natural /= gcd;
    MySqrt {
        numerator_natural,
        numerator_irrational,
        denominator,
    }
}

fn print_degree_two_real(left_part: Rational, sqrt_delta: MySqrt, imaginary: bool) {
    let mut a = left_part.numerator;
    let mut b = sqrt_delta.numerator_natural;
    if b < 0 {
        b = -b;
    }
    let c = sqrt_delta.numerator_irrational;
    let mut d = left_part.denominator;
    if d < 0 {
        a = -a;
        d = -d;
    }

    let mut string = String::from("");
    if a != 0 {
        string = format!("{a} ");
    }

    if b != 0 || c != 0 {
        if imaginary == true {
            if b != 1 && c != 1 {
                string = format!("{string}± ({b} * sqrt({c}))i");
            } else if b != 1 {
                string = format!("{string}± {b}i");
            } else if c != 1 {
                string = format!("{string}± sqrt({c})i");
            } else {
                string = format!("{string}± i");
            }
        } else {
            if b != 1 && c != 1 {
                string = format!("{string}± {b} * sqrt({c})");
            } else if b != 1 {
                string = format!("{string}± {b}");
            } else if c != 1 {
                string = format!("{string}± sqrt({c})");
            } else {
                string = format!("{string}± 1");
            }
        }
    }


    let len = std::cmp::max(string.chars().count(), format!("{d}").len()) + 2;

    if b != 0 || c != 0 {
        if imaginary == true {
            println!("Discriminant is less than zero. The two complex solutions are:");
        } else {
            println!("Discriminant is more than zero. The two solutions are:");
        }
    } else {
        println!("Discriminant is zero. The solution is:");
    }
    if d != 1 {
        println!("{}{string}", " ".repeat(4 + (len - string.chars().count()) / 2));
        println!("x = {}", "-".repeat(len));
        println!("{}{d}", " ".repeat(4 + (len - d.to_string().len()) / 2));
        let mut right = (sqrt_delta.numerator_natural as f64 * approx_sqrt(sqrt_delta.numerator_irrational)) / left_part.denominator as f64;
        let left = left_part.numerator as f64 / left_part.denominator as f64;
        if right < 0. {
            right *= -1.;
        }
        println!("Approximations:");
        if imaginary == true {
            println!("{} - {}i", left, right);
            println!("{} + {}i", left, right);
        } else {
            println!("{}", left - right);
            println!("{}", left + right);

        }
    } else if string.len() > 0 {
        println!("x = {string}");
    } else {
        println!("x = 0");
    }
}

fn degree_two_real(a: Rational, b: Rational, c: Rational) {
    let mut delta = b * b - 4 * a * c;
    let mut imaginary = false;
    if delta.numerator < 0 {
        imaginary = true;
        delta = -1 * delta;
    }
    let mut sqrt_delta = sqrt(delta, a);
    let mut left_part = (-1 * b) / (2 * a);

    let tmp = left_part.denominator;

    left_part.denominator *= sqrt_delta.denominator;
    left_part.numerator *= sqrt_delta.denominator;
    sqrt_delta.numerator_natural *= tmp;
    sqrt_delta.denominator *= tmp;

    print_degree_two_real(left_part, sqrt_delta, imaginary);
}

fn degree_two(equation: Vec<GaussianRational>) {
    println!("Polynomial degree: 2");
    let (a, b, c) = (equation[2], equation[1], equation[0]);
    if a.is_real() && b.is_real() && c.is_real() {
        degree_two_real(a.real(), b.real(), c.real());
    } else {
        degree_two_complex(a, b, c);
    }
}

pub fn solve(mut equation: Vec<GaussianRational>) -> Vec<GaussianRational> {
    for i in (1..equation.len()).rev() {
        if equation[i] == zero() {
            equation.pop();
        } else {
            break;
        }
    }
    match equation.len() {
        0 => unreachable!(),
        1 => degree_zero(equation),
        2 => degree_one(equation),
        3 => degree_two(equation),
        _ => println!("Polynomial of degree greater than 2 detected, I can't solve that !"),
    }
    vec![]
}