?? mremez1.c
字號(hào):
extern float ad[66],x[66],y[66],grid[1045],des[1045],wt[1045],alpha[66];
extern int iext[66];
extern float pi2,dev;
extern int nfcns,ngrid,niter;
/*-------------------------------------------------------------------*/
float d(int k,int n,int m)
{float z,q;
int l,j;
z=1.;
q=x[k];
for(l=1;l<=m;l++)
{for(j=l;j<=n;j+=m)
if((j-k)!=0)
z*=2.*(q-x[j]);
}
z=1./z;
return(z);
}
/*------------------------------------------------------------------*/
float gee(int k,int n)
{float p,xf,d1,c;
int j;
p=0.0;
xf=grid[k];
xf=cos(pi2*xf);
d1=0.;
for(j=1;j<=n;j++)
{c=xf-x[j];
c=ad[j]/c;
d1+=c;
p+=c*y[j];
}
p/=d1;
return(p);
}
/*---------------------------------------------------------------------
Subroutine mremez1.c
in chapter 8
---------------------------------------------------------------------*/
void mremez1()
{
float a[66],p[66],q[66];
float devl,dnum,dden,dtemp,fsh,gtemp,aa,bb,ft,xt,xe;
float ynz,comp,y1,err,delf,xt1;
int itrmax,nz,nzz,jxt,jet,k,l,nu,jchnge;
int nut,klow,k1,knz,kup,luck,nut1;
int j,nzzmj,nzmj,kn,nm1,kkk,cn,jm1,nf1j;
int flag;
comp=0.0;
itrmax=25;
devl=-1.;
nz=nfcns+1;
nzz=nfcns+2;
niter=0;
while(1)
{
iext[nzz]=ngrid+1;
niter++;
if(niter>itrmax)break;
for(j=1;j<=nz;j++)
{jxt=iext[j];
dtemp=grid[jxt];
dtemp=cos(dtemp*pi2);
x[j]=dtemp;
}
jet=(nfcns-1)/15+1;
for (j=1;j<=nz;j++)
ad[j]=d(j,nz,jet);
dnum=0.0;
dden=0.0;
k=1;
for(j=1;j<=nz;j++)
{
l=iext[j];
dtemp=ad[j]*des[l];
dnum+=dtemp;
dtemp=k*ad[j]/wt[l];
dden+=dtemp;
k=-k;
}
dev=dnum/dden;
nu=1;
if(dev>0.)
nu=-1;
dev*=-nu;
k=nu;
for(j=1;j<=nz;j++)
{
l=iext[j];
dtemp=k*dev/wt[l];
y[j]=des[l]+dtemp;
k=-k;
}
if(dev<=devl)break;
devl=dev;
jchnge=0;
k1=iext[1];
knz=iext[nz];
klow=0;
nut=-nu;
j=1;
while(1)
{
while(1)
{
if(j==nzz)
ynz=comp;
if(j>=nzz)break;
kup=iext[j+1];
l=iext[j]+1;
nut=-nut;
if(j==2)
y1=comp;
comp=dev;
if(l>=kup)break;
err=gee(l,nz);
err=(err-des[l])*wt[l];
dtemp=nut*err-comp;
if(dtemp<=0.)break;
comp=nut*err;
while(1)
{
l++;
if(l>=kup)break;
err=gee(l,nz);
err=(err-des[l])*wt[l];
dtemp=nut*err-comp;
if(dtemp<=0.)break;
comp=nut*err;
}
iext[j]=l-1;
j++;
klow=l-1;
jchnge++;
}
if(j<nzz)
{
l--;
flag=0;
while(1)
{
l--;
if(l<=klow)break;
err=gee(l,nz);
err=(err-des[l])*wt[l];
dtemp=nut*err-comp;
if(dtemp>0.)break;
if(jchnge>0.)
{
klow=iext[j];
j++;
flag=1;
break;
}
}
if(flag==1)continue;
flag=0;
if(l>klow)
{
comp=nut*err;
while(1)
{
l--;
if(l<=klow)break;
err=gee(l,nz);
err=(err-des[l])*wt[l];
dtemp=nut*err-comp;
if(dtemp<=0.)break;
comp=nut*err;
}
klow=iext[j];
iext[j]=l+1;
j++;
jchnge++;
flag=1;
}
if(flag==1)continue;
l=iext[j]+1;
flag=0;
if(jchnge>0.)
{
iext[j]=l-1;
j++;
klow=l-1;
jchnge++;
flag=1;
}
if(flag!=1)
{
while(1)
{
l++;
if(l>=kup)break;
err=gee(l,nz);
err=(err-des[l])*wt[l];
dtemp=nut*err-comp;
if(dtemp>0.)
{
comp=nut*err;
while(1)
{
l++;
if(l>=kup)break;
err=gee(l,nz);
err=(err-des[l])*wt[l];
dtemp=nut*err-comp;
if(dtemp<=0.)break;
comp=nut*err;
}
iext[j]=l-1;
jchnge++;
break;
}
}
klow=iext[j];
j++;
}
continue;
}
if(j>nzz)
{
flag=0;
if(luck>9)
{
kn=iext[nzz];
for(j=1;j<=nfcns;j++)
iext[j]=iext[j+1];
iext[nz]=kn;
flag=2;
break;
}
if(comp>y1)
y1=comp;
k1=iext[nzz];
}
else
{
flag=0;
if(k1>iext[1])
k1=iext[1];
if(knz<iext[nz])
knz=iext[nz];
nut1=nut;
nut=-nu;
l=0;
kup=k1;
comp=ynz*1.00001;
luck=1;
while(1)
{
l++;
if(l>=kup)break;
err=gee(l,nz);
err=(err-des[l])*wt[l];
dtemp=nut*err-comp;
flag=0;
if(dtemp>0.)
{
comp=nut*err;
j=nzz;
while(1)
{
l++;
if(l>=kup)break;
err=gee(l,nz);
err=(err-des[l])*wt[l];
dtemp=nut*err-comp;
if(dtemp<=0.)break;
comp=nut*err;
}
iext[j]=l-1;
j++;
klow=l-1;
jchnge++;
flag=1;
break;
}
}
if(flag==1)continue;
luck=6;
}
l=ngrid+1;
klow=knz;
nut=-nut1;
comp=y1*1.00001;
while(1)
{
l--;
if(l<=klow)break;
err=gee(l,nz);
err=(err-des[l])*wt[l];
dtemp=nut*err-comp;
flag=0;
if(dtemp>0.)
{
j=nzz;
comp=nut*err;
luck+=10;
while(1)
{
l--;
if(l<=klow)break;
err=gee(l,nz);
err=(err-des[l])*wt[l];
dtemp=nut*err-comp;
if(dtemp<=0.)break;
comp=nut*err;
}
klow=iext[j];
iext[j]=l+1;
j++;
jchnge++;
flag=1;
break;
}
}
if(flag!=1)break;
}
if(flag==2)continue;
if(luck!=6)
{
for(j=1;j<=nfcns;j++)
{
nzzmj=nzz-j;
nzmj=nz-j;
iext[nzzmj]=iext[nzmj];
}
iext[1]=k1;
continue;
}
if(jchnge>0)continue;
}
nm1=nfcns-1;
fsh=1.e-06;
gtemp=grid[1];
x[nzz]=-2.;
cn=2*nfcns-1;
delf=1./cn;
l=1;
kkk=0;
if(grid[1]<=.01&&grid[ngrid]>0.49)
kkk=1;
if(nfcns<=3)
kkk=1;
if(kkk!=1)
{
dtemp=cos(pi2*grid[1]);
dnum=cos(pi2*grid[ngrid]);
aa=2./(dtemp-dnum);
bb=-(dtemp+dnum)/(dtemp-dnum);
}
for(j=1;j<=nfcns;j++)
{
ft=j-1;
ft*=delf;
xt=cos(pi2*ft);
if(kkk!=1)
{
xt=(xt-bb)/aa;
xt1=sqrt(1.-xt*xt);
ft=atan2(xt1,xt)/pi2;
}
while(1)
{
xe=x[l];
if(xt>xe)break;
if((xe-xt)<fsh)break;
l++;
}
if((xe-xt)<fsh||(xt-xe)<fsh)
a[j]=y[l];
if(xt>xe)
{
grid[1]=ft;
a[j]=gee(1,nz);
}
if(l>1)l--;
}
grid[1]=gtemp;
dden=pi2/cn;
for(j=1;j<=nfcns;j++)
{
dtemp=0.;
dnum=j-1;
dnum*=dden;
if(nm1>=1)
for(k=1;k<=nm1;k++)
dtemp+=a[k+1]*cos(dnum*k);
dtemp=dtemp*2.+a[1];
alpha[j]=dtemp;
}
for(j=2;j<=nfcns;j++)
alpha[j]*=2./cn;
alpha[1]/=cn;
if(kkk!=1)
{
p[1]=2.*alpha[nfcns]*bb+alpha[nm1];
p[2]=2.*alpha[nfcns]*aa;
q[1]=alpha[nfcns-2]-alpha[nfcns];
for(j=2;j<=nm1;j++)
{
if(j>=nm1)
{
aa*=.5;
bb*=.5;
}
p[j+1]=0.;
for(k=1;k<=j;k++)
{
a[k]=p[k];
p[k]=2.*bb*a[k];
}
p[2]+=a[1]*2.*aa;
jm1=j-1;
for(k=1;k<=jm1;k++)
p[k]+=q[k]+aa*a[k+1];
for(k=3;k<=j+1;k++)
p[k]+=aa*a[k-1];
if(j==nm1) continue;
for(k=1;k<=j;k++)
q[k]=-a[k];
nf1j=nfcns-j-1;
q[1]+=alpha[nf1j];
}
for(j=1;j<=nfcns;j++)
alpha[j]=p[j];
}
if(nfcns>3)
return;
alpha[nfcns+1]=0.0;
alpha[nfcns+2]=0.0;
return;
}
?? 快捷鍵說(shuō)明
復(fù)制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號(hào)
Ctrl + =
減小字號(hào)
Ctrl + -