Prev Next cppad_poly.cpp

CppAD Speed: Second Derivative of a Polynomial

link_poly
Routine that computes the second derivative of a polynomial using CppAD:
# include <cppad/cppad.hpp>
# include <cppad/speed/uniform_01.hpp>

bool link_poly(
     size_t                     size     , 
     size_t                     repeat   , 
     bool                       retape   ,
     CppAD::vector<double>     &a        ,  // coefficients of polynomial
     CppAD::vector<double>     &z        ,  // polynomial argument value
     CppAD::vector<double>     &ddp      )  // second derivative w.r.t z  
{
     // -----------------------------------------------------
     // setup
     typedef CppAD::AD<double>     ADScalar; 
     typedef CppAD::vector<ADScalar> ADVector; 

     size_t i;      // temporary index
     size_t m = 1;  // number of dependent variables
     size_t n = 1;  // number of independent variables
     ADVector Z(n); // AD domain space vector
     ADVector P(m); // AD range space vector

     // choose the polynomial coefficients
     CppAD::uniform_01(size, a);

     // AD copy of the polynomial coefficients
     ADVector A(size);
     for(i = 0; i < size; i++)
          A[i] = a[i];

     // forward mode first and second differentials
     CppAD::vector<double> p(1), dp(1), dz(1), ddz(1);
     dz[0]  = 1.;
     ddz[0] = 0.;

     CppAD::ADFun<double> f;

     if( retape ) while(repeat--)
     {
          // choose an argument value
          CppAD::uniform_01(1, z);
          Z[0] = z[0];

          // declare independent variables
          Independent(Z);

          // AD computation of the function value 
          P[0] = CppAD::Poly(0, A, Z[0]);

          // create function object f : A -> detA
          f.Dependent(Z, P);

          // pre-allocate memory for three forward mode calculations
          f.capacity_taylor(3);

          // get the next argument value
          CppAD::uniform_01(1, z);

          // evaluate the polynomial at the new argument value
          p = f.Forward(0, z);

          // evaluate first order Taylor coefficient
          dp = f.Forward(1, dz);

          // second derivative is twice second order Taylor coef
          ddp     = f.Forward(2, ddz);
          ddp[0] *= 2.;
     }
     else
     {
          // choose an argument value
          CppAD::uniform_01(1, z);
          Z[0] = z[0];

          // declare independent variables
          Independent(Z);

          // AD computation of the function value 
          P[0] = CppAD::Poly(0, A, Z[0]);

          // create function object f : A -> detA
          f.Dependent(Z, P);

          while(repeat--)
          {    // sufficient memory is allocated by second repetition

               // get the next argument value
               CppAD::uniform_01(1, z);

               // evaluate the polynomial at the new argument value
               p = f.Forward(0, z);

               // evaluate first order Taylor coefficient
               dp = f.Forward(1, dz);

               // second derivative is twice second order Taylor coef
               ddp     = f.Forward(2, ddz);
               ddp[0] *= 2.;
          }
     }
     return true;
}

Input File: speed/cppad/poly.cpp