tpsa.cc
/*
* Copyright (C) 2007 by Lingyun Yang
* lingyun (.dot.] yang@gmail.com
*
* This program is free software; you can redistribute it and/or modify it
* under the terms of the GNU General Public License as published by the
* Free Software Foundation; either version 2 of the License, or (at your
* option) any later version.
*/
//
// enscript tpsa.cc -Ecpp --color --language=html -o ../tpsa.html
//
//! \brief A TPSA lib
//! \file tpsa.cc
//! \author Lingyun Yang, lingyun ([dot].) yang@gmail.com
//!
//!
#ifndef TPSA_HPP
#define TPSA_HPP
#include <utility>
#include <iostream>
#include <cstdlib>
namespace TPSA {
//! maximum size of pre-generated binomial coefficient
const size_t C_SIZE = 13;
//! Pre generated binomial coefficient
//! \[ C_{ij}=i!/(j!(i-j)!) \]
static const int C[13][13] = {
{1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{1, 2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{1, 3, 3, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{1, 4, 6, 4, 1, 0, 0, 0, 0, 0, 0, 0, 0},
{1, 5, 10, 10, 5, 1, 0, 0, 0, 0, 0, 0, 0},
{1, 6, 15, 20, 15, 6, 1, 0, 0, 0, 0, 0, 0},
{1, 7, 21, 35, 35, 21, 7, 1, 0, 0, 0, 0, 0},
{1, 8, 28, 56, 70, 56, 28, 8, 1, 0, 0, 0, 0},
{1, 9, 36, 84, 126, 126, 84, 36, 9, 1, 0, 0, 0},
{1, 10, 45, 120, 210, 252, 210, 120, 45, 10, 1, 0, 0},
{1, 11, 55, 165, 330, 462, 462, 330, 165, 55, 11, 1, 0},
{1, 12, 66, 220, 495, 792, 924, 792, 495, 220, 66, 12, 1}
};
//! \brief Vector for Truncated Power Series Algebra.
//!
//! To calculate the derivative of a function $f=x*x-1$ at $x=1$ up to N order,
//! declare f as a Vector : Vector f(N), x(N); x=1.0; f = x*x - 1;
//! then f.df(i) is the ith order of derivative at $x=1.0$
//!
template<class T> class Vector{
T* v;
//! dimension of vector
size_t size;
//! order of the vector, = size - 1.
size_t order;
//! memoty copy
void copy_from(const Vector& src) {
memcpy(v, src.v, src.size*sizeof(T));
}
public:
//! a vector with order n
explicit Vector(const int n) {
if (n >= C_SIZE) {
std::cerr << "Vector maximum dimension reached ( MAX= "
<< C_SIZE << " ) "
<< " " << n
<< std::endl;
std::exit(-1);
}
v = new T[n+1];
for (size_t i = 0; i < n; ++i) {
v[i] = T(0);
}
size = n+1;
order = n;
}
//! a vector with value x, order n
Vector(const T& x, const size_t n) {
v = new T[n+1];
v[0] = x;
v[1] = 1;
size = n + 1;
order = n;
}
Vector(const Vector& rhs) {
v = new T[rhs.size];
for (size_t i = 0; i < rhs.size; ++i) {
v[i] = rhs.v[i];
}
size = rhs.size;
order = rhs.order;
}
~Vector() {
delete [] v;
}
//! ith value in the vector
T operator() (size_t i) const { return v[i]; }
//! ith order derivative
T df(size_t i) { return v[i]; };
Vector& operator= (const Vector& rhs) {
if (rhs.v == v)
return *this;
#ifdef TPSA_CHECK_BOUNDARY
if (v->size == rhs.v->size)
#endif
copy_from(rhs);
}
Vector& operator= (const T& rhs) {
v[0] = rhs;
v[1] = T(1);
for (size_t i = 2; i < size; ++i)
v[i] = T(0);
}
Vector& operator+= (const Vector<T>& rhs) {
for (size_t i = 0; i < size; ++i)
v[i] += rhs.v[i];
}
Vector& operator-= (const Vector<T>& rhs) {
for (size_t i = 0; i < size; ++i)
v[i] -= rhs.v[i];
}
Vector& operator*= (const Vector<T>& rhs) {
T *t = new T[size];
for (size_t i = 0; i < size; ++i) {
t[i] = T(0);
for (size_t k = 0; k <= i; ++k) {
t[i] += C[i][k]*v[i]*rhs.v[i-k];
}
}
for (size_t i = 0; i < size; ++i)
v[i] = t[i];
delete [] t;
}
Vector& operator+= (const T& rhs) {
v[0] += rhs;
}
Vector& operator-= (const T& rhs) {
v[0] -= rhs;
}
Vector& operator*= (const T& rhs) {
v[0] *= rhs;
}
Vector& operator/= (const T& rhs) {
v[0] /= rhs;
}
Vector operator+(const Vector<T>& rhs) const{
Vector<T> ret(rhs.order);
for (size_t i = 0; i < size; ++i) {
ret.v[i] = v[i] + rhs.v[i];
}
return ret;
}
Vector operator-(const Vector<T>& rhs) const{
Vector<T> ret(rhs.order);
for (size_t i = 0; i < size; ++i) {
ret.v[i] = v[i] - rhs.v[i];
}
return ret;
}
Vector operator*(const Vector<T>& rhs) const{
Vector<T> ret(rhs.order);
for (size_t i = 0; i < size; ++i) {
ret.v[i] = T(0);
for (size_t k = 0; k <= i; ++k) {
ret.v[i] += C[i][k]*v[k]*rhs.v[i-k];
}
}
return ret;
}
//! \brief Get dimension as a pair
size_t len() const {
return size;
}
//! \brief I/O
template<class T2> friend std::ostream& operator<<(std::ostream& os, const Vector<T2>& V);
template<typename T2> friend Vector<T2> operator+(const T2&x, const Vector<T2>& v);
template<typename T2> friend Vector<T2> operator-(const T2&x, const Vector<T2>& v);
template<typename T2> friend Vector<T2> operator*(const T2&x, const Vector<T2>& v);
template<typename T2> friend Vector<T2> operator/(const T2&x, const Vector<T2>& v);
};
template<class T>
Vector<T> operator+(const T& x, const Vector<T>& v) {
Vector<T> t(v.order);
t.v[0] = x + v(0);
for (size_t i = 1; i < v.size; ++i) {
t.v[i] = v.v[i];
}
return t;
}
template<class T>
Vector<T> operator-(const T& x, const Vector<T>& v) {
Vector<T> t(v.order);
t.v[0] = x - v(0);
for (size_t i = 1; i < v.size; ++i) {
t.v[i] = -v.v[i];
}
return t;
}
template<class T>
Vector<T> operator*(const T& x, const Vector<T>& v) {
Vector<T> t(v.order);
T b(x);
for (size_t i = 0; i < v.size; ++i) {
t.v[i] = x*v.v[i];
}
return t;
}
template<class T>
Vector<T> operator/(const T& x, const Vector<T>& v) {
Vector<T> t(v.order);
T b;
t.v[0] = x/v.v[0];
for (size_t i = 1; i < v.size; ++i) {
b = T(0);
for (size_t j = 0; j < i; ++j) {
b -= C[i][j]*v.v[i-j]*t.v[j];
}
t.v[i] = b/v.v[0];
}
return t;
}
template<class T>
std::ostream& operator<< (std::ostream& os, const Vector<T>& V)
{
os << "(";
for (size_t i = 0; i < V.size-1; ++i) {
os << V.v[i] << ' ';
}
os << V.v[V.size-1] << ")";
return os;
}
}
#endif
#ifdef TEST
using namespace TPSA;
using namespace std;
int main(int argc, char** argv) {
Vector<double> x(3), y(3), z(3);
x = 2.0;
cout << "x=" << x << endl;
y = 3.0;
x = y;
x += 3.0;
x -= 3.0;
x *= 3.0;
x /= 3.1;
x = 2.0;
cout << "x=" << x << endl;
z = 1.0/x;
cout << "z=1.0/x=" << z << endl;
x = 1.0;
z = x*x*x+2.0*x;
z = x*x*x-2.0*x;
cout << "z=x^3+2x=" << z << endl;
//x += y;
//x -= y;
//x *= y;
//x /= y;
//z = x + y;
//z = x - y;
//z = x * y;
//z = x / y;
//z = exp(x);
z = 3.0+3.0+x;
z = 3.0 -x;
cout << z << endl;
return 0;
}
#endif
Generated by GNU enscript 1.6.4.