jiri/LA_library
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

*** empty log message ***

Browse Source
This commit is contained in:
jiri 2009-11-05 09:55:44 +00:00
parent ea6547f9f6
commit e090cc5712
1 changed files with 6 additions and 4 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

10
matexp.h
View File

@ -176,7 +176,7 @@ return int(ceil(log(n)/log2-log(.75)));
template<class T, class C> template<class T, class C>
NRVec<C> exp_aux(const T &x, int &power,int maxpower= -1, int maxtaylor= -1, C prescale=1.) NRVec<C> exp_aux(const T &x, int &power,int maxpower= -1, int maxtaylor= -1, typename LA_traits<T>::elementtype prescale=1.)
{ {
//should better be computed by mathematica to have accurate last digits, chebyshev instead, see exp in glibc //should better be computed by mathematica to have accurate last digits, chebyshev instead, see exp in glibc
static double exptaylor[]={ static double exptaylor[]={
@ -209,7 +209,7 @@ double scale=exp(-log(2.)*power);
//find how long taylor expansion will be necessary //find how long taylor expansion will be necessary
const double precision=1e-14; //decreasing brings nothing const double precision=1e-14; //further decreasing brings nothing
double s,t; double s,t;
s=mnorm*scale; s=mnorm*scale;
int n=0; int n=0;
@ -230,6 +230,8 @@ for(i=0,t=1.;i<=n;i++)
taylor2[i]=exptaylor[i]*t; taylor2[i]=exptaylor[i]*t;
t*=scale; t*=scale;
} }
//cout <<"TEST power, scale "<<power<<" "<<scale<<endl;
//cout <<"TEST taylor2 "<<taylor2<<endl;
return taylor2; return taylor2;
} }
@ -242,7 +244,7 @@ const T exp(const T &x, bool horner=true, int maxpower= -1, int maxtaylor= -1 )
int power; int power;
//prepare the polynom of and effectively scale T //prepare the polynom of and effectively scale T
NRVec<typename LA_traits<T>::elementtype> taylor2=exp_aux<T,typename LA_traits<T>::elementtype>(x,power,maxpower,maxtaylor); NRVec<typename LA_traits<T>::normtype> taylor2=exp_aux<T,typename LA_traits<T>::normtype>(x,power,maxpower,maxtaylor);
T r= horner?polynom0(x,taylor2):polynom(x,taylor2); T r= horner?polynom0(x,taylor2):polynom(x,taylor2);
@ -267,7 +269,7 @@ if(mat.nrows()!=mat.ncols()||(unsigned int) mat.nrows() != (unsigned int)rhs.siz
int power; int power;
//prepare the polynom of and effectively scale the matrix //prepare the polynom of and effectively scale the matrix
NRVec<typename LA_traits<V>::elementtype> taylor2=exp_aux<M,typename LA_traits<V>::elementtype>(mat,power,maxpower,maxtaylor,scale); NRVec<typename LA_traits<V>::normtype> taylor2=exp_aux<M,typename LA_traits<V>::normtype>(mat,power,maxpower,maxtaylor,scale);
V tmp; V tmp;