C语言间接平差程序

教材《误差理论与测量平差基础》第二版武汉大学出版社P108页的例7—1的运行结果：源程序：＃define N 5 /*N是观测值个数*/#define T 3 /*T是必要观测数*/＃include<stdio.h>#include<math。

h>float Nbb［T］［T］,Nb[T］[T]，W[T］［1],x［T]［1］；main（）｛float D（float a［T]［N］，float b［N][N］，float c[N][T]）;float K（float a[T][N］,float b［N］[N］，float c［N］[1])；float G(float a[T］[T］）；float F（float ca［T-1］[T—1］);float DM(float a[1][N］，float b［N］[N] ，float c［N］[1]）；int i,j,m,n;float B［N]［T］,BT［T］［N],V［N]［1］,VT［1］[N］，P［N］［N],C[N］[1]，Bx[N][1],f，g，h，x1; printf（"请输入V的系数B［N]［T]：\n”);for（i=0；i<N;i++）for(j=0;j〈T;j++)scanf（”％8f”，＆B［i］[j］)；printf（"请输入观测值的权阵P［N］[N］:\n");for（i=0；i<N；i++)for(j=0;j<N；j++)scanf（"%8f”,＆P[i]［j])；printf（"请输入常数C［N][1］:\n”）；for(i=0;i〈N；i++）for（j=0;j<1；j++)scanf(”％8f”，＆C[i][j］)；for(i=0；i〈N；i++）for（j=0;j<T；j++）BT[j]［i]=B［i］［j];g=D(BT, P， B);h=K（BT， P， C）；f=G(Nbb);for（i=0；i〈T;i++）for（j=0；j〈1;j++）{x[i］［j］=Nb［i］［0]*W［0］[j]；for(m=1；m〈T；m++)x［i］［j]+=（Nb[i］［m］＊W［m］［j］);｝for（i=0;i〈T;i++）x［i］［0]=x［i］[0]/f;for（i=0;i〈N;i++)for（j=0;j〈1；j++）{Bx［i][j]=B[i］[0］＊x［0］[j]；for(m=1;m〈T；m++)Bx［i］[j]+=(B［i］［m]＊x［m][j］);｝for(i=0;i〈N；i++)V[i］[0]=（Bx[i][0]-C[i][0]）;for(i=0；i<N；i++)for(j=0；j〈1;j++)VT[j][i]=V[i]［j]；x1=DM（VT，P,V)；x1=x1/（N-T）；printf("参数x［T］［1]=\n");for(i=0；i〈T;i++）printf("%15f"，x［i］[0]);printf（"\n"）；printf("改正数V［N][1］=\n”)；for（i=0；i<N；i++)printf（"%15f”,V[i］［0］）；printf（"\n单位权中误差x1=%15f”,sqrt(x1）); printf("\n协因数阵Qxx[T］[T]:\n"）；for(i=0;i<T;i++）｛for（j=0；j〈T；j++)printf（”%15f”，Nb［i］[j］/f);printf（”\n")；｝}float G（float a[T］［T］){int i，j，m，n;float c[T—1］[T—1］,y=0；for(i=0;i<T;i++）for（j=0;j<T;j++）{for（m=0；m<T;m++)for（n=0；n〈T；n++）{if(m〈i＆＆n<j）c[m][n］=a[m］［n］;if（m>i&&n〈j)c[m-1］[n]=a［m］[n］;if（m<i＆&n>j）c［m］[n—1]=a[m］［n]；if（m>i&&n>j）c[m-1]［n-1]=a［m]［n］；｝if（（i+j)%2==0)Nb[j］[i]=F（c);elseNb[j］［i］=(-1)＊F(c）;｝for(m=0；m<T；m++)y+=(a［0][m］＊Nb[m][0］）；return (y);}float F(float ca[T—1］[T—1]){int i,j,m,n，s，t，k=1;float f=1，c，x，sn;for （i=0，j=0；i<T—1&＆j〈T—1;i++，j++）｛if （ca［i］[j］==0){for （m=i；ca[m］［j]==0；m++）;if (m==T—1）｛sn=0;return （sn）;｝elsefor （n=j;n<T—1;n++)｛c=ca[i］[n］；ca［i］[n］=ca［m］［n];ca[m］［n］=c；｝k＊=（-1)；｝for (s=T—2；s>i；s--）｛x=ca[s][j];for (t=j;t〈T-1;t++）ca［s][t］—=ca[i][t]*（x/ca［i][j])；}｝for （i=0;i<T-1;i++)f*=ca［i]［i];sn=k＊f;return (sn);}float D（float a［T］[N］,float b[N][N] ,float c［N］［T])｛int i,j，m；float d［T］[N];for（i=0;i<T；i++)for（j=0;j〈N;j++）｛d[i][j］=a［i］[0]*b［0][j］;for(m=1；m〈N;m++）d[i］[j]+=(a［i][m］*b[m］[j]）；｝for（i=0;i〈T；i++）for（j=0;j〈T；j++)｛Nbb［i]［j］=d[i][0］*c［0］［j］；for(m=1；m〈N;m++)Nbb［i][j］+=（d[i］[m]＊c［m][j])；}return （Nbb［0］［0］）；｝float K(float a[T］[N］,float b[N]［N],float c［N]［1］）｛int i，j,m;float d［T]［N］;for（i=0；i<T;i++）for(j=0;j<N;j++）{d[i］［j］=a［i]［0］＊b［0]［j]；for(m=1;m<N；m++)d[i]［j]+=（a[i］[m］＊b[m］[j])；｝for（i=0;i<T;i++）for(j=0；j〈1；j++){W[i]［j]=d[i]［0]*c［0］[j]；for(m=1;m〈N；m++)W［i］[j］+=(d［i]［m］*c[m］[j］)；}return （W［0］[0］)；｝float DM（float a[1]［N],float b[N]［N] ，float c[N][1]）{int i，j,m;float d[1][N］，x;for（i=0;i〈1;i++）for(j=0;j<N；j++）｛d[i]［j］=a［i］[0]＊b[0］[j]；for（m=1；m<N;m++）d［i][j］+=（a［i］[m］*b［m］［j]);}for(i=0；i<1;i++)for（j=0；j<1;j++)｛x=d[i］[0］＊c［0］[j]；for（m=1；m<N;m++）x+=(d［i][m］*c[m][j]);｝return (x）；}程序说明：1) 用该程序前，根据具体情况输入N和T；2）该程序中的（N-T）自由度(即多余观测数）必须大于等于2，不然程序运行时会出错，原因在于求行列式的逆时，有语句for (s=R-2；s>i；s—-)，R=1时s=-1。

间接平差法是一种测量平差方法，通常用于解决线性系统中的超定方程问题，即有多余观测的情况。

这种方法通过解最小二乘问题来找到最佳的参数估计值。

在间接平差法中，待估参数和已知参数是通过最小二乘目标函数（通常是误差项的平方和）进行连接的。

具体步骤如下：
1. 列出误差方程：误差方程是观测值与计算出的观测值初值之间的差值。

2. 计算出V后与观测值求和，得到最终的平差值。

此外，间接平差法还可以用于解决GNSS SPP、摄影测量解算、光束法平差、控制网平差等测绘问题。

在解决这些问题时，通常会使用非线性最小二乘解法中的其他方法，如梯度下降法（最速下降法）来获得最佳的参数估计值。

总之，间接平差法是一种广泛应用的测量平差方法，通过最小化目标函数来求解线性系统中的超定方程问题。

它被广泛应用于各种测绘解算问题，为参数估计提供了有效的解决方案。

测量平差程序设计1.角度（度分秒）到弧度AngleToRadian#define PI 3.14159265double AngleToRadian（double angle）{int D，M；double S,radian,degree, angle,MS；D=int（angle+0.3）；MS=angle-D；M=int（（MS）*100+0.3）；S=（MS*100-M）*100;degree=D+M/60.0+S/3600.0;radian=degree*PI/180.0;return radian;}注意:防止数据溢出,要加个微小量,例如0.3.2.弧度换角度(度分秒) RadianToAngle#define PI 3.14159265double RadianToAngle(double radian){int D,M;double S,radian,degree,MS,angle;degree=radian*180/PI;D=int(degree);MS=degree-D;M=int(MS*60);S=(MS*60-M)*60;angle=D+M/100.0+S/10000.0;return angle;}3.已知两点求坐标方位角Azimuth#include <math.h>double Azimuth(double xi,double yi,double xj,double yj) {double Dx,Dy,S,T;Dx=xj-xi;Dy=yj-yi;S=sqrt(Dx*Dx+Dy*Dy);if(S<1e-10) return 0;T=asin(Dy/S);if(Dx<0) T=PI-T;if(Dx>0&&(Dy<0)||T<0) T=2*PI+T;return T;}4.开辟二维数组的动态空间的宏#include <malloc.h>#define NewArray2D(type,A,i,n,m){A=(type**)malloc(n*sizeof(type*));\for(i=0;i<m;i++)\A[i]=(type*)malloc(m*sizeof(type));\}5.释放开辟的二维数组的空间#define FreeSpace(A,i,m){for(i=0;i<m;i++)\free(A[ i]);\free(A);\}注意:释放空间与开辟空间相反,释放空间是先释放列,后释放行.6.矩阵求转置transformmatrixvoid transformmatrix（double **A，double **B，int i,int j）{int m,n;for(m=0;m<=i;m++)for(n=0;n<=j;n++){B[n][m]=A[m][n]:}}7.矩阵相乘(mulmatrix)void mulmatrix(double **A,double **B,double **C,int i,int j,int k) {int m,n,p;for(m=0;m<i;m++)for(n=0;n<j;n++){C[m][n]=0;for(p=0;p<k;p++){C[m][n]+=A[m][p]*B[p][n]:}}}8.矩阵求逆(countermatrix)#include <math.h>void countermatrix(double **T, double **s, double **r, double **Q,double **N, double **rt,int n){for(i=0;i<n;i++){s=N[i][i];for(k=0;k<i;k++){s-=T[k][i]*T[k][i];}T[i][i]=sqrt(s)for(j=i+1;j<n;j++){s=N[i][j];for(k=0;k<i;k++){s-=T[k][i]*T[k][j];}T[i][j]=s/T[i][i];}}for(i=0;i<n;i++)for(j=0;j<n;j++){T[i][j]=0;}for(i=n-1;i>=0;i++){r[i][i]=1/T[i][i];for(j=i+1;j<n;j++){s=0;for(k=i;k<j-1;k++){s-=r[i][k]*T[k][j];}r[i][j]=s/T[i][i];}}for(i=0;i<n;i++)for(j=0;j<n;j++){r[i][j]=0;}transformmatrix(r,rt,n,n)mulmatrix(r,rt,Q,n,n)}9.平差主程序之读入数据typedef struct POINT{char name[8];double x,y;int type;}POINT;typedef struct READV ALUE{POINT *begin;POINT *end;double value;}READV ALUE;POINT *GETPOINT(char *name,POINT *pPoint,int nPoint) {int i;for(i=0;i<nPoint;i++){if (strcmp(pPoint[i].name,name)==0)return (pPoint+i)}for(i=0;i<nPoint;i++){if(pPoint[i]=NULL)strcmp(pPoint[i].name,name);pPoint[i].type=0;return(pPoint+i);}}double AngleToRadian（double angle）{int D，M；double S,radian,degree, angle,MS；D=int（angle+0.3）；MS=angle-D；M=int（（MS）*100+0.3）；S=（MS*100-M）*100;degree=D+M/60.0+S/3600.0;radian=degree*PI/180.0;return radian;}main(){POINT *pPoint=NULL;READV ALUE *pDirect=NULL;READV ALUE *pDistance=NULL;int nPoint,nKnownPoint,nDirect,nDistance,i;double mo,mf,ms;char begin[8],end[8];FILE *fp=0;fp=fopen(“c:\\dat\\t1.txt”,”r”)fscanf(fp,”%d,%d,%d,%d\n”,&nPoint,&nKnowPoint,&nDirect,&nDistance) if(nPoint>0)pPoint=(POINT*)malloc(nDirect*sizeof(POINT));if(nDirect>0)pDirect=(READV ALUE*)malloc(nDirect*sizeof(READV ALUE));if(nDistance>0)pDistance=(READV ALUE*)malloc(nDistance*sizeof(RAADV ALUE));fscanf(fp,”%lf,%lf,%lf\n”,&mo,&mf,&ms);for(i=0;i<nKnownPoint;i++){fscanf(fp,”%s,%lf,%lf\n”,pPoint[i].name,&pPoint[i].x,&pPoint[i].y);type=1;}for( ;i<nPoint;i++){pPoint[i].name=NULL;pPoint[i].x=0;pPoint[i].y=0;pPoint[i].type=0;}for(i=0;i<nDirect;i++){fscanf(fp,”%s,%s,%lf\n”,begin,end,&pDirect[i].value);pDirect[i].begin=GetPoint(begin,pPoint,nPoint);pDirect[i].end=GetPoint(end,pPoint,nPoint);}for(i=0;i<nDistance;i++){fscanf(fp,”%s,%s,%lf\n”,begin,end,&pDistance[i].value);pDistance[i].begin=GetPoint(begin,pPoint,nPoint);pDistance[i].end=GetPoint(end,pPoint,nPoint);}fclose(fp);}10.角度检验（checkangle）#include <math.h>int checkangle(double angle){int M,S;double MS;if(angle>=0&&angle<360){MS=angle-(int)(angle);if(M<6){S=(int)(MS*1000);if(S%10<6){return 1;}}}return 0;}11.前方交会#define PI=3014159265/***此处调用程序角度换弧度AngleToRadian***/Qianfang(double XE, double YE, double XF, double YF, doubleDEG, double DEF, double DFG, double DFE, double *DFE, double *DFG){double C,A,B;C=DGE-DGF;A=DEF-DEG;B=DFG-DFE;if((C<-PI&&C>-2*PI)||(C>0&&C<PI){XG=(XE/tan(B)+XF/tan(A)-YE+YF)/(1/tan(A)+ 1/tan(B);YG=(YE/tan(B)+YF/tan(A)+XE-XF)/ (1/tan(A)+ 1/tan(B);}if((C>-PI&&C<0)||(C>PI&&C<2*PI)){XG=(XE/tan(B)+XF/tan(A)+YE-YF)/(1/tan(A)+ 1/tan(B);YG=(YE/tan(B)+YF/tan(A)-XE+XF)/ (1/tan(A)+ 1/tan(B);}}12.坐标概算全方向法子函数取出观测方向GetAllDirectint GetAllDirect(char *name,int nDirect,READV ALUE *pDirect, READV ALUE *pStation) {int i,nCount=0;for(i=0;i<nDirect;i++)if(strcmp(pDirect[i].begin->name,name)==0)){pStation[nCount].begin=p(pDirect[nCount].begin;pStation[nCount].end=p(pDirect[nCount].end;pStation[nCount].value=p(pDirect[nCount].value;nCount++;}return nCount;}坐标概算全方向法子程序实现流程(coordinate)coordinate (入口参数设置){READV ALUE pStation[50],pObject[50];int nCount,i,j,k,m,n,p,nobject;for(i=0;i<nPoint;i++){nCount=GetAllDirect(pPoint[i].name,nDirect,pStation)if((nCount>1)||( nCount=1)){for(j=0;j<nCount;j++){if(pStation[j].end->type==1){for(k=0;k<nCount;k++){if(pStation[k].end->type==0)nobject=GetAllDirect(pStation[j].end->name,nDirect,pDirect,pobject)m=-1;n=-1;for(p=0;p<nobject;p++){if(strcmp(pobject[p].end->name,pPoint[i].name)==0){m=p;}if(strcmp(pobject[p].end->name,pStation[k].end->name)==0){n=p;}if(m>=0&&n>=0){pPoint[i]=pStation[k].end-pStation[j].end;pStation[j].end=pObject[m].value-pObject[n].value;{Xe=pPoint[i].x;Ye=pPoint[i].y;Xf=pStation[j].end->x;Yf=pStation[j].end->y;Lef=pStation[j].value;Leg=pStation[k].value;Lfe=pObject[m].value;Lfg=pObject[n].value;Qianfang(Xe,Xf,Ye,Yf,Lef,Leg,Lfe,Lfg,*Xg,*Yg;)pStation[k].end->x=*xg;pStation[k].end->y=*yg;pStation[k].end.type=2;}}}}}}}}13.坐标增量法(calcoordinate)子函数由端点名称得边长值的函数GetDistancedouble GetDistance(char *begin,char *end,int nDistance,READV ALUE *pDistance){int i;for(i=0;i<nDistance;i++){if((strcmp(pDistance[i].begin->name,begin)==0&&strcmp(pDistance[i].end->name,end==0)| |(strcmp(pDistance[i].begin->name,end)==0&&strcmp(pDistance[i].end,begin)==0))return pDistance[i].value;}return -1;}/***函数取出观测方向GetAllDirect***/void calcoordinate(int nDirect,READV ALUE *pDirect,int nDistace,READV ALUE *pDistance,int nPoint,POINT *pPoint){int nPoint,nCount,nDirect,nDistance;int m=-1,i,j,k;double x1,y1,x2,y2,A0,A,S,dx,dy;READV ALUE*pDirect=NULL;READV ALUE pStation[50];for(i=0;i<nPoint;i++){if(pPoint[i].type>0){nCount=GetAllDirect(pPoint[i].name,nDirect,pDirect,pStation[50]);for(j=0;j<nCount;j++){if(pStation[j].end->type>0)m=j;if(m!=-1){for(k=0;k<nCount;k++){if(pStation[k].end->type==0){x1=pPoint[i].x;y1=pPoint[i].y;x2=pStation[j].end->x;y2=pStation[j].end->y;A0=Bearing(x1,y1,x2,y2);A=A0-(DMSToRAD(pStation[m].value)-DMSToRAD(pStation[k].value));if(A<0)A=A+2*PI;if(A>2*PI)A=A-2*PI;S=GetDistance(pPoint[i],pStation[k].end,nDistance,pDistance);if(S<0)continue;else{dx=S*cos(A);dy=S*sin(A);pStation[k].end->x=pPoint[i].x+dx;pStation[k].end->y=pPoint[i].y+dy;pStation[k].end->type=2;}}}}}}}}14.高斯正反算高斯正算：#include <math.h>#include <stdio.h>#define PI 3.14159265double DMSToRAD（double dDMS）{int L1,L2；double T,L3；L1=(int)（dDMS+0.3）；L2=(int)((dDMS-L1)*100+0.3)；L3=((dDMS-L1)*100-L2)*100；T=(L1+L2/60.0+L3/3600.0)*PI/180.0;return T;}void PreGausePositive(double B,double L,double L0, double a, double b, double *N, double *l, double *c, double *t, double *X,double *B1){double a0,a2,a4,a6,a8,m0,m2,m4,m6,m8;double e,e1;e=(sqrt(a*a-b*b))/a;e1=(sqrt(a*a-b*b))/b;B1=DMSToRAD(B);t=tanB1;c=sqrt(e1*e1*cosB1*cos*B1);l=L-L0;N=a/(sqrt(1-e*e*sinB1*sinB1));m0=a*(1-e*e);m2=3/2*e*e*m0;m4=5/4*e*e*m2;m6=7/6*e*e*m4;m8=9/8*e*e*m6;a0=m0+m2/2+3*m4/8+5*m6/16+35*m8/128;a2=m2/2+m4/2+15*m6/32+7/16*m8;a4=m4/8+3*m6/16+7*m8/32;a6=m6/32+m8/16;a8=m8/128;X=a0*B1-a2*(sin(2*B1))/2+a4*(sin(4*B1))/4-a6*(sin(6*B1))/6+a8*(sin(8*B1))/8;}Void BLToXY(double *x,double *y,double N,double l,double c,double t,double B1,double X)x=X+N*l*l*t*cosB1*cosB1*((3+l*l*cosB1*cosB1*(5-t*t+9*c*c+4*c*c*c*c)/4+l*l*\cosB1*cosB1*(61-58*t*t+t*t*t*t)/30))/6;y=N*l*cosB1(1+l*l*cosB1*((1+c*c-t*t)+l*l*cosB1*cosB1(5-18*t*t+t*t*t*t+14*c*c\-58*t*t*c*c)));}高斯反算void XYToBL(double x,double y,double L0,double a,double b,double q,double *B,\double *L){double Bf,c,t,y,N,e1,ee=(sqrt(a*a-b*b))/(a*a);e1=(sqrt(a*a-b*b))/(b*b);for(Bf=0;;){t=tanBf;c=e1*e1*cosBf;N=a/(sqrt(1-e*e*sinBf*sinBf));B=Bf-(1+c*c)*t*y*y/(2*N*N)*(1-y*y)/(12*N*N)*(15+3*t*t+c*c-9*t*t*c*c)-y*y/ (30*N*N)*(61+90*t*t+45*t*t*t*t);if(fabs(B-Bf)<q)break;Bf=B;}L=L0+y/(N*cosBf)*(1-y*y/(6*N*N))*((1+2*t*t+c*c)-y*y/(20*N*N)*(5+6*c*c)+28*t* t+24*t*t*t*t+8*t*t*c*c)-y*y/(42*N*N)*(61+662*t*t+1320*t*t*t*t+720*t*t*t*t*t*t));B=RADToDMS(B);L=RADToDMS(L);}附录资料：不需要的可以自行删除SHA算法的实现C语言程序：#include <iostream>#include <vector> //定义vector数组#include <string> //记录消息using namespace std;const int NUM = 8; //一个字由32比特(或者8个16进制数)const int BIT = 512; //消息认证码要以512比特一组//字常量string H0 = "67452301";string H1 = "EFCDAB89";string H2 = "98BADCFE";string H3 = "10325476";string H4 = "C3D2E1F0";//定义SHA1(安全哈希算法)类class SHA1{public://将一个字符串形式的字转化为vector数组vector<int> hex_into_dec(string word);//将vector转化为string字符串形式string num_into_message(vector<int> A);//两个字X和Y的逻辑"和"vector<int> word_AND(vector<int> A,vector<int> B);//两个字X和Y的逻辑"或"vector<int> word_OR(vector<int> A,vector<int> B);//两个字X和Y的逻辑"异或"vector<int> word_XOR(vector<int> A,vector<int> B);//两个字X和Y的逻辑"补"vector<int> word_COMPLEMENT(vector<int> A);//两个字X和Y的摸2^32整数加vector<int> word_ADD(vector<int> A,vector<int> B);//将字X循环左移s个位置vector<int> ROTL(vector<int> A,int s);//SHA-1的填充方案,我们设定msg由ASCII码组成vector<vector<int> > SHA_1_PAD(string msg);//将SHA-1压成以字为单位vector<vector<vector<int> > > compress(vector<vector<int> > result);//定义ft函数，每个ft函数都有B,C,D三个字作为输入，并产生一个字作为输出vector<int> Ft(int t,vector<int> B,vector<int> C,vector<int> D);//定义字常数Kvector<int> K(int t);//开始进行SHA-1(安全Hash算法)的加密vector<vector<int> > SHA_1(string msg);};//将vector转化为string字符串形式string SHA1::num_into_message(vector<int> A){int i;string msg = "";for(i = 0;i < A.size();i++){if(A[i] >= 0 && A[i] <= 9)msg += '0' + A[i];else if(A[i] >= 10 && A[i] <= 15)msg += 'A' + (A[i] - 10);}return msg;}//将一个字符串形式的字转化为vector数组vector<int> SHA1::hex_into_dec(string word){int i;vector<int> result(NUM,0);for(i = 0;i < NUM;i++){if(word[i] >= '0' && word[i] <= '9'){result[i] = word[i] - '0';}else if(word[i] >= 'A' && word[i] <= 'F'){result[i] = 10 + word[i] - 'A';}}return result;}//两个字X和Y的逻辑"和"vector<int> SHA1::word_AND(vector<int> A,vector<int> B) {vector<int> result(NUM,0);int i;for(i = 0;i < NUM;i++){result[i] = A[i] & B[i];}return result;}//两个字X和Y的逻辑"或"vector<int> SHA1::word_OR(vector<int> A,vector<int> B) {vector<int> result(NUM,0);int i;for(i = 0;i < NUM;i++){result[i] = A[i] | B[i];}return result;}//两个字X和Y的逻辑"异或"vector<int> SHA1::word_XOR(vector<int> A,vector<int> B) {vector<int> result(NUM,0);int i;for(i = 0;i < NUM;i++){result[i] = A[i] ^ B[i];}return result;}//两个字X和Y的逻辑"补"vector<int> SHA1::word_COMPLEMENT(vector<int> A) {vector<int> result(NUM,0);int i;for(i = 0;i < NUM;i++){result[i] = 15 - A[i];}return result;}//两个字X和Y的摸2^32整数加vector<int> SHA1::word_ADD(vector<int> A,vector<int> B) {vector<int> result(NUM,0);int i;for(i = NUM - 1;i >= 0;i--){result[i] = A[i] + B[i];if(i != 0){int temp = result[i] / 16;result[i-1] += temp;}result[i] %= 16;}return result;}//将字X循环左移s个位置vector<int> SHA1::ROTL(vector<int> A,int s){vector<int> result = A;vector<int> temp(NUM,0);int i,j;for(i = 0;i < s;i++){for(j = NUM - 1;j >= 0;j--){if(result[j] / 8 >= 1){temp[j] = 1;result[j] <<= 1;result[j] %= 16;if(j < NUM - 1)result[j] += temp[j + 1];}else if(result[j] / 8 == 0){temp[j] = 0;result[j] <<= 1;result[j] %= 16;}}result[NUM - 1] += temp[0];}return result;}//SHA-1的填充方案,我们设定msg由ASCII码组成vector<vector<int> > SHA1::SHA_1_PAD(string msg){int len = msg.length();int bit_num = len * 8;int i,j;int num,lest = bit_num % 512;if(lest != 0) //看消息长度是否超过512字节，我们需要将它补成512的倍数num = bit_num / 512 + 1;elsenum = bit_num / 512;//首先我们以8位字节为一组保存到vector里面，512比特为一组，即一组里面有64位元素vector<vector<int> > result;result.resize(num);for(i = 0;i < num;i++){result[i].resize(64);}for(i = 0;i < num;i++){for(j = 0;j < 64 && i * 64 + j < len;j++){result[i][j] = msg[i * 64 + j];}}//下面开始为未够512比特的消息分组进行补长度操作if(lest != 0){int x = num - 1,last_len = lest / 8;result[x][last_len] = 128; //先补一个"1"for(i = last_len + 1;i < 56;i++){result[x][i] = 0;}int last_l = lest;j = 63;while(j >= 56){result[x][j] = last_l % 128;last_l /= 128;j--;}}return result;}//将SHA-1压成以字为单位(三维数组有点复杂)vector<vector<vector<int> > > SHA1::compress(vector<vector<int> > result) {vector<vector<int> > rr;rr.resize(result.size());int i,j;for(i = 0;i < rr.size();i++){rr[i].resize(128);}for(i = 0;i < result.size();i++){for(j = 0;j < result[i].size();j++){rr[i][2 * j] = result[i][j] / 16;rr[i][2 * j + 1] = result[i][j] % 16;}}vector<vector<vector<int> > > rrr;rrr.resize(result.size());for(i = 0;i < rrr.size();i++){rrr[i].resize(16);}for(i = 0;i < rrr.size();i++){for(j = 0;j < 16;j++){rrr[i][j].resize(8);}}for(i = 0;i < rr.size();i++){for(j = 0;j < rr[i].size();j++){rrr[i][j / 8][j % 8] = rr[i][j];}}return rrr;}//定义ft函数，每个ft函数都有B,C,D三个字作为输入，并产生一个字作为输出vector<int> SHA1::Ft(int t,vector<int> B,vector<int> C,vector<int> D){vector<int> result;if(t >= 0 && t <= 19){vector<int> a1 = word_AND(B,C);vector<int> a2 = word_AND(word_COMPLEMENT(B),D);result = word_OR(a1,a2);}else if((t >= 20 && t <= 39) || (t >= 60 && t <= 79)){vector<int> a1 = word_XOR(B,C);result = word_XOR(a1,D);}else if(t >= 40 && t <= 59){vector<int> a1 = word_AND(B,C);vector<int> a2 = word_AND(B,D);vector<int> a3 = word_AND(C,D);vector<int> a4 = word_OR(a1,a2);result = word_OR(a4,a3);}return result;}//定义字常数Kvector<int> SHA1::K(int t){vector<int> result;if(t >= 0 && t <= 19){result = hex_into_dec("5A827999");}else if(t >= 20 && t <= 39){result = hex_into_dec("6ED9EBA1");}else if(t >= 40 && t <= 59){result = hex_into_dec("8F1BBCDC");}else if(t >= 60 && t <= 79){result = hex_into_dec("CA62C1D6");}return result;}//开始进行SHA-1(安全Hash算法)的加密vector<vector<int> > SHA1::SHA_1(string msg){vector<int> h0 = hex_into_dec(H0);vector<int> h1 = hex_into_dec(H1);vector<int> h2 = hex_into_dec(H2);vector<int> h3 = hex_into_dec(H3);vector<int> h4 = hex_into_dec(H4);vector<vector<int> > result1 = SHA_1_PAD(msg);vector<vector<vector<int> > > result2 = compress(result1);int n = result2.size();int i,j;for(i = 0;i < n;i++){vector<vector<int> > W;W.resize(80);for(j = 0;j < 16;j++){W[j] = result2[i][j];}for(j = 16;j < 80;j++){vector<int> a1 = word_XOR(W[j-3],W[j-8]);vector<int> a2 = word_XOR(a1,W[j-14]);vector<int> a3 = word_XOR(a2,W[j-16]);W[j] = ROTL(a3,1);}//将string转化为vector数组vector<int> A = hex_into_dec(H0);vector<int> B = hex_into_dec(H1);vector<int> C = hex_into_dec(H2);vector<int> D = hex_into_dec(H3);vector<int> E = hex_into_dec(H4);for(j = 0;j < 80;j++){vector<int> a1 = ROTL(A,5);vector<int> a2 = Ft(j,B,C,D);vector<int> a3 = word_ADD(a1,a2);vector<int> a4 = word_ADD(a3,E);vector<int> a5 = word_ADD(a4,W[j]);vector<int> temp = word_ADD(a5,K(j));E = D;D = C;C = ROTL(B,30);B = A;A = temp;}h0 = word_ADD(h0,A);h1 = word_ADD(h1,B);h2 = word_ADD(h2,C);h3 = word_ADD(h3,D);h4 = word_ADD(h4,E);}//返回结果(H0||H1||H2||H3||H4)vector<vector<int> > result;result.push_back(h0);result.push_back(h1);result.push_back(h2);result.push_back(h3);result.push_back(h4);return result;}int main(){SHA1 sha1; //定义SHA1算法类string message = "cryptographyisthepracticeandstudyoftechniquesforsecurecommunicationinthepresenceofthirdpartiesmoregenerallyitisaboutconstructingandanalyzingprotocolsthatovercome theinfluenceofadversariesandwhicharerelatedtovariousaspectsininformationsecuritysu chasdataconfidentialitydataintegrityauthenticationandnonrepudiationmoderncryptogra phyintersectsthedisciplinesofmathematicscomputerscienceandelectricalengineeringap plicationsofcryptographyincludeATMcardscomputerpasswordsandelectroniccommerc e";vector<vector<int> > result;result = sha1.SHA_1(message);cout << "消息为:" << endl << message << endl;cout << "利用填充方案SHA-1-PAD给出对消息的填充,得出SHA-1(x)得:" << endl;int i;for(i = 0;i < result.size();i++){cout << sha1.num_into_message(result[i]);}cout << endl;return 0;}程序运行结果：- 21 -。