测绘程序设计课程实习报告

1、测绘程序设计课程实习报告一、实习目的巩固之前所学的c+类定义、文件操作等各方面知识，掌握矩阵计算的基本算法，熟悉并熟练运用平差模型，掌握测量平差的基本原理和公式，并处理测量观测数据。二、实习内容该课程的实习主要为c语言及c+编程基础，矩阵操作算法，条件平差，间接平差，以及在水准网平差中的应用。三、实习过程（1） 实习1（文件的基本操作以及角度、弧度的转化）1. 准备C+文件操作函数包含于头文件<fstream><ofstream><ifstream>，一般情况下使用<fstream>即可。在计算过程中需要文件输出格式控制，函数包含于<iom

2、ainip>。c+的流操作函数包含于<iostream>，相应的流函数在命名空间“std”。角度转弧度的过程中，首先判断正负，对符号进行不同情况的处理。将度数的分，秒位分别除以60、3600，然后把得到的结果与度位绝对值相加，判断符号。弧度转角度的过程中，同样判断正负，取小数位乘以60取整数得到相应角度的分，再取剩余的小数乘以60即得到相应角度的秒。2. 代码#include <iostream>#include <math.h>#include <fstream>#include <string>using namespace

3、 std;const double pi = 3.14159;int main()fstream file("a.txt"); /定义文件操作对象string aa;file>>aa;cout<<aa<<endl; /输出文件中的内容int a = 1;/定义重复输入使用的变量while(a) /限制循环，重复输入float a1,a2,a3;cout<<endl;double deg;cout<<"输入度"<<endl;cin>>a1;cout<<"

4、;输入分"<<endl;cin>>a2;cout<<"输入秒"<<endl;cin>>a3;cout<<"角度为 "<<a1<<"° "<<a2<<" "<<a3<<""<<endl;if(a1 >= 0) /判断输入度数的正负deg = a1 + a2/60 + a3/3600;elsedeg = a1 - a2/

5、60 - a3/3600;double radian;radian = deg/180*pi;cout<<"转换得弧度为 "<<radian<<endl;cout<<endl;cout<<"输入弧度： "<<endl;cin>>radian;bool tm = 0;deg = radian*180/pi;if(deg<0)tm = 1;deg = -deg;a1 = int(deg);a2 = int( (deg-a1)*60 );a3 = ( (deg-a1-a2

6、/60)*3600 );if(tm = 1) a1 = -a1; tm = 0;cout<<"转换后角度为 "<<a1<<"° "<<a2<<" "<<a3<<""<<endl;return 0;3. 运行结果 输入度20输入分55输入秒32角度为 20° 55 32转换得弧度为 0.36522输入弧度：3.14转化后角度为 179° 54 32.0387Press any key to c

7、ontinue4. 出现的问题角度转换过程中需要考虑角度的正负号问题，角度转化可以添加是否超过59的判断。（2） 实习2（概率计算）1. 准备在误差检验、粗差探测等平差计算中经常需要计算某些特定分布函数的反函数值，为满足计算需要，使用CProbability类对相关计算进行封装。2. 代码及运行结果与书上相同。（3） 实习3（矩阵基本运算）1. 准备矩阵的运算有加法、减法、乘法、转置、求逆，封装在一个类中。为了方便之后的计算与使用，类中还添加了从文件中读入矩阵、输出矩阵，向控制台输入输出矩阵的函数。代入运算的矩阵均用动态的一维数组表示。先定义double型指针，在确定矩阵大小后动态申请空间，r

8、ow行col列的矩阵定义为：double *mat = new doublerow*col; 第i行j列的元素为mati*col+j。矩阵的加法、减法只需遍历两个矩阵的每个元素，相加之后放到结果矩阵中（如有必要先判断相加的两矩阵行数与列数是否相等）。乘法计算中设相成的两个矩阵mat1、mat2的行列分别为row1、col,col、col2，则所求结果矩阵第i行j列的元素 用三个循环即可计算矩阵中所有的元素。矩阵转置将结果矩阵中的i行j列元素赋值为需要计算的矩阵中的j行i列个元素即可。矩阵求逆，由于平差计算中需要求逆的矩阵一般为实对称正定矩阵，为了提高大型矩阵的计算效率，可以使用对实对称正定矩阵

