Polynomial Manipulations
************************

In Octave, a polynomial is represented by its coefficients (arranged in descending order).  For example, a vector C of length N+1 corresponds to the following polynomial of order N

     p(x) = C(1) x^N + ... + C(N) x + C(N+1).

Evaluating Polynomials
======================

The value of a polynomial represented by the vector C can be evaluated at the point X very easily, as the following example shows:

     N = length (c) - 1;
     val = dot (x.^(N:-1:0), c);

While the above example shows how easy it is to compute the value of a polynomial, it isn't the most stable algorithm.  With larger polynomials you should use more elegant algorithms, such as Horner's Method, which is exactly what the Octave function polyval does.

In the case where X is a square matrix, the polynomial given by C is still well-defined.  As when X is a scalar the obvious implementation is easily expressed in Octave, but also in this case more elegant algorithms perform better.  The polyvalm function provides such an algorithm.

 -- Y =: polyval (P, X)
 -- Y =: polyval (P, X, [], MU)
 -- [Y, DY] =: polyval (P, X, S)
 -- [Y, DY] =: polyval (P, X, S, MU)

     Evaluate the polynomial P at the specified values of X.

     If X is a vector or matrix, the polynomial is evaluated for each of the elements of X.

     When MU is present, evaluate the polynomial for (X - MU(1)) / MU(2).

     In addition to evaluating the polynomial, the second output represents the prediction interval, Y +/- DY, which contains at least 50% of the future predictions.  To calculate the prediction interval, the structured variable S, originating from polyfit, must be supplied.

     polyvalm, polyaffine, polyfit, roots, poly

 -- Y =: polyvalm (C, X)
     Evaluate a polynomial in the matrix sense.

     polyvalm (C, X) will evaluate the polynomial in the matrix sense, i.e., matrix multiplication is used instead of element by element multiplication as used in polyval.

     The argument X must be a square matrix. polyval, roots, poly

Finding Roots
=============

Octave can find the roots of a given polynomial.  This is done by computing the companion matrix of the polynomial (see the compan function for a definition), and then finding its eigenvalues.

 -- R =: roots (C)

     Compute the roots of the polynomial C.

     For a vector C with N components, return the roots of the polynomial

          c(1) * x^(N-1) + ... + c(N-1) * x + c(N)

     As an example, the following code finds the roots of the quadratic polynomial

          p(x) = x^2 - 5.

          c = [1, 0, -5];
          roots (c)
          =>  2.2361
          => -2.2361

     Note that the true result is +/- sqrt(5) which is roughly +/- 2.2361. poly, compan, fzero

 -- Z =: polyeig (C0, C1, ..., CL)
 -- [V, Z] =: polyeig (C0, C1, ..., CL)

     Solve the polynomial eigenvalue problem of degree L.

     Given an NxN matrix polynomial

     C(S) = C0 + C1 S + ... + CL S^L

     polyeig solves the eigenvalue problem

     (C0 + C1 Z + ... + CL Z^L) V = 0.

     Note that the eigenvalues Z are the zeros of the matrix polynomial. Z is a row vector with N*L elements.  V is a matrix (N x N*L) with columns that correspond to the eigenvectors.

     eig, eigs, compan

 -- A =: compan (C)
     Compute the companion matrix corresponding to polynomial coefficient vector C.

     The companion matrix is

               _                                                        _
              |  -c(2)/c(1)   -c(3)/c(1)  ...  -c(N)/c(1)  -c(N+1)/c(1)  |
              |       1            0      ...       0             0      |
              |       0            1      ...       0             0      |
          A = |       .            .      .         .             .      |
              |       .            .       .        .             .      |
              |       .            .        .       .             .      |
              |_      0            0      ...       1             0     _|

     The eigenvalues of the companion matrix are equal to the roots of the polynomial. roots, poly, eig

 -- [MULTP, IDXP] =: mpoles (P)
 -- [MULTP, IDXP] =: mpoles (P, TOL)
 -- [MULTP, IDXP] =: mpoles (P, TOL, REORDER)
     Identify unique poles in P and their associated multiplicity.

     By default, the output is ordered from the pole with the largest magnitude to the smallest magnitude.

     Two poles are considered to be multiples if the difference between them is less than the relative tolerance TOL.

          abs (P1 - P0) / abs (P0) < TOL

     If the pole is 0 then no scaling is done and TOL is interpreted as an absolute tolerance.  The default value for TOL is 0.001.

     If the optional parameter REORDER is false/zero, poles are not sorted.

     The output MULTP is a vector specifying the multiplicity of the poles. MULTP(n) refers to the multiplicity of the Nth pole P(IDXP(n)).

     For example:

          p = [2 3 1 1 2];
          [m, n] = mpoles (p)
             => m = [1; 1; 2; 1; 2]
             => n = [2; 5; 1; 4; 3]
             => p(n) = [3, 2, 2, 1, 1]

     residue, poly, roots, conv, deconv

