poly: add divide function, example and tests

examples/poly_operations/main.v

@@ -5,19 +5,29 @@ import vsl.poly
 fn main() {
 	// Addition
 	// Degree is not modified unless highest coefficients cancel each other out
-	poly_1 := [1.0, 12, 3]
-	poly_2 := [3.0, 2, 7]
+	poly_1 := [1.0, 12, 3] // 1 + 12x + 3x^2
+	poly_2 := [3.0, 2, 7] // 3 + 2x + 7x^2
 	result_add := poly.add(poly_1, poly_2)
-	println('Addition result: ${result_add}')
+	println('Addition result: ${result_add}') // Expected: [4.0, 14.0, 10.0] (4 + 14x + 10x^2)
 
 	// Subtraction
 	// Degree is not modified unless highest coefficients cancel each other out
 	result_sub := poly.subtract(poly_1, poly_2)
-	println('Subtraction result: ${result_sub}')
+	println('Subtraction result: ${result_sub}') // Expected: [-2.0, 10.0, -4.0] (-2 + 10x - 4x^2)
 
 	// Multiplication
 	// with given degree n and m for poly_1 and poly_2
 	// resulting polynomial will be of degree n + m
 	result_mult := poly.multiply(poly_1, poly_2)
-	println('Multplication result: ${result_mult}')
+	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
+	// Result includes the quotient and the remainder
+	// To get the real remainder, divide it by the divisor.
+	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)
+	println('Division remainder: ${remainder}') // Expected remainder: [-14.0]
+	// Real remainder: -14 / (-2 + x)
 }

poly/poly.v

@@ -306,3 +306,36 @@ pub fn multiply(a []f64, b []f64) []f64 {
 	}
 	return c
 }
+
+pub fn divide(dividend []f64, divisor []f64) ([]f64, []f64) {
+	if divisor.len == 0 {
+		panic('divisor cannot be an empty polynomial')
+	}
+	if dividend.len == 0 {
+		return []f64{len: 0}, []f64{len: 0}
+	}
+
+	mut quotient := []f64{len: dividend.len - divisor.len + 1, init: 0.0}
+	mut remainder := dividend.clone()
+
+	divisor_degree := divisor.len - 1
+	divisor_lead_coeff := divisor[divisor_degree]
+
+	for remainder.len >= divisor.len {
+		remainder_degree := remainder.len - 1
+		lead_coeff := remainder[remainder_degree]
+
+		quotient_term := lead_coeff / divisor_lead_coeff
+		quotient_idx := remainder_degree - divisor_degree
+		quotient[quotient_idx] = quotient_term
+
+		for i in 0 .. divisor.len {
+			remainder[quotient_idx + i] -= quotient_term * divisor[i]
+		}
+
+		for remainder.len > 0 && remainder[remainder.len - 1] == 0.0 {
+			remainder = remainder[0..remainder.len - 1].clone()
+		}
+	}
+	return quotient, remainder
+}

poly/poly_test.v

@@ -21,7 +21,9 @@ fn test_sorted_3_() {
 	b := 7.0
 	c := -8.0
 	x, y, z := sorted_3_(a, b, c)
-	assert x == -8.0 && y == 5.0 && z == 7.0
+	assert x == -8.0
+	assert y == 5.0
+	assert z == 7.0
 }
 
 fn test_add() {
@@ -39,6 +41,13 @@ fn test_subtract() {
 fn test_multiply() {
 	a := [9.0, -1, 5]
 	b := [0.0, -1, 7]
-
 	assert multiply(a, b) == [0.0, -9, 64, -12, 35, 0]
 }
+
+fn test_division() {
+	dividend := [6.0, -5.0, 1.0]
+	divisor := [2.0, 1.0]
+	quotient, remainder := divide(dividend, divisor)
+	assert quotient == [-7.0, 1]
+	assert remainder == [20.0]
+}