9、求逆的方法。设C为n阶对称正定矩阵，y、x为n维向量，其关系式为： y=C·x （3-1）确定了上的一个映像，如能写出逆关系： x=B·y （3-2）则B为C的逆阵，即B=现将式（3-1）写成 . . （3-3） . .因C对称正定，必有0,用除式（3-3）第一个方程的两端，解出，把和的位置交换，并将代入其他各式得 . . （3-4） . . 事实上式（3-4）可以改写为 . + . + . （3-5） . +如果对式（3-5）中的变量按如下规则重新编号 （3-6） 经n次变换后恢复原状，采用变量循环重新编号法的计算公式如下：对于 由于变量循环重新编号法求逆均在下三角阵（包

10、括对角元素）进行，因而运行速快，在大型矩阵的计算中可以很好的适用。矩阵运算的函数封装在类mat中，函数在mat.h中定义，在mat.cpp中实现。2. 代码文件mat.h：#include <fstream>#include <iostream>#include<math.h>using namespace std;class matpublic:void add(int row,int col,double m,double n,double result);/加法void sub(int row,int col,double m,double n,dou

11、ble result);/减法double* multiply(int row, int col, int col2, double m,double n);/乘法void trans(int row, int col, double m, double result);/转置bool inverse(double a, int n);/求逆void fileIn(ifstream &in, int size, double *mat);/从文件中读入void (ofstream out, int row, int col, double *mat);/向文件输出void outPut

12、(int row, int col, double *mat);/向控制台输出void inPut(int row, int col, double *mat);/从控制台输入int ij(int i, int j)/下三角矩阵元素判断return (i>=j)? i*(i+1)/2+j : j*(j+1)/2+i; ;文件mat.cpp：#include "mat.h"void mat:add(int row,int col,double m,double n,double result)for(int i=0; i<row*col; i+)resulti =

13、 mi+ni;void mat:sub(int row,int col,double m,double n,double result)for(int i=0; i<row*col; i+)resulti = mi-ni;double* mat:multiply(int row, int col, int col2, double m,double n)int i,j,k;double *result = new doublerow*col2;for(i=0; i<row; i+)for(j=0; j<col2; j+)resulti*col2+j = 0;for(k=0;

14、k<col; k+)resulti*col2+j += mi*col+k * nk*col2+j; return result;void mat:trans(int row, int col, double m, double result)int i,j;double *temp = new doublerow*col;for(i=0; i<col; i+)for(j=0; j<row; j+)resulti*row+j = mj*col+i;bool mat:inverse(double mat, int n)int i;int j;int k;double *c = n

15、ew doublen;for(k=0; k<n; k+)double c11 = mat0;if(c11 + 1.0 = 1.0)delete c;return false;for(i=1; i<n; i+)double ai0 = mati*(i+1)/2;if(i<=n-k-1)ci = -ai0/c11;elseci = ai0/c11;for(j=1; j<=i; j+)mat(i-1)*i/2 + j -1 = mati*(i+1)/2+j + ai0*cj;for(i=1; i<n; i+)mat(n-1)*n/2+i-1 = ci;matn*(n+1

16、)/2-1 = 1.0/c11;delete c;return true;void mat:fileIn(ifstream &in, int size, double *mat)for(int i=0; i<size; i+)in>>mati;void mat:(ofstream out, int row, int col, double *mat)for(int i=0; i<row; i+)for(int j=0; j<col; j+)out<<mati*col+j<<" "<<endl;out&

17、lt;<endl;void mat:outPut(int row, int col, double *mat)int i,j;for(i=0; i<row; i+)for(j=0; j<col; j+)cout<<mati*col+j<<" "cout<<endl;void mat:inPut(int row, int col, double *mat)int i,j;for(i=0; i<row; i+)for(j=0; j<col; j+)cin>>mati*col+j;测试文件main.cp

18、p:#include "mat.h"mat g_mat;void main()double *mat1;double *mat2;double *mat3;double *mat4;double *mat_result;mat1 = new double2*3;mat2 = new double2*3;mat3 = new double3*1;mat4 = new double2*2;cout<<"从文件中读取"<<endl;ifstream in("in.txt");g_mat.fileIn(in,2*3,m