Products of Polynomials
=======================

 -- Y =: conv (A, B)
 -- Y =: conv (A, B, SHAPE)
     Convolve two vectors A and B.

     When A and B are the coefficient vectors of two polynomials, the convolution represents the coefficient vector of the product polynomial.

     The size of the result is determined by the optional SHAPE argument which takes the following values

     SHAPE = "full"
          Return the full convolution.  (default) The result is a vector with length equal to length (A) + length (B) - 1.

     SHAPE = "same"
          Return the central part of the convolution with the same size as A.

     SHAPE = "valid"
          Return only the parts which do not include zero-padded edges. The size of the result is max (size (A) - size (B) + 1, 0).

     deconv, conv2, convn, fftconv

 -- C =: convn (A, B)
 -- C =: convn (A, B, SHAPE)
     Return the n-D convolution of A and B.

     The size of the result is determined by the optional SHAPE argument which takes the following values

     SHAPE = "full"
          Return the full convolution.  (default)

     SHAPE = "same"
          Return central part of the convolution with the same size as A. The central part of the convolution begins at the indices floor ([size(B)/2] + 1).

     SHAPE = "valid"
          Return only the parts which do not include zero-padded edges. The size of the result is max (size (A) - size (B) + 1, 0).

     conv2, conv

 -- B =: deconv (Y, A)
 -- [B, R] =: deconv (Y, A)
     Deconvolve two vectors (polynomial division).

     [B, R] = deconv (Y, A) solves for B and R such that Y = conv (A, B) + R.

     If Y and A are polynomial coefficient vectors, B will contain the coefficients of the polynomial quotient and R will be a remainder polynomial of lowest order. conv, residue

 -- C =: conv2 (A, B)
 -- C =: conv2 (V1, V2, M)
 -- C =: conv2 (..., SHAPE)
     Return the 2-D convolution of A and B.

     The size of the result is determined by the optional SHAPE argument which takes the following values

     SHAPE = "full"
          Return the full convolution.  (default)

     SHAPE = "same"
          Return the central part of the convolution with the same size as A. The central part of the convolution begins at the indices floor ([size(B)/2] + 1).

     SHAPE = "valid"
          Return only the parts which do not include zero-padded edges. The size of the result is max (size (A) - size (B) + 1, 0).

     When the third argument is a matrix, return the convolution of the matrix M by the vector V1 in the column direction and by the vector V2 in the row direction. conv, convn

 -- Q =: polygcd (B, A)
 -- Q =: polygcd (B, A, TOL)

     Find the greatest common divisor of two polynomials.

     This is equivalent to the polynomial found by multiplying together all the common roots.  Together with deconv, you can reduce a ratio of two polynomials.

     The tolerance TOL defaults to sqrt (eps).

     Caution: This is a numerically unstable algorithm and should not be used on large polynomials.

     Example code:

          polygcd (poly (1:8), poly (3:12)) - poly (3:8)
          => [ 0, 0, 0, 0, 0, 0, 0 ]
          deconv (poly (1:8), polygcd (poly (1:8), poly (3:12))) - poly (1:2)
          => [ 0, 0, 0 ]

     poly, roots, conv, deconv, residue

 -- [R, P, K, E] =: residue (B, A)
 -- [B, A] =: residue (R, P, K)
 -- [B, A] =: residue (R, P, K, E)
     The first calling form computes the partial fraction expansion for the quotient of the polynomials, B and A.

     The quotient is defined as

          B(s)    M       r(m)        N
          ---- = SUM ------------- + SUM k(i)*s^(N-i)
          A(s)   m=1 (s-p(m))^e(m)   i=1

     where M is the number of poles (the length of the R, P, and E), the K vector is a polynomial of order N-1 representing the direct contribution, and the E vector specifies the multiplicity of the m-th residue's pole.

     For example,

          b = [1, 1, 1];
          a = [1, -5, 8, -4];
          [r, p, k, e] = residue (b, a)
             => r = [-2; 7; 3]
             => p = [2; 2; 1]
             => k = [](0x0)
             => e = [1; 2; 1]

     which represents the following partial fraction expansion

                  s^2 + s + 1       -2        7        3
             ------------------- = ----- + ------- + -----
             s^3 - 5s^2 + 8s - 4   (s-2)   (s-2)^2   (s-1)

     The second calling form performs the inverse operation and computes the reconstituted quotient of polynomials, B(s)/A(s), from the partial fraction expansion; represented by the residues, poles, and a direct polynomial specified by R, P and K, and the pole multiplicity E.

     If the multiplicity, E, is not explicitly specified the multiplicity is determined by the function mpoles.

     For example:

          r = [-2; 7; 3];
          p = [2; 2; 1];
          k = [1, 0];
          [b, a] = residue (r, p, k)
             => b = [1, -5, 9, -3, 1]
             => a = [1, -5, 8, -4]

          where mpoles is used to determine e = [1; 2; 1]

     Alternatively the multiplicity may be defined explicitly, for example,

          r = [7; 3; -2];
          p = [2; 1; 2];
          k = [1, 0];
          e = [2; 1; 1];
          [b, a] = residue (r, p, k, e)
             => b = [1, -5, 9, -3, 1]
             => a = [1, -5, 8, -4]

     which represents the following partial fraction expansion

           -2        7        3         s^4 - 5s^3 + 9s^2 - 3s + 1
          ----- + ------- + ----- + s = --------------------------
          (s-2)   (s-2)^2   (s-1)          s^3 - 5s^2 + 8s - 4

     mpoles, poly, roots, conv, deconv

