/************************************************************************
 * ms5x5.c: ５×５魔方陣の総数を求める．
 ------------------------------------------------------------------------
 * 奥村晴彦「アルゴリズム辞典」(日本評論社) に載っている全ての４方陣を
 * 求めるプログラムをそのまま５方陣に書き換えたものです．
 * 10000個毎に解を表示します．
 *
 * Num[m33 = i] == Num[m33 = 26-i] だから，
 * m33 = 1, ..., 13 に対して求まれば，総数がわかる．
 *
 * 所要時間の目安： 約２時間３０分 / Pentium-M 2GHz
 -----------------------------------------------------------------------*/
#include <stdio.h>

#define B(m)        if (ok[m]) { ok[m] = 0;
#define E(m)        ok[m] = 1; }
#define foreach(m)  for (m = 1; m <= 25; m++)
#define FORMAT     " %2d %2d %2d %2d %2d\n"

int main()
{
    int m11, m12, m13, m14, m15,
        m21, m22, m23, m24, m25,
        m31, m32, m33, m34, m35,
        m41, m42, m43, m44, m45,
        m51, m52, m53, m54, m55;
    int i, count = 0;
    int Num[14] = {0,};
    int ok34[90], *ok=ok34+34; /* ok[-34..55] */

    for (i =-34; i <=  0; i++) ok[i] = 0;
    for (i =  1; i <= 25; i++) ok[i] = 1;
    for (i = 26; i <= 55; i++) ok[i] = 0;

    for (m33 = 1; m33 <= 13; m33++) B(m33)
    foreach(m11) B(m11)
    for (m55 = m11+1; m55 <= 25; m55++) B(m55)
    for (m15 = m11+1; m15 <= 25; m15++) B(m15)
    for (m51 = m15+1; m51 <= 25; m51++) B(m51)
        foreach(m22) B(m22) m44 = 65-m11-m22-m33-m55; B(m44)
        foreach(m24) B(m24) m42 = 65-m15-m24-m33-m51; B(m42)
        foreach(m12) B(m12)
        foreach(m13) B(m13) m14 = 65-m11-m12-m13-m15; B(m14) 
        foreach(m52) B(m52) m32 = 65-m12-m22-m42-m52; B(m32)
        foreach(m53) B(m53) m54 = 65-m51-m52-m53-m55; B(m54)
                            m34 = 65-m14-m24-m44-m54; B(m34)
        foreach(m23) B(m23) m43 = 65-m13-m23-m33-m53; B(m43)
        foreach(m21) B(m21) m25 = 65-m21-m22-m23-m24; B(m25)
        foreach(m31) B(m31) m35 = 65-m31-m32-m33-m34; B(m35)
                            m41 = 65-m11-m21-m31-m51; B(m41)
                            m45 = 65-m41-m42-m43-m44;
                            if (ok[m45]) {
                                if ((++count % 10000) == 0) {
                                    printf("No.%9d\n", count);
                                    printf(FORMAT,m11,m12,m13,m14,m15);
                                    printf(FORMAT,m21,m22,m23,m24,m25);
                                    printf(FORMAT,m31,m32,m33,m34,m35);
                                    printf(FORMAT,m41,m42,m43,m44,m45);
                                    printf(FORMAT,m51,m52,m53,m54,m55);
                                    printf("\n");
                                }
                            }
        E(m41) E(m35) E(m31)
        E(m25) E(m21)
        E(m43) E(m23)
        E(m34) E(m54) E(m53)
        E(m32) E(m52)
        E(m14) E(m13)
        E(m12)
        E(m42) E(m24)
        E(m44) E(m22)
    E(m51) E(m15) E(m55) 
    printf("Sub-Total to [m33=%2d, m11=%2d] = %9d\n\n", m33, m11, count);
    E(m11) 
    Num[m33]=count-Num[0];
    Num[0]=count;
    printf("Total of Magic5[m33=%2d] = %9d\n\n", m33, Num[m33]);
    E(m33)
    printf("Done!  Total = %9d\n\n", 2*Num[0]-Num[13]);
    getchar();
    return 0;
}