fix divide function

examples/poly_operations/main.v

@@ -21,7 +21,7 @@ fn main() {
 	result_mult := poly.multiply(poly_1, poly_2)
 	println('Multplication result: ${result_mult}') // Expected: [3.0, 38.0, 40.0, 90.0, 21.0] (3 + 38x + 400x^2 + 90x^3 + 21x^4)
 
-	// Division TODO:Fix this
+	// Division
 	// Result includes the quotient and the remainder
 	// To get the real remainder, divide it by the divisor.
 
@@ -30,8 +30,8 @@ fn main() {
 	// poly_divisor := [-2.0, 1.0] // -2 + x
 
 	// REVERSED WAY
-	poly_dividend := [1.0, -4.0, -4.0, 2.0] // 2 - 4x - 4x^2 + x^3
-	poly_divisor := [1.0, -2.0] // -2 + x
+	poly_dividend := [2.0, -4.0, -4.0, 1.0] // 2 - 4x - 4x^2 + x^3
+	poly_divisor := [-2.0, 1.0] // -2 + x
 
 	quotient, remainder := poly.divide(poly_dividend, poly_divisor)
 	println('Division quotient: ${quotient}') // Expected quotient: [-8.0, -2.0, 1.0] (-8 - 2x + x^2)

poly/poly.v

@@ -215,17 +215,17 @@ pub fn balance_companion_matrix(cm [][]f64) [][]f64 {
 			if col_norm == 0.0 || row_norm == 0.0 {
 				continue
 			}
-			mut g := row_norm / poly.radix
+			mut g := row_norm / radix
 			mut f := 1.0
 			s := col_norm + row_norm
 			for col_norm < g {
-				f *= poly.radix
-				col_norm *= poly.radix2
+				f *= radix
+				col_norm *= radix2
 			}
-			g = row_norm * poly.radix
+			g = row_norm * radix
 			for col_norm > g {
-				f /= poly.radix
-				col_norm /= poly.radix2
+				f /= radix
+				col_norm /= radix2
 			}
 			if (row_norm + col_norm) < 0.95 * s * f {
 				not_converged = true
@@ -290,25 +290,34 @@ pub fn multiply(a []f64, b []f64) []f64 {
 // Output: (q, r) where q is the quotient and r is the remainder
 // such that a(x) = b(x) * q(x) + r(x) and degree(r) < degree(b)
 pub fn divide(a []f64, b []f64) ([]f64, []f64) {
-	mut quotient := []f64{}
+	if b.len == 0 {
+		panic('divisor cannot be an empty polynomial')
+	}
+	if a.len == 0 {
+		return []f64{len: 0}, []f64{len: 0}
+	}
+
+	mut quotient := []f64{len: a.len - b.len + 1, init: 0.0}
 	mut remainder := a.clone()
-	b_lead_coef := b[0]
+
+	b_degree := b.len - 1
+	b_lead_coeff := b[b_degree]
 
 	for remainder.len >= b.len {
-		lead_coef := remainder[0] / b_lead_coef
-		quotient << lead_coef
+		remainder_degree := remainder.len - 1
+		lead_coeff := remainder[remainder_degree]
+
+		quotient_term := lead_coeff / b_lead_coeff
+		quotient_idx := remainder_degree - b_degree
+		quotient[quotient_idx] = quotient_term
+
 		for i in 0 .. b.len {
-			remainder[i] -= lead_coef * b[i]
-		}
-		remainder = unsafe { remainder[1..] }
-		for remainder.len > 0 && math.abs(remainder[0]) < 1e-10 {
-			remainder = unsafe { remainder[1..] }
+			remainder[quotient_idx + i] -= quotient_term * b[i]
 		}
-	}
 
-	if remainder.len == 0 {
-		remainder = []f64{}
+		for remainder.len > 0 && remainder[remainder.len - 1] == 0.0 {
+			remainder = remainder[0..remainder.len - 1].clone()
+		}
 	}
-
 	return quotient, remainder
 }

poly/poly_test.v

@@ -78,10 +78,15 @@ fn test_multiply() {
 }
 
 fn test_divide() {
-	// (x^2 + 2x + 1) / (x + 1) = (x + 1)
-	a := [1.0, 2.0, 1.0]
-	b := [1.0, 1.0]
-	quotient, remainder := divide(a, b)
-	assert quotient == [1.0, 1.0]
-	assert remainder == [] // The empty set indicates that two polynomials divide each other exactly (without remainder).
+	a := [0.0, 0.0, 1.0, 1.0]
+	b := [-2.0, 1.0, 1]
+	mut quotient, mut remainder := divide(a, b)
+	assert quotient == [0.0, 1.0]
+	assert remainder == [0.0, 2.0]
+
+	c := [5.0, -11.0, -7.0, 4.0]
+	d := [5.0, 4.0]
+	quotient, remainder = divide(c, d)
+	assert quotient == [1.0, -3.0, 1]
+	assert remainder == []
 }