#include<stdio.h>
#include<stdlib.h>
#include<math.h>
void PrintM2(FILE *fp,double M[],int n,int t,char *fmt,char *title,bool IsLabel)
{
int i,j;
if(title) fprintf(fp,"\n%s:",title);
int index=0;
for(i=0;i<n;i++)
{
if(IsLabel) fprintf(fp,"\n%3d",i+1);
else fprintf(fp,"\n");
for(j=0;j<=i;j++)
{
if(j>0&&j%t==0) fprintf(fp,"\n");
fprintf(fp,fmt,M[index++]);
}
}
fprintf(fp,"\n");
}
bool inverse(double a[],int n)
{
int i,j,k;
double *a0=new double[n];
for(k=0;k<n;k++)
{
double a00=a[0];
if(a00+1.0==1.0)
{
delete[]a0;
return false;
}
for(i=1;i<n;i++)
{
double ai0=a[i*(i+1)/2];
if(i<=n-k-1) a0[i]=-ai0/a00;
else a0[i]=ai0/a00;
for(j=1;j<=i;j++)
a[(i-1)*i/2+j-1]=a[i*(i+1)/2+j]+ai0*a0[j];
}
for(i=1;i<n;i++)
a[(n-1)*n/2+i-1]=a0[i];
a[n*(n+1)/2-1]=1.0/a00;
}
delete []a0;
return true;
}
int ij(int i,int j)
{
return (i>=j)?i*(i+1)/2+j:j*(j+1)/2+i;
}
double Adjust(int n,int t,double A[],double L[],double P[],double X[],double N[],double V[])
{
int i,j,k;
double *U=new double[t];
for(i=0;i<t*(t+1)/2;i++)
N[i]=0.0;
for(j=0;j<t;j++)
U[j]=0.0;
for(k=0;k<n;k++)
{
for(i=0;i<t;i++)
{
U[i]+=A[k*t+i]*P[k]*L[k];
for(j=0;j<=i;j++)
N[ij(i,j)]+=A[k*t+i]*P[k]*A[k*t+j];
}
}
if(inverse(N,t)==false)
{
delete []U;
return -1.0;
}
for(i=0;i<t;i++)
{
X[i]=0.0;
for(j=0;j<t;j++)
X[i]+=N[ij(i,j)]*U[j];
}
delete []U;
double pvv=0;
for(i=0;i<n;i++)
{
double vi=-L[i];
for(j=0;j<t;j++)
vi+=A[i*t+j]*X[j];
V[i]=vi;
pvv+=vi*vi*P[i];
}
if(n==t)
return 0.0;
return sqrt(pvv/(n-t));
}
void main()
{
int i,j,n,t;
double *A=new double[6*3];
double *L=new double[n];
double *P=new double[n];
double *X=new double[t];
double *N=new double[t*(t+1)/2];
double *V=new double[n];
for(i=0;i<n;i++)
{
for(j=0;j<t;j++)
{fscanf(fp,"%lf",A+i*t+j);}
fscanf(fp,"%lf",L+i);
}
for(i=0;i<n;i++)
{
fscanf(fp,"%lf",P+i);
}
fclose(fp);
fp=fopen("result.txt","w");
fprintf(fp,"===参数平差（观测值独立）计算===\n");
fprintf(fp,"\n误差方程：\n");
for(i=0;i<n;i++)
{
for(j=0;j<t;j++)
fprintf(fp,"%7.3lf",A[i*t+j]);
fprintf(fp,"%7.3lf\n",L[i]);
}
double m=Adjust(n,t,A,L,P,X,N,V);
fprintf(fp,"参数平差值及其中误差：\n");
for(i=0;i<t;i++)
{
fprintf(fp,"%4d%10.4lf%10.4lf\n",i+1,X[i],sqrt(N[ij(i,i)])*m);
}
fprintf(fp,"\n最小二乘残差（V）:\n");
for(i=0;i<n;i++)
fprintf(fp,"%5d%10.3lf\n",i+1,V[i]);
fprintf(fp,"单位权中误差:u=%.3lf",m);
PrintM2(fp,N,n,t,"%9.6lf","参数平差值的权逆阵",false);
fclose(fp);
}
