vlang · suleyman-kaya · Jul 29, 2024 · Jul 29, 2024 · Jul 29, 2024 · Jul 29, 2024
diff --git a/poly/poly.v b/poly/poly.v
@@ -6,6 +6,11 @@ import vsl.errors
 const radix = 2
 const radix2 = (radix * radix)
 
+/* Evaluates a polynomial P(x) = a_n * x^n + a_{n-1} * x^{n-1} + ... + a_1 * x + a_0
+   using Horner's method: P(x) = (...((a_n * x + a_{n-1}) * x + a_{n-2}) * x + ... + a_1) * x + a_0
+   Input: c = [a_0, a_1, ..., a_n], x
+   Output: P(x)
+*/
 pub fn eval(c []f64, x f64) f64 {
 	if c.len == 0 {
 		errors.vsl_panic('coeficients can not be empty', .efailed)
@@ -18,6 +23,10 @@ pub fn eval(c []f64, x f64) f64 {
 	return ans
 }
 
+/* Evaluates a polynomial P(x) and its derivatives P'(x), P''(x), ..., P^(k)(x)
+   Input: c = [a_0, a_1, ..., a_n] representing P(x), x, and lenres (k+1)
+   Output: [P(x), P'(x), P''(x), ..., P^(k)(x)]
+*/
 pub fn eval_derivs(c []f64, x f64, lenres int) []f64 {
 	mut res := []f64{}
 	lenc := c.len
@@ -49,7 +58,12 @@ pub fn eval_derivs(c []f64, x f64, lenres int) []f64 {
 	return res
 }
 
-pub fn solve_quadratic(a f64, b f64, c f64) []f64 { // Handle linear case
+/* Solves the quadratic equation ax^2 + bx + c = 0
+   using the quadratic formula: x = (-b ± √(b^2 - 4ac)) / (2a)
+   Input: a, b, c
+   Output: Array of real roots (if any)
+*/
+pub fn solve_quadratic(a f64, b f64, c f64) []f64 {
 	if a == 0 {
 		if b == 0 {
 			return []
@@ -76,6 +90,11 @@ pub fn solve_quadratic(a f64, b f64, c f64) []f64 { // Handle linear case
 	}
 }
 
+/* Solves the cubic equation x^3 + ax^2 + bx + c = 0
+   using Cardano's formula and trigonometric solution
+   Input: a, b, c
+   Output: Array of real roots
+*/
 pub fn solve_cubic(a f64, b f64, c f64) []f64 {
 	q_ := (a * a - 3.0 * b)
 	r_ := (2.0 * a * a * a - 9.0 * a * b + 27.0 * c)
@@ -88,15 +107,6 @@ pub fn solve_cubic(a f64, b f64, c f64) []f64 {
 	if r == 0.0 && q == 0.0 {
 		return [-a / 3.0, -a / 3.0, -a / 3.0]
 	} else if cr2 == cq3 {
-		/*
-		this test is actually r2 == q3, written in a form suitable
-                        for exact computation with integers
-		*/
-		/*
-		Due to finite precision some double roots may be missed, and
-                        considered to be a pair of complex roots z = x +/- epsilon i
-                        close to the real axis.
-		*/
 		sqrt_q := math.sqrt(q)
 		if r > 0.0 {
 			return [-2.0 * sqrt_q - a / 3.0, sqrt_q - a / 3.0, sqrt_q - a / 3.0]
@@ -121,11 +131,19 @@ pub fn solve_cubic(a f64, b f64, c f64) []f64 {
 	}
 }
 
+/* Swaps two numbers: f(a, b) = (b, a)
+   Input: a, b
+   Output: (b, a)
+*/
 @[inline]
 fn swap_(a f64, b f64) (f64, f64) {
 	return b, a
 }
 
+/* Sorts three numbers in ascending order: f(x, y, z) = (min(x,y,z), median(x,y,z), max(x,y,z))
+   Input: x, y, z
+   Output: (min, median, max)
+*/
 @[inline]
 fn sorted_3_(x_ f64, y_ f64, z_ f64) (f64, f64, f64) {
 	mut x := x_
@@ -143,6 +161,16 @@ fn sorted_3_(x_ f64, y_ f64, z_ f64) (f64, f64, f64) {
 	return x, y, z
 }
 
+/* Creates a companion matrix for the polynomial P(x) = a_n * x^n + a_{n-1} * x^{n-1} + ... + a_1 * x + a_0
+   The companion matrix C is defined as:
+   [0 0 0 ... 0 -a_0/a_n]
+   [1 0 0 ... 0 -a_1/a_n]
+   [0 1 0 ... 0 -a_2/a_n]
+   [. . . ... . ........]
+   [0 0 0 ... 1 -a_{n-1}/a_n]
+   Input: a = [a_0, a_1, ..., a_n]
+   Output: Companion matrix C
+*/
 pub fn companion_matrix(a []f64) [][]f64 {
 	nc := a.len - 1
 	mut cm := [][]f64{len: nc, init: []f64{len: nc}}
@@ -161,6 +189,11 @@ pub fn companion_matrix(a []f64) [][]f64 {
 	return cm
 }
 
+/* Balances a companion matrix C to improve numerical stability
+   Uses an iterative scaling process to make the row and column norms as close to each other as possible
+   Input: Companion matrix C
+   Output: Balanced matrix B such that D^(-1)CD = B, where D is a diagonal matrix
+*/
 pub fn balance_companion_matrix(cm [][]f64) [][]f64 {
 	nc := cm.len
 	mut m := cm.clone()
@@ -169,15 +202,15 @@ pub fn balance_companion_matrix(cm [][]f64) [][]f64 {
 	mut col_norm := 0.0
 	for not_converged {
 		not_converged = false
-		for i := 0; i < nc; i++ { // column norm, excluding the diagonal
+		for i := 0; i < nc; i++ {
 			if i != nc - 1 {
 				col_norm = math.abs(m[i + 1][i])
 			} else {
 				col_norm = 0.0
 				for j := 0; j < nc - 1; j++ {
 					col_norm += math.abs(m[j][nc - 1])
 				}
-			} // row norm, excluding the diagonal
+			}
 			if i == 0 {
 				row_norm = math.abs(m[0][nc - 1])
 			} else if i == nc - 1 {
@@ -222,86 +255,70 @@ pub fn balance_companion_matrix(cm [][]f64) [][]f64 {
 	return m
 }
 
-// Arithmetic operations on polynomials
-//
-// In the following descriptions a, b, c are polynomials of degree
-// na, nb, nc respectively.
-//
-// Each polynomial is represented by an array containing its
-// coefficients, together with a separately declared integer equal
-// to the degree of the polynomial.  The coefficients appear in
-// ascending order; that is,
-//
-// a(x)  =  a[0]  +  a[1] * x  +  a[2] * x^2   +  ...  +  a[na] * x^na  .
-//
-//
-//
-// sum = eval( a, x )    Evaluate polynomial a(t) at t = x.
-// c = add( a, b )   c = a + b, nc = max(na, nb)
-// c = sub( a, b )   c = a - b, nc = max(na, nb)
-// c = mul( a, b )   c = a * b, nc = na+nb
-//
-//
-// Division:
-//
-// c = div( a, b )       c = a / b, nc = MAXPOL
-//
-// returns i = the degree of the first nonzero coefficient of a.
-// The computed quotient c must be divided by x^i.  An error message
-// is printed if a is identically zero.
-//
-//
-// Change of variables:
-// If a and b are polynomials, and t = a(x), then
-//     c(t) = b(a(x))
-// is a polynomial found by substituting a(x) for t.  The
-// subroutine call for this is
-//
+/* Adds two polynomials: (a_n * x^n + ... + a_0) + (b_m * x^m + ... + b_0)
+   Input: a = [a_0, ..., a_n], b = [b_0, ..., b_m]
+   Output: [a_0 + b_0, a_1 + b_1, ..., max(a_k, b_k), ...]
+*/
 pub fn add(a []f64, b []f64) []f64 {
-	na := a.len
-	nb := b.len
-	nc := int(math.max(na, nb))
-	mut c := []f64{len: nc}
-	for i := 0; i < nc; i++ {
-		if i > na {
-			c[i] = b[i]
-		} else if i > nb {
-			c[i] = a[i]
-		} else {
-			c[i] = a[i] + b[i]
-		}
+	mut result := []f64{len: math.max(a.len, b.len)}
+	for i in 0 .. result.len {
+		result[i] = if i < a.len { a[i] } else { 0.0 } + if i < b.len { b[i] } else { 0.0 }
 	}
-	return c
+	return result
 }
 
-pub fn substract(a []f64, b []f64) []f64 {
-	na := a.len
-	nb := b.len
-	nc := int(math.max(na, nb))
-	mut c := []f64{len: nc}
-	for i := 0; i < nc; i++ {
-		if i > na {
-			c[i] = -b[i]
-		} else if i > nb {
-			c[i] = a[i]
-		} else {
-			c[i] = a[i] - b[i]
-		}
+/* Subtracts two polynomials: (a_n * x^n + ... + a_0) - (b_m * x^m + ... + b_0)
+   Input: a = [a_0, ..., a_n], b = [b_0, ..., b_m]
+   Output: [a_0 - b_0, a_1 - b_1, ..., a_k - b_k, ...]
+*/
+pub fn subtract(a []f64, b []f64) []f64 {
+	mut result := []f64{len: math.max(a.len, b.len)}
+	for i in 0 .. result.len {
+		result[i] = if i < a.len { a[i] } else { 0.0 } - if i < b.len { b[i] } else { 0.0 }
 	}
-	return c
+	return result
 }
 
+/* Multiplies two polynomials: (a_n * x^n + ... + a_0) * (b_m * x^m + ... + b_0)
+   Input: a = [a_0, ..., a_n], b = [b_0, ..., b_m]
+   Output: [c_0, c_1, ..., c_{n+m}] where c_k = ∑_{i+j=k} a_i * b_j
+*/
 pub fn multiply(a []f64, b []f64) []f64 {
-	na := a.len
-	nb := b.len
-	nc := na + nb
-	mut c := []f64{len: nc}
-	for i := 0; i < na; i++ {
-		x := a[i]
-		for j := 0; j < nb; j++ {
-			k := i + j
-			c[k] += x * b[j]
+	mut result := []f64{len: a.len + b.len - 1}
+	for i in 0 .. a.len {
+		for j in 0 .. b.len {
+			result[i + j] += a[i] * b[j]
+		}
+	}
+	return result
+}
+
+/* Divides two polynomials: (a_n * x^n + ... + a_0) / (b_m * x^m + ... + b_0)
+   Uses polynomial long division algorithm
+   Input: a = [a_0, ..., a_n], b = [b_0, ..., b_m]
+   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{}
+	mut remainder := a.clone()
+	b_lead_coef := b[0]
+
+	for remainder.len >= b.len {
+		lead_coef := remainder[0] / b_lead_coef
+		quotient << lead_coef
+		for i in 0 .. b.len {
+			remainder[i] -= lead_coef * b[i]
+		}
+		remainder = remainder[1..]
+		for remainder.len > 0 && math.abs(remainder[0]) < 1e-10 {
+			remainder = remainder[1..]
 		}
 	}
-	return c
+
+	if remainder.len == 0 {
+		remainder = []f64{}
+	}
+
+	return quotient, remainder
 }
diff --git a/poly/poly_test.v b/poly/poly_test.v
@@ -1,5 +1,7 @@
 module poly
 
+import math
+
 fn test_eval() {
 	// ans = 2
 	// ans = 4.0 + 4 * 2 = 12
@@ -16,29 +18,38 @@ 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_add() {
-	a := [6.0, 777, -3]
-	b := [1.0, -755, -4]
-	assert add(a, b) == [7.0, 22, -7]
+	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]
 }
 
-fn test_substract() {
-	a := [6.0, 777, -3]
-	b := [1.0, -755, -4]
-	assert substract(a, b) == [5.0, 1532, 1]
+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]
 }
 
-// fn test_multiply(){
-// 	a := [9.0, -1, 5]
-// 	b := [0.0, -1, 7]
+fn test_multiply() {
+	// (2+3x+4x^2) * (-3x+2x^2) = (-6x -5x^2 -6x^3 + 8x^4)
+	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]
+}
 
-// 	assert multiply(a, b) == [0.0, -9, 73, -12, 35, 0]
-// }
+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)
+	println('Divide quotient: ${quotient}')
+	println('Divide remainder: ${remainder}')
+	assert quotient == [1.0, 1.0]
+	assert remainder == [] // boş küme bu iki polinomun tam bölündüğünü(kalansız bölündüğünü) ifade eder
+}