using System;
using System.Collections.Generic;
using System.Text;
using System.Diagnostics;
using System.Collections;
namespace PolyLib
{
/// <summary>
/// By hanzzoid, 3 Jul 2007. Taken from http://www.codeproject.com/KB/cs/PolyLib.aspx
/// </summary>
public class Polynomial
{
#region Fields
/// <summary>
/// Coefficients a_0,...,a_n of a polynomial p, such that
/// p(x) = a_0 + a_1*x + a_2*x^2 + ... + a_n*x^n.
/// </summary>
public Complex[] Coefficients;
#endregion
#region constructors
/// <summary>
/// Inits zero polynomial p = 0.
/// </summary>
public Polynomial()
{
Coefficients = new Complex[1];
Coefficients[0] = Complex.Zero;
}
/// <summary>
/// Inits polynomial from given complex coefficient array.
/// </summary>
/// <param name="coeffs"></param>
public Polynomial(params Complex[] coeffs)
{
if (coeffs == null || coeffs.Length < 1)
{
Coefficients = new Complex[1];
Coefficients[0] = Complex.Zero;
}
else
{
Coefficients = (Complex[])coeffs.Clone();
}
}
/// <summary>
/// Inits polynomial from given real coefficient array.
/// </summary>
/// <param name="coeffs"></param>
public Polynomial(params double[] coeffs)
{
if (coeffs == null || coeffs.Length < 1)
{
Coefficients = new Complex[1];
Coefficients[0] = Complex.Zero;
}
else
{
Coefficients = new Complex[coeffs.Length];
for (int i = 0; i < coeffs.Length; i++)
Coefficients[i] = new Complex(coeffs[i]);
}
}
/// <summary>
/// Inits constant polynomial.
/// </summary>
/// <param name="coeffs"></param>
public Polynomial(Complex c)
{
Coefficients = new Complex[1];
//if (c == null)
// Coefficients[0] = Complex.Zero;
//else
Coefficients[0] = c;
}
/// <summary>
/// Inits constant polynomial.
/// </summary>
/// <param name="coeffs"></param>
public Polynomial(double c)
{
Coefficients = new Complex[1];
Coefficients[0] = new Complex(c);
}
/// <summary>
/// Inits polynomial from string like "2x^2 + 4x + (2+2i)"
/// </summary>
/// <param name="p"></param>
public Polynomial(string p)
{
throw new NotImplementedException();
}
#endregion
#region dynamics
/// <summary>
/// Computes value of the differentiated polynomial at x.
/// </summary>
/// <param name="p"></param>
/// <param name="x"></param>
/// <returns></returns>
public Complex Differentiate(Complex x)
{
Complex[] buf = new Complex[Degree];
for (int i = 0; i < buf.Length; i++)
buf[i] = (i + 1) * Coefficients[i + 1];
return (new Polynomial(buf)).Evaluate(x);
}
/// <summary>
/// Computes the definite integral within the borders a and b.
/// </summary>
/// <param name="a">Left integration border.</param>
/// <param name="b">Right integration border.</param>
/// <returns></returns>
public Complex Integrate(Complex a, Complex b)
{
Complex[] buf = new Complex[Degree + 2];
buf[0] = Complex.Zero; // this value can be arbitrary, in fact
for (int i = 1; i < buf.Length; i++)
buf[i] = Coefficients[i - 1] / i;
Polynomial p = new Polynomial(buf);
return (p.Evaluate(b) - p.Evaluate(a));
}
/// <summary>
/// Degree of the polynomial.
/// </summary>
public int Degree
{
get
{
return Coefficients.Length - 1;
}
}
/// <summary>
/// Checks if given polynomial is zero.
/// </summary>
/// <returns></returns>
public bool IsZero()
{
for (int i = 0; i < Coefficients.Length; i++)
if (Coefficients[i] != 0) return false;
return true;
}
/// <summary>
/// Evaluates polynomial by using the horner scheme.
/// </summary>
/// <param name="x"></param>
/// <returns></returns>
public Complex Evaluate(Complex x)
{
Complex buf = Coefficients[Degree];
for (int i = Degree - 1; i >= 0; i--)
{
buf = Coefficients[i] + x * buf;
}
return buf;
}
/// <summary>
/// Normalizes the polynomial, e.i. divides each coefficient by the
/// coefficient of a_n the greatest term if a_n != 1.
/// </summary>
public void Normalize()
{
this.Clean();
if (Coefficients[Degree] != Complex.One)
for (int k = 0; k <= Degree; k++)
Coefficients[k] /= Coefficients[Degree];
}
/// <summary>
/// Removes unnecessary zero terms.
/// </summary>
public void Clean()
{
int i;
for (i = Degree; i >= 0 && Coefficients[i] == 0; i--) ;
Complex[] coeffs = new Complex[i + 1];
for (int k = 0; k <= i; k++)
coeffs[k] = Coefficients[k];
Coefficients = (Complex[])coeffs.Clone();
}
/// <summary>
/// Factorizes polynomial to its linear factors.
/// </summary>
/// <returns></returns>
public FactorizedPolynomial Factorize()
{
// this is to be returned
FactorizedPolynomial p = new FactorizedPolynomial();
// cannot factorize polynomial of degree 0 or 1
if (this.Degree <= 1)
{
p.Factor = new Polynomial[] { this };
p.Power = new int[] { 1 };
return p;
}
Complex[] roots = Roots(this);
//ArrayList rootlist = new ArrayList();
//foreach (Complex z in roots) rootlist.Add(z);
//roots = null; // don't need you anymore
//rootlist.Sort();
//// number of different roots
//int num = 1; // ... at least one
//// ...or more?
//for (int i = 1; i < rootlist.Count; i++)
// if (rootlist[i] != rootlist[i - 1]) num++;
Polynomial[] factor = new Polynomial[roots.Length];
int[] power = new int[roots.Length];
//factor[0] = new Polynomial( new Complex[]{ -(Complex)rootlist[0] * Coefficients[Degree],
// Coefficients[Degree] } );
//power[0] = 1;
//num = 1;
//len = 0;
//for (int i = 1; i < rootlist.Count; i++)
//{
// len++;
// if (rootlist[i] != rootlist[i - 1])
// {
// factor[num] = new Polynomial(new Complex[] { -(Complex)rootlist[i], Complex.One });
// power[num] = len;
// num++;
// len = 0;
// }
//}
power[0] = 1;
factor[0] = new Polynomial(new Complex[] { -Coefficients[Degree] * (Complex)roots[0], Coefficients[Degree] });
for (int i = 1; i < roots.Length; i++)
{
power[i] = 1;
factor[i] = new Polynomial(new Complex[] { -(Complex)roots[i], Complex.One });
}
p.Factor = factor;
p.Power = power;
return p;
}
/// <summary>
/// Computes the roots of polynomial via Weierstrass iteration.
/// </summary>
/// <returns></returns>
public Complex[] Roots()
{
double tolerance = 1e-12;
int max_iterations = 30;
Polynomial q = Normalize(this);
//Polynomial q = p;
Complex[] z = new Complex[q.Degree]; // approx. for roots
Complex[] w = new Complex[q.Degree]; // Weierstra� corrections
// init z
for (int k = 0; k < q.Degree; k++)
//z[k] = (new Complex(.4, .9)) ^ k;
z[k] = Complex.Exp(2 * Math.PI * Complex.I * k / q.Degree);
for (int iter = 0; iter < max_iterations
&& MaxValue(q, z) > tolerance; iter++)
for (int i = 0; i < 10; i++)
{
for (int k = 0; k < q.Degree; k++)
w[k] = q.Evaluate(z[k]) / WeierNull(z, k);
for (int k = 0; k < q.Degree; k++)
z[k] -= w[k];
}
// clean...
for (int k = 0; k < q.Degree; k++)
{
z[k].Re = Math.Round(z[k].Re, 12);
z[k].Im = Math.Round(z[k].Im, 12);
}
return z;
}
/// <summary>
/// Computes the roots of polynomial p via Weierstrass iteration.
/// </summary>
/// <param name="p">Polynomial to compute the roots of.</param>
/// <param name="tolerance">Computation precision; e.g. 1e-12 denotes 12 exact digits.</param>
/// <param name="max_iterations">Maximum number of iterations; this value is used to bound
/// the computation effort if desired precision is hard to achieve.</param>
/// <returns></returns>
public Complex[] Roots(double tolerance, int max_iterations)
{
Polynomial q = Normalize(this);
Complex[] z = new Complex[q.Degree]; // approx. for roots
Complex[] w = new Complex[q.Degree]; // Weierstra� corrections
// init z
for (int k = 0; k < q.Degree; k++)
//z[k] = (new Complex(.4, .9)) ^ k;
z[k] = Complex.Exp(2 * Math.PI * Complex.I * k / q.Degree);
for (int iter = 0; iter < max_iterations
&& MaxValue(q, z) > tolerance; iter++)
for (int i = 0; i < 10; i++)
{
for (int k = 0; k < q.Degree; k++)
w[k] = q.Evaluate(z[k]) / WeierNull(z, k);
for (int k = 0; k < q.Degree; k++)
z[k] -= w[k];
}
// clean...
for (int k = 0; k < q.Degree; k++)
{
z[k].Re = Math.Round(z[k].Re, 12);
z[k].Im = Math.Round(z[k].Im, 12);
}
return z;
}
#endregion
#region statics
/// <summary>
/// Expands factorized polynomial p_1(x)^(k_1)*...*p_r(x)^(k_r) to its normal form a_0 + a_1 x + ... + a_n x^n.
/// </summary>
/// <param name="p"></param>
/// <returns></returns>
public static Polynomial Expand(FactorizedPolynomial p)
{
Polynomial q = new Polynomial(new Complex[] { Complex.One });
for (int i = 0; i < p.Factor.Length; i++)
{
for (int j = 0; j < p.Power[i]; j++)
q *= p.Factor[i];
q.Clean();
}
// clean...
for (int k = 0; k <= q.Degree; k++)
{
q.Coefficients[k].Re = Math.Round(q.Coefficients[k].Re, 12);
q.Coefficients[k].Im = Math.Round(q.Coefficients[k].Im, 12);
}
return q;
}
/// <summary>
/// Evaluates factorized polynomial p at point x.
/// </summary>
/// <param name="p"></param>
/// <returns></returns>
public static Complex Evaluate(FactorizedPolynomial p, Complex x)
{
Complex z = Complex.One;
for (int i = 0; i < p.Factor.Length; i++)
{
z *= Complex.Pow(p.Factor[i].Evaluate(x), p.Power[i]);
}
return z;
}
/// <summary>
/// Removes unncessary leading zeros.
/// </summary>
/// <param name="p"></param>
/// <returns></returns>
public static Polynomial Clean(Polynomial p)
{
int i;
for (i = p.Degree; i >= 0 && p.Coefficients[i] == 0; i--) ;
Complex[] coeffs = new Complex[i + 1];
for (int k = 0; k <= i; k++)
coeffs[k] = p.Coefficients[k];
return new Polynomial(coeffs);
}
/// <summary>
/// Normalizes the polynomial, e.i. divides each coefficient by the
/// coefficient of a_n the greatest term if a_n != 1.
/// </summary>
public static Polynomial Normalize(Polynomial p)
{
Polynomial q = Clean(p);
if (q.Coefficients[q.Degree] != Complex.One)
for (int k = 0; k <= q.Degree; k++)
q.Coefficients[k] /= q.Coefficients[q.Degree];
return q;
}
/// <summary>
/// Computes the roots of polynomial p via Weierstrass iteration.
/// </summary>
/// <param name="p"></param>
/// <returns></returns>
public static Complex[] Roots(Polynomial p)
{
double tolerance = 1e-12;
int max_iterations = 30;
Polynomial q = Normalize(p);
//Polynomial q = p;
Complex[] z = new Complex[q.Degree]; // approx. for roots
Complex[] w = new Complex[q.Degree]; // Weierstra� corrections
// init z
for (int k = 0; k < q.Degree; k++)
//z[k] = (new Complex(.4, .9)) ^ k;
z[k] = Complex.Exp(2 * Math.PI * Complex.I * k / q.Degree);
for (int iter = 0; iter < max_iterations
&& MaxValue(q, z) > tolerance; iter++)
for (int i = 0; i < 10; i++)
{
for (int k = 0; k < q.Degree; k++)
w[k] = q.Evaluate(z[k]) / WeierNull(z, k);
for (int k = 0; k < q.Degree; k++)
z[k] -= w[k];
}
// clean...
for (int k = 0; k < q.Degree; k++)
{
z[k].Re = Math.Round(z[k].Re, 12);
z[k].Im = Math.Round(z[k].Im, 12);
}
return z;
}
/// <summary>
/// Computes the roots of polynomial p via Weierstrass iteration.
/// </summary>
/// <param name="p">Polynomial to compute the roots of.</param>
/// <param name="tolerance">Computation precision; e.g. 1e-12 denotes 12 exact digits.</param>
/// <param name="max_iterations">Maximum number of iterations; this value is used to bound
/// the computation effort if desired pecision is hard to achieve.</param>
/// <returns></returns>
public static Complex[] Roots(Polynomial p, double tolerance, int max_iterations)
{
Polynomial q = Normalize(p);
Complex[] z = new Complex[q.Degree]; // approx. for roots
Complex[] w = new Complex[q.Degree]; // Weierstra� corrections
// init z
for (int k = 0; k < q.Degree; k++)
//z[k] = (new Complex(.4, .9)) ^ k;
z[k] = Complex.Exp(2 * Math.PI * Complex.I * k / q.Degree);
for (int iter = 0; iter < max_iterations
&& MaxValue(q, z) > tolerance; iter++)
for (int i = 0; i < 10; i++)
{
for (int k = 0; k < q.Degree; k++)
w[k] = q.Evaluate(z[k]) / WeierNull(z, k);
for (int k = 0; k < q.Degree; k++)
z[k] -= w[k];
}
// clean...
for (int k = 0; k < q.Degree; k++)
{
z[k].Re = Math.Round(z[k].Re, 12);
z[k].Im = Math.Round(z[k].Im, 12);
}
return z;
}
/// <summary>
/// Computes the greatest value |p(z_k)|.
/// </summary>
/// <param name="p"></param>
/// <param name="z"></param>
/// <returns></returns>
public static double MaxValue(Polynomial p, Complex[] z)
{
double buf = 0;
for (int i = 0; i < z.Length; i++)
{
if(Complex.Abs(p.Evaluate(z[i])) > buf)
buf = Complex.Abs(p.Evaluate(z[i]));
}
return buf;
}
/// <summary>
/// For g(x) = (x-z_0)*...*(x-z_n), this method returns
/// g'(z_k) = \prod_{j != k} (z_k - z_j).
/// </summary>
/// <param name="z"></param>
/// <returns></returns>
private static Complex WeierNull(Complex[] z, int k)
{
if (k < 0 || k >= z.Length)
throw new ArgumentOutOfRangeException();
Complex buf = Complex.One;
for (int j = 0; j < z.Length; j++)
if (j != k) buf *= (z[k] - z[j]);
return buf;
}
/// <summary>
/// Differentiates given polynomial p.
/// </summary>
/// <param name="p"></param>
/// <returns></returns>
public static Polynomial Derivative(Polynomial p)
{
Complex[] buf = new Complex[p.Degree];
for (int i = 0; i < buf.Length; i++)
buf[i] = (i + 1) * p.Coefficients[i + 1];
return new Polynomial(buf);
}
/// <summary>
/// Integrates given polynomial p.
/// </summary>
/// <param name="p"></param>
/// <returns></returns>
public static Polynomial Integral(Polynomial p)
{
Complex[] buf = new Complex[p.Degree + 2];
buf[0] = Complex.Zero; // this value can be arbitrary, in fact
for (int i = 1; i < buf.Length; i++)
buf[i] = p.Coefficients[i - 1] / i;
return new Polynomial(buf);
}
/// <summary>
/// Computes the monomial x^degree.
/// </summary>
/// <param name="degree"></param>
/// <returns></returns>
public static Polynomial Monomial(int degree)
{
if (degree == 0) return new Polynomial(1);
Complex[] coeffs = new Complex[degree + 1];
for (int i = 0; i < degree; i++)
coeffs[i] = Complex.Zero;
coeffs[degree] = Complex.One;
return new Polynomial(coeffs);
}
public static Polynomial[] GetStandardBase(int dim)
{
if(dim < 1)
throw new ArgumentException("Dimension expected to be greater than zero.");
Polynomial[] buf = new Polynomial[dim];
for (int i = 0; i < dim; i++)
buf[i] = Monomial(i);
return buf;
}
#endregion
#region overrides & operators
public static Polynomial operator +(Polynomial p, Polynomial q)
{
int degree = Math.Max(p.Degree, q.Degree);
Complex[] coeffs = new Complex[degree + 1];
for (int i = 0; i <= degree; i++)
{
if (i > p.Degree) coeffs[i] = q.Coefficients[i];
else if (i > q.Degree) coeffs[i] = p.Coefficients[i];
else coeffs[i] = p.Coefficients[i] + q.Coefficients[i];
}
return new Polynomial(coeffs);
}
public static Polynomial operator -(Polynomial p, Polynomial q)
{
return p + (-q);
}
public static Polynomial operator -(Polynomial p)
{
Complex[] coeffs = new Complex[p.Degree + 1];
for (int i = 0; i < coeffs.Length; i++)
coeffs[i] = -p.Coefficients[i];
return new Polynomial(coeffs);
}
public static Polynomial operator *(Complex d, Polynomial p)
{
Complex[] coeffs = new Complex[p.Degree + 1];
for (int i = 0; i < coeffs.Length; i++)
coeffs[i] = d * p.Coefficients[i];
return new Polynomial(coeffs);
}
public static Polynomial operator *(Polynomial p, Complex d)
{
Complex[] coeffs = new Complex[p.Degree + 1];
for (int i = 0; i < coeffs.Length; i++)
coeffs[i] = d * p.Coefficients[i];
return new Polynomial(coeffs);
}
public static Polynomial operator *(double d, Polynomial p)
{
Complex[] coeffs = new Complex[p.Degree + 1];
for (int i = 0; i < coeffs.Length; i++)
coeffs[i] = d * p.Coefficients[i];
return new Polynomial(coeffs);
}
public static Polynomial operator *(Polynomial p, double d)
{
Complex[] coeffs = new Complex[p.Degree + 1];
for (int i = 0; i < coeffs.Length; i++)
coeffs[i] = d * p.Coefficients[i];
return new Polynomial(coeffs);
}
public static Polynomial operator /(Polynomial p, Complex d)
{
Complex[] coeffs = new Complex[p.Degree + 1];
for (int i = 0; i < coeffs.Length; i++)
coeffs[i] = p.Coefficients[i] / d;
return new Polynomial(coeffs);
}
public static Polynomial operator /(Polynomial p, double d)
{
Complex[] coeffs = new Complex[p.Degree + 1];
for (int i = 0; i < coeffs.Length; i++)
coeffs[i] = p.Coefficients[i] / d;
return new Polynomial(coeffs);
}
public static Polynomial operator *(Polynomial p, Polynomial q)
{
int degree = p.Degree + q.Degree;
Polynomial r = new Polynomial();
for (int i = 0; i <= p.Degree; i++)
for (int j = 0; j <= q.Degree; j++)
r += (p.Coefficients[i] * q.Coefficients[j]) * Monomial(i + j);
return r;
}
public static Polynomial operator ^(Polynomial p, uint k)
{
if (k == 0)
return Monomial(0);
else if (k == 1)
return p;
else
return p * (p ^ (k - 1));
}
public override string ToString()
{
if (this.IsZero()) return "0";
else
{
string s = "";
for (int i = 0; i < Degree + 1; i++)
{
if (Coefficients[i] != Complex.Zero)
{
if (Coefficients[i] == Complex.I)
s += "i";
else if (Coefficients[i] != Complex.One)
{
if (Coefficients[i].IsReal() && Coefficients[i].Re > 0)
s += Coefficients[i].ToString();
else
s += "(" + Coefficients[i].ToString() + ")";
}
else if (/*Coefficients[i] == Complex.One && */i == 0)
s += 1;
if (i == 1)
s += "x";
else if (i > 1)
s += "x^" + i.ToString();
}
if (i < Degree && Coefficients[i + 1] != 0 && s.Length > 0)
s += " + ";
}
return s;
}
}
public string ToString(string format)
{
if (this.IsZero()) return "0";
else
{
string s = "";
for (int i = 0; i < Degree + 1; i++)
{
if (Coefficients[i] != Complex.Zero)
{
if (Coefficients[i] == Complex.I)
s += "i";
else if (Coefficients[i] != Complex.One)
{
if (Coefficients[i].IsReal() && Coefficients[i].Re > 0)
s += Coefficients[i].ToString(format);
else
s += "(" + Coefficients[i].ToString(format) + ")";
}
else if (/*Coefficients[i] == Complex.One && */i == 0)
s += 1;
if (i == 1)
s += "x";
else if (i > 1)
s += "x^" + i.ToString(format);
}
if (i < Degree && Coefficients[i + 1] != 0 && s.Length > 0)
s += " + ";
}
return s;
}
}
public override bool Equals(object obj)
{
return (this.ToString() == ((Polynomial)obj).ToString());
}
#endregion
#region structs
/// <summary>
/// Factorized polynomial p := set of polynomials p_1,...,p_k and their corresponding
/// powers n_1,...,n_k, such that p = (p_1)^(n_1)*...*(p_k)^(n_k).
/// </summary>
public struct FactorizedPolynomial
{
/// <summary>
/// Set of factors the polynomial consists of.
/// </summary>
public Polynomial[] Factor;
/// <summary>
/// Set of powers, where Factor[i] is lifted
/// to Power[i].
/// </summary>
public int[] Power;
}
#endregion
}
}