Laurent series algorithm
From aemwiki
- 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); }