CppAD: A C++ Algorithmic Differentiation Package 20110419
poly.hpp
Go to the documentation of this file.
00001 /* $Id$ */
00002 # ifndef CPPAD_POLY_INCLUDED
00003 # define CPPAD_POLY_INCLUDED
00004 
00005 /* --------------------------------------------------------------------------
00006 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-07 Bradley M. Bell
00007 
00008 CppAD is distributed under multiple licenses. This distribution is under
00009 the terms of the 
00010                     Common Public License Version 1.0.
00011 
00012 A copy of this license is included in the COPYING file of this distribution.
00013 Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
00014 -------------------------------------------------------------------------- */
00015 /*
00016 $begin Poly$$
00017 $spell
00018         cppad.hpp
00019         CppAD
00020         namespace
00021         cstddef
00022         ifndef
00023         endif
00024         deg
00025         const
00026         std
00027         da
00028 $$
00029 
00030 $index Poly$$
00031 $index polynomial$$
00032 $index derivative, polynomial template$$
00033 $index template, polynomial derivative$$
00034 
00035 $section Evaluate a Polynomial or its Derivative$$
00036 
00037 $head Syntax$$
00038 $code # include <cppad/poly.hpp>$$
00039 $pre
00040 $$
00041 $syntax%%p% = Poly(%k%, %a%, %z%)%$$
00042 
00043 
00044 $head Description$$
00045 Computes the $th k$$ derivative of the polynomial 
00046 $latex \[
00047         P(z) = a_0 + a_1 z^1 + \cdots + a_d z^d
00048 \] $$
00049 If $italic k$$ is equal to zero, the return value is $latex P(z)$$.
00050 
00051 $head Include$$
00052 The file $code cppad/poly.hpp$$ is included by $code cppad/cppad.hpp$$
00053 but it can also be included separately with out the rest of 
00054 the $code CppAD$$ routines.
00055 Including this file defines
00056 $code Poly$$ within the $code CppAD$$ namespace.
00057 
00058 $head k$$
00059 The argument $italic k$$ has prototype
00060 $syntax%
00061         size_t %k%
00062 %$$
00063 It specifies the order of the derivative to calculate.
00064 
00065 $head a$$
00066 The argument $italic a$$ has prototype
00067 $syntax%
00068         const %Vector% &%a%
00069 %$$
00070 (see $xref/Poly/Vector/Vector/$$ below).
00071 It specifies the vector corresponding to the polynomial $latex P(z)$$.
00072 
00073 $head z$$
00074 The argument $italic z$$ has prototype
00075 $syntax%
00076         const %Type% &%z%
00077 %$$
00078 (see $italic Type$$ below).
00079 It specifies the point at which to evaluate the polynomial
00080 
00081 $head p$$
00082 The result $italic p$$  has prototype
00083 $syntax%
00084         %Type% %p%
00085 %$$
00086 (see $xref/Poly/Type/Type/$$ below)
00087 and it is equal to the $th k$$ derivative of $latex P(z)$$; i.e., 
00088 $latex \[
00089 p = \frac{k !}{0 !} a_k 
00090   + \frac{(k+1) !}{1 !} a_{k+1} z^1 
00091   + \ldots
00092   + \frac{d !}{(d - k) !} a_d z^{d - k}
00093 \]
00094 $$
00095 If $latex k > d$$, $syntax%%p% = %Type%(0)%$$.
00096 
00097 $head Type$$
00098 The type $italic Type$$ is determined by the argument $italic z$$.
00099 It is assumed that
00100 multiplication and addition of $italic Type$$ objects
00101 are commutative.
00102 
00103 $subhead Operations$$
00104 The following operations must be supported where
00105 $italic x$$ and $italic y$$ are objects of type $italic Type$$
00106 and $italic i$$ is an $code int$$:
00107 $table
00108 $syntax%%x%  = %i%$$   $cnext assignment     $rnext
00109 $syntax%%x%  = %y%$$   $cnext assignment     $rnext
00110 $syntax%%x% *= %y%$$   $cnext multiplication computed assignment $rnext
00111 $syntax%%x% += %y%$$   $cnext addition computed assignment
00112 
00113 $tend
00114 
00115 
00116 $head Vector$$
00117 The type $italic Vector$$ must be a $xref/SimpleVector/$$ class with
00118 $xref/SimpleVector/Elements of Specified Type/elements of type/$$
00119 $italic Type$$.
00120 The routine $xref/CheckSimpleVector/$$ will generate an error message
00121 if this is not the case.
00122 
00123 $head Operation Sequence$$
00124 The $italic Type$$ operation sequence used to calculate $italic p$$ is 
00125 $xref/glossary/Operation/Independent/independent/1/$$
00126 of $italic z$$ and the elements of $italic a$$
00127 (it does depend on the size of the vector $italic a$$).
00128 
00129 
00130 $children%
00131         example/poly.cpp%
00132         omh/poly_hpp.omh
00133 %$$
00134 
00135 $head Example$$
00136 The file
00137 $xref/Poly.cpp/$$
00138 contains an example and test of this routine.
00139 It returns true if it succeeds and false otherwise.
00140 
00141 $head Source$$
00142 The file $cref/poly.hpp/$$ contains the 
00143 current source code that implements these specifications.
00144 
00145 $end
00146 ------------------------------------------------------------------------------
00147 */
00148 // BEGIN PROGRAM
00149 # include <cstddef>  // used to defined size_t
00150 # include <cppad/check_simple_vector.hpp>
00151 
00152 namespace CppAD {    // BEGIN CppAD namespace
00153 
00154 template <class Type, class Vector>
00155 Type Poly(size_t k, const Vector &a, const Type &z)
00156 {       size_t i;
00157         size_t d = a.size() - 1;
00158 
00159         Type tmp;
00160 
00161         // check Vector is Simple Vector class with Type elements
00162         CheckSimpleVector<Type, Vector>();
00163 
00164         // case where derivative order greater than degree of polynomial
00165         if( k > d )
00166         {       tmp = 0;
00167                 return tmp;
00168         }
00169         // case where we are evaluating a derivative
00170         if( k > 0 )
00171         {       // initialize factor as (k-1) !
00172                 size_t factor = 1;
00173                 for(i = 2; i < k; i++)
00174                         factor *= i;
00175 
00176                 // set b to coefficient vector corresponding to derivative
00177                 Vector b(d - k + 1);
00178                 for(i = k; i <= d; i++)
00179                 {       factor   *= i;
00180                         tmp       = factor;
00181                         b[i - k]  = a[i] * tmp; 
00182                         factor   /= (i - k + 1);
00183                 }
00184                 // value of derivative polynomial
00185                 return Poly(0, b, z);
00186         }
00187         // case where we are evaluating the original polynomial
00188         Type sum = a[d];
00189         i        = d;
00190         while(i > 0)
00191         {       sum *= z;
00192                 sum += a[--i];
00193         }
00194         return sum;
00195 }
00196 } // END CppAD namespace
00197 // END PROGRAM
00198 # endif