ch7.h

将C语言的常用程序集移植到VC开发环境的源代码

                for (i=0; i<=n-1; i++)
                  { dd=0.0;
                    for (j=0; j<=n-1; j++)
                      dd=dd+s02[j]*p[i*n+j];
                    s[i*10+8]=dd;
                  }
                nt=n;
                for (i=0; i<=n-1; i++)
                  { y[i*8]=y[i*8]+aa[0]*s[i*10+8];
                    y[i*8+1]=y[i*8+1]-s[i*10+8];
                    er[i]=er[i]+s[i*10+8];
                    if (fabs(s[i*10+8])<=(bnd*ym[i]))
                      nt=nt-1;
                  }
              }
          }
      }
    if (nt>0)
      { t0=td;
        if ((h>(hmin*1.00001))||(iw>=0))
          { if (iw!=0) rm=0.25*rm;
            iw=1; irt1=2;
            rr=fabs(hmin/hd);
            if (rm<rr) rm=rr;
            rr=fabs(hmax/hd);
            if (rm>rr) rm=rr;
            r=1.0;
            for (j=1; j<=m-1; j++)
              { r=r*rm;
                for (i=0; i<=n-1; i++)
                  y[i*8+j]=s[i*10+j]*r;
              }
            h=hd*rm;
            for (i=0; i<=n-1; i++)
              y[i*8]=s[i*10];
            idb=m;
            goto l80;
          }
        for (i=0; i<=n-1; i++)
          for (j=0; j<=m-1; j++)
            y[i*8+j]=s[i*10+j];
        h=hd; nq=nqd; jt=nq;
        free(d); free(p); free(s); free(s02);
        free(ym); free(er); free(yy);
        free(iis); free(jjs); free(y); return(-3);
      }
    dd=0.0;
    for (i=0; i<=n-1; i++)
      dd=dd+(er[i]/ym[i])*(er[i]/ym[i]);
    iw=0;
    if (dd<=e)
      { if (m>=3)
          for (j=2; j<=m-1; j++)
            for (i=0; i<=n-1; i++)
              y[i*8+j]=y[i*8+j]+aa[j]*er[i];
        kf=1; hw=h;
        if (idb>1)
          { idb=idb-1;
            if (idb<=1)
              for (i=0; i<=n-1; i++)
                s[i*10+9]=er[i];
            for (i=0; i<=n-1; i++)
	      if (ym[i]<fabs(y[i*8])) ym[i]=fabs(y[i*8]);
            jt=nq;
            goto l20;
          }
      }
    if (dd>e)
      { kf=kf-2;
        if (h<=(hmin*1.00001))
          { free(d); free(p); free(s); free(s02);
            free(ym); free(er); free(yy);
            free(iis); free(jjs); free(y);
            hw=h; jt=nq; return(-1);
          }
        t0=td;
        if (kf<=-5)
          { if (nq==1)
              { for (i=0; i<=n-1; i++)
                  for (j=0; j<=m-1; j++)
                    y[i*8+j]=s[i*10+j];
                h=hd; nq=nqd; jt=nq;
		free(d); free(p); free(s); free(s02);
		free(ym); free(er); free(yy);
		free(iis); free(jjs); free(y); return(-4);
              }
            for (i=0; i<=n-1; i++) yy[i]=y[i*8];
            ggearf(t0,yy,n,d);
            r=h/hd;
            for (i=0; i<=n-1; i++)
              { y[i*8]=s[i*10];
                s[i*10+1]=hd*d[i];
                y[i*8+1]=s[i*10+1]*r;
              }
            nq=1; kf=1; goto l60;
          }
      }
    pr2=log(dd/e); pr2=enq2*pr2; pr2=exp(pr2);
    pr2=1.2*pr2;
    pr3=1.0e+20;
    if (nq<7)
      if (kf>-1)
        { dd=0.0;
          for (i=0; i<=n-1; i++)
            { pr3=(er[i]-s[i*10+9])/ym[i];
              dd=dd+pr3*pr3;
            }
          pr3=log(dd/eup); pr3=enq3*pr3;
          pr3=exp(pr3); pr3=1.4*pr3;
        }
    pr1=1.0e+20;
    if (nq>1)
      { dd=0.0;
        for (i=0; i<=n-1; i++)
          { pr1=y[i*8+m-1]/ym[i];
            dd=dd+pr1*pr1;
          }
        pr1=log(dd/edwn); pr1=enq1*pr1;
        pr1=exp(pr1); pr1=1.3*pr1;
      }
    if (pr2<=pr3)
      { if (pr2>pr1)
          { r=1.0e+04;
            if (pr1>1.0e-04) r=1.0/pr1;
            nqw=nq-1;
          }
        else
          { nqw=nq; r=1.0e+04;
            if (pr2>1.0e-04) r=1.0/pr2;
          }
      }
    else
      { if (pr3<pr1)
          { r=1.0e+04;
            if (pr3>1.0e-04) r=1.0/pr3;
            nqw=nq+1;
          }
        else
          { r=1.0e+04;
            if (pr1>1.0e-04) r=1.0/pr1;
            nqw=nq-1;
          }
      }
    idb=10;
    if (kf==1)
      if (r<1.1)
        { for (i=0; i<=n-1; i++)
            if (ym[i]<fabs(y[i*8])) ym[i]=fabs(y[i*8]);
          jt=nq; goto l20;
        }
    if (nqw>nq)
      for (i=0; i<=n-1; i++)
        y[i*8+nqw]=er[i]*aa[m-1]/(m+0.0);
    m=nqw+1;
    if (kf==1)
      { irt=2; rr=hmax/fabs(h);
        if (r>rr) r=rr;
        h=h*r; hw=h;
        if (nq==nqw)
          { r1=1.0;
            for (j=1; j<=m-1; j++)
              { r1=r1*r;
                for (i=0; i<=n-1; i++)
                  y[i*8+j]=y[i*8+j]*r1;
              }
            idb=m;
            for (i=0; i<=n-1; i++)
              if (ym[i]<fabs(y[i*8])) ym[i]=fabs(y[i*8]);
            jt=nq; goto l20;
          }
        nq=nqw;
        goto l60;
      }
    rm=rm*r; irt1=3;
    rr=fabs(hmin/hd);
    if (rm<rr) rm=rr;
    rr=fabs(hmax/hd);
    if (rm>rr) rm=rr;
    r=1.0;
    for (j=1; j<=m-1; j++)
      { r=r*rm;
        for (i=0; i<=n-1; i++)
          y[i*8+j]=s[i*10+j]*r;
      }
    h=hd*rm;
    for (i=0; i<=n-1; i++)
      y[i*8]=s[i*10];
    idb=m;
    if (nqw==nq) goto l80;
    nq=nqw; goto l60;
}
/////////////////////////////////////////////////////////////
void gdfte(double a,double b,double ya,double yb,int n,double y[],
		void (*gdftef)(double x,double z[]))
{
	//extern void gdftef();
    //extern int atrde();//#include "ch1.h"
    int i,k,nn,m1;//,j  : unreferenced local variable
    double z[4],h,x,*g,*d;
    g=(double*)malloc(6*n*sizeof(double));
    d=(double*)malloc(2*n*sizeof(double));
    h=(b-a)/(n-1.0); nn=2*n-1;
    g[0]=1.0; g[1]=0.0;
    y[0]=ya; y[n-1]=yb;
    g[3*n-3]=1.0; g[3*n-4]=0.0;
    for (i=2; i<=n-1; i++)
      { x=a+(i-1)*h;
        gdftef(x,z);
        k=3*(i-1)-1;
        g[k]=z[0]-h*z[1]/2.0;
        g[k+1]=h*h*z[2]-2.0*z[0];
        g[k+2]=z[0]+h*z[1]/2.0;
        y[i-1]=h*h*z[3];
      }
    m1=3*n-2;
    atrde(g,n,m1,y);
    h=h/2.0;
    g[0]=1.0; g[1]=0.0;
    d[0]=ya; d[nn-1]=yb;
    g[3*nn-3]=1.0; g[3*nn-4]=0.0;
    for (i=2; i<=nn-1; i++)
      { x=a+(i-1)*h;
        gdftef(x,z);
        k=3*(i-1)-1;
        g[k]=z[0]-h*z[1]/2.0;
        g[k+1]=h*h*z[2]-2.0*z[0];
        g[k+2]=z[0]+h*z[1]/2.0;
        d[i-1]=h*h*z[3];
      }
    m1=3*nn-2;
    atrde(g,nn,m1,d);
    for (i=2; i<=n-1; i++)
      { k=i+i-1;
        y[i-1]=(4.0*d[k-1]-y[i-1])/3.0;
      }
    free(g); free(d);
    return;
}
#endif
/////////////////////////////////////////////////////////////
/*
void gelr1f(double t,double y[],int n,double d[])
{ 
	t=t; n=n;
    d[0]=y[1]; d[1]=-y[0]; d[2]=-y[2];
    return;
}

void gelr2f(double t,double y[],int n,double d[])
{
	t=t; n=n;
    d[0]=y[1]; d[1]=-y[0]; d[2]=-y[2];
    return;
}

void gwityf(double t,double y[],int n,double d[])
{ 
	t=t; n=n;
    d[0]=y[1]; d[1]=-y[0]; d[2]=-y[2];
    return;
}

void grkt1f(double t,double y[],int n,double d[])
{
	t=t; n=n;
    d[0]=y[1]; d[1]=-y[0]; d[2]=-y[2];
    return;
}

void grkt2f(double t,double y[],int n,double d[])
{ 
	t=t; n=n;
    d[0]=y[1]; d[1]=-y[0];
    return;
}

void ggil1f(double t,double y[],int n,double d[])
{ 
	t=t; n=n;
    d[0]=y[1]; d[1]=-y[0]; d[2]=-y[2];
    return;
}

void ggil2f(double t,double y[],int n,double d[])
{
	t=t; n=n;
    d[0]=y[1]; d[1]=-y[0]; d[2]=-y[2];
    return;
}

void gmrsnf(double t,double y[],int n,double d[])
{ 
	double q;
    n=n;
    q=60.0*(0.06+t*(t-0.6));
    d[0]=q*y[1]; d[1]=-q*y[0];
    return;
}

void gpbs1f(double t,double y[],int n,double d[])
{
	t=t; n=n;
    d[0]=-y[1]; d[1]=y[0];
    return;
}

void gpbs2f(double t,double y[],int n,double d[])
{ 
	t=t; n=n;
    d[0]=-y[1]; d[1]=y[0];
    return;
}

void ggjfqf(double t,double y[],int n,double d[])
{
	t=t; n=n;
    d[0]=-y[1]; d[1]=y[0];
    return;
}

void gadmsf(double t,double y[],int n,double d[])
{ 
	t=t; n=n;
    d[0]=y[1]; d[1]=-y[0]; d[2]=-y[2];
    return;
}

void ghamgf(double t,double y[],int n,double d[])
{
	t=t; n=n;
    d[0]=y[1]; d[1]=-y[0]; d[2]=y[2];
    return;
}

void gtnr1f(double t,double y[],int n,double d[])
{
	t=t; n=n;
    d[0]=-21.0*y[0]+19.0*y[1]-20.0*y[2];
    d[1]=19.0*y[0]-21.0*y[1]+20.0*y[2];
    d[2]=40.0*y[0]-40.0*y[1]-40.0*y[2];
    return;
}
void gtnr2f(double t,double y[],int n,double d[])
{
	t=t; n=n;
    d[0]=-21.0*y[0]+19.0*y[1]-20.0*y[2];
    d[1]=19.0*y[0]-21.0*y[1]+20.0*y[2];
    d[2]=40.0*y[0]-40.0*y[1]-40.0*y[2];
    return;
}
//--------------
void ggearf(double t,double y[],int n,double d[])
{
	t=t; n=n;
    d[0]=-21.0*y[0]+19.0*y[1]-20.0*y[2];
    d[1]=19.0*y[0]-21.0*y[1]+20.0*y[2];
    d[2]=40.0*y[0]-40.0*y[1]-40.0*y[2];
    return;
}

void ggears(double t,double y[],int n,double *p)// int n; double t,y[],p[3][3];
{ 
	t=t; n=n; y[0]=y[0];
	p[0]=-21.0; p[1]=19.0; p[2]=-20.0;
    p[3]=19.0; p[4]=-21.0; p[5]=20.0;
    p[6]=40.0; p[7]=-40.0; p[8]=-40.0;
	//p[0][0]=-21.0; p[0][1]=19.0; p[0][2]=-20.0;
    //p[1][0]=19.0; p[1][1]=-21.0; p[1][2]=20.0;
    //p[2][0]=40.0; p[2][1]=-40.0; p[2][2]=-40.0;
	
    return;
}
//--------------
void gdftef(double x,double z[])
{
	z[0]=-1.0; z[1]=0.0;
    z[2]=2.0/(x*x); z[3]=1.0/x;
    return;
}
/**/