Derivatives / Integrals / Transforms
====================================

Octave comes with functions for computing the derivative and the integral of a polynomial.  The functions polyder and polyint both return new polynomials describing the result.  As an example we'll compute the definite integral of p(x) = x^2 + 1 from 0 to 3.

     c = [1, 0, 1];
     integral = polyint (c);
     area = polyval (integral, 3) - polyval (integral, 0)
     => 12

 -- K =: polyder (P)
 -- K =: polyder (A, B)
 -- [Q, D] =: polyder (B, A)
     Return the coefficients of the derivative of the polynomial whose coefficients are given by the vector P.

     If a pair of polynomials is given, return the derivative of the product A*B.

     If two inputs and two outputs are given, return the derivative of the polynomial quotient B/A.  The quotient numerator is in Q and the denominator in D. polyint, polyval, polyreduce

 -- Q =: polyint (P)
 -- Q =: polyint (P, K)
     Return the coefficients of the integral of the polynomial whose coefficients are represented by the vector P.

     The variable K is the constant of integration, which by default is set to zero. polyder, polyval

 -- G =: polyaffine (F, MU)
     Return the coefficients of the polynomial vector F after an affine transformation.

     If F is the vector representing the polynomial f(x), then G = polyaffine (F, MU) is the vector representing:

          g(x) = f( (x - MU(1)) / MU(2) )

     polyval, polyfit

Polynomial Interpolation
========================

