PolynomialConvolution.cs

Go to the documentation of this file.

using System;

namespace FFT
{
   ///
   /// Convolution based on polynomial evaluation and interpolation. 
   ///
   /// The vectors \f$ \mathbf{a} = [a_0,\ldots,a_{n-1}] \f$ and 
   /// \f$ \mathbf{b} = [b_0,\ldots,b_{m-1}] \f$ 
   /// are treated as cooeficients of two polynomials:
   /// \f$ A(x) = a_{n-1}*x^{n-1} + \ldots + a_1*x + a_0 \f$ and 
   /// \f$ B(x) = b_{m-1}*x^{m-1} + \ldots + b_1*x + b_0 \f$.
   /// The convolution is calculated in the following way 
   /// (#Convolution(Complex[], Complex[])):
   ///  -# generate the roots of the \f$(n+m)\f$-th Chebyshev polynomial
   ///     of the first kind (#GenerateDomain(Complex[], Complex[]))
   ///     \f[ \mathbf{x}=[x_0,\ldots,x_{n+m-1}], \quad\text{where}\quad x_k=\frac{\cos[(2k+1)\pi]}{2(n+m-1)+2}; \f]
   ///  -# evaluate the two polynomials on the generated domain 
   ///     (#Evaluate(Complex[], Complex[]))
   ///     \f[ \mathbf{y}^A = [A(x_0),\ldots,A(x_{n+m-1})]\quad\text{and}\quad\mathbf{y}^B = [B(x_0),\ldots,B(x_{n+m-1})] \f]
   ///  -# pointwise multiply the values (#PointwiseMultiply(Complex[], Complex[]))
   ///     \f[ \mathbf{y} = [A(x_0)*B(x_0),\ldots,A(x_{n+m-1})*B(x_{n+m-1})] \f]
   ///  -# use \f$ \mathbf{x} \f$ and \f$ \mathbf{y} \f$ to interpolate a polynomial 
   ///     (#Interpolate(Complex[], Complex[]))
   ///     \f[ C(x) = c_{n+m-1}* x^{m+n-1} + \ldots + c_1 * x + c_0; \f]
   /// 
   /// The convolution of \f$ \mathbf{a} \f$ and \f$ \mathbf{b} \f$ is 
   /// \f$ \mathbf{c}=[a_0,\ldots,c_{n+m-1}] \f$.  
   ///
   public class PolynomialConvolution : IConvolution
   {
      ///
      /// Computes the convolution.
      ///
      public virtual Complex[] Convolution(Complex[] a, Complex[] b)
      {
         Complex[] x  = GenerateDomain(a.Length + b.Length - 1);
         Complex[] ya = Evaluate(a,x);
         Complex[] yb = Evaluate(b,x);
         Complex[] y  = PointwiseMultiply(ya,yb);
         Complex[] c  = Interpolate(x,y);
         
         return c;
      }
      
      ///
      /// Generates domain for evaluation and interpolation. 
      ///
      /// Uses the roots of the Chebyshev polynomial of the first kind of degree <tt>n</tt> 
      /// on [0,1].
      /// 
      /// \param n the degree of the Chebyshev polynomial 
      /// \returns roots of canonical Chebyshev polynomial of first kind of degree  
      ///          <tt>n</tt>.
      ///
      protected internal Complex[] GenerateDomain(int n)
      {
         Complex[] res = new Complex[n];
         
         //Roots of the Chebyshev Polynoials of the first kind
         for (int i = 0; i < n; i++) {
            res[i] = Math.Cos((2*i+1)*Math.PI/(2*n+2));
         }
         
         return res;
      }
      
      ///
      /// Newton Interpolation.
      /// 
      /// Calculates the coefficients of a polynomial that takes values <tt>y</tt> 
      /// at points <tt>x</tt>.
      ///
      protected internal Complex[] Interpolate(Complex[] x, Complex[] y)
      {
         if (x.Length != y.Length) {
            throw new ArgumentException("The input vectors are of different length.");
         }
         
         int n = x.Length;
         
         Complex[] a; // resulting polynomial
         Complex[] d; // divided differences (denoted here [y_i,...,y_{i+k}]) 
         Complex[] p; // polynomial p(x)=(x-x_0)*(x-x_1)*...*(x-x_k) 
         
         // Initial values:
         a = new Complex[n];            // a(x) = 0 
         d = (Complex[]) y.Clone();     // d[i] = [y_i]  
         p = new Complex[n]; p[0] = 1;  // p(x) = 1 
         
         
         for (int k = 0; k < n; k++) {
            // Add next Newton combination to the solution
            // a(x) = a(x) + [y_0,...,y_k]*p(x) 
            for (int i = 0; i < n; i++) {
               a[i] = a[i] + d[0] * p[i];
            }
            
            // Calculate divided differences for the next step
            // Before: d[i] = [y_i,...,y_{i+k}]
            // Note: if k = n, no iteration is preformed.
            for (int i = 0; i < n-(k+1); i++) {
               d[i] = (d[i+1] - d[i])/(x[k+i+1] - x[i]);
            }
            // After: d[i] = [y_i,...,y_{i+k+1}] 
            
            // Calculate the p for next iteration.
            // Multiply p(x) by (x-x_k)
            // Before: p(x) = (x-x_0)*...*(x-x_k) 
            // Note: Execution during the last iteration of the main loop is redundant. 
            Complex c = 0; //carry value
            for (int i = 0; i < n; i++) {
               Complex tmp = p[i];
               p[i] = - x[k]*p[i] + c;
               c = tmp;
            }
            // After: p(x) = (x-x_0)*...*(x_{k+1})
         }
         
         return a;
      }
      
      ///
      /// Evaluates a polynomial at a sequence of points. 
      ///
      /// \param a the coefficients of the polynomial
      /// \param x the points of evaluation
      /// \returns vector of values the polynomial takes at each point in <tt>x</tt>
      ///
      protected internal Complex[] Evaluate(Complex[] a, Complex[] x)
      {
         Complex[] y = new Complex[x.Length];
         for (int i = 0; i < x.Length; i++) {
            y[i] = 0;
            for (int j = a.Length - 1; j >= 0; j--) {
               y[i] *= x[i];
               y[i] += a[j];
            }
         }
         
         return y;
      }

      ///
      /// Pointwise multiplication of two vectors. 
      ///
      protected internal Complex[] PointwiseMultiply(Complex[] a, Complex[] b)
      {
         if (a.Length != b.Length) {
            throw new ArgumentException("Both tables must have the same size.");
         }
           
         Complex[] res = new Complex[a.Length];
         
         for (int i = 0; i < a.Length; i++) {
            res[i] = a[i] * b[i];
         }
         
         return res;
      }

      ///
      /// \copybrief FFT::IConvolution::ToString()
      ///
      /// This method always returns <tt>"Pointwise"</tt>.
      public override string ToString()
      {
         return "Polynomial";
      }
      
   } //class PolynomialConvolution
   
} //namespace FFT

---