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