From 4c41c7518457de1041aef0225e722085180c6544 Mon Sep 17 00:00:00 2001
From: PottierLoic <loic.pottier@live.fr>
Date: Mon, 15 Jul 2024 22:28:04 +0200
Subject: [PATCH 01/11] poly: fix eval and sorted_3_ + fix tests

---
 poly/poly.v      | 18 ++++++++++--------
 poly/poly_test.v | 31 ++++++++++++++++---------------
 2 files changed, 26 insertions(+), 23 deletions(-)

diff --git a/poly/poly.v b/poly/poly.v
index 30e842e2f..ddd87d4e6 100644
--- a/poly/poly.v
+++ b/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
 }
diff --git a/poly/poly_test.v b/poly/poly_test.v
index 8db4f0e1d..a644e9af1 100644
--- a/poly/poly_test.v
+++ b/poly/poly_test.v
@@ -2,11 +2,11 @@ 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
+	assert eval(cof, 4) == 56
 }
 
 fn test_swap() {
@@ -16,13 +16,13 @@ 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 y == 5.0
+}
 
 fn test_add() {
 	a := [6.0, 777, -3]
@@ -36,9 +36,10 @@ fn test_substract() {
 	assert substract(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]
-// }

From e73652674b0ffa334094683a9b6f6e1ad1d2d101 Mon Sep 17 00:00:00 2001
From: PottierLoic <loic.pottier@live.fr>
Date: Mon, 15 Jul 2024 22:44:24 +0200
Subject: [PATCH 02/11] poly: fmt + complete test

---
 poly/poly.v      |  6 +++---
 poly/poly_test.v | 13 ++++++-------
 2 files changed, 9 insertions(+), 10 deletions(-)

diff --git a/poly/poly.v b/poly/poly.v
index ddd87d4e6..d53f8e5fd 100644
--- a/poly/poly.v
+++ b/poly/poly.v
@@ -11,8 +11,8 @@ pub fn eval(c []f64, x f64) f64 {
 	}
 	len := c.len
 	mut ans := 0.0
-	mut i := len-1
-	for i >= 0{
+	mut i := len - 1
+	for i >= 0 {
 		ans = c[i] + x * ans
 		i--
 	}
@@ -136,7 +136,7 @@ fn sorted_3_(x_ f64, y_ f64, z_ f64) (f64, f64, f64) {
 		x, y = swap_(x, y)
 	}
 	if x > z {
-		 x, z = swap_(x, z)
+		x, z = swap_(x, z)
 	}
 	if y > z {
 		y, z = swap_(y, z)
diff --git a/poly/poly_test.v b/poly/poly_test.v
index a644e9af1..068edda4c 100644
--- a/poly/poly_test.v
+++ b/poly/poly_test.v
@@ -16,12 +16,12 @@ fn test_swap() {
 	assert a == 202.0 && b == 101.0
 }
 
-fn test_sorted_3_(){
+fn test_sorted_3_() {
 	a := 5.0
 	b := 7.0
-  c := -8.0
-  x, y, z := sorted_3_(a, b, c)
-  assert y == 5.0
+	c := -8.0
+	x, y, z := sorted_3_(a, b, c)
+	assert x == -8.0 && y == 5.0 && z == 7.0
 }
 
 fn test_add() {
@@ -36,10 +36,9 @@ fn test_substract() {
 	assert substract(a, b) == [5.0, 1532, 1]
 }
 
-fn test_multiply(){
-  a := [9.0, -1, 5]
+fn test_multiply() {
+	a := [9.0, -1, 5]
 	b := [0.0, -1, 7]
 
 	assert multiply(a, b) == [0.0, -9, 64, -12, 35, 0]
 }
-

From 6e57a00b5dfdd599b19a41e6b7a52a6449059e56 Mon Sep 17 00:00:00 2001
From: PottierLoic <loic.pottier@live.fr>
Date: Mon, 15 Jul 2024 23:34:59 +0200
Subject: [PATCH 03/11] poly: add examples for eval and basic operations

---
 examples/poly_eval/README.md       | 17 +++++++++++++++++
 examples/poly_eval/main.v          | 11 +++++++++++
 examples/poly_operations/README.md | 17 +++++++++++++++++
 examples/poly_operations/main.v    | 23 +++++++++++++++++++++++
 4 files changed, 68 insertions(+)
 create mode 100644 examples/poly_eval/README.md
 create mode 100644 examples/poly_eval/main.v
 create mode 100644 examples/poly_operations/README.md
 create mode 100644 examples/poly_operations/main.v

diff --git a/examples/poly_eval/README.md b/examples/poly_eval/README.md
new file mode 100644
index 000000000..f4162db0b
--- /dev/null
+++ b/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! 😊
diff --git a/examples/poly_eval/main.v b/examples/poly_eval/main.v
new file mode 100644
index 000000000..546f11179
--- /dev/null
+++ b/examples/poly_eval/main.v
@@ -0,0 +1,11 @@
+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('Evalutated value: ${result}')
+}
diff --git a/examples/poly_operations/README.md b/examples/poly_operations/README.md
new file mode 100644
index 000000000..79e818955
--- /dev/null
+++ b/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! 😊
diff --git a/examples/poly_operations/main.v b/examples/poly_operations/main.v
new file mode 100644
index 000000000..adf818861
--- /dev/null
+++ b/examples/poly_operations/main.v
@@ -0,0 +1,23 @@
+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]
+	poly_2 := [3.0, 2, 7]
+	result_add := poly.add(poly_1, poly_2)
+	println('Addition result: ${result_add}')
+
+	// Substraction
+	// Degree is not modified unless highest coefficients cancel each other out
+	result_sub := poly.substract(poly_1, poly_2)
+	println('Substraction result: ${result_sub}')
+
+	// 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}')
+}

From cc2d76c0f1c3f7c5ae626923fb74955e2697059e Mon Sep 17 00:00:00 2001
From: PottierLoic <loic.pottier@live.fr>
Date: Mon, 15 Jul 2024 23:45:28 +0200
Subject: [PATCH 04/11] poly: fix typo

---
 examples/poly_operations/main.v | 6 +++---
 poly/poly.v                     | 2 +-
 poly/poly_test.v                | 4 ++--
 3 files changed, 6 insertions(+), 6 deletions(-)

diff --git a/examples/poly_operations/main.v b/examples/poly_operations/main.v
index adf818861..c7ebb399d 100644
--- a/examples/poly_operations/main.v
+++ b/examples/poly_operations/main.v
@@ -10,10 +10,10 @@ fn main() {
 	result_add := poly.add(poly_1, poly_2)
 	println('Addition result: ${result_add}')
 
-	// Substraction
+	// Subtraction
 	// Degree is not modified unless highest coefficients cancel each other out
-	result_sub := poly.substract(poly_1, poly_2)
-	println('Substraction result: ${result_sub}')
+	result_sub := poly.subtract(poly_1, poly_2)
+	println('Subtraction result: ${result_sub}')
 
 	// Multiplication
 	// with given degree n and m for poly_1 and poly_2
diff --git a/poly/poly.v b/poly/poly.v
index d53f8e5fd..e050a4290 100644
--- a/poly/poly.v
+++ b/poly/poly.v
@@ -275,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))
diff --git a/poly/poly_test.v b/poly/poly_test.v
index 068edda4c..0d33c8ce1 100644
--- a/poly/poly_test.v
+++ b/poly/poly_test.v
@@ -30,10 +30,10 @@ fn test_add() {
 	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() {

From 9277af11f76f372d69674a181891c5fec09696a6 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Lo=C3=AFc=20Pottier?= <loic.pottier@live.fr>
Date: Tue, 16 Jul 2024 15:47:46 +0200
Subject: [PATCH 05/11] poly: add divide function, example and tests

---
 examples/poly_operations/main.v | 20 +++++++++++++++-----
 poly/poly.v                     | 33 +++++++++++++++++++++++++++++++++
 poly/poly_test.v                | 13 +++++++++++--
 3 files changed, 59 insertions(+), 7 deletions(-)

diff --git a/examples/poly_operations/main.v b/examples/poly_operations/main.v
index c7ebb399d..cbba520e5 100644
--- a/examples/poly_operations/main.v
+++ b/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)
 }
diff --git a/poly/poly.v b/poly/poly.v
index e050a4290..31fcce2e3 100644
--- a/poly/poly.v
+++ b/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
+}
diff --git a/poly/poly_test.v b/poly/poly_test.v
index 0d33c8ce1..99d4cf566 100644
--- a/poly/poly_test.v
+++ b/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]
+}

