1 文本格式

using System;
using System.Text;

namespace Legalsoft.Truffer
{
    /// <summary>
    /// operations on polynomials
    /// </summary>
    public class Poly
    {
        /// <summary>
        /// polynomial c[0]+c[1]x+c[2]x^2+ ... + c[n-2]x^n-2 + c[n-1]x^n-1
        /// </summary>
        private double[] c { get; set; }

        /// <summary>
        /// Construct polynomial
        /// </summary>
        /// <param name="cc"></param>
        public Poly(double[] cc)
        {
            this.c = cc;
        }

        public double poly(double x)
        {
            int j;
            double p = c[j = c.Length - 1];
            while (j > 0)
            {
                p = p * x + c[--j];
            }
            return p;
        }

        public String toString()
        {
            StringBuilder sb = new StringBuilder(32);
            int j = c.Length - 1;
            sb.Append(String.Format("%fx^%d", c[j], j));
            j--;
            for (; j != 0; j--)
            {
                sb.Append(String.Format("%+fx^%d", c[j], j));
            }
            sb.Append(String.Format("%+f ", c[0]));
            return sb.ToString().Substring(0);
        }

        /// <summary>
        /// Build Polynomial from roots 
        /// </summary>
        /// <param name="z"></param>
        /// <returns></returns>
        /// <exception cref="Exception"></exception>
        public static Poly buildFromRoots(Complex[] z)
        {
            for (int i = 0; i < z.Length; i++)
            {
                bool found = false;
                for (int j = 0; j < z.Length; j++)
                {
                    //if (z[i].re == z[j].re && z[i].im == -z[j].im)
                    if (Math.Abs(z[i].re - z[j].re) <= float.Epsilon && Math.Abs(z[i].im - (-z[j].im)) <= float.Epsilon)
                    {
                        found = true;
                        break;
                    }
                }
                if (!found)
                {
                    throw new Exception("Roots must be conjugate");
                }
            }

            Complex[] c = new Complex[z.Length + 1];
            c[0] = z[0].neg();
            c[1] = new Complex(1, 0);
            for (int i = 1; i < z.Length; i++)
            {
                Complex d = c[0];
                c[0] = c[0].mul(z[i].neg());
                for (int j = 1; j < i + 1; j++)
                {
                    Complex dd = c[j];
                    c[j] = d.sub(z[i].mul(c[j]));
                    d = dd;
                }
                c[i + 1] = d;
            }
            double[] cc = new double[c.Length];
            for (int i = 0; i < cc.Length; i++)
            {
                cc[i] = c[i].re;
            }
            return new Poly(cc);
        }

        /// <summary>
        /// Build Polynomial from roots 
        /// </summary>
        /// <param name="z"></param>
        /// <returns></returns>
        public static Poly buildFromRoots(double[] z)
        {
            double[] c = new double[z.Length + 1];
            c[0] = -z[0]; c[1] = 1;
            for (int i = 1; i < z.Length; i++)
            {
                double d = c[0];
                c[0] *= -z[i];
                for (int j = 1; j < i + 1; j++)
                {
                    double dd = c[j];
                    c[j] = d - z[i] * c[j];
                    d = dd;
                }
                c[i + 1] = d;
            }

            return new Poly(c);
        }

        /// <summary>
        /// Given the coefficients of a polynomial of degree nc as an array c[0..nc] of
        /// size nc+1 (with c[0] being the constant term), and given a value x, this
        /// routine fills an output array pd of size nd+1 with the value of the
        /// polynomial evaluated at x in pd[0], and the first nd derivatives at x in
        /// pd[1..nd].
        /// </summary>
        /// <param name="c"></param>
        /// <param name="x"></param>
        /// <param name="pd"></param>
        public static void ddpoly(double[] c, double x, double[] pd)
        {
            int nc = c.Length - 1;
            int nd = pd.Length - 1;
            double cnst = 1.0;
            pd[0] = c[nc];
            for (int j = 1; j < nd + 1; j++)
            {
                pd[j] = 0.0;
            }
            for (int i = nc - 1; i >= 0; i--)
            {
                int nnd = (nd < (nc - i) ? nd : nc - i);
                for (int j = nnd; j > 0; j--)
                {
                    pd[j] = pd[j] * x + pd[j - 1];
                }
                pd[0] = pd[0] * x + c[i];
            }
            for (int i = 2; i < nd + 1; i++)
            {
                cnst *= i;
                pd[i] *= cnst;
            }
        }