Octave comes with good support for various kinds of interpolation, most of which are described in Interpolation.  One simple alternative to the functions described in the aforementioned chapter, is to fit a single polynomial, or a piecewise polynomial (spline) to some given data points.  To avoid a highly fluctuating polynomial, one most often wants to fit a low-order polynomial to data.  This usually means that it is necessary to fit the polynomial in a least-squares sense, which just is what the polyfit function does.

 -- P =: polyfit (X, Y, N)
 -- [P, S] =: polyfit (X, Y, N)
 -- [P, S, MU] =: polyfit (X, Y, N)
     Return the coefficients of a polynomial P(X) of degree N that minimizes the least-squares-error of the fit to the points [X(:), Y(:)].

     N is typically an integer ≥ 0 specifying the degree of the approximating polynomial.  If N is a logical vector, it is used as a mask to selectively force the corresponding polynomial coefficients to be used or ignored.

     The polynomial coefficients are returned in the row vector P.  The output P may be directly used with polyval to estimate values using the fitted polynomial.

     The optional output S is a structure containing the following fields:

     yf
          The values of the polynomial for each value of X.

     X
          The Vandermonde matrix used to compute the polynomial coefficients.

     R
          Triangular factor R from the QR decomposition.

     C
          The unscaled covariance matrix, formally equal to the inverse of X'*X, but computed in a way minimizing roundoff error propagation.

     df
          The degrees of freedom.

     normr
          The norm of the residuals.

     The second output may be used by polyval to calculate the statistical error limits of the predicted values.  In particular, the standard deviation of P coefficients is given by

     sqrt (diag (S.C)/S.DF) * S.NORMR.

     When the third output, MU, is present the original data is centered and scaled which can improve the numerical stability of the fit.  The coefficients P are associated with a polynomial in

     XHAT = (X - MU(1)) / MU(2)

     where MU(1) = mean (X), and MU(2) = std (X).

     Example 1 : logical N and integer N

          f = @(x) x.^2 + 5;   # data-generating function
          x = 0:5;
          y = f (x);
          ## Fit data to polynomial A*x^3 + B*x^1
          p = polyfit (x, y, logical ([1, 0, 1, 0]))
          => p = [ 0.0680, 0, 4.2444, 0 ]
          ## Fit data to polynomial using all terms up to x^3
          p = polyfit (x, y, 3)
          => p = [ -4.9608e-17, 1.0000e+00, -3.3813e-15, 5.0000e+00 ]

     polyval, polyaffine, roots, vander, zscore

In situations where a single polynomial isn't good enough, a solution is to use several polynomials pieced together.  The function splinefit fits a piecewise polynomial (spline) to a set of data.

 -- PP =: splinefit (X, Y, BREAKS)
 -- PP =: splinefit (X, Y, P)
 -- PP =: splinefit (..., "periodic", PERIODIC)
 -- PP =: splinefit (..., "robust", ROBUST)
 -- PP =: splinefit (..., "beta", BETA)
 -- PP =: splinefit (..., "order", ORDER)
 -- PP =: splinefit (..., "constraints", CONSTRAINTS)

     Fit a piecewise cubic spline with breaks (knots) BREAKS to the noisy data, X and Y.

     X is a vector, and Y is a vector or N-D array.  If Y is an N-D array, then X(j) is matched to Y(:,...,:,j).

     P is a positive integer defining the number of intervals along X, and P+1 is the number of breaks.  The number of points in each interval differ by no more than 1.

     The optional property PERIODIC is a logical value which specifies whether a periodic boundary condition is applied to the spline.  The length of the period is max (BREAKS) - min (BREAKS). The default value is false.

     The optional property ROBUST is a logical value which specifies if robust fitting is to be applied to reduce the influence of outlying data points.  Three iterations of weighted least squares are performed. Weights are computed from previous residuals.  The sensitivity of outlier identification is controlled by the property BETA.  The value of BETA is restricted to the range, 0 < BETA < 1.  The default value is BETA = 1/2.  Values close to 0 give all data equal weighting.  Increasing values of BETA reduce the influence of outlying data.  Values close to unity may cause instability or rank deficiency.

     The fitted spline is returned as a piecewise polynomial, PP, and may be evaluated using ppval.

     The splines are constructed of polynomials with degree ORDER. The default is a cubic, ORDER=3.  A spline with P pieces has P+ORDER degrees of freedom.  With periodic boundary conditions the degrees of freedom are reduced to P.

     The optional property, CONSTRAINTS, is a structure specifying linear constraints on the fit.  The structure has three fields, "xc", "yc", and "cc".

     "xc"
          Vector of the x-locations of the constraints.

     "yc"
          Constraining values at the locations XC. The default is an array of zeros.

     "cc"
          Coefficients (matrix).  The default is an array of ones.  The number of rows is limited to the order of the piecewise polynomials, ORDER.

     Constraints are linear combinations of derivatives of order 0 to ORDER-1 according to

          cc(1,j) * y(xc(j)) + cc(2,j) * y'(xc(j)) + ... = yc(:,...,:,j).

     interp1, unmkpp, ppval, spline, pchip, ppder, ppint, ppjumps