19、at1);cout<<"矩阵mat1为："<<endl;g_mat.outPut(2,3,mat1);cout<<endl;g_mat.fileIn(in,2*3,mat2);cout<<"矩阵mat2为："<<endl;g_mat.outPut(2,3,mat1);cout<<endl;g_mat.fileIn(in,3*1,mat3);cout<<"矩阵mat3为："<<endl;g_mat.outPut(3,1,mat1);cout&

20、lt;<endl;g_mat.fileIn(in,2*2,mat4);cout<<"矩阵mat4为："<<endl;g_mat.outPut(2,2,mat1);cout<<endl;in.close();cout<<"mat1 + mat2 = "<<endl;mat_result = new double2*3;g_mat.add (2, 3, mat1, mat2, mat_result);g_mat.fileIn(in, 2*3, mat_result);g_mat.outPut(

21、2, 3, mat_result);cout<<"mat1 - mat2 = "<<endl;g_mat.sub (2, 3, mat1, mat2, mat_result);g_mat.fileIn(in, 2*3, mat_result);g_mat.outPut(2, 3, mat_result);cout<<"mat1 * mat3 = "<<endl;delete mat_result;mat_result = g_mat.multiply (2, 3, 1, mat1, mat3);g_mat.

22、fileIn(in, 2*3, mat_result);g_mat.outPut(2, 1, mat_result);cout<<"inv(mat4) = "<<endl;g_mat.inverse (mat4, 2);g_mat.outPut (2, 2, mat4);cout<<"转置 mat3 为"<<endl;delete mat_result;mat_result = new double1*3;g_mat.trans (3, 1, mat3, mat_result);g_mat.outPut(1

23、, 3, mat_result);delete mat1;delete mat2;delete mat3;delete mat_result;3. 运行结果测试文件“in.txt”： 1 2 34 5 62 3 45 6 71231 34 5 测试结果： 从文件中读取矩阵mat1为：1 2 34 5 6矩阵mat2为：1 2 34 5 6矩阵mat3为：123矩阵mat4为：1 23 4mat1 + mat2 =3 5 79 11 13mat1 - mat2 =-1 -1 -1-1 -1 -1mat1 * mat3 =1432inv(mat4) =-0.8 0.6-0.2 5转置 mat3 为

24、1 2 3Press any key to continue4. 出现的问题矩阵的输入与输出都要与一维数组相对应，在用动态指针操作内存容易出现数组越界的问题，直接造成程序崩溃，可以在关键的代码上使用trycatch排错。定义的类可以将函数体放在类中形成内联函数以增加程序的效率。函数除了使用类外还可以考虑使用一个命名空间对函数进行封装，或是在类中封装矩阵的其他元素，重载操作符元素（在简单的运算中比较方便，在较为复杂的运算中不太适用，比如用下三角存放的矩阵）。（4） 实习3（间接平差）1. 准备间接平差法（参数平差法）是通过选定t个与观测值有一定关系的独立未知量作为参数，将每个观测值都分别表达成这

25、t个参数的函数，建立函数模型，按最小二乘原理，用求自由极值的方法解出参数的最或然值，从而求得各观测值的平差值。设误差方程 V = Bx-l 其中V为观测值改正数，A为系数矩阵，x为参数向量，L为观测值向量，则x的最小二乘解 其中P为观测值权阵单位权中误差在间接平差的计算函数中，需要输入已知参数：n（观测值数）t（参数个数）B（误差方程系数矩阵，大小n*t）L（观测值向量，大小n*1）P（观测值权阵简化为对角数组）需要求解的参数：X（参数平差值向量，大小t*1）V（观测值平差值向量，大小n*1）u（单位权中误差）计算过程：1) 通过B、P、L计算、2) 求的逆矩阵3) 计算平差值 4) 将X代入

26、计算 V = Bx-l 5) 代入计算 计算函数void adjust(int n, int t, double V, double B, double X, double L, double P,double &u)2. 代码文件“mat.cpp” “mat.h”同上；文件“main.cpp”:#include"mat.h"#include<iomanip>mat mat_c;/*V = BX-L;V: n*1B: n*tX: t*1L: n*1P: n*1已知：B、L、P求：X、V、*/void adjust(int n, int t, double

