xqinghua.c

来自「该程序是用vc开发的对动态数组进行管理的DLL」· C语言 代码 · 共 2,502 行 · 第 1/5 页

  #include "xqinghua.h"
  #include "math.h"
  #include "stdio.h"
  int aband(b,d,n,l,il,m)
  int n,l,il,m;
  double b[],d[];
  { int ls,k,i,j,is,u,v;
    double p,t;
    if (il!=(2*l+1))
      { printf("fail\n"); return(-2);}
    ls=l;
    for (k=0;k<=n-2;k++)
      { p=0.0;
        for (i=k;i<=ls;i++)
          { t=fabs(b[i*il]);
            if (t>p) {p=t; is=i;}
          }
        if (p+1.0==1.0)
          { printf("fail\n"); return(0);}
        for (j=0;j<=m-1;j++)
          { u=k*m+j; v=is*m+j;
            t=d[u]; d[u]=d[v]; d[v]=t;
          }
        for (j=0;j<=il-1;j++)
          { u=k*il+j; v=is*il+j;
            t=b[u]; b[u]=b[v]; b[v]=t;
          }
        for (j=0;j<=m-1;j++)
          { u=k*m+j; d[u]=d[u]/b[k*il];}
        for (j=1;j<=il-1;j++)
          { u=k*il+j; b[u]=b[u]/b[k*il];}
        for (i=k+1;i<=ls;i++)
          { t=b[i*il];
            for (j=0;j<=m-1;j++)
              { u=i*m+j; v=k*m+j;
                d[u]=d[u]-t*d[v];
              }
            for (j=1;j<=il-1;j++)
              { u=i*il+j; v=k*il+j;
                b[u-1]=b[u]-t*b[v];
              }
            u=i*il+il-1; b[u]=0.0;
          }
        if (ls!=(n-1)) ls=ls+1;
      }
    p=b[(n-1)*il];
    if (fabs(p)+1.0==1.0)
      { printf("fail\n"); return(0);}
    for (j=0;j<=m-1;j++)
      { u=(n-1)*m+j; d[u]=d[u]/p;}
    ls=1;
    for (i=n-2;i>=0;i--)
      { for (k=0;k<=m-1;k++)
          { u=i*m+k;
            for (j=1;j<=ls;j++)
              { v=i*il+j; is=(i+j)*m+k;
                d[u]=d[u]-b[v]*d[is];
              }
          }
        if (ls!=(il-1)) ls=ls+1;
      }
    return(2);
  }



/*********** Example Start **********/
/*

  #include "aband.c"
  #include "stdio.h"
  main()
  { int i;
    static double b[8][5]={ {3.0,-4.0,1.0,0.0,0.0},
    {-2.0,-5.0,6.0,1.0,0.0},{1.0,3.0,-1.0,2.0,-3.0},
    {2.0,5.0,-5.0,6.0,-1.0},{-3.0,1.0,-1.0,2.0,-5.0},
    {6.0,1.0,-3.0,2.0,-9.0},{-4.0,1.0,-1.0,2.0,0.0},
    {5.0,1.0,-7.0,0.0,0.0}};
    static double d[8][3]={ {13.0,29.0,-13.0},
    {-6.0,17.0,-21.0},{-31.0,-6.0,4.0},{64.0,3.0,16.0},
    {-20.0,1.0,-5.0},{-22.0,-41.0,56.0},{-29.0,10.0,-21.0},
    {7.0,-24.0,20.0}};
    if (aband(b,d,8,2,5,3)>0)
      for (i=0;i<=7;i++)
      printf("x(%d)=%13.7e,  %13.7e,  %13.7e\n",
            i,d[i][0],d[i][1],d[i][2]);
  }




*/
/*********** Example  End **********/