The number of BREAKS (or knots) used to construct the piecewise polynomial is a significant factor in suppressing the noise present in the input data, X and Y.  This is demonstrated by the example below.

     x = 2 * pi * rand (1, 200);
     y = sin (x) + sin (2 * x) + 0.2 * randn (size (x));
     ## Uniform breaks
     breaks = linspace (0, 2 * pi, 41); % 41 breaks, 40 pieces
     pp1 = splinefit (x, y, breaks);
     ## Breaks interpolated from data
     pp2 = splinefit (x, y, 10);  % 11 breaks, 10 pieces
     ## Plot
     xx = linspace (0, 2 * pi, 400);
     y1 = ppval (pp1, xx);
     y2 = ppval (pp2, xx);
     plot (x, y, ".", xx, [y1; y2])
     axis tight
     ylim auto
     legend ({"data", "41 breaks, 40 pieces", "11 breaks, 10 pieces"})

The piecewise polynomial fit, provided by splinefit, has continuous derivatives up to the ORDER-1.  For example, a cubic fit has continuous first and second derivatives.  This is demonstrated by the code

     ## Data (200 points)
     x = 2 * pi * rand (1, 200);
     y = sin (x) + sin (2 * x) + 0.1 * randn (size (x));
     ## Piecewise constant
     pp1 = splinefit (x, y, 8, "order", 0);
     ## Piecewise linear
     pp2 = splinefit (x, y, 8, "order", 1);
     ## Piecewise quadratic
     pp3 = splinefit (x, y, 8, "order", 2);
     ## Piecewise cubic
     pp4 = splinefit (x, y, 8, "order", 3);
     ## Piecewise quartic
     pp5 = splinefit (x, y, 8, "order", 4);
     ## Plot
     xx = linspace (0, 2 * pi, 400);
     y1 = ppval (pp1, xx);
     y2 = ppval (pp2, xx);
     y3 = ppval (pp3, xx);
     y4 = ppval (pp4, xx);
     y5 = ppval (pp5, xx);
     plot (x, y, ".", xx, [y1; y2; y3; y4; y5])
     axis tight
     ylim auto
     legend ({"data", "order 0", "order 1", "order 2", "order 3", "order 4"})

When the underlying function to provide a fit to is periodic, splinefit is able to apply the boundary conditions needed to manifest a periodic fit. This is demonstrated by the code below.

     ## Data (100 points)
     x = 2 * pi * [0, (rand (1, 98)), 1];
     y = sin (x) - cos (2 * x) + 0.2 * randn (size (x));
     ## No constraints
     pp1 = splinefit (x, y, 10, "order", 5);
     ## Periodic boundaries
     pp2 = splinefit (x, y, 10, "order", 5, "periodic", true);
     ## Plot
     xx = linspace (0, 2 * pi, 400);
     y1 = ppval (pp1, xx);
     y2 = ppval (pp2, xx);
     plot (x, y, ".", xx, [y1; y2])
     axis tight
     ylim auto
     legend ({"data", "no constraints", "periodic"})

More complex constraints may be added as well.  For example, the code below illustrates a periodic fit with values that have been clamped at the endpoints, and a second periodic fit which is hinged at the endpoints.

     ## Data (200 points)
     x = 2 * pi * rand (1, 200);
     y = sin (2 * x) + 0.1 * randn (size (x));
     ## Breaks
     breaks = linspace (0, 2 * pi, 10);
     ## Clamped endpoints, y = y' = 0
     xc = [0, 0, 2*pi, 2*pi];
     cc = [(eye (2)), (eye (2))];
     con = struct ("xc", xc, "cc", cc);
     pp1 = splinefit (x, y, breaks, "constraints", con);
     ## Hinged periodic endpoints, y = 0
     con = struct ("xc", 0);
     pp2 = splinefit (x, y, breaks, "constraints", con, "periodic", true);
     ## Plot
     xx = linspace (0, 2 * pi, 400);
     y1 = ppval (pp1, xx);
     y2 = ppval (pp2, xx);
     plot (x, y, ".", xx, [y1; y2])
     axis tight
     ylim auto
     legend ({"data", "clamped", "hinged periodic"})

