Polynomial Arithmetics

Polynomial Arithmetic

In SepalSolver, polynomials are not just static formulas; they are dynamic objects that can be added, subtracted, multiplied, and divided. By representing polynomials as an array of doubles in descending order, we can leverage classical signal processing algorithms, such as convolution, to perform these operations with high numerical efficiency.

1. Addition and Subtraction (PolyAdd / PolySub)

Adding or subtracting polynomials involves summing the coefficients of corresponding powers. Since our arrays are in descending order, we must carefully align the “tails” of the arrays (the constant terms) before performing the operation, especially if the polynomials have different degrees.

double[] p1 = [1.0, 2.0, 3.0]; // x^2 + 2x + 3
double[] p2 = [4.0, 5.0];      // 4x + 5

// Result of PolyAdd: { 1.0, 6.0, 8.0 } -> x^2 + 6x + 8
double[] sum = Polyadd(p1, p2);
Console.WriteLine($"Sum:\n = [{string.Join(", ", sum)}]");

// Result of PolySub: { 1.0, -2.0, -2.0 } -> x^2 - 2x - 2
double[] diff = Polysub(p1, p2);
Console.WriteLine($"Diff:\n = [{string.Join(", ", diff)}]");

Ouput

Sum:
 = [1, 6, 8]
Diff:
 = [1, -2, -2]

2. Multiplication (Conv)

Multiplying two polynomials is mathematically equivalent to the convolution of their coefficient arrays. If you multiply a polynomial of degree \(N\) by one of degree \(M\), the resulting degree will be \(N + M\). We use a sliding-window approach to compute each term of the resulting array.

double[] p1 = [1.0, 1.0]; // x + 1
double[] p2 = [1.0, -1.0]; // x - 1

//Conv(p1, p2) results in { 1.0, 0.0, -1.0 } -> x^2 - 1
double[] product = Conv(p1, p2);
Console.WriteLine($"Product:\n = [{string.Join(", ", product)}]");

Ouput

Product:
 = [1, 0, -1]

3. Division and Remainder (Deconv)

Polynomial division is implemented via deconvolution (synthetic division). Dividing \(P_1\) by \(P_2\) yields two results: the Quotient \((Q)\) and the Remainder \((R)\). This is vital for operations like partial fraction expansion or checking for factorability.

double[] dividend = [1.0, 0.0, -1.0]; // x^2 - 1
double[] divisor = [1.0, 1.0];        // x + 1

// Deconv returns Quotient [1.0, -1.0] (x - 1) and Remainder { 0.0 }
(double[] Q, double[] R) = Deconv(dividend, divisor);

Console.WriteLine($"Quotient:\n = [{string.Join(", ", Q)}]");
Console.WriteLine($"Remainder:\n = [{string.Join(", ", R)}]");

Ouput

Quotient:
 = [1, -1]
Remainder:
 = [0, 0, 0]

Examples

Example 1 : Combining Structural Models

Suppose you are combining two transfer functions in a control system. You need to add the characteristic polynomials of two sub-components. Using PolyAdd, the solver automatically handles the degree mismatch between a quadratic and a linear system.

double[] modelA = { 1.0, -5.0, 6.0 }; // s^2 - 5s + 6
double[] modelB = { 1.0, 2.0 };       // s + 2
double[] combined = Polyadd(modelA, modelB);

// Result: s^2 - 4s + 8
Console.WriteLine($"New Degree: {combined.Length - 1}");
Console.WriteLine($"Combined Coefficients: [{string.Join(", ", combined)}]");

Ouput

New Degree: 2
Combined Coefficients: [1, -4, 8]

In root-finding verification, you might want to multiply factors \((x - r_1)(x - r_2)\) to see if you recover the original polynomial. Conv makes this/// expansion trivial and numerically stable.

double[] factor1 = [1.0, -2.0]; // (x - 2)
double[] factor2 = { 1.0, -3.0 }; // (x - 3)
double[] expanded = Conv(factor1, factor2);

// Result: { 1.0, -5.0, 6.0 } -> x^2 - 5x + 6
Console.WriteLine($"Expanded = [{string.Join(", ", expanded)}]");

Ouput

Expanded = [1, -5, 6]

Example 3 : Finding Remainders in Signal Filtering

When dividing polynomials to check for divisibility or performing long division in a filter design, the remainder tells us if one polynomial is a perfect factor of another.

double[] P = [1.0, 2.0, 1.0]; // x^2 + 2x + 1
double[] D = { 1.0, 0.0 };      // x
var (Q, R) = Deconv(P, D);
// Q is {1.0, 2.0} (x + 2), R is {1.0} (remainder 1)
Console.WriteLine($"Remainder is zero? {R.All(val => val == 0)}");

Ouput

Remainder is zero? False

Numerical Note: Array Sizing

During PolyAdd and PolySub, SepalSolver creates a result array equal to the size of the largest input. During Conv, the result size is always p1.Length + p2.Length - 1. Understanding these allocation rules helps in managing memory when performing recursive polynomial arithmetic.