The following text is written by Walter Trump and is coming from his web site ,Thanks for his authorization to use his text.

How many magic squares are there?
Results of historical and computer enumeration

Order semi-magic
3 9   1   1   0   0  
4 68 688   880   48   48   0  
5 579 043 051 200   275 405 224   48 544   3 600   16  
6 9.4597 (13) ·1022   1.775399 (42) ·1019   0   0   0  
7 4.2848 (17) ·1038   3.79809 (50) ·1034   1.125151 (51) ·1018   1.21 (12) ·1017   20 190 684  
8 1.0804 (13) ·1059   5.2210 (70) ·1054   2.5228 (14) ·1027   C8 + ?   4.677 (17) ·1015  
9 2.8997 (69) ·1084   7.8448 (38) ·1079   7.28 (15) ·1040   81·E9 + ?   1.363 (21) ·1024  
10 1.477 (29) ·10115   2.4160 (35) ·10110   0   0   0  

Variations of a square by means of rotations and reflections are not counted.
Statistical notation: 1.2345 (25) ·109 means, that the (unknown) correct number is in the interval (1.2345 ±0.0025) ·109 with a probability of 99%.
Ultramagic squares are associative (centrally symmetrical) and pandiagonal.

B3 = 1, the so-called Lo Shu is unique.
B4 was found by the Frenchman Bernard Frénicle de Bessy in 1693. First analytical proof by Kathleen Ollerenshaw and Herman Bondi (1982).
A4, C5 and E5 could be found on the former website of Mutsumi Suzuki.
B5 was calculated in 1973 by Richard Schroeppel (computer-program), published in January, 1976 in Scientific American.
A5 was calculated by myself in march 2000 using a common PC. Suzuki published the result on his website. I could confirm the result by using other methods.
D5 is equal to the number of regular panmagic squares. They may be generated using latin squares, as Leonhard Euler pointed out already in the 18th century.
D6 and D10 were proved by A.H. Frost (1878) and more elegantly by C. Planck (1919).
C6 and C10 are also equal to 0, because each associative (symmetrical) magic square of even order can be transformed into a pandiagonal magic square.
B6 was first estimated (with proper result) by Karl Pinn and Christian Wieczerkowski and published in May, 1998. They used a method called 'Parallel Tempering Monte Carlo' and got: 1.7745(16)·1019
They also estimated B7 with lower accuracy: (3.760 ±0.052)·1034
My own estimates (see table) match these results. I used another method that could be called 'Monte Carlo Backtracking'.
All estimates in the columns B, C and D are found with the same method, that is more equal to the approach of Schroeppel than the one of Pinn and Wieczerkowski. There should be no severe systematical errors, because all results could be confirmed by another but more inaccurate method.
E7 was calculated in May 2001 by myself. Special transformations made it possible to consider only two positions of the integers 1, 25 and 49. With advanced equations and a heuristic backtracking algorithm the calculation time could be reduced to a few days. All ultramagic squares of order 7 have been saved and are available for further research. For more details see:Ultramagic Squares of Order 7
D7 (November 2001) was a big surprise. There are 38,102,400 regular pandiagonal magic squares of order 7. Albert L. Candy found 640,120,320 irregular ones. From E7 can be created nearly one billion more. But who would have believed that there are more than 1017 such squares? This estimate is very difficult, because the probability that a normal order-7 square is pandiagonal equals 1 : 3·1017.
D8 ist greater than C8 because each associative magic square of order 8 can be transformed into a pandiagonal one and there are examples of additional pandiagonal squares that could not be derived from an associative square.
E8 and E9 were estimated in March 2002. In the case of E8 I could find 64 transformations and several equations with only 4 variables.
D9 is greater than 81·E9, because each ultramagic square of order 9 can be transformed by cyclic permutation of rows and columns into 80 other pandiagonal magic squares that are not associative.

Walter Trump, Nürnberg, (c) 2001-11-01 (last modified: 2003-05-16)