27、 V, double B, double X, double L, double P,double &u)int i,j,k;double *BTPB = new double(t+1)*t/2;double *BTPB2 = new double(t+1)*t/2;for(i=0; i<(t+1)*t/2; i+)BTPB2i = 0;cout<<"BTPB2"<<endl;for(i=0; i<t; i+)for(j=0; j<=i; j+)for(k=0; k<n; k+)BTPB2i*(i+1)/2+j +=

28、Bk*t+i*Pk*Bk*t+j;cout<<BTPB2i*(i+1)/2+j;cout<<endl;mat_c.inverse(BTPB2,t);double *BTPL = new doublet*1;for(i=0; i<t; i+)BTPLi = 0;for(k=0; k<n; k+)BTPLi += Bk*t+i*Pk*Lk;for(i=0; i<t; i+)Xi = 0;for(j=0; j<t; j+)Xi += BTPB2mat_c.ij(i,j)*BTPLj;/X = mat_c.multiply(t,t,1,BTPB2,BTP

29、L); 下三角矩阵不能用一般乘法double *BX;BX = mat_c.multiply(n,t,1,B,X);mat_c.sub(n,1,BX,L,V);/ V = BX-L u = 0;for(i=0; i<n; i+)u += Vi*Pi*Vi;u = sqrt(u/(n-t);void main()double* B = new double3*6;double* L = new double6*1;double* P = new double6;double* X = new double3;double* V = new double6;double u = 0;ifst

30、ream file("");mat_c.fileIn(file,3*6,B);mat_c.fileIn();mat_c.fileIn();int n,t;file>>n;file>>t;cout<<setiosflags(ios:fixed) << setprecision(6) ;adjust(n,t,V,B,X,L,P,u);cout<<u<<endl;mat_c.outPut(6,1,V);double H = 0.023,1.114,1.142,0.079,0.099,1.21;for(int

31、i=0; i<6; i+)cout<<"改正后H"<<i<<":"<<Hi+Vi<<endl;Hi = Hi+Vi;cout<<"A: "<<10+H3<<endl;cout<<"B: "<<10+H4<<endl;cout<<"D: "<<10+H5<<endl;测试文件数据：-1 1 0 0 -1 1-1 0 11 0

32、00 1 00 0 10.0030.0030.0110002221116 33. 运行结果BTPB25.000000-2.0000005.000000-2.000000-2.0000005.0000000.0042430.0010000.001000-0.003000-0.004000-0.0000000.004000改正后H0:0.024000改正后H1:1.115000改正后H2:1.139000改正后H3:0.075000改正后H4:0.099000改正后H5:1.214000A: 10.075000B: 10.099000D: 11.214000Press any key to con

33、tinue4. 出现的问题由于double类型计算精度问题，输出时容易出现很趋近于0的小数，需要控制输出流格式，包含头文件<iomainip>，使用cout<<setiosflags(ios:fixed) << setprecision(6) ;最后的数字为控制小数输出的位数。为了提高计算效率，使用下三角矩阵存储，使用向量存储权阵，使用这种方法时，之前所用的矩阵类中乘法不能使用，需要直接用for循环进行计算。（5） 实习4（条件平差）1. 准备条件平差是根据各观测元素间所构成的几何条件以及起始数据间的强制条件，按最小二乘法的原理求得各观测值的最或然值，以消除

34、由于多次观测产生的矛盾的平差方法。条件平差模型为 其中V为观测值改正数，A为系数矩阵，W为条件方程自由项向量，则由最小二乘原理： 其中P为观测值权阵，K为拉格朗日算法中联系数向量，解为：改正数向量V的解为：单位权中误差在条件平差的计算函数中，需要输入已知参数：n（观测值数）r（条件方程个数）A（条件方程系数矩阵，大小r*n）W（条件方程自由项向量，大小r*1）P（观测值权阵简化为对角数组）需要求解的参数：V（观测值平差值向量，大小n*1）u（单位权中误差）计算过程：1) 通过A、P计算2) 求的逆矩阵3) 计算平联系数向量 4) 计算单位权中误差 2. 代码文件“mat.cpp” “mat.h

35、”同上；文件“main.cpp”:#include <iostream>#include <cmath>#include <iomanip>#include "mat.h"using namespace std;mat mat_cul;void condition(int n, int r, double A, double W, double P, double V,double &u)int i,j,k;double *AQAT = new doubler*(r+1)/2;/法方程系数阵for(i=0; i<r; i+)

