add readme

README.md

@@ -0,0 +1,59 @@
+# Symbolic solver for polynomials with rational coefficients #
+
+## Program output ##
+
+Using `cargo run --release`. On my laptop, the last example times out on debug settings.
+
+	Equation: (x + 50) * 10 - 150 = 100
+	Parsed: (- (- (* (+ x 50) 10) 150) 100)
+	Factorized normal form: 10(25 + 1x)
+	Solutions: { x = -25 }
+
+	Equation: (x - 2) * (x + 2) = 0
+	Parsed: (- (* (- x 2) (+ x 2)) 0)
+	Factorized normal form: (-2 + 1x)(2 + 1x)
+	Solutions: { x = 2, x = -2 }
+
+	Equation: x ^ 2 = 4
+	Parsed: (- (^ x 2) 4)
+	Factorized normal form: (-4 + 0x + 1x^2)
+	Solutions: { x = 2, x = -2 }
+
+	Equation: x ^ 2 - 2 = 0
+	Parsed: (- (- (^ x 2) 2) 0)
+	Factorized normal form: (-2 + 0x + 1x^2)
+	Solutions: { x = 2 ^ (1/2), x = -2 ^ (1/2) }
+
+	Equation: x ^ 2 = 2 * x + 15
+	Parsed: (- (^ x 2) (+ (* 2 x) 15))
+	Factorized normal form: (-15 + -2x + 1x^2)
+	Solutions: { x = 5, x = -3 }
+
+	Equation: (x ^ 2 - 2 * x - 15) * (x + 5) = 0
+	Parsed: (- (* (- (- (^ x 2) (* 2 x)) 15) (+ x 5)) 0)
+	Factorized normal form: (-15 + -2x + 1x^2)(5 + 1x)
+	Solutions: { x = 5, x = -3, x = -5 }
+
+	Equation: x ^ 3 + 3 * x ^ 2 - 25 * x - 75 = 0
+	Parsed: (- (- (- (+ (^ x 3) (* 3 (^ x 2))) (* 25 x)) 75) 0)
+	Factorized normal form: (3 + 1x)(-25 + 0x + 1x^2)
+	Solutions: { x = -3, x = 5, x = -5 }
+
+	Equation: x ^ 3 - 91 * x  - 90 = 0
+	Parsed: (- (- (- (^ x 3) (* 91 x)) 90) 0)
+	Factorized normal form: (-90 + -91x + 0x^2 + 1x^3)
+	Guessing rational solutions: 10, -9, -1
+	Verified guessed solutions!
+	Solutions: { x = 10, x = -9, x = -1 }
+
+## Notes ##
+
+This uses an egraph with rewrite rules representing the usual rules of arithmetic, and constant folding with rational numbers. This alone is enough to solve linear equations, by just starting with the difference of the two sides and minimizing AST size. However, it can't suffice for quadratic equations, as the egraph can't create values "out of this air". For example, to solve `x^2 = 2` the egraph must at least have "seen" the expression `2^(1/2)` once.
+
+So instead, the egraph is used to bring the expression into a standard form: a product of monic polynomials which is factored as much as possible. Since this is difficult to model with a cost function, I wrote a custom extractor (in `normal_form.rs`).
+
+Then we can use the usual quadratic formula to find solutions. The obtained solutions are afterwards simplified by another egraph, with a special "square root of a square integer" rule. This is pretty minimal, a more complete solution would also be able to deal with square roots of rational numbers etc.
+
+Finally, the unfactorized cubic equation from the examples gets factored into a linear and a quadratic factor by the egraph. This is kind of "magic" though (the right factor just seems to randomly come up elsewhere), and for general cubics we can't expect that to happen, even if they have integer solutions. The additional last equation is an example for this. So we use another strategy to deal with cubic equations with rational solutions: they are solved numerically by the general cubic formula, the numerical solutions are then approximated by rational numbers using continued fractions, and finally the result is verified by using an egraph to see if the obtained factization is equivalent to the original polynomial. This required to increase node and iteration limits, and takes a few seconds.
+
+This has essentially only been tested on the 8 examples above. So I would expect there to still be a few bugs. The algorithm is also spectacularly inefficient, considering that the equality of rational polynomials is almost trivial to compute.

