dscz.c

求解差值真的范围场馆哟个自行车不是你沙发是但书

#include <stdio.h>
#include <math.h>
#include <malloc.h>
double dscz(double x[],double y[],int x0,int n)
{
 int i,j,m;
 double *x1,*a,y0;
 x1=calloc(n*n,sizeof(double));
 a=calloc(n,sizeof(double));
 for(i=0;i<n;i++)
     {m=x[i];
      for(j=0;j<n;j++)
        *(x1+i*n+j)=m^j;}
 ddqjf(x1,y,a,n);
 y0=*(a+n-1);
 for(i=n-2;i>=0;i--)
     y0=x0*y0+*(a+i);
 free(x1);
 free(a);
 return y0;
}


ddqjf(double *a,double b[],double x[],int n)     /*迭代求精法*/
{
 double *a0,*b0,*r,*e,amin,alphe,epsmax;
 int i,j,k,kmax;
 a0=calloc(n*n,sizeof(double));
 b0=calloc(n,sizeof(double));
 r=calloc(n,sizeof(double));
 e=calloc(n,sizeof(double));
 kmax=1000;
 amin=fabs(*a);
 for(i=0;i<n;i++)
    {*(b0+i)=b[i];
     for(j=0;j<n;j++)
         {*(a0+i*n+j)=*(a+i*n+j);
         if(amin>fabs(*(a0+i*n+j))) amin=fabs(*(a+i*n+j));
         }
     }
 alphe=0.2*amin*amin;
 if(alphe>0.002)
    alphe=0.0;
 for(i=0;i<n;i++)
    *(a+i*n+j)=(*(a0+i*n+j))+alphe;
 gaoi(a,b,x,n);
 k=0;
 do
 {
  k=k+1;
  epsmax=0.0;
  for(i=0;i<n;i++)
     for(j=0;j<n;j++)
        {if(i==j) *(a+i*n+j)=(*(a0+i*n+j))+alphe;
         else     *(a+i*n+j)=*(a0+i*n+j);
        }
  rk(a0,b0,x,r,n);
  gaoi(a,r,e,n);
  for(i=0;i<n;i++)
    {x[i]=x[i]+e[i];
     if(epsmax<fabs(e[i])) epsmax=fabs(e[i]);
    }
 }while(epsmax>1.0e-5&&k<kmax);
 free(a0);
 free(b0);
 free(r);
 free(e);
 if(k==kmax)
 printf("ddqjf k=%d",kmax);
}

rk(double *a0,double b0[],double x[],double r[],int n)
{
 int i,j;
 double t;
 for(i=0;i<n;i++)
    {t=0;
     for(j=0;j<n;j++)
         t+=(*(a0+i*n+j))*x[j];
     r[i]=b0[i]-t;
    }
}

gaoi(double *a,double r[],double x[],int n)      /*普通线性方程求解*/
{
 int i,j,k;
 double me,c;
 for(k=0;k<n-1;k++)
    {fini(a,r,n,k);
     for(i=k+1;i<n;i++)
     {c=(*(a+i*n+k))/(*(a+k*n+k));r[i]=r[i]-r[k]*c;
     for(j=k+1;j<n;j++)
        *(a+i*n+j)=(*(a+i*n+j))-(*(a+k*n+j))*(*(a+i*n+k))/(*(a+k*n+k));
     }
    }
 x[n-1]=r[n-1]/(*(a+(n-1)*n+n-1));
 for( k=n-2;k>=0;k--)
    {me=0;
     for(j=k+1;j<n;j++)
         me=me+(*(a+k*n+j))*x[j];
     x[k]=(r[k]-me)/(*(a+k*n+k));
    }
}

fini(double *a,double b[],int n,int k)
{
 int r=k;
 double amax=fabs(*(a+k*n+k));
 int i,j,t;
 for(i=k;i<n;i++)
 {
 if(fabs(*(a+i*n+k))>amax)
 {
 amax=fabs(*(a+i*n+k));
 r=i;
 }
 }
 if(amax==0)
 {
 printf("Can not get the answer.\n");
 exit(0);
 }
 if(r!=k)
 {
 for(j=k;j<n;j++)
 {
 t=*(a+k*n+j);
 *(a+k*n+j)=*(a+r*n+j);
 *(a+r*n+j)=t;
 }
 t=b[k];
 b[k]=b[r];
 b[r]=t;
 }
}