From 969539365bfdd2fe08a24de0d6c1cfc348d8f18f Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Lo=C3=AFc=20Pottier?= <loic.pottier@live.fr>
Date: Fri, 19 Jul 2024 16:31:10 +0200
Subject: [PATCH 06/11] poly: add tests and examples

---
 examples/poly_eval/main.v |  7 ++++++-
 poly/poly_test.v          | 22 +++++++++++++++++-----
 2 files changed, 23 insertions(+), 6 deletions(-)

diff --git a/examples/poly_eval/main.v b/examples/poly_eval/main.v
index 546f11179..a82a51205 100644
--- a/examples/poly_eval/main.v
+++ b/examples/poly_eval/main.v
@@ -7,5 +7,10 @@ fn main() {
 	// with x = 4, result is 2 * 4^2 + 5 * 4 + 4 = 56
 	coef := [4.0, 5, 2]
 	result := poly.eval(coef, 4)
-	println('Evalutated value: ${result}')
+	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}')
 }
diff --git a/poly/poly_test.v b/poly/poly_test.v
index 99d4cf566..13033b166 100644
--- a/poly/poly_test.v
+++ b/poly/poly_test.v
@@ -9,6 +9,21 @@ fn test_eval() {
 	assert eval(cof, 4) == 56
 }
 