src/cubic.rs

@@ -12,7 +12,7 @@ pub fn approximate_rational_cubic(b: &Rational, c: &Rational, d: &Rational, limi
 	numerical_sols.into_iter().map(|x|rational_approx(x, limit)).collect()
 }
 
-// assuming leading coefficient is not 0
+// assuming leading coefficient is not 0, and also that we don't have a triple root
 pub fn solve_cubic_numerically(a: f64, b: f64, c: f64, d: f64) -> Vec<f64> {
 	assert_ne!(a, 0.0);
 
@@ -30,7 +30,7 @@ pub fn solve_cubic_numerically(a: f64, b: f64, c: f64, d: f64) -> Vec<f64> {
 }
 
 fn solve_depressed_cubic_numerically(p: f64, q: f64) -> Vec<f64> {
-	let disc = 4.0 * p * p * p + 27.0 * q * q;
+	let disc = p * p * p / 27.0 + q * q / 4.0;
 
 	if disc < 0.0 {
 		let r = 2.0 * (-p / 3.0).sqrt();
@@ -42,10 +42,14 @@ fn solve_depressed_cubic_numerically(p: f64, q: f64) -> Vec<f64> {
 			 r * (phi + 4.0 * PI / 3.0).cos()]
 	} else {
 		// not implemented at the moment
-		vec![]
+		let u = - q / 2.0 + disc.sqrt();
+		let v = - q / 2.0 - disc.sqrt();
+		vec![u.powf(1.0/3.0) + v.powf(1.0/3.0)]
 	}
 }
 
+// pretty inefficient implementation of continued fraction approximation
+// but doesn't matter here
 pub fn rational_approx(x: f64, limit: u64) -> Rational {
 	if x < 0.0 {
 		let Rational { num, denom } = rational_approx(-x, limit);
@@ -55,7 +59,7 @@ pub fn rational_approx(x: f64, limit: u64) -> Rational {
 	let mut num = 0;
 	let mut denom = 0;
 	for l in 0 .. 10 {
-		let (p,q) = rational_approx_level(x, l);
+		let Some((p,q)) = rational_approx_level(x, l) else { break; };
 
 		if q > limit {
 			break;
@@ -66,17 +70,22 @@ pub fn rational_approx(x: f64, limit: u64) -> Rational {
 
 	Rational{
 		num: num as i64,
-		denom: denom
+		denom
 	}
 }
 
-fn rational_approx_level(x: f64, level: usize) -> (u64, u64) {
+fn rational_approx_level(x: f64, level: usize) -> Option<(u64, u64)> {
+	// we expect very big numbers or infinity if the previous iteration was exact
+	if x > 1e9 {
+		return None;
+	}
+
 	if level == 0 {
-		(x as u64, 1)
+		Some((x as u64, 1))
 	} else {
-		let (p,q) = rational_approx_level(1.0 / (x - x.floor()), level-1);
+		let (p,q) = rational_approx_level(1.0 / (x - x.floor()), level-1)?;
 		let floorx = x as u64;
 
-		(q + floorx * p, p)
+		Some((q + floorx * p, p))
 	}
 }

src/normal_form.rs

@@ -141,7 +141,7 @@ fn extract_polynomial(egraph: &EGraph, stats: &HashMap<Id, PolyStats>, id: Id) -
 
 						assert!(leading_deg >= remainder.len());
 
-						remainder.resize(leading_deg, RATIONAL_ZERO.clone());
+						remainder.resize(leading_deg, RATIONAL_ZERO);
 						remainder.push(leading_coeff);
 						return Some(remainder);
 					}