The
following text is written by Christian Boyer and is coming from his web site
www.multimagie.com/indexengl.htm
,
Thanks for his authorization to use his text.
Multimagic formulas
How do we know what should be the theoretical values of sums of rows, columns and diagonals of the multimagic squares?
Start with the single magic square. It is known, and easy to demonstrate (…except for an 8year old child like Gauss who is supposed to have demonstrated at this age at school…), that the sum of the first integer numbers from 1 to N is:
N(N+1)/2
A nth order magic square contains all the numbers from 1 to n^{2}. So, the sum of the first integer numbers from 1 to n^{2} is:
n^{2}(n^{2}+1)/2
Since this square has n rows, it is sufficient to divide by n in order to know the magic sum S1 of each row, column or diagonal, so:
S1 = n(n²+1)/2
For the bimagic square, the sum of the squares of the first integers from 1 to N is:
N(N+1)(2N+1)/6
Replacing N by n^{2}, then dividing by n, we get the sum S2 of each row, column or diagonal, so:
S2 =
n(n²+1)(2n²+1)/6 = S1*(2n²+1)/3
For the trimagic square, the sum of the cubes:
S3 = n * S1²
For the tetramagic, Fermat gave as early as 1636, in a letter sent to Roberval, the solution for the sum of the integers from 1 to N to the 4thpower. Here is the Fermat's formula after having replaced N by n^{2}, then dividing by n:
S4 = ((4n^{2} + 2)*n*S1^{2}  S2) / 5
Replacing (4n^{2} +2) * S1 by 6S2, we can write an easier formula than Fermat:
S4 = S2 * (6n * S1  1) / 5
For the pentamagic (sum of 5thpower):
S5 = (3n*S2² 
S3) / 2
If we prefer to develop Sp, it gives the following generous table, which permits us to calculate the sums of the future hexa, hepta or octomagic squares:
S1 = (1/2) n^3 + (1/2) n
S2 = (1/3) n^5 + (1/2) n^3 + (1/6) n
S3 = (1/4) n^7 + (1/2) n^5 + (1/4) n^3
S4 = (1/5) n^9 + (1/2) n^7 + (1/3) n^5  (1/30) n
S5 = (1/6) n^11 + (1/2) n^9 + (5/12) n^7  (1/12) n^3
S6 = (1/7) n^13 + (1/2) n^11 + (1/2) n^9  (1/6) n^5+(1/42)n
S7 = (1/8)n^15+(1/2)n^13+(7/12)n^11(7/24)n^7+(1/12)n^3
S8 = (1/9)n^17+(1/2)n^15+(2/3)n^13(7/15)n^9+(2/9)n^5(1/30)n 
In order to calculate Sp when the nth order pmultimagic square contains the integers from 0 to n^{2}1 (instead of 1 to n^{2}), it is sufficient to put a  sign (instead of the + red ) in the 2nd column above. You may remark that the first column of Sp is equal to (1/(p+1)) n^(2p+1). This feature has been known since the middle of the XVIIth century, thanks to Pascal who has dedicated a Treatise on the subject. Jacques Bernoulli later took an interest in this question of the sum of integer powers. After his Ars conjectandi published after his death at the beginning of the XVIIIth century, the numbers used in this development are now called Bernoulli numbers (1/6, 1/30, 1/42, 1/30, 5/66, …).
Here are the formulas for some norders of pmultimagic squares, containing numbers from 1 to n^{2} . Here we will find the sums of the Pfeffermann's 8thorder bimagic square, Pfeffermann's 9thorder bimagic square, and William Benson's 32th order trimagic square:
Sp 
3rdorder 
4thorder 
5thorder 
6thorder 
7thorder 
8thorder 
9thorder 
... 
32thorder 
S1 
15 
34 
65 
111 
175 
260 
369 
... 
16400 
S2 
95 
374 
1105 
2701 
5775 
11180 
20049 
... 
11201200 
S3 
675 
4624 
21125 
73926 
214375 
540800 
1225449 
... 
8606720000 
And here are the formulas of some norders of pmultimagic squares, containing numbers from
0 to n²1. Here we will find the sums of our record 512thorder tetramagic and 1024thorder pentamagic squares:
Sp 
32ndorder 
512thorder 
1024thorder 
S1 
16368 
67108608 
536870400 
S2 
11168432 
11728056920832 
375299432076800 
S3 
8573165568 
2305825417061203968 
295147342229667840000 
S4 
7019705733392 
483565716171561366524160 
247587417561640996243120640 
S5 
5987221633671168 
105636341097042573844228866048 
216345083469423421673932062720000 
Bibliography
For more details about sum of powers, see the excellent chapter 14 Sommation des puissances numériques, book II of the Théorie des Nombres by Edouard Lucas , Librairie Blanchard, Paris.
