幻方，是把自然数：1,2,3,.....,n*n，填入n*n的方阵里，使得方阵每行、每列、两条主对角线的和等于 1/2*n(n*n+1).  (n>=3)
n称为幻方的阶。
例如下面是一个3阶幻方：

4 9 2
3 5 7
8 1 6


以前我总想用pascal编一个幻方的程序，但对于偶数阶幻方中的单偶阶(4k+2)幻方，总是束手无策。查了好多网页，要么没有详细介绍，不知道采用哪个神仙招法，要么是个人网页，已经过期，被取消.
(现在的网站个人网页好多都要收费，没有交费就删删删删，查了好多都是无效页，
我%$#@@#%^&&***#@#!)

最后没办法，跑到台湾网站去找，直到看了一篇论文，才知道编出4k+2阶幻方的具体方法。总算完成了心愿！ (^_^)
下面是我编的程序，能打印出100阶以内的幻方。之所以只限于100阶是怕溢出和变量超过TP7所限制的64K.由于用到数组，再多阶的幻方，怕是要用到Delphi来编了！(算法可以一样)

========================================
program MagicSquare;
{幻方构造程序}

{uses crt; }
const  maxn=100;  {能打印出的幻方最大阶数 }
var
   a:array[1..maxn,1..maxn] of integer;
   n:integer;
   i,j:integer;
   ch:char;
procedure prn(n:integer);      {打印幻方到屏幕 }
var
   i,j:integer;
   sum:longint;
begin
   for i:=1 to n do
     begin
       for j:=1 to n do write(a[j,i]:4,' ');
        writeln;
     end;
   writeln;
   sum:=n;
   sum:=sum*(sum*sum+1) div 2;
   writeln('Magic Square n= ',n,'   Sum= ',sum);

end;

procedure prntofile(n:integer);   {打印幻方到文件}
var
   i,j,m:integer;
   sum:longint;
   f:text;
begin
    assign(f,'HFTEMP.txt');
    rewrite(f);

   for i:=1 to n do
     begin
       m:=0;
       for j:=1 to n do
         begin
          write(f,a[j,i]:6);
          m:=m+1;
          if m=40 then
             begin
                 m:=0;
                 writeln(f);
             end;
         end;
        writeln(f); writeln(f);
     end;
   sum:=n;
   sum:=sum*(sum*sum+1) div 2;
   writeln(f,'Magic Square n= ',n,'   Sum= ',sum);
   writeln(f,'Program by j.t.chang');
   close(f);
   writeln('Save to file: HFTEMP.TXT');
end;

procedure  oddhf(n:integer);        {奇数阶幻方的构造}
var
   x,k,p,xx,yy:integer;
begin
    for x:=1 to n*n do
      begin
         k:=(x-1) div n + (n+3) div 2 +x-1;
         yy:=k-(k-1) div n *n;
         p:=(n+1) div 2 +x -1 - (x-1) div n;
         xx:=n+1-p+ (p-1) div n *n;
         a[xx,yy]:=x;
      end;
end;

procedure DevenHF(n:integer);    {双偶阶(4k)幻方构造}
var
   i,j,k,L,m:integer;
begin
   for i:=1 to n do
     for j:=1 to n do
      begin
          if j mod 4>1 then m:=1
            else m:=0;
           k:=n-i-(n-2*i+1)*m;
          if i mod 4>1 then m:=1
            else m:=0;
          L:=n-j+1-(n-2*j+1)*m;
          a[i,j]:=k*N+L;
      end;
end;

procedure SevenHF(n:integer);    {单偶阶(4k+2)幻方构造 }
var
    i,j,k:integer;
begin
    DevenHF(n-2);               {采用先构造4k阶幻方}
    k:=(n-2) div 4;
    for i:=n-1 downto 2 do
         for j:=2 to n-1 do a[j,i]:=a[j-1,i-1]+8*k+2;  {扩成4k+2阶}
    for i:=1 to n-2 do
      begin
         a[i,1]:=0;
         a[1,i]:=0;
      end;
     {填外圈}
     a[1,1]:=1;  a[n,n]:=n*n;
     a[n,1]:=4;  a[1,n]:=n*n-3;
     a[n-1,1]:=10; a[2,n]:=3;
     a[3,n]:=5; a[4,n]:=7; a[n,n-2]:=2;  a[n,n-1]:=9;
     a[1,2]:=6; a[1,3]:=8;
     for i:=4 to k+2 do a[1,i]:=i+7;
     for j:=5 to k+3 do a[j,n]:=k+5+j;
     for i:=2*k+2 to 3*k do a[n,i]:=7+i;
     for j:=2*k+3 to 4*K do a[j,1]:=k+5+j;
     for i:=3*k+1 to 4*k-1 do a[n,i]:=2*k+5+i;
     for j:=k+4 to 2*k+2 do a[j,n]:=5*k+1+j;
     for i:=k+3 to 2*k+1 do a[1,i]:=6*k+1+i;

     for i:=2 to n-1 do
      begin
       if a[i,1]=0 then a[i,1]:=n*n+1-a[i,n]
         else a[i,n]:=n*n+1-a[i,1];
       if a[1,i]=0 then a[1,i]:=n*n+1-a[n,i]
         else a[n,i]:=n*n+1-a[1,i];
      end;
end;

(*************************)
begin

{     clrscr; }
     for i:=1 to  maxn do
       for j:=1 to maxn do a[i,j]:=0;

     write('Enter n:');
     readln(n);
     if( n<=2 ) or (n>maxn) then exit;
     if odd(n) then  oddHF(n)
      else if n mod 4=0 then  DevenHF(n)
        else  SevenHF(n);
     prn(n);
     write('Save to file  (y/n)? ');
     readln(ch);
     if (ch='y') or  (ch='Y') then  prntofile(n);
     writeln('Program by j.t.chang');
end.
========================================