From c76b09797402eea9f39ad50f1a1791e3c53ab07c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=BCleyman=20Kaya?= <111174999+suleyman-kaya@users.noreply.github.com> Date: Mon, 5 Aug 2024 03:44:20 +0300 Subject: [PATCH] poly: edit multiply and add divide functions (#215) * edit multiply function, add divison, and fix comment lines * update readme * fix readme * Update README.md --- poly/README.md | 151 ++++++++++++++++++++++++++++++----------- poly/poly.v | 173 ++++++++++++++++++++++++----------------------- poly/poly_test.v | 51 ++++++++------ 3 files changed, 231 insertions(+), 144 deletions(-) diff --git a/poly/README.md b/poly/README.md index 2fe043dcb..0e6855662 100644 --- a/poly/README.md +++ b/poly/README.md @@ -1,13 +1,16 @@ # Polynomials -This chapter describes functions for evaluating and solving polynomials. -There are routines for finding real and complex roots of quadratic and -cubic equations using analytic methods. An iterative polynomial solver -is also available for finding the roots of general polynomials with real -coefficients (of any order). The functions are declared in the module `vsl.poly`. +This chapter describes functions for evaluating and solving polynomials. There are routines +for finding real and complex roots of quadratic and cubic equations using analytic methods. +An iterative polynomial solver is also available for finding the roots of general polynomials +with real coefficients (of any order). The functions are declared in the module `vsl.poly`. ## Polynomial Evaluation +```v ignore +fn eval(c []f64, x f64) f64 +``` + The functions described here evaluate the polynomial ```console @@ -16,20 +19,13 @@ P(x) = c[0] + c[1] x + c[2] x^2 + . . . + c[len-1] x^(len-1) using Horner's method for stability. -```v ignore -fn eval(c []f64, x f64) f64 -``` - -This function evaluates a polynomial with real coefficients for the real variable `x`. - ```v ignore fn eval_derivs(c []f64, x f64, lenres u64) []f64 ``` -This function evaluates a polynomial and its derivatives storing the -results in the array `res` of size `lenres`. The output array -contains the values of `d^k P(x)/d x^k` for the specified value of -`x` starting with `k = 0`. +This function evaluates a polynomial and its derivatives, storing the results in the array +`res` of size `lenres`. The output array contains the values of `d^k P(x)/d x^k` for the +specified value of `x`, starting with `k = 0`. ## Quadratic Equations @@ -43,22 +39,17 @@ This function finds the real roots of the quadratic equation, a x^2 + b x + c = 0 ``` -The number of real roots (either zero, one or two) is returned, and -their locations are are returned as `[ x0, x1 ]`. If no real roots -are found then `[]` is returned. If one real root -is found (i.e. if `a=0`) then it is are returned as `[ x0 ]`. When two -real roots are found they are are returned as `[ x0, x1 ]` in -ascending order. The case of coincident roots is not considered -special. For example `(x-1)^2=0` will have two roots, which happen -to have exactly equal values. +The number of real roots (either zero, one or two) is returned, and their locations are +returned as `[ x0, x1 ]`. If no real roots are found then `[]` is returned. If one real root +is found (i.e. if `a=0`) then it is returned as `[ x0 ]`. When two real roots are found they +are returned as `[ x0, x1 ]` in ascending order. The case of coincident roots is not considered +special. For example `(x-1)^2=0` will have two roots, which happen to have exactly equal values. -The number of roots found depends on the sign of the discriminant -`b^2 - 4 a c`. This will be subject to rounding and cancellation -errors when computed in double precision, and will also be subject to -errors if the coefficients of the polynomial are inexact. These errors -may cause a discrete change in the number of roots. However, for -polynomials with small integer coefficients the discriminant can always -be computed exactly. +The number of roots found depends on the sign of the discriminant `b^2 - 4 a c`. This will +be subject to rounding and cancellation errors when computed in double precision, and will +also be subject to errors if the coefficients of the polynomial are inexact. These errors may +cause a discrete change in the number of roots. However, for polynomials with small integer +coefficients the discriminant can always be computed exactly. ## Cubic Equations @@ -72,13 +63,93 @@ This function finds the real roots of the cubic equation, x^3 + a x^2 + b x + c = 0 ``` -with a leading coefficient of unity. The number of real roots (either -one or three) is returned, and their locations are returned as `[ x0, x1, x2 ]`. -If one real root is found then only `[ x0 ]` -is returned. When three real roots are found they are returned as -`[ x0, x1, x2 ]` in ascending order. The case of -coincident roots is not considered special. For example, the equation -`(x-1)^3=0` will have three roots with exactly equal values. As -in the quadratic case, finite precision may cause equal or -closely-spaced real roots to move off the real axis into the complex -plane, leading to a discrete change in the number of real roots. +with a leading coefficient of unity. The number of real roots (either one or three) is +returned, and their locations are returned as `[ x0, x1, x2 ]`. If one real root is found +then only `[ x0 ]` is returned. When three real roots are found they are returned as +`[ x0, x1, x2 ]` in ascending order. The case of coincident roots is not considered special. +For example, the equation `(x-1)^3=0` will have three roots with exactly equal values. As +in the quadratic case, finite precision may cause equal or closely-spaced real roots to move +off the real axis into the complex plane, leading to a discrete change in the number of real roots. + +## Companion Matrix + +```v ignore +fn companion_matrix(a []f64) [][]f64 +``` + +Creates a companion matrix for the polynomial + +```console +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] +``` + +## Balanced Companion Matrix + +```v ignore +fn balance_companion_matrix(cm [][]f64) [][]f64 +``` + +Balances a companion matrix `C` to improve numerical stability. It uses an iterative scaling +process to make the row and column norms as close to each other as possible. The output is +a balanced matrix `B` such that `D^(-1)CD = B`, where `D` is a diagonal matrix. + +## Polynomial Operations + +```v ignore +fn add(a []f64, b []f64) []f64 +``` + +Adds two polynomials: + +```console +(a_n * x^n + ... + a_0) + (b_m * x^m + ... + b_0) +``` + +Returns the result as `[a_0 + b_0, a_1 + b_1, ..., a_k + b_k ...]`. + +```v ignore +fn subtract(a []f64, b []f64) []f64 +``` + +Subtracts two polynomials: + +```console +(a_n * x^n + ... + a_0) - (b_m * x^m + ... + b_0) +``` + +Returns the result as `[a_0 - b_0, a_1 - b_1, ..., a_k - b_k, ...]`. + +```v ignore +fn multiply(a []f64, b []f64) []f64 +``` + +Multiplies two polynomials: + +```console +(a_n * x^n + ... + a_0) * (b_m * x^m + ... + b_0) +``` + +Returns the result as `[c_0, c_1, ..., c_{n+m}]` where `c_k = ∑_{i+j=k} a_i * b_j`. + +```v ignore +fn divide(a []f64, b []f64) ([]f64, []f64) +``` + +Divides two polynomials: + +```console +(a_n * x^n + ... + a_0) / (b_m * x^m + ... + b_0) +``` + +Uses polynomial long division algorithm. Returns `(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)`. diff --git a/poly/poly.v b/poly/poly.v index 41926d9ab..401ce60e6 100644 --- a/poly/poly.v +++ b/poly/poly.v @@ -6,6 +6,10 @@ import vsl.errors const radix = 2 const radix2 = (radix * radix) +// eval is a function that 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 +22,9 @@ pub fn eval(c []f64, x f64) f64 { return ans } +// eval_derivs 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 +56,11 @@ 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 +// solve_quadratic 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 +87,10 @@ pub fn solve_quadratic(a f64, b f64, c f64) []f64 { // Handle linear case } } +// solve_cubic 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 +103,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 +127,17 @@ pub fn solve_cubic(a f64, b f64, c f64) []f64 { } } +// swap_ 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 } +// sorted_3_ 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 +155,15 @@ fn sorted_3_(x_ f64, y_ f64, z_ f64) (f64, f64, f64) { return x, y, z } +// companion_matrix 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 +182,10 @@ pub fn companion_matrix(a []f64) [][]f64 { return cm } +// balance_companion_matrix 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,7 +194,7 @@ 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 { @@ -177,7 +202,7 @@ pub fn balance_companion_matrix(cm [][]f64) [][]f64 { 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 +247,66 @@ 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 -// +// add 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] - } +// subtract 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 } +// multiply 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 +} + +// divide 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 index 8db4f0e1d..b9c25cdb2 100644 --- 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 == [] // The empty set indicates that two polynomials divide each other exactly (without remainder). +}