The splinefit function also provides the convenience of a ROBUST fitting, where the effect of outlying data is reduced.  In the example below, three different fits are provided.  Two with differing levels of outlier suppression and a third illustrating the non-robust solution.

     ## Data
     x = linspace (0, 2*pi, 200);
     y = sin (x) + sin (2 * x) + 0.05 * randn (size (x));
     ## Add outliers
     x = [x, linspace(0,2*pi,60)];
     y = [y, -ones(1,60)];
     ## Fit splines with hinged conditions
     con = struct ("xc", [0, 2*pi]);
     ## Robust fitting, beta = 0.25
     pp1 = splinefit (x, y, 8, "constraints", con, "beta", 0.25);
     ## Robust fitting, beta = 0.75
     pp2 = splinefit (x, y, 8, "constraints", con, "beta", 0.75);
     ## No robust fitting
     pp3 = splinefit (x, y, 8, "constraints", con);
     ## Plot
     xx = linspace (0, 2*pi, 400);
     y1 = ppval (pp1, xx);
     y2 = ppval (pp2, xx);
     y3 = ppval (pp3, xx);
     plot (x, y, ".", xx, [y1; y2; y3])
     legend ({"data with outliers","robust, beta = 0.25", ...
              "robust, beta = 0.75", "no robust fitting"})
     axis tight
     ylim auto

A very specific form of polynomial interpretation is the Padé approximant. For control systems, a continuous-time delay can be modeled very simply with the approximant.

 -- [NUM, DEN] =: padecoef (T)
 -- [NUM, DEN] =: padecoef (T, N)
     Compute the Nth-order Padé approximant of the continuous-time delay T in transfer function form.

     The Padé approximant of exp (-sT) is defined by the following equation

                       Pn(s)
          exp (-sT) ~ -------
                       Qn(s)

     Where both Pn(s) and Qn(s) are Nth-order rational functions defined by the following expressions

                   N    (2N - k)!N!        k
          Pn(s) = SUM --------------- (-sT)
                  k=0 (2N)!k!(N - k)!

          Qn(s) = Pn(-s)

     The inputs T and N must be non-negative numeric scalars.  If N is unspecified it defaults to 1.

     The output row vectors NUM and DEN contain the numerator and denominator coefficients in descending powers of s.  Both are Nth-order polynomials.

     For example:

          t = 0.1;
          n = 4;
          [num, den] = padecoef (t, n)
          => num =

                1.0000e-04  -2.0000e-02   1.8000e+00  -8.4000e+01   1.6800e+03

          => den =

                1.0000e-04   2.0000e-02   1.8000e+00   8.4000e+01   1.6800e+03

The function, ppval, evaluates the piecewise polynomials, created by mkpp or other means, and unmkpp returns detailed information about the piecewise polynomial.