        /// <summary>
        /// Given the coefficients of a polynomial of degree nc as an array c[0..nc] of
        /// size nc+1 (with c[0] being the constant term), and given a value x, this
        /// routine fills an output array pd of size nd+1 with the value of the
        /// polynomial evaluated at x in pd[0], and the first nd derivatives at x in
        /// pd[1..nd].
        /// </summary>
        /// <param name="u"></param>
        /// <param name="v"></param>
        /// <param name="q"></param>
        /// <param name="r"></param>
        /// <exception cref="Exception"></exception>
        public static void poldiv(double[] u, double[] v, double[] q, double[] r)
        {
            int n = u.Length - 1;
            int nv = v.Length - 1;
            //while (nv >= 0 && v[nv] == 0.0)
            while (nv >= 0 && Math.Abs(v[nv]) <= float.Epsilon)
            {
                nv--;
            }
            if (nv < 0)
            {
                throw new Exception("poldiv divide by zero polynomial");
            }

            //r = u;
            r = Globals.CopyFrom(u);
            //q.assign(u.Length, 0.0);
            for (int k = n - nv; k >= 0; k--)
            {
                q[k] = r[nv + k] / v[nv];
                for (int j = nv + k - 1; j >= k; j--)
                {
                    r[j] -= q[k] * v[j - k];
                }
            }
            for (int j = nv; j <= n; j++)
            {
                r[j] = 0.0;
            }
        }

    }
}
 

2 代码格式

using System;
using System.Text;

namespace Legalsoft.Truffer
{
    /// <summary>
    /// operations on polynomials
    /// </summary>
    public class Poly
    {
        /// <summary>
        /// polynomial c[0]+c[1]x+c[2]x^2+ ... + c[n-2]x^n-2 + c[n-1]x^n-1
        /// </summary>
        private double[] c { get; set; }

        /// <summary>
        /// Construct polynomial
        /// </summary>
        /// <param name="cc"></param>
        public Poly(double[] cc)
        {
            this.c = cc;
        }

        public double poly(double x)
        {
            int j;
            double p = c[j = c.Length - 1];
            while (j > 0)
            {
                p = p * x + c[--j];
            }
            return p;
        }

        public String toString()
        {
            StringBuilder sb = new StringBuilder(32);
            int j = c.Length - 1;
            sb.Append(String.Format("%fx^%d", c[j], j));
            j--;
            for (; j != 0; j--)
            {
                sb.Append(String.Format("%+fx^%d", c[j], j));
            }
            sb.Append(String.Format("%+f ", c[0]));
            return sb.ToString().Substring(0);
        }

        /// <summary>
        /// Build Polynomial from roots 
        /// </summary>
        /// <param name="z"></param>
        /// <returns></returns>
        /// <exception cref="Exception"></exception>
        public static Poly buildFromRoots(Complex[] z)
        {
            for (int i = 0; i < z.Length; i++)
            {
                bool found = false;
                for (int j = 0; j < z.Length; j++)
                {
                    //if (z[i].re == z[j].re && z[i].im == -z[j].im)
                    if (Math.Abs(z[i].re - z[j].re) <= float.Epsilon && Math.Abs(z[i].im - (-z[j].im)) <= float.Epsilon)
                    {
                        found = true;
                        break;
                    }
                }
                if (!found)
                {
                    throw new Exception("Roots must be conjugate");
                }
            }

            Complex[] c = new Complex[z.Length + 1];
            c[0] = z[0].neg();
            c[1] = new Complex(1, 0);
            for (int i = 1; i < z.Length; i++)
            {
                Complex d = c[0];
                c[0] = c[0].mul(z[i].neg());
                for (int j = 1; j < i + 1; j++)
                {
                    Complex dd = c[j];
                    c[j] = d.sub(z[i].mul(c[j]));
                    d = dd;
                }
                c[i + 1] = d;
            }
            double[] cc = new double[c.Length];
            for (int i = 0; i < cc.Length; i++)
            {
                cc[i] = c[i].re;
            }
            return new Poly(cc);
        }

