55 lines
976 B
C
55 lines
976 B
C
|
|
#include <math.h>
|
|
#include "nrutil.h"
|
|
#define BIG 1.0e30
|
|
|
|
void pade(cof,n,resid)
|
|
double cof[];
|
|
float *resid;
|
|
int n;
|
|
{
|
|
void lubksb(),ludcmp(),mprove();
|
|
int j,k,*indx;
|
|
float d,rr,rrold,sum,**q,**qlu,*x,*y,*z;
|
|
|
|
indx=ivector(1,n);
|
|
q=matrix(1,n,1,n);
|
|
qlu=matrix(1,n,1,n);
|
|
x=vector(1,n);
|
|
y=vector(1,n);
|
|
z=vector(1,n);
|
|
for (j=1;j<=n;j++) {
|
|
y[j]=x[j]=cof[n+j];
|
|
for (k=1;k<=n;k++) {
|
|
q[j][k]=cof[j-k+n];
|
|
qlu[j][k]=q[j][k];
|
|
}
|
|
}
|
|
ludcmp(qlu,n,indx,&d);
|
|
lubksb(qlu,n,indx,x);
|
|
rr=BIG;
|
|
do {
|
|
rrold=rr;
|
|
for (j=1;j<=n;j++) z[j]=x[j];
|
|
mprove(q,qlu,n,indx,y,x);
|
|
for (rr=0.0,j=1;j<=n;j++)
|
|
rr += SQR(z[j]-x[j]);
|
|
} while (rr < rrold);
|
|
*resid=sqrt(rrold);
|
|
for (k=1;k<=n;k++) {
|
|
for (sum=cof[k],j=1;j<=k;j++) sum -= z[j]*cof[k-j];
|
|
y[k]=sum;
|
|
}
|
|
for (j=1;j<=n;j++) {
|
|
cof[j]=y[j];
|
|
cof[j+n] = -z[j];
|
|
}
|
|
free_vector(z,1,n);
|
|
free_vector(y,1,n);
|
|
free_vector(x,1,n);
|
|
free_matrix(qlu,1,n,1,n);
|
|
free_matrix(q,1,n,1,n);
|
|
free_ivector(indx,1,n);
|
|
}
|
|
#undef BIG
|