+fn test_eval_derivs() {
+	cof := [4.0, 3.0, 2.0]
+	x := 1.0
+	lenres := 3
+	assert eval_derivs(cof, 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() {
 	mut a := 101.0
 	mut b := 202.0
@@ -21,9 +36,7 @@ fn test_sorted_3_() {
 	b := 7.0
 	c := -8.0
 	x, y, z := sorted_3_(a, b, c)
-	assert x == -8.0
-	assert y == 5.0
-	assert z == 7.0
+	assert x == -8.0 && y == 5.0 && z == 7.0
 }
 
 fn test_add() {
@@ -48,6 +61,5 @@ 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]
+	assert quotient == [-7.0, 1] && remainder == [20.0]
 }

From 8c08e420f123b74abfde93a92a3769573d728f0b Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Lo=C3=AFc=20Pottier?= <loic.pottier@live.fr>
Date: Tue, 23 Jul 2024 11:30:22 +0200
Subject: [PATCH 07/11] noise: add tests/examples for companion matrix

---
 examples/poly_matrix/README.md | 16 ++++++++++++++++
 examples/poly_matrix/main.v    | 27 +++++++++++++++++++++++++++
 poly/poly_test.v               | 27 +++++++++++++++++++++++----
 3 files changed, 66 insertions(+), 4 deletions(-)
 create mode 100644 examples/poly_matrix/README.md
 create mode 100644 examples/poly_matrix/main.v

diff --git a/examples/poly_matrix/README.md b/examples/poly_matrix/README.md
new file mode 100644
index 000000000..76b0e48d9
--- /dev/null
+++ b/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! 😊
diff --git a/examples/poly_matrix/main.v b/examples/poly_matrix/main.v
new file mode 100644
index 000000000..fb4dd9818
--- /dev/null
+++ b/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)
+	}
+}
diff --git a/poly/poly_test.v b/poly/poly_test.v
index 13033b166..1bde71a5c 100644
--- a/poly/poly_test.v
+++ b/poly/poly_test.v
@@ -5,15 +5,15 @@ fn test_eval() {
 	// ans = 5 + 4 * 2 = 13
 	// ans = 4 + 4 * 13 = 56
 	x := 4
-	cof := [4.0, 5, 2]
-	assert eval(cof, 4) == 56
+	coef := [4.0, 5, 2]
+	assert eval(coef, 4) == 56
 }
 
 fn test_eval_derivs() {
-	cof := [4.0, 3.0, 2.0]
+	coef := [4.0, 3.0, 2.0]
 	x := 1.0
 	lenres := 3
-	assert eval_derivs(cof, x, lenres) == [9.0, 7.0, 4.0]
+	assert eval_derivs(coef, x, lenres) == [9.0, 7.0, 4.0]
 }
 
 fn test_solve_quadratic() {
@@ -39,6 +39,25 @@ fn test_sorted_3_() {
 	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]

From eae3bf6e035e9aa14e57e21def4a9b4d33ad0cde Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Lo=C3=AFc=20Pottier?= <loic.pottier@live.fr>
Date: Tue, 23 Jul 2024 11:31:37 +0200
Subject: [PATCH 08/11] noise: fmt

---
 examples/poly_matrix/main.v | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/examples/poly_matrix/main.v b/examples/poly_matrix/main.v
index fb4dd9818..23a966e14 100644
--- a/examples/poly_matrix/main.v
+++ b/examples/poly_matrix/main.v
@@ -3,24 +3,24 @@ module main
 import vsl.poly
 
 fn main() {
-  // Polynomial coefficients for p(x) = 2x^3 -4x^2 + 3x - 5
+	// 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
+	// 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:")
+	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:")
+	balanced_matrix := poly.balance_companion_matrix(comp_matrix)
+	println('Balanced companion matrix:')
 	for row in balanced_matrix {
 		println(row)
 	}

From 2f05eb11a6b0a1617ca4cb9a2d1066d9ebd75ce6 Mon Sep 17 00:00:00 2001
From: PottierLoic <loic.pottier@live.fr>
Date: Tue, 10 Sep 2024 19:23:06 +0200
Subject: [PATCH 09/11] fix(poly): merge @suleyman-kaya changes

Merging with @suleyman-kaya changes
removed println and useless math import in poly_test.v
---
 examples/poly_operations/main.v |  2 +-
 poly/poly.v                     | 33 ---------------------------------
 poly/poly_test.v                |  7 -------
 3 files changed, 1 insertion(+), 41 deletions(-)

diff --git a/examples/poly_operations/main.v b/examples/poly_operations/main.v
index cbba520e5..47b9fb6e4 100644
--- a/examples/poly_operations/main.v
+++ b/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
+	// Division TODO:Fix this
 	// 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
diff --git a/poly/poly.v b/poly/poly.v
index 697c8de86..da10bf14c 100644
--- a/poly/poly.v
+++ b/poly/poly.v
@@ -312,36 +312,3 @@ pub fn divide(a []f64, b []f64) ([]f64, []f64) {
 
 	return quotient, remainder
 }
-
-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
-}
diff --git a/poly/poly_test.v b/poly/poly_test.v
index 13573b3c0..d28a2fb40 100644
--- a/poly/poly_test.v
+++ b/poly/poly_test.v
@@ -1,7 +1,5 @@
 module poly
 
-import math
-
 fn test_eval() {
 	x := 4
 	coef := [4.0, 5, 2]
@@ -61,7 +59,6 @@ fn test_add() {
 	a := [1.0, 2.0, 3.0]
 	b := [6.0, 20.0, -10.0]
 	result := add(a, b)
-	println('Add result: ${result}')
 	assert result == [7.0, 22.0, -7.0]
 }
 
@@ -69,7 +66,6 @@ fn test_subtract() {
 	a := [6.0, 1532.0, -4.0]
 	b := [1.0, -1.0, -5.0]
 	result := subtract(a, b)
-	println('Subtract result: ${result}')
 	assert result == [5.0, 1533.0, 1.0]
 }
 
@@ -78,7 +74,6 @@ fn test_multiply() {
 	a := [2.0, 3.0, 4.0]
 	b := [0.0, -3.0, 2.0]
 	result := multiply(a, b)
-	println('Multiply result: ${result}')
 	assert result == [0.0, -6.0, -5.0, -6.0, 8.0]
 }
 
@@ -87,8 +82,6 @@ fn test_divide() {
 	a := [1.0, 2.0, 1.0]
 	b := [1.0, 1.0]
 	quotient, remainder := divide(a, b)
-	println('Divide quotient: ${quotient}')
-	println('Divide remainder: ${remainder}')
 	assert quotient == [1.0, 1.0]
 	assert remainder == [] // The empty set indicates that two polynomials divide each other exactly (without remainder).
 }

From d0b4fd567f51634a8863e3d2bac5a13516c2f103 Mon Sep 17 00:00:00 2001
From: PottierLoic <loic.pottier@live.fr>
Date: Tue, 10 Sep 2024 19:40:58 +0200
Subject: [PATCH 10/11] run fmt

---
 examples/poly_operations/main.v | 11 +++++++++--
 1 file changed, 9 insertions(+), 2 deletions(-)

diff --git a/examples/poly_operations/main.v b/examples/poly_operations/main.v
index 47b9fb6e4..dd3d552b2 100644
--- a/examples/poly_operations/main.v
+++ b/examples/poly_operations/main.v
@@ -24,8 +24,15 @@ fn main() {
 	// Division TODO:Fix this
 	// 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
+
+	// OLD WAY
+	// poly_dividend := [2.0, -4.0, -4.0, 1.0] // 2 - 4x - 4x^2 + x^3
+	// 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
+
 	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]

From 6ebfb97404f308e6c554d5eb6e91270720f4dd25 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Lo=C3=AFc=20Pottier?= <loic.pottier@live.fr>
Date: Fri, 13 Sep 2024 18:17:13 +0200
Subject: [PATCH 11/11] fix divide function

---
 examples/poly_operations/main.v |  6 ++---
 poly/poly.v                     | 47 ++++++++++++++++++++-------------
 poly/poly_test.v                | 17 +++++++-----
 3 files changed, 42 insertions(+), 28 deletions(-)

diff --git a/examples/poly_operations/main.v b/examples/poly_operations/main.v
index dd3d552b2..8985948d5 100644
--- a/examples/poly_operations/main.v
+++ b/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)
diff --git a/poly/poly.v b/poly/poly.v
index da10bf14c..fdaaa70f0 100644
--- a/poly/poly.v
+++ b/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
 }
diff --git a/poly/poly_test.v b/poly/poly_test.v
index d28a2fb40..82571b6eb 100644
--- a/poly/poly_test.v
+++ b/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 == []
 }