        /// <summary>
        /// Build Polynomial from roots 
        /// </summary>
        /// <param name="z"></param>
        /// <returns></returns>
        public static Poly buildFromRoots(double[] z)
        {
            double[] c = new double[z.Length + 1];
            c[0] = -z[0]; c[1] = 1;
            for (int i = 1; i < z.Length; i++)
            {
                double d = c[0];
                c[0] *= -z[i];
                for (int j = 1; j < i + 1; j++)
                {
                    double dd = c[j];
                    c[j] = d - z[i] * c[j];
                    d = dd;
                }
                c[i + 1] = d;
            }

            return new Poly(c);
        }

        /// <summary>
        /// Given the coefficients of a polynomial of degree nc as an array c[0..nc] of
        /// size nc+1 (with c[0] being the constant term), and given a value x, this
        /// routine fills an output array pd of size nd+1 with the value of the
        /// polynomial evaluated at x in pd[0], and the first nd derivatives at x in
        /// pd[1..nd].
        /// </summary>
        /// <param name="c"></param>
        /// <param name="x"></param>
        /// <param name="pd"></param>
        public static void ddpoly(double[] c, double x, double[] pd)
        {
            int nc = c.Length - 1;
            int nd = pd.Length - 1;
            double cnst = 1.0;
            pd[0] = c[nc];
            for (int j = 1; j < nd + 1; j++)
            {
                pd[j] = 0.0;
            }
            for (int i = nc - 1; i >= 0; i--)
            {
                int nnd = (nd < (nc - i) ? nd : nc - i);
                for (int j = nnd; j > 0; j--)
                {
                    pd[j] = pd[j] * x + pd[j - 1];
                }
                pd[0] = pd[0] * x + c[i];
            }
            for (int i = 2; i < nd + 1; i++)
            {
                cnst *= i;
                pd[i] *= cnst;
            }
        }

        /// <summary>
        /// Given the coefficients of a polynomial of degree nc as an array c[0..nc] of
        /// size nc+1 (with c[0] being the constant term), and given a value x, this
        /// routine fills an output array pd of size nd+1 with the value of the
        /// polynomial evaluated at x in pd[0], and the first nd derivatives at x in
        /// pd[1..nd].
        /// </summary>
        /// <param name="u"></param>
        /// <param name="v"></param>
        /// <param name="q"></param>
        /// <param name="r"></param>
        /// <exception cref="Exception"></exception>
        public static void poldiv(double[] u, double[] v, double[] q, double[] r)
        {
            int n = u.Length - 1;
            int nv = v.Length - 1;
            //while (nv >= 0 && v[nv] == 0.0)
            while (nv >= 0 && Math.Abs(v[nv]) <= float.Epsilon)
            {
                nv--;
            }
            if (nv < 0)
            {
                throw new Exception("poldiv divide by zero polynomial");
            }

            //r = u;
            r = Globals.CopyFrom(u);
            //q.assign(u.Length, 0.0);
            for (int k = n - nv; k >= 0; k--)
            {
                q[k] = r[nv + k] / v[nv];
                for (int j = nv + k - 1; j >= k; j--)
                {
                    r[j] -= q[k] * v[j - k];
                }
            }
            for (int j = nv; j <= n; j++)
            {
                r[j] = 0.0;
            }
        }

    }
}