/*********** Result  Start **********/
/*
x(0)= 3.000000e+00,   5.000000e+00,  -4.070908e-16
x(1)=-1.000000e+00,  -3.000000e+00,   3.000000e+00
x(2)= 1.908873e-16,   2.000000e+00,  -1.000000e+00
x(3)=-5.000000e+00,  -2.220446e-15,  -1.332268e-15
x(4)= 7.000000e+00,  -9.688431e-16,   2.000000e+00
x(5)= 1.000000e+00,   1.000000e+00,  -3.000000e+00
x(6)= 2.000000e+00,  -1.000000e+00,  -2.385245e-17
x(7)=-9.288423e-16,   4.000000e+00,  -5.000000e+00
*/
/*********** Result  End **********/


  #include "stdlib.h"
  #include "math.h"
  #include "stdio.h"
  int abint(a,n,b,eps,x)
  int n;
  double a[],b[],x[],eps;
  { int i,j,k,kk;
    double *p,*r,*e,q,qq;
    extern int agaus();
    extern void brmul();
    p=malloc(n*n*sizeof(double));
    r=malloc(n*sizeof(double));
    e=malloc(n*sizeof(double));
    i=60;
    for (k=0; k<=n-1; k++)
    for (j=0; j<=n-1; j++)
      p[k*n+j]=a[k*n+j];
    for (k=0; k<=n-1; k++) x[k]=b[k];
    kk=agaus(p,x,n);
    if (kk==0)
      { free(p); free(r); free(e); return(kk); }
    q=1.0+eps;
    while (q>=eps)
      { if (i==0)
          { free(p); free(r); free(e); return(i); }
        i=i-1;
        brmul(a,x,n,n,1,e);
        for ( k=0; k<=n-1; k++) r[k]=b[k]-e[k];
        for ( k=0; k<=n-1; k++)
        for ( j=0; j<=n-1; j++)
           p[k*n+j]=a[k*n+j];
        kk=agaus(p,r,n);
        if (kk==0)
          { free(p); free(r); free(e); return(kk); }
        q=0.0;
        for ( k=0; k<=n-1; k++)
          { qq=fabs(r[k])/(1.0+fabs(x[k]+r[k]));
            if (qq>q) q=qq;
          }
        for ( k=0; k<=n-1; k++) x[k]=x[k]+r[k];
      }
    free(p); free(r); free(e); return(1);
  }



/*********** Example Start **********/
/*

  #include "stdio.h"
  #include "abint.c"
  #include "brmul.c"
  #include "agaus.c"
  main()
  { int k;
    double eps=0.000001,x[4];
    static double a[4][4]=
                 {{3.4336,-0.5238,0.67105,-0.15272},
                  {-0.5238,3.28326,-0.73051,-0.2689},
                  {0.67105,-0.73051,4.02612,0.01835},
                  {-0.15272,-0.2689,0.01835,2.75702}};
    static double b[4]={-1.0,1.5,2.5,-2.0};
    if (abint(a,4,b,eps,x)!=0)
      { printf("\n");
        for (k=0; k<=3; k++)
        printf("x[%d]=%13.7e\n",k,x[k]);
        printf("\n");
      }
  }




*/
/*********** Example  End **********/

/*********** Result  Start **********/
/*

x[0]=-3.977180e-01
x[1]= 5.100536e-01
x[2]= 7.829837e-01
x[3]=-7.029161e-01




*/
/*********** Result  End **********/


  #include "stdlib.h"
  #include "stdio.h"
  int acgas(ar,ai,n,br,bi)
  int n;
  double ar[],ai[],br[],bi[];
  { int *js,l,k,i,j,is,u,v;
    double p,q,s,d;
    js=malloc(n*sizeof(int));
    for (k=0;k<=n-2;k++)
      { d=0.0;
        for (i=k;i<=n-1;i++)
        for (j=k;j<=n-1;j++)
          { u=i*n+j;
            p=ar[u]*ar[u]+ai[u]*ai[u];
            if (p>d) {d=p;js[k]=j;is=i;}
          }
        if (d+1.0==1.0)
          { free(js); printf("err**fail\n");
            return(0);
          }
        if (is!=k)
          { for (j=k;j<=n-1;j++)
              { u=k*n+j; v=is*n+j;
                p=ar[u]; ar[u]=ar[v]; ar[v]=p;
                p=ai[u]; ai[u]=ai[v]; ai[v]=p;
              }
            p=br[k]; br[k]=br[is]; br[is]=p;
            p=bi[k]; bi[k]=bi[is]; bi[is]=p;
          }
        if (js[k]!=k)
          for (i=0;i<=n-1;i++)
            { u=i*n+k; v=i*n+js[k];
              p=ar[u]; ar[u]=ar[v]; ar[v]=p;
              p=ai[u]; ai[u]=ai[v]; ai[v]=p;
            }
        v=k*n+k;
        for (j=k+1;j<=n-1;j++)
          { u=k*n+j;
            p=ar[u]*ar[v]; q=-ai[u]*ai[v];
            s=(ar[v]-ai[v])*(ar[u]+ai[u]);
            ar[u]=(p-q)/d; ai[u]=(s-p-q)/d;
          }
        p=br[k]*ar[v]; q=-bi[k]*ai[v];
        s=(ar[v]-ai[v])*(br[k]+bi[k]);
        br[k]=(p-q)/d; bi[k]=(s-p-q)/d;
        for (i=k+1;i<=n-1;i++)
          { u=i*n+k;
            for (j=k+1;j<=n-1;j++)
              { v=k*n+j; l=i*n+j;
                p=ar[u]*ar[v]; q=ai[u]*ai[v];
                s=(ar[u]+ai[u])*(ar[v]+ai[v]);
                ar[l]=ar[l]-p+q;
                ai[l]=ai[l]-s+p+q;
              }
            p=ar[u]*br[k]; q=ai[u]*bi[k];
            s=(ar[u]+ai[u])*(br[k]+bi[k]);
            br[i]=br[i]-p+q; bi[i]=bi[i]-s+p+q;
          }
      }
    u=(n-1)*n+n-1;
    d=ar[u]*ar[u]+ai[u]*ai[u];
    if (d+1.0==1.0)
      { free(js); printf("err**fail\n");
        return(0);
      }
    p=ar[u]*br[n-1]; q=-ai[u]*bi[n-1];
    s=(ar[u]-ai[u])*(br[n-1]+bi[n-1]);
    br[n-1]=(p-q)/d; bi[n-1]=(s-p-q)/d;
    for (i=n-2;i>=0;i--)
    for (j=i+1;j<=n-1;j++)
      { u=i*n+j;
        p=ar[u]*br[j]; q=ai[u]*bi[j];
        s=(ar[u]+ai[u])*(br[j]+bi[j]);
        br[i]=br[i]-p+q;
        bi[i]=bi[i]-s+p+q;
      }
    js[n-1]=n-1;
    for (k=n-1;k>=0;k--)
      if (js[k]!=k)
        { p=br[k]; br[k]=br[js[k]]; br[js[k]]=p;
          p=bi[k]; bi[k]=bi[js[k]]; bi[js[k]]=p;
        }
    free(js);
    return(1);
  }