The following example shows how to combine two linear functions and a quadratic into one function.  Each of these functions is expressed on adjoined intervals.

     x = [-2, -1, 1, 2];
     p = [ 0,  1, 0;
           1, -2, 1;
           0, -1, 1 ];
     pp = mkpp (x, p);
     xi = linspace (-2, 2, 50);
     yi = ppval (pp, xi);
     plot (xi, yi);

 -- PP =: mkpp (BREAKS, COEFS)
 -- PP =: mkpp (BREAKS, COEFS, D)

     Construct a piecewise polynomial (pp) structure from sample points BREAKS and coefficients COEFS.

     BREAKS must be a vector of strictly increasing values.  The number of intervals is given by NI = length (BREAKS) - 1.

     When M is the polynomial order COEFS must be of size: NI-by-(M + 1).

     The i-th row of COEFS, COEFS(I,:), contains the coefficients for the polynomial over the I-th interval, ordered from highest (M) to lowest (0) degree.

     COEFS may also be a multi-dimensional array, specifying a vector-valued or array-valued polynomial.  In that case the polynomial order M is defined by the length of the last dimension of COEFS. The size of first dimension(s) are given by the scalar or vector D. If D is not given it is set to 1.  In this case P(R, I, :) contains the coefficients for the R-th polynomial defined on interval I.  In any case COEFS is reshaped to a 2-D matrix of size [NI*prod(D) M].

     Programming Note: ppval evaluates polynomials at XI - BREAKS(i), i.e., it subtracts the lower endpoint of the current interval from XI.  This must be taken into account when creating piecewise polynomials objects with mkpp. unmkpp, ppval, spline, pchip, ppder, ppint, ppjumps

 -- [X, P, N, K, D] =: unmkpp (PP)

     Extract the components of a piecewise polynomial structure PP.

     This function is the inverse of mkpp: it extracts the inputs to mkpp needed to create the piecewise polynomial structure PP. The code below makes this relation explicit:

          [breaks, coefs, numinter, order, dim] = unmkpp (pp);
          pp2  = mkpp (breaks, coefs, dim);

     The piecewise polynomial structure pp2 obtained in this way, is identical to the original pp.  The same can be obtained by directly accessing the fields of the structure pp.

     The components are:

     X
          Sample points or breaks.

     P
          Polynomial coefficients for points in sample interval. P(I, :) contains the coefficients for the polynomial over interval I ordered from highest to lowest degree. If D > 1, then P is a matrix of size [N*prod(D) M], where the I + (1:D) rows are the coefficients of all the D polynomials in the interval I.

     N
          Number of polynomial pieces or intervals, N = length (X) - 1.

     K
          Order of the polynomial plus 1.

     D
          Number of polynomials defined for each interval.

     mkpp, ppval, spline, pchip

 -- YI =: ppval (PP, XI)
     Evaluate the piecewise polynomial structure PP at the points XI.

     If PP describes a scalar polynomial function, the result is an array of the same shape as XI.  Otherwise, the size of the result is [pp.dim, length(XI)] if XI is a vector, or [pp.dim, size(XI)] if it is a multi-dimensional array. mkpp, unmkpp, spline, pchip

 -- ppd =: ppder (pp)
 -- ppd =: ppder (pp, m)
     Compute the piecewise M-th derivative of a piecewise polynomial struct PP.

     If M is omitted the first derivative is calculated. mkpp, ppval, ppint

 -- PPI =: ppint (PP)
 -- PPI =: ppint (PP, C)
     Compute the integral of the piecewise polynomial struct PP.

     C, if given, is the constant of integration. mkpp, ppval, ppder

 -- JUMPS =: ppjumps (PP)
     Evaluate the boundary jumps of a piecewise polynomial.

     If there are n intervals, and the dimensionality of PP is d, the resulting array has dimensions [d, n-1]. mkpp

Miscellaneous Functions
=======================

 -- Y =: poly (A)
 -- Y =: poly (X)
     If A is a square N-by-N matrix, poly (A) is the row vector of the coefficients of det (z * eye (N) - A), the characteristic polynomial of A.

     For example, the following code finds the eigenvalues of A which are the roots of poly (A).

          roots (poly (eye (3)))
              => 1.00001 + 0.00001i
                 1.00001 - 0.00001i
                 0.99999 + 0.00000i

     In fact, all three eigenvalues are exactly 1 which emphasizes that for numerical performance the eig function should be used to compute eigenvalues.

     If X is a vector, poly (X) is a vector of the coefficients of the polynomial whose roots are the elements of X. That is, if C is a polynomial, then the elements of D = roots (poly (C)) are contained in C.  The vectors C and D are not identical, however, due to sorting and numerical errors. roots, eig

 -- polyout: (C)
 -- polyout: (C, X)
 -- STR =: polyout (...)
     Display a formatted version of the polynomial C.

     The formatted polynomial

          c(x) = c(1) * x^n + ... + c(n) x + c(n+1)

     is returned as a string or written to the screen if nargout is zero.

     The second argument X specifies the variable name to use for each term and defaults to the string "s". polyreduce

 -- P =: polyreduce (C)
     Reduce a polynomial coefficient vector to a minimum number of terms by stripping off any leading zeros. polyout