// Rascal, the Advanced Scientific CALculator
// Copyright (C) 2001, Sebastian Ritterbusch (Rascal@Ritterbusch.de)
//
// 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.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more detauls.
//
// You should have received a copy of the GNU General Public License
// along with this program; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
//
#ifndef TAYLOR_HPP_INCLUDED
#define TAYLOR_HPP_INCLUDED
#include <iostream.h>
template <class T>
class mtaylor
{
private:
int n;
T *a;
public:
mtaylor(void) : n(0),a(0) {}
mtaylor(T x,int iN=1) :a(new T[iN])
{
if(a)
n=iN;
else
n=0;
int i;
if(n)
a[0]=x;
for(i=1;i<n;i++)
a[i]=0;
}
mtaylor(T x,T y,int iN=2) : a(new T[iN])
{
if(a)
n=iN;
else
n=0;
int i;
if(n)
a[0]=x;
if(n>1)
a[1]=y;
for(i=2;i<n;i++)
a[i]=0;
}
mtaylor(const mtaylor<T> & b) : a(new T[b.n])
{
if(a)
n=b.n;
else
n=0;
int i;
for(i=0;i<n;i++)
a[i]=b.a[i];
}
~mtaylor(void) { if(a) delete [] a; }
mtaylor<T> & operator =(const mtaylor<T> & b)
{
T * new_a=new T[b.n];
n=new_a?b.n:0;
int i;
for(i=0;i<n;i++)
new_a[i]=b.a[i];
delete [] a;
a=new_a;
return *this;
}
void minsize(int nn)
{
if(nn<=n)
return;
T * new_a=new T[nn];
int i;
for(i=0;i<n;i++)
new_a[i]=a[i];
for(;i<nn;i++)
new_a[i]=0;
delete [] a;
a=new_a;
n=new_a?nn:0;
}
friend inline mtaylor<T> & operator +=(mtaylor<T> & a,const mtaylor<T> & b) { int i; a.minsize(b.n); for(i=0;i<a.n && i<b.n;i++) a.a[i]=a.a[i]+b.a[i]; return a;}
friend inline mtaylor<T> & operator +=(mtaylor<T> & a,const T & b) { if(a.n) a.a[0]=a.a[0]+b; return a;}
friend inline mtaylor<T> operator +(mtaylor<T> a,const mtaylor<T> & b) { return a+=b; }
friend inline mtaylor<T> operator +(mtaylor<T> a,const T & b) { return a+=b; }
friend inline mtaylor<T> operator +(const T & b,mtaylor<T> a) { return a+=b; }
friend inline mtaylor<T> operator -(mtaylor<T> a) { int i; for(i=0;i<a.n;i++) a.a[i]=-a.a[i]; return a;}
friend inline mtaylor<T> & operator -=(mtaylor<T> & a,const mtaylor<T> & b) { int i; a.minsize(b.n);for(i=0;i<a.n && i<b.n;i++) a.a[i]=a.a[i]-b.a[i]; return a;}
friend inline mtaylor<T> & operator -=(mtaylor<T> & a,const T & b) { if(a.n) a.a[0]=a.a[0]-b; return a;}
friend inline mtaylor<T> operator -(mtaylor<T> a,const mtaylor<T> & b) { return a-=b; }
friend inline mtaylor<T> operator -(mtaylor<T> a,const T & b) { return a-=b; }
friend inline mtaylor<T> operator -(const T & b,mtaylor<T> a) { return -(a-=b); }
friend inline mtaylor<T> & operator *=(mtaylor<T> & a,const mtaylor<T> & b)
{
int i,j;
a.minsize(b.n);
for(i=a.n-1;i>=0;i--)
{
a.a[i]=a.a[i]*b.a[0];
for(j=1;j<=i && j<b.n;j++)
a.a[i]=a.a[i]+a.a[i-j]*b.a[j];
}
return a;
}
friend inline mtaylor<T> & operator *=(mtaylor<T> & a,const T & b) { int i;for(i=0;i<a.n;i++) a.a[i]=a.a[i]*b; return a;}
friend inline mtaylor<T> operator *(mtaylor<T> a,const mtaylor<T> & b) { return a*=b; }
friend inline mtaylor<T> operator *(mtaylor<T> a,const T & b) { return a*=b; }
friend inline mtaylor<T> operator *(const T & b,mtaylor<T> a) { return a*=b; }
friend inline mtaylor<T> & operator /=(mtaylor<T> & a,const mtaylor<T> & b)
{
int i,j;
a.minsize(b.n);
for(i=0;i<a.n;i++)
{
for(j=1;j<=i && j<b.n;j++)
a.a[i]=a.a[i]-a.a[i-j]*b.a[j];
a.a[i]=a.a[i]/b.a[0];
}
return a;
}
friend inline mtaylor<T> & operator /=(mtaylor<T> & a,const T & b) { int i; for(i=0;i<a.n;i++) a.a[i]=a.a[i]/b; return a;}
friend inline mtaylor<T> operator /(mtaylor<T> a,const mtaylor<T> & b) { return a/=b; }
friend inline mtaylor<T> operator /(mtaylor<T> a,const T & b) { return a/=b; }
friend inline mtaylor<T> operator /(const T & b,mtaylor<T> a) { return mtaylor<T>(b,a.n)/a; }
friend inline ostream & operator <<(ostream & o,const mtaylor<T> & a) { int i; for(i=0;i<a.n-1;i++) o << a.a[i] << ' '; return (o << a.a[a.n-1]); }
inline T operator ()(int i) const { return a[i]; }
inline T & operator ()(int i) { return a[i]; }
inline int dim(void) const { return n; }
friend inline int operator==(const mtaylor<T> & a,const mtaylor<T> & b)
{
return (a.n!=0) && (b.n!=0) && (a.a[0]==b.a[0]);
}
friend inline int operator!=(const mtaylor<T> & a,const mtaylor<T> & b)
{
return (a.n==0) || (b.n==0) || (a.a[0]!=b.a[0]);
}
friend inline mtaylor<T> sqrt(const mtaylor<T> & a)
{
if(a.n<=0)
return a; // nonsense
mtaylor<T> result(sqrt(a.a[0]),a.n);
int k,j;
T factor=T(1)/(T(2)*result.a[0]);
for(k=1;k<a.n;k++)
{
result.a[k]=a.a[k];
if(k&1) // odd k
{
for(j=1;j<=(k-1)/2;j++)
result.a[k]=result.a[k]-T(2)*result.a[j]*result.a[k-j];
} else // even k
{
result.a[k]=result.a[k]-result.a[k/2]*result.a[k/2];
for(j=1;j<k/2;j++)
result.a[k]=result.a[k]-T(2)*result.a[j]*result.a[k-j];
}
result.a[k]=result.a[k]*factor;
}
return result;
}
friend inline mtaylor<T> sqr(const mtaylor<T> & a)
{
if(a.n<=0)
return a; // nonsense
mtaylor<T> result(sqr(a.a[0]),a.n);
int k,j;
for(k=1;k<a.n;k++)
{
if(k&1) // odd k
{
result.a[k]=a.a[0]*a.a[k];
for(j=1;j<=(k-1)/2;j++)
result.a[k]=result.a[k]+a.a[j]*a.a[k-j];
result.a[k]=T(2)*result.a[k];
} else // even k
{
result.a[k]=sqr(a.a[k/2]);
for(j=0;j<k/2;j++)
result.a[k]=result.a[k]+T(2)*a.a[j]*a.a[k-j];
}
}
return result;
}
friend inline mtaylor<T> exp(const mtaylor<T> & a)
{
if(a.n<=0)
return a; // nonsense
mtaylor<T> result(exp(a.a[0]),a.n);
int k,j;
for(k=1;k<a.n;k++)
{
result.a[k]=result.a[0]*T(k)*a.a[k];
for(j=1;j<k;j++)
result.a[k]=result.a[k]+result.a[j]*T(k-j)*a.a[k-j];
result.a[k]=result.a[k]/T(k);
}
return result;
}
friend inline mtaylor<T> sin(const mtaylor<T> & a)
{
if(a.n<=0)
return a; // nonsense
mtaylor<T> sresult(sin(a.a[0]),a.n);
mtaylor<T> cresult(cos(a.a[0]),a.n);
int k,j;
for(k=1;k<a.n;k++)
{
sresult.a[k]=cresult.a[0]*a.a[k]*T(k);
for(j=1;j<k;j++)
sresult.a[k]=sresult.a[k]+cresult.a[j]*T(k-j)*a.a[k-j];
sresult.a[k]=sresult.a[k]/T(k);
if(k<a.n-1)
{
cresult.a[k]=sresult.a[0]*a.a[k]*T(k);
for(j=1;j<k;j++)
cresult.a[k]=cresult.a[k]+sresult.a[j]*T(k-j)*a.a[k-j];
cresult.a[k]=cresult.a[k]/T(-k);
}
}
return sresult;
}
friend inline mtaylor<T> cos(const mtaylor<T> & a)
{
if(a.n<=0)
return a; // nonsense
mtaylor<T> sresult(sin(a.a[0]),a.n);
mtaylor<T> cresult(cos(a.a[0]),a.n);
int k,j;
for(k=1;k<a.n;k++)
{
if(k<a.n-1)
{
sresult.a[k]=cresult.a[0]*a.a[k]*T(k);
for(j=1;j<k;j++)
sresult.a[k]=sresult.a[k]+cresult.a[j]*T(k-j)*a.a[k-j];
sresult.a[k]=sresult.a[k]/T(k);
}
cresult.a[k]=sresult.a[0]*a.a[k]*T(k);
for(j=1;j<k;j++)
cresult.a[k]=cresult.a[k]+sresult.a[j]*T(k-j)*a.a[k-j];
cresult.a[k]=cresult.a[k]/T(-k);
}
return cresult;
}
friend inline mtaylor<T> tan(const mtaylor<T> & a)
{
if(a.n<=0)
return a; // nonsense
mtaylor<T> sresult(sin(a.a[0]),a.n);
mtaylor<T> cresult(cos(a.a[0]),a.n);
int k,j;
for(k=1;k<a.n;k++)
{
sresult.a[k]=cresult.a[0]*a.a[k]*T(k);
for(j=1;j<k;j++)
sresult.a[k]=sresult.a[k]+cresult.a[j]*T(k-j)*a.a[k-j];
sresult.a[k]=sresult.a[k]/T(k);
cresult.a[k]=sresult.a[0]*a.a[k]*T(k);
for(j=1;j<k;j++)
cresult.a[k]=cresult.a[k]+sresult.a[j]*T(k-j)*a.a[k-j];
cresult.a[k]=cresult.a[k]/T(-k);
}
return sresult/cresult;
}
friend inline mtaylor<T> cot(const mtaylor<T> & a)
{
if(a.n<=0)
return a; // nonsense
mtaylor<T> sresult(sin(a.a[0]),a.n);
mtaylor<T> cresult(cos(a.a[0]),a.n);
int k,j;
for(k=1;k<a.n;k++)
{
sresult.a[k]=cresult.a[0]*a.a[k]*T(k);
for(j=1;j<k;j++)
sresult.a[k]=sresult.a[k]+cresult.a[j]*T(k-j)*a.a[k-j];
sresult.a[k]=sresult.a[k]/T(k);
cresult.a[k]=sresult.a[0]*a.a[k]*T(k);
for(j=1;j<k;j++)
cresult.a[k]=cresult.a[k]+sresult.a[j]*T(k-j)*a.a[k-j];
cresult.a[k]=cresult.a[k]/T(-k);
}
return cresult/sresult;
}
friend inline mtaylor<T> sinh(const mtaylor<T> & a)
{
if(a.n<=0)
return a; // nonsense
mtaylor<T> sresult(sinh(a.a[0]),a.n);
mtaylor<T> cresult(cosh(a.a[0]),a.n);
int k,j;
for(k=1;k<a.n;k++)
{
sresult.a[k]=cresult.a[0]*a.a[k]*T(k);
for(j=1;j<k;j++)
sresult.a[k]=sresult.a[k]+cresult.a[j]*T(k-j)*a.a[k-j];
sresult.a[k]=sresult.a[k]/T(k);
if(k<a.n-1)
{
cresult.a[k]=sresult.a[0]*a.a[k]*T(k);
for(j=1;j<k;j++)
cresult.a[k]=cresult.a[k]+sresult.a[j]*T(k-j)*a.a[k-j];
cresult.a[k]=cresult.a[k]/T(k);
}
}
return sresult;
}
friend inline mtaylor<T> cosh(const mtaylor<T> & a)
{
if(a.n<=0)
return a; // nonsense
mtaylor<T> sresult(sinh(a.a[0]),a.n);
mtaylor<T> cresult(cosh(a.a[0]),a.n);
int k,j;
for(k=1;k<a.n;k++)
{
if(k<a.n-1)
{
sresult.a[k]=cresult.a[0]*a.a[k]*T(k);
for(j=1;j<k;j++)
sresult.a[k]=sresult.a[k]+cresult.a[j]*T(k-j)*a.a[k-j];
sresult.a[k]=sresult.a[k]/T(k);
}
cresult.a[k]=sresult.a[0]*a.a[k]*T(k);
for(j=1;j<k;j++)
cresult.a[k]=cresult.a[k]+sresult.a[j]*T(k-j)*a.a[k-j];
cresult.a[k]=cresult.a[k]/T(k);
}
return cresult;
}
friend inline mtaylor<T> tanh(const mtaylor<T> & a)
{
if(a.n<=0)
return a; // nonsense
mtaylor<T> sresult(sinh(a.a[0]),a.n);
mtaylor<T> cresult(cosh(a.a[0]),a.n);
int k,j;
for(k=1;k<a.n;k++)
{
sresult.a[k]=cresult.a[0]*a.a[k]*T(k);
for(j=1;j<k;j++)
sresult.a[k]=sresult.a[k]+cresult.a[j]*T(k-j)*a.a[k-j];
sresult.a[k]=sresult.a[k]/T(k);
cresult.a[k]=sresult.a[0]*a.a[k]*T(k);
for(j=1;j<k;j++)
cresult.a[k]=cresult.a[k]+sresult.a[j]*T(k-j)*a.a[k-j];
cresult.a[k]=cresult.a[k]/T(k);
}
return sresult/cresult;
}
friend inline mtaylor<T> coth(const mtaylor<T> & a)
{
if(a.n<=0)
return a; // nonsense
mtaylor<T> sresult(sinh(a.a[0]),a.n);
mtaylor<T> cresult(cosh(a.a[0]),a.n);
int k,j;
for(k=1;k<a.n;k++)
{
sresult.a[k]=cresult.a[0]*a.a[k]*T(k);
for(j=1;j<k;j++)
sresult.a[k]=sresult.a[k]+cresult.a[j]*T(k-j)*a.a[k-j];
sresult.a[k]=sresult.a[k]/T(k);
cresult.a[k]=sresult.a[0]*a.a[k]*T(k);
for(j=1;j<k;j++)
cresult.a[k]=cresult.a[k]+sresult.a[j]*T(k-j)*a.a[k-j];
cresult.a[k]=cresult.a[k]/T(k);
}
return cresult/sresult;
}
friend inline mtaylor<T> log(const mtaylor<T> & a)
{
if(a.n<=0)
return a; // nonsense
mtaylor<T> result(log(a.a[0]),a.n);
int k,j;
for(k=1;k<a.n;k++)
{
result.a[k]=a.a[k];
for(j=1;j<k;j++)
{
result.a[k]=result.a[k]-T(j)*result.a[j]*a.a[k-j]/T(k);
}
result.a[k]=result.a[k]/a.a[0];
}
return result;
}
friend inline mtaylor<T> asin(const mtaylor<T> & a)
{
if(a.n<=0)
return a; // nonsense
mtaylor<T> result(asin(a.a[0]),a.n);
mtaylor<T> invdiff(sqrt(T(1)-a*a));
int k,j;
for(k=1;k<a.n;k++)
{
result.a[k]=a.a[k];
for(j=1;j<k;j++)
{
result.a[k]=result.a[k]-T(j)*result.a[j]*invdiff.a[k-j]/T(k);
}
result.a[k]=result.a[k]/invdiff.a[0];
}
return result;
}
friend inline mtaylor<T> acos(const mtaylor<T> & a)
{
if(a.n<=0)
return a; // nonsense
mtaylor<T> result(acos(a.a[0]),a.n);
mtaylor<T> invdiff(-sqrt(T(1)-a*a));
int k,j;
for(k=1;k<a.n;k++)
{
result.a[k]=a.a[k];
for(j=1;j<k;j++)
{
result.a[k]=result.a[k]-T(j)*result.a[j]*invdiff.a[k-j]/T(k);
}
result.a[k]=result.a[k]/invdiff.a[0];
}
return result;
}
friend inline mtaylor<T> atan(const mtaylor<T> & a)
{
if(a.n<=0)
return a; // nonsense
mtaylor<T> result(atan(a.a[0]),a.n);
mtaylor<T> invdiff(T(1)+a*a);
int k,j;
for(k=1;k<a.n;k++)
{
result.a[k]=a.a[k];
for(j=1;j<k;j++)
{
result.a[k]=result.a[k]-T(j)*result.a[j]*invdiff.a[k-j]/T(k);
}
result.a[k]=result.a[k]/invdiff.a[0];
}
return result;
}
friend inline mtaylor<T> acot(const mtaylor<T> & a)
{
if(a.n<=0)
return a; // nonsense
mtaylor<T> result(acot(a.a[0]),a.n);
mtaylor<T> invdiff(T(-1)-a*a);
int k,j;
for(k=1;k<a.n;k++)
{
result.a[k]=a.a[k];
for(j=1;j<k;j++)
{
result.a[k]=result.a[k]-T(j)*result.a[j]*invdiff.a[k-j]/T(k);
}
result.a[k]=result.a[k]/invdiff.a[0];
}
return result;
}
friend inline mtaylor<T> asinh(const mtaylor<T> & a)
{
return log(a+sqrt(a*a+T(1)));
}
friend inline mtaylor<T> acosh(const mtaylor<T> & a)
{
return log(a+sqrt(a*a-T(1)));
}
friend inline mtaylor<T> atanh(const mtaylor<T> & a)
{
return log((T(1)+a)/(T(1)-a))/T(2);
}
friend inline mtaylor<T> acoth(const mtaylor<T> & a)
{
return log((a+T(1))/(a-T(1)))/T(2);
}
friend inline mtaylor<T> pow(const mtaylor<T> & a,const T &b)
{
// if(a.n>1)
// a.a[1]=b*a.a[1]*pow(a.a[0],b-1);
// a.a[0]=pow(a.a[0],b);
// return a;
return exp(b*log(a));
}
friend inline mtaylor<T> power(const mtaylor<T> & a,int b)
{
// if(a.n>1)
// a.a[1]=b*a.a[1]*pow(a.a[0],b-1);
// a.a[0]=pow(a.a[0],b);
// return a;
if(b<0)
return power(1/a,-b);
if(b==1)
return a;
if(b==0)
return mtaylor<T>(1);
if(b&1)
return a*power(a,b-1);
else
return sqr(power(a,b/2));
}
};
#endif
syntax highlighted by Code2HTML, v. 0.9.1