/*********** Example Start **********/
/*

  #include "stdio.h"
  #include "acgas.c"
  main()
  { int i;
    static double ar[4][4]={ {1.0,3.0,2.0,13.0},
                             {7.0,2.0,1.0,-2.0},
                             {9.0,15.0,3.0,-2.0},
                             {-2.0,-2.0,11.0,5.0}};
    static double ai[4][4]={ {3.0,-2.0,1.0,6.0},
                             {-2.0,7.0,5.0,8.0},
                             {9.0,-3.0,15.0,1.0},
                             {-2.0,-2.0,7.0,6.0}};
    static double br[4]={2.0,7.0,3.0,9.0};
    static double bi[4]={1.0,2.0,-2.0,3.0};
    if (acgas(ar,ai,4,br,bi)!=0)
      for (i=0;i<=3;i++)
        printf("b(%d)=%13.7e  +j %13.7e\n",i,br[i],bi[i]);
  }




*/
/*********** Example  End **********/

/*********** Result  Start **********/
/*
b(0)= 6.782330e-02  +j  7.078231e-02
b(1)=-1.623413e-01  +j -7.612936e-01
b(2)= 5.985239e-01  +j -4.371312e-01
b(3)= 2.464562e-01  +j  1.139959e-01



*/
/*********** Result  End **********/


  #include "math.h"
  #include "stdio.h"
  int achol(a,n,m,d)
  int n,m;
  double a[],d[];
  { int i,j,k,u,v;
    if ((a[0]+1.0==1.0)||(a[0]<0.0))
      { printf("fail\n"); return(-2);}
    a[0]=sqrt(a[0]);
    for (j=1; j<=n-1; j++) a[j]=a[j]/a[0];
    for (i=1; i<=n-1; i++)
      { u=i*n+i;
        for (j=1; j<=i; j++)
          { v=(j-1)*n+i;
            a[u]=a[u]-a[v]*a[v];
          }
        if ((a[u]+1.0==1.0)||(a[u]<0.0))
          { printf("fail\n"); return(-2);}
        a[u]=sqrt(a[u]);
        if (i!=(n-1))
          { for (j=i+1; j<=n-1; j++)
              { v=i*n+j;
                for (k=1; k<=i; k++)
                  a[v]=a[v]-a[(k-1)*n+i]*a[(k-1)*n+j];
                a[v]=a[v]/a[u];
              }
          }
      }
    for (j=0; j<=m-1; j++)
      { d[j]=d[j]/a[0];
        for (i=1; i<=n-1; i++)
          { u=i*n+i; v=i*m+j;
            for (k=1; k<=i; k++)
              d[v]=d[v]-a[(k-1)*n+i]*d[(k-1)*m+j];
            d[v]=d[v]/a[u];
          }
      }
    for (j=0; j<=m-1; j++)
      { u=(n-1)*m+j;
        d[u]=d[u]/a[n*n-1];
        for (k=n-1; k>=1; k--)
          { u=(k-1)*m+j;
            for (i=k; i<=n-1; i++)
              { v=(k-1)*n+i;
                d[u]=d[u]-a[v]*d[i*m+j];
              }
            v=(k-1)*n+k-1;
            d[u]=d[u]/a[v];
          }
      }
    return(2);
  }