36、for(j=0; j<=i; j+)AQATmat_cul.ij(i,j) = 0;for(k=0; k<n; k+)AQATmat_cul.ij(i,j) += Ai*n+k*Aj*n+k/Pk; /系数阵求逆if(mat_cul.inverse(AQAT,r) = false)return ;double *K = new double r;for(i=0; i<r; i+)Ki = 0.0;for(j=0; j<r; j+)Ki -= AQATmat_cul.ij(i,j)*Wj;/计算Vdouble pvv = 0;for(i=0; i<n; i+)dou

37、ble vi =0.0;for(int j=0; j<r; j+)vi += Aj*n+i*Kj;vi /= Pi;Vi = vi;pvv+= vi*vi*Pi;delete AQAT;delete K;u = sqrt(pvv/r);void main()double V6;double A = 1,0,0,1,-1,0, 0,0,-1,-1,0,1, 0,1,0,0,1,-1;double W = 0.003,-0.011,0.003;double P = 2,2,2,1,1,1;double H = 0.023,1.114,1.142,0.079,0.099,1.21;double

38、 u;condition(6,3,A,W,P,V,u);cout<<"中误差： "<<u<<endl;cout<<setiosflags(ios:fixed) << setprecision(4) ;int i;for(i=0; i<6; i+)cout<<Vi<<endl;for(i=0; i<6; i+)cout<<"改正后H"<<i<<":"<<Hi+Vi<<endl;Hi

39、= Hi+Vi;cout<<"A: "<<10+H3<<endl;cout<<"B: "<<10+H4<<endl;cout<<"D: "<<10+H5<<endl;3. 运行结果中误差： 0.004242640.00100.0010-0.0030-0.0040-0.00000.0040改正后H0:0.0240改正后H1:1.1150改正后H2:1.1390改正后H3:0.0750改正后H4:0.0990改正后H5:1.2140

40、A: 10.0750B: 10.0990D: 11.2140Press any key to continue4. 出现的问题和之前间接平差的问题类似，主要为输出流控制与矩阵计算方面的问题。（6） 实习5（水准网平差）1. 准备利用前面的间接平差函数，通过已知的数据求出平差计算所需要的n（观测值数）t（参数个数）B（误差方程系数矩阵，大小n*t）L（观测值向量，大小n*1）P（观测值权阵简化为对角数组）水准网平差用类Clevel进行封装，类中主要包含计算所需的矩阵和数据，以及进行计算的函数。通过类的封装在建立对象时输入水准网数据所在文件名，并计算后输出相应的结果。类中的函数成员如下：1) 构造

41、函数Clevel(char *) 输入文件的名称，并调用其它的计算函数，输出平差结果。2) void inputdata(char *);/从文件中输入读入计算所需的数据，存到相应的数组和变量中。首先从文件中读取观测值n、总点数pNum、已知点数p_k、精度m。从中得到t = pNum-p_k。然后从文件中读取起点号、终点号以及相应的高程值、距离，通过n判断循环次数。3) void calculataHBPL();/通过数据计算H、B、L、P。计算近似高程值。首先将H中已知点以后的高程数据赋初始值-9999.9，以判断是否已经得到H近似值（没赋近似值的高程值小于-9999）。然后判断数据中的水

42、准路线，若已知高程值差的两个高程点中一个为未知点，一个为已知点，则可得出未知点的近似高程值，直至所有高程值都已经求出时跳出循环。计算系数矩阵B。先将矩阵中所有数据赋值0。分别判断观测水准路线的高差点，如果第一个高程点不是已知点和控制点则将系数矩阵B中对应位置赋值-1，如果第二个高程点不是已知点和控制点则将系数矩阵B中对应位置赋值1。计算观测值向量L。观测值向量中每一个元素 Li=hi-(Hk2-Hk1); 式中k1、k2为观测值高程点对应高程向量H中的点号。计算权阵P。权阵中每个元素与距离成反比，通过Pi=m/Si;求出每个权阵。4) void adjust()；函数体内容与上面间接平差函数的

43、内容基本相同，通过上面求得的B、L、P求的平差后的X、V。由于函数里运算的变量命名与之前间接平差相同，只是把从参数输入数据改为直接调用类成员数据。5) void outputResult()；通过间接平差计算出的X、V，输出平差结果。改正后的观测值hi = Hi+Vi；改正后的高程值i = Hi+Xi；并将中误差，改正数等输出。2. 代码文件mat.cpp mat.h 同上。文件Clevel.h：#include <iostream>#include <fstream>#include <string>#include <iomanip>#inc

