with(linalg): digit1:=proc(n);(n-1) mod 5;end: digit2:=proc(n);iquo(n-1,5);end: magic:=proc(i,j) 1+((digit1(i)+digit2(i)+2*digit1(j)+2*digit2(j)) mod 5)+ 5*((digit1(i)+2*digit2(i)+3*digit1(j)+digit2(j)) mod 5)+ 25*((digit1(i)+3*digit2(i)+2*digit1(j)+digit2(j)) mod 5)+ 125*((digit1(i)+4*digit2(i)+3*digit1(j)+2*digit2(j)) mod 5); end: m:=matrix(25,25,magic); # m is the MAGIC SQUARE # now we'll check some properties: one:=vector(25,1): print(`row sums:`); multiply(m,one); # output: columnvector of row sums print(`column sums:`); multiply(transpose(one),m); # output: row vector of column sums print(`diagonal sums:`); sum(m[i,i],i=1..25); sum(m[i,26-i],i=1..25); # output: diagonal sums mm:=matrix(25,25,(i,j)->m[i,j]^2): # mm is the matrix obtained by squaring each entry of m print(`after squaring each entry we get:`); print(`row sums:`); multiply(mm,one); # output: columnvector of row sums print(`column sums:`); multiply(transpose(one),mm); # output: row vector of column sums print(`diagonal sums:`); sum(mm[i,i],i=1..25); sum(mm[i,26-i],i=1..25); # output: diagonal sums