- Language
- C++
- Purpose
- Calculates Laurent Series: <math>\underset{n=0}{\overset{N-1}{\sum}}a_nZ^{-n}</math>
- Author
- James R. Craig, Igor Jankovic, Randal Barnes
#include <math.h>
#include <complex>
typedef complex<double> cmplex;
/***********************************************************************
Laurent
************************************************************************
Calculates Laurent Expansion about Z=0.0 using Horner's rule (Knuth's rule
for the real version)
Laurent() = sum{0..N-1} a_n Z^{-n}
LaurentRe() = Laurent with real coefficients only
Author: James R. Craig
-----------------------------------------------------------------------*/
cmplex Laurent(const cmplex &Z, const int N, const cmplex * const &a )
{
static cmplex w,sum;
w=1.0/Z;
sum=a[N-1];
for (int i=N-2; i>=0; i--){ sum=(sum*w)+a[i]; }
return sum;
}
//-----------------------------------------------------------------------
cmplex Laurent(const cmplex &Z, const int N, const double * const &a)
{
//Algorithm provided by Randal J. Barnes, Univ. of Minnesota
static double r,s,c,d,e;
static double wr,wi,den;
den=1.0/(Z.real()*Z.real()+Z.imag()*Z.imag());
wr= Z.real()*den;
wi=-Z.imag()*den;
r=2.0*wr;
s=wr*wr+wi*wi;
c=a[N-1];
d=a[N-2];
for (int i=N-3; i>=0; i--){
e=c;
c=d+r*e;
d=a[i]-s*e;
}
return -(c*cmplex(wr,wi)+d);
}