50 lines
989 B
C++
50 lines
989 B
C++
#include <cmath>
|
|
#include <complex>
|
|
#include <limits>
|
|
#include "nr.h"
|
|
using namespace std;
|
|
|
|
void NR::laguer(Vec_I_CPLX_DP &a, complex<DP> &x, int &its)
|
|
{
|
|
const int MR=8,MT=10,MAXIT=MT*MR;
|
|
const DP EPS=numeric_limits<DP>::epsilon();
|
|
static const DP frac[MR+1]=
|
|
{0.0,0.5,0.25,0.75,0.13,0.38,0.62,0.88,1.0};
|
|
int iter,j;
|
|
DP abx,abp,abm,err;
|
|
complex<DP> dx,x1,b,d,f,g,h,sq,gp,gm,g2;
|
|
|
|
int m=a.size()-1;
|
|
for (iter=1;iter<=MAXIT;iter++) {
|
|
its=iter;
|
|
b=a[m];
|
|
err=abs(b);
|
|
d=f=0.0;
|
|
abx=abs(x);
|
|
for (j=m-1;j>=0;j--) {
|
|
f=x*f+d;
|
|
d=x*d+b;
|
|
b=x*b+a[j];
|
|
err=abs(b)+abx*err;
|
|
}
|
|
err *= EPS;
|
|
if (abs(b) <= err) return;
|
|
g=d/b;
|
|
g2=g*g;
|
|
h=g2-2.0*f/b;
|
|
sq=sqrt(DP(m-1)*(DP(m)*h-g2));
|
|
gp=g+sq;
|
|
gm=g-sq;
|
|
abp=abs(gp);
|
|
abm=abs(gm);
|
|
if (abp < abm) gp=gm;
|
|
dx=MAX(abp,abm) > 0.0 ? DP(m)/gp : polar(1+abx,DP(iter));
|
|
x1=x-dx;
|
|
if (x == x1) return;
|
|
if (iter % MT != 0) x=x1;
|
|
else x -= frac[iter/MT]*dx;
|
|
}
|
|
nrerror("too many iterations in laguer");
|
|
return;
|
|
}
|