#include <iostream.h>
#include <math.h>
#include "feed.h"
// ----------- constructors * destructors * equality -----------
// Default constructor
CFeed::CFeed(unsigned long N, double f, double *df, double *ddf)
{
unsigned long i, j; // counters
m_N = N; // store values of size
m_f = f; // and function
m_df = m_ddf =0; // zero pointers
if(m_N == 0) return; // zero length class
m_df = new double [m_N]; // allocate memory for 1st
m_ddf = new double [m_N*m_N]; // order partials and second
// order
for(i = 0; i < m_N; i++) // copy derivatives
{
if(df != 0)
m_df[i] = df[i]; // copy first order partials
else
m_df[i] = 0;
for(j = 0; j < m_N; j++)
{
if(ddf != 0) m_ddf[i*m_N+j] = ddf[i*N+j];
else m_ddf[i*m_N+j] = 0; // copy 2nd order partials
}
}
}
// destructor
CFeed::~CFeed(void)
{
m_f = 0;
m_N = 0; // zero length
if(m_df != 0) delete [] m_df; m_df = 0; // first partials
if(m_ddf != 0) delete [] m_ddf; m_ddf = 0; // second partials
}
// equal operator
CFeed & CFeed::operator=(const CFeed &feed)
{
unsigned long i,j; // counters
if(this == &feed) return(*this); // itself
if(m_df != 0) delete [] m_df; m_df = 0; // free memory
if(m_ddf != 0) delete [] m_ddf; m_ddf = 0; // and zero
// pointers
m_N = feed.m_N; // store values of size
m_f = feed.m_f; // and function
if(m_N == 0) return(*this); // done; zero length class
m_df = new double [m_N]; // allocate memory for 1st and
m_ddf = new double [m_N*m_N]; // 2nd order partials
for(i = 0; i < m_N; i++)
{
m_df[i] = feed.m_df[i]; // copy first order partials
for(j = 0; j < m_N; j++) // copy second order partials
m_ddf[i*m_N+j] = feed.m_ddf[i*feed.m_N+j];
}
return(*this); // return class
}
// ----------------- Multiplication operators -------------------
// multiplication operator
CFeed CFeed::operator*(const CFeed &feed)
{
CFeed product(*this); // initialize product to 'this'
product *= feed; // and multiply by 'feed'
return(product); // return result
}
// multiplication equal operator
CFeed & CFeed::operator*=(const CFeed &feed)
{
unsigned long i,j; // counters
CFeed THIS(*this); // make copy of 'this'
if(m_N < feed.m_N) // 'this' smaller so set 'this' to
*this = feed; // 'feed' to increase its size
for(i = 0; i < m_N; i++)
{
for(j = 0; j < m_N; j++) // multiply 2nd partials:
{ // chain rule
m_ddf[i*m_N+j] = 0;
if((i < feed.m_N) && (j < feed.m_N))
m_ddf[i*m_N+j] += THIS.m_f*feed.m_ddf[i*feed.m_N+j];
if((i < THIS.m_N) && (j < feed.m_N))
m_ddf[i*m_N+j] += THIS.m_df[i]*feed.m_df[j];
if((j < THIS.m_N) && (i < feed.m_N))
m_ddf[i*m_N+j] += THIS.m_df[j]*feed.m_df[i];
if((i < THIS.m_N) && (j < THIS.m_N))
m_ddf[i*m_N+j] += THIS.m_ddf[i*THIS.m_N+j]*feed.m_f;
}
}
for(j = 0; j < m_N; j++) // multiply
{
m_df[j] = 0;
if( j < feed.m_N)
m_df[j] += THIS.m_f*feed.m_df[j]; // 1st partials:
// chain rule
if( j < THIS.m_N)
m_df[j] += feed.m_f*THIS.m_df[j];
}
m_f = feed.m_f*THIS.m_f; // multiply bases
return(*this);
}
// multiplication operator
CFeed CFeed::operator*(const double c) // real constant
{
CFeed product(*this); // initialize product to 'this'
product *= c; // and multiply by constant
return(product); // return result
}
// multiplication equal operator
CFeed & CFeed::operator*=(const double c) // real constant
{
unsigned long i,j; // counters
m_f = c*m_f; // multiply base
for(i = 0; i < m_N; i++)
{
m_df[i] = c*m_df[i]; // 1st partials
for(j = 0; j < m_N; j++)
m_ddf[i*m_N+j] = c*m_ddf[i*m_N+j]; // and second
// partials
}
return(*this);
}
// --------------------------------------------------------------
// call operator() * composite function chain rule
CFeed CFeed::operator()(const CFeed &feed)
{
CFeed composite(feed.m_N); // define composite
unsigned long i,j; // counters
if(m_N == 0) return(*this); // constant
for(i = 0; i < composite.m_N; i++)
{
for(j = 0; j < composite.m_N; j++) // 2nd partials
composite.m_ddf[i*composite.m_N+j] =
feed.m_df[i]*feed.m_df[j]*m_ddf[0] +
feed.m_ddf[i*feed.m_N+j]*m_df[0];
}
for(j = 0; j < composite.m_N; j++) // 1st partials
composite.m_df[j] = feed.m_df[j]*m_df[0];
composite.m_f = m_f; // base
return(composite);
}
// -------------------- I/O operators and functions ------------
//
double CFeed::operator()(void) {return(m_f);} // base value
double CFeed::operator()(unsigned long i)
{
if(i > m_N) return(0); // band index
else if(i > 0) return(m_df[i-1]); // first partial
else return(0); // bad index
}
double CFeed::operator()(unsigned long i, unsigned long j)
{
if((i > m_N) || (j > m_N))
return(0); // band index
else if((i > 0) && (j > 0))
return(m_ddf[(i-1)*m_N+(j-1)]); // second partial
else
return(0); // bad index
}
// ------------------- NON-Member operators ---------------------
CFeed operator*(const double c, const CFeed &feed)
{
CFeed product(feed); // initialize product to 'feed'
product *= c; // multiply by constant
return(product); // return result
}
// ----------------- Math functions (friends) -------------------
// exponential fuction
CFeed exp(const CFeed &feed)
{
double df[1] = {exp(feed.m_f)}; // first derivative
double ddf[1] = {exp(feed.m_f)}; // second derivative
CFeed e(1,exp(feed.m_f),df,ddf); // e - feed class
return(e(feed)); // chain rule; return
}
// log function
CFeed log(const CFeed &feed)
{
double x = feed.m_f;
double df[1] = {1.0/x}; // first derivative
double ddf[1] = {-1.0/(x*x)}; // second derivative
CFeed ln(1,log(x),df,ddf); // ln - feed class
return(ln(feed)); // log; return
}