/*********** Example Start **********/
/*

  #include "stdio.h"
  #include "achol.c"
  main()
  { int i;
    static double a[4][4]={ {5.0,7.0,6.0,5.0},
    {7.0,10.0,8.0,7.0},{6.0,8.0,10.0,9.0},{5.0,7.0,9.0,10.0}};
    static double d[4][2]={ {23.0,92.0},{32.0,128.0},
                            {33.0,132.0},{31.0,124.0}};
    if (achol(a,4,2,d)>0)
      for (i=0; i<=3; i++)
        printf("x(%d)=%13.7e,   %13.7e\n",i,d[i][0],d[i][1]);
  }




*/
/*********** Example  End **********/

/*********** Result  Start **********/
/*
x(0)= 1.000000e+00,    4.000000e+00
x(1)= 1.000000e+00,    4.000000e+00
x(2)= 1.000000e+00,    4.000000e+00
x(3)= 1.000000e+00,    4.000000e+00



*/
/*********** Result  End **********/


  #include "stdlib.h"
  #include "stdio.h"
  int acjdn(ar,ai,br,bi,n,m)
  int n,m;
  double ar[],ai[],br[],bi[];
  { int *js,l,k,i,j,is,u,v;
    double p,q,s,d;
    js=malloc(n*sizeof(int));
    for (k=0;k<=n-1;k++)
      { d=0.0;
        for (i=k;i<=n-1;i++)
        for (j=k;j<=n-1;j++)
          { u=i*n+j;
            p=ar[u]*ar[u]+ai[u]*ai[u];
            if (p>d) {d=p;js[k]=j;is=i;}
          }
        if (d+1.0==1.0)
          { free(js); printf("err**fail\n");
            return(0);
          }
        if (is!=k)
          { for (j=k;j<=n-1;j++)
              { u=k*n+j; v=is*n+j;
                p=ar[u]; ar[u]=ar[v]; ar[v]=p;
                p=ai[u]; ai[u]=ai[v]; ai[v]=p;
              }
            for (j=0;j<=m-1;j++)
              { u=k*m+j; v=is*m+j;
                p=br[u]; br[u]=br[v]; br[v]=p;
                p=bi[u]; bi[u]=bi[v]; bi[v]=p;
              }
          }
        if (js[k]!=k)
          for (i=0;i<=n-1;i++)
            { u=i*n+k; v=i*n+js[k];
              p=ar[u]; ar[u]=ar[v]; ar[v]=p;
              p=ai[u]; ai[u]=ai[v]; ai[v]=p;
            }
        v=k*n+k;
        for (j=k+1;j<=n-1;j++)
          { u=k*n+j;
            p=ar[u]*ar[v]; q=-ai[u]*ai[v];
            s=(ar[v]-ai[v])*(ar[u]+ai[u]);
            ar[u]=(p-q)/d; ai[u]=(s-p-q)/d;
          }
        for (j=0;j<=m-1;j++)
          { u=k*m+j;
            p=br[u]*ar[v]; q=-bi[u]*ai[v];
            s=(ar[v]-ai[v])*(br[u]+bi[u]);
            br[u]=(p-q)/d; bi[u]=(s-p-q)/d;
          }
        for (i=0;i<=n-1;i++)
          if (i!=k)
            { u=i*n+k;
              for (j=k+1;j<=n-1;j++)
                { v=k*n+j; l=i*n+j;
                  p=ar[u]*ar[v]; q=ai[u]*ai[v];
                  s=(ar[u]+ai[u])*(ar[v]+ai[v]);
                  ar[l]=ar[l]-p+q;
                  ai[l]=ai[l]-s+p+q;
                }
              for (j=0;j<=m-1;j++)
                { l=i*m+j; v=k*m+j;
                  p=ar[u]*br[v]; q=ai[u]*bi[v];
                  s=(ar[u]+ai[u])*(br[v]+bi[v]);
                  br[l]=br[l]-p+q; bi[l]=bi[l]-s+p+q;
                }
            }
      }
    for (k=n-1;k>=0;k--)
      if (js[k]!=k)
        for (j=0;j<=m-1;j++)
          { u=k*m+j;v=js[k]*m+j;
            p=br[u]; br[u]=br[v]; br[v]=p;
            p=bi[u]; bi[u]=bi[v]; bi[v]=p;
          }
    free(js);
    return(1);
  }



/*********** Example Start **********/