Laurent series algorithm

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);
}