poly: improvements, tests and examples

examples/poly_eval/README.md

@@ -0,0 +1,17 @@
+# Example - poly_eval 📘
+
+This example demonstrates the usage of the V Scientific Library for evaluating polynomial at given
+value of x.
+
+## Instructions
+
+1. Ensure you have the V compiler installed. You can download it from [here](https://vlang.io).
+2. Ensure you have the VSL installed. You can do it following the [installation guide](https://github.com/vlang/vsl?tab=readme-ov-file#-installation)!
+3. Navigate to this directory.
+4. Run the example using the following command:
+
+```sh
+v run main.v
+```
+
+Enjoy exploring the capabilities of the V Scientific Library! 😊

examples/poly_eval/main.v

@@ -0,0 +1,16 @@
+module main
+
+import vsl.poly
+
+fn main() {
+	// represent the polynomial 2 * x^2 + 5 * x + 4
+	// with x = 4, result is 2 * 4^2 + 5 * 4 + 4 = 56
+	coef := [4.0, 5, 2]
+	result := poly.eval(coef, 4)
+	println('Evaluated value: ${result}')
+
+	// base polynomial: 2 * x^2 + 5 * x + 4 = 56 with x = 4
+	// 1st derivative: 4 * x + 5 = 21 with x = 4
+	result_derivs := poly.eval_derivs(coef, 4, 2)
+	println('Evaluated values: ${result_derivs}')
+}

examples/poly_matrix/README.md

@@ -0,0 +1,16 @@
+# Example - poly_matrix 📘
+
+This example demonstrates the usage of the V Scientific Library for constructing companion matrix.
+
+## Instructions
+
+1. Ensure you have the V compiler installed. You can download it from [here](https://vlang.io).
+2. Ensure you have the VSL installed. You can do it following the [installation guide](https://github.com/vlang/vsl?tab=readme-ov-file#-installation)!
+3. Navigate to this directory.
+4. Run the example using the following command:
+
+```sh
+v run main.v
+```
+
+Enjoy exploring the capabilities of the V Scientific Library! 😊

examples/poly_matrix/main.v

@@ -0,0 +1,27 @@
+module main
+
+import vsl.poly
+
+fn main() {
+	// Polynomial coefficients for p(x) = 2x^3 -4x^2 + 3x - 5
+	coef := [2.0, -4.0, 3.0, -5.0]
+
+	// For the polynomial p(x) = a_n x^n + a_{n-1} x^{n-1} + ... + a_1 x + a_0
+	// The companion matrix C will be:
+	// [ 0    0    0   -a_0/a_n ]
+	// [ 1    0    0   -a_1/a_n ]
+	// [ 0    1    0   -a_2/a_n ]
+	// [ 0    0    1   -a_3/a_n ]
+	comp_matrix := poly.companion_matrix(coef)
+	println('Companion matrix:')
+	for row in comp_matrix {
+		println(row)
+	}
+
+	// Balancing matrix if needed
+	balanced_matrix := poly.balance_companion_matrix(comp_matrix)
+	println('Balanced companion matrix:')
+	for row in balanced_matrix {
+		println(row)
+	}
+}

examples/poly_operations/README.md

@@ -0,0 +1,17 @@
+# Example - poly_eval 📘
+
+This example demonstrates the usage of the V Scientific Library for additioning, substracting and
+multiplying polynomials.
+
+## Instructions
+
+1. Ensure you have the V compiler installed. You can download it from [here](https://vlang.io).
+2. Ensure you have the VSL installed. You can do it following the [installation guide](https://github.com/vlang/vsl?tab=readme-ov-file#-installation)!
+3. Navigate to this directory.
+4. Run the example using the following command:
+
+```sh
+v run main.v
+```
+
+Enjoy exploring the capabilities of the V Scientific Library! 😊

examples/poly_operations/main.v

@@ -0,0 +1,33 @@
+module main
+
+import vsl.poly
+
+fn main() {
+	// Addition
+	// Degree is not modified unless highest coefficients cancel each other out
+	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}') // 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}') // 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}') // 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

@@ -10,9 +10,11 @@ pub fn eval(c []f64, x f64) f64 {
 		panic('coeficients can not be empty')
 	}
 	len := c.len
-	mut ans := c[len - 1]
-	for e in c[..len - 1] {
-		ans = e + x * ans
+	mut ans := 0.0
+	mut i := len - 1
+	for i >= 0 {
+		ans = c[i] + x * ans
+		i--
 	}
 	return ans
 }
@@ -131,13 +133,13 @@ fn sorted_3_(x_ f64, y_ f64, z_ f64) (f64, f64, f64) {
 	mut y := y_
 	mut z := z_
 	if x > y {
-		y, x = swap_(x, y)
+		x, y = swap_(x, y)
 	}
-	if y > z {
-		z, y = swap_(y, z)
+	if x > z {
+		x, z = swap_(x, z)
 	}
-	if x > y {
-		y, x = swap_(x, y)
+	if y > z {
+		y, z = swap_(y, z)
 	}
 	return x, y, z
 }
@@ -273,7 +275,7 @@ pub fn add(a []f64, b []f64) []f64 {
 	return c
 }
 
-pub fn substract(a []f64, b []f64) []f64 {
+pub fn subtract(a []f64, b []f64) []f64 {
 	na := a.len
 	nb := b.len
 	nc := int(math.max(na, nb))
@@ -304,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

@@ -2,11 +2,26 @@ module poly
 
 fn test_eval() {
 	// ans = 2
-	// ans = 4.0 + 4 * 2 = 12
-	// ans = 5 + 4 * 12 = 53
+	// ans = 5 + 4 * 2 = 13
+	// ans = 4 + 4 * 13 = 56
 	x := 4
-	cof := [4.0, 5, 2]
-	assert eval(cof, 4) == 53
+	coef := [4.0, 5, 2]
+	assert eval(coef, 4) == 56
+}
+
+fn test_eval_derivs() {
+	coef := [4.0, 3.0, 2.0]
+	x := 1.0
+	lenres := 3
+	assert eval_derivs(coef, x, lenres) == [9.0, 7.0, 4.0]
+}
+
+fn test_solve_quadratic() {
+	a := 1.0
+	b := -3.0
+	c := 2.0
+	roots := solve_quadratic(a, b, c)
+	assert roots == [1.0, 2.0]
 }
 
 fn test_swap() {
@@ -16,29 +31,54 @@ fn test_swap() {
 	assert a == 202.0 && b == 101.0
 }
 
-// fn test_sorted_3_(){
-// 	a := 5.0
-// 	b := 7.0
-// 	c := -8.0
-// 	x, y, z := sorted_3_(a, b, c)
-// 	assert y == 7.0
-// }
+fn test_sorted_3_() {
+	a := 5.0
+	b := 7.0
+	c := -8.0
+	x, y, z := sorted_3_(a, b, c)
+	assert x == -8.0 && y == 5.0 && z == 7.0
+}
+
+fn test_companion_matrix() {
+	coef := [2.0, -4.0, 3.0, -5.0]
+	assert companion_matrix(coef) == [
+		[0.0, 0.0, 0.4],
+		[1.0, 0.0, -0.8],
+		[0.0, 1.0, 0.6],
+	]
+}
+
+fn test_balance_companion_matrix() {
+	coef := [2.0, -4.0, 3.0, -5.0]
+	matrix := companion_matrix(coef)
+	assert balance_companion_matrix(matrix) == [
+		[0.0, 0.0, 0.8],
+		[0.5, 0.0, -0.8],
+		[0.0, 1.0, 0.6],
+	]
+}
 
 fn test_add() {
 	a := [6.0, 777, -3]
 	b := [1.0, -755, -4]
 	assert add(a, b) == [7.0, 22, -7]
 }
 
-fn test_substract() {
+fn test_subtract() {
 	a := [6.0, 777, -3]
 	b := [1.0, -755, -4]
-	assert substract(a, b) == [5.0, 1532, 1]
+	assert subtract(a, b) == [5.0, 1532, 1]
 }
 
-// fn test_multiply(){
-// 	a := [9.0, -1, 5]
-// 	b := [0.0, -1, 7]
+fn test_multiply() {
+	a := [9.0, -1, 5]
+	b := [0.0, -1, 7]
+	assert multiply(a, b) == [0.0, -9, 64, -12, 35, 0]
+}
 
-// 	assert multiply(a, b) == [0.0, -9, 73, -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] && remainder == [20.0]
+}