44、lude <math.h>#include "mat.h"using namespace std; class Clevelpublic:int n; /总观测值数int p_k; /已知点数int t; /未知点数int pNum; /高程点总数double m; /单位权string *Pname; /点名地址数组int *startP; /起点号int *endP; /终点号double *h; /观测值double *H; /高程值double *S; /距离长度double *P; /权double *B;double *L;double *X;dou

45、ble *V;double u; /中误差Clevel(char *);void inputdata(char *);int GetPNum(string name);void adjust();void outputResult();void calculataHBPL();mat mat_c;/矩阵计算类private:void calH();void calB();void calP();void calL();文件Clevel.cpp:#include "Clevel.h"Clevel:Clevel(char* )inputdata();calculataHBPL(

46、);adjust();outputResult();/*V = BX-L;V: n*1B: n*tX: t*1L: n*1P: n*1已知：B、L、P求：X、V、*/void Clevel:adjust()int i,j,k;double *BTPB = new double(t+1)*t/2;double *BTPB2 = new double(t+1)*t/2;/ 下三角求BTPB2 for(i=0; i<(t+1)*t/2; i+)BTPB2i = 0;for(i=0; i<t; i+)for(j=0; j<=i; j+)for(k=0; k<n; k+)BTPB

47、2i*(i+1)/2+j += Bk*t+i*Pk*Bk*t+j;/ BTPB2求逆mat_c.inverse(BTPB2,t);double *BTPL = new doublet*1;for(i=0; i<t; i+)BTPLi = 0;for(k=0; k<n; k+)BTPLi += Bk*t+i*Pk*Lk;for(i=0; i<t; i+)Xi = 0;for(j=0; j<t; j+)Xi += BTPB2mat_c.ij(i,j)*BTPLj;/X = mat_c.multiply(t,t,1,BTPB2,BTPL); 下三角矩阵不能用一般乘法doubl

48、e *BX;BX = mat_c.multiply(n,t,1,B,X);mat_c.sub(n,1,BX,L,V);/ V = BX-L u = 0;for(i=0; i<n; i+)u += Vi*Pi*Vi;u = sqrt(u/(n-t);int Clevel:GetPNum(string name)/点名数字化int i;for(i=0; i<pNum; i+)if(Pnamei != "zz")if(name = Pnamei)return i;else Pnamei = name;return i;return -1;void Clevel:inp

49、utdata(char *)/输入原始数据fstream in();in>>n>>pNum>>p_k>>m;t = pNum-p_k;startP = new intn;/起点号endP = new intn;/终点号L = new doublen;/观测值H = new doublepNum;/高程值S = new doublen;/距离长度P = new doublen*n;/权 B = new doublet*n;h = new doublen;Pname = new stringpNum;X = new doublet;V = new d

50、oublen;int i;for(i=0;i<pNum;i+)Pnamei="zz"/Getstationnumber函数根据Pname是否为空确定Pnamei是否为点名string buffer;/临时数组，保存从文件中读到的点名/读取已知高程数据for(i=0;i<p_k;i+)in>>buffer;int bb = GetPNum(buffer);in>>Hbb;/cout<<Heighty;/读取观测值数据for( i=0;i<n;i+)in>>buffer;/读取起点名startPi = GetPN

51、um(buffer);/获取起点对应的点号in>>buffer;/读取终点名endPi = GetPNum(buffer);in>>hi>>Si;/读取观测值高差及对应的距离in.close();/cout<<"1"/高程近似值计算void Clevel:calH()int i;for(i=p_k; i<pNum; i+)Hi=-9999.9;int Hcount=0;/计算出近似高程的点数while(1)for(i=0;i<n;i+)int k1=startPi;int k2=endPi;if(Hk1>-9999.0 && Hk2<-9999.0)Hk2=Hk1+hi;Hcount+;else if(Hk1<-9999.0 && Hk2>-9999.0)Hk1=Hk2-hi;Hcount+;if(Hcount = t)break;void Clevel:calB()int i,j;for(i=0; i<n*(pNum-p_k); i+)Bi = 0;for(i=0;i<n;i+)for(j=0;j