Cマガ電脳クラブ(第076回) 素数魔方陣

問題

今回は魔方陣の問題だ。3×3魔方陣とは、相違なる9つの数を3×3に配置し、
縦3本、横3本、対角2本、計8本の直線上にある3数の合計(これを「定和」という)がすべて同じになる、
数の配列の事である。

ここでは、この3×3魔方陣にちょっと趣向を凝らし、
各マスの数を素数で構成した「素数魔方陣」を作ってみた(Fig.1)。

この例では、図中の8本の破線上の数の合計はどこも177で、
この177を定和に持つ素数魔方陣は回転・対称の配置を除くとこれひとつしかない。
そこでひとつの定和に対してふたつの素数魔方陣が作れるもののなかで、定和が最小のものを見つけてほしい。

Fig.1 素数魔方陣

ソース

using System;
using System.Collections.Generic;
using System.Linq;

class Program
{
    const int TeiwaLimit = 400; //定和の上限
    static int[] SosuuArr;

    //エラトステネスの篩
    static void Eratosthenes()
    {
        var IsSosuuArr = new System.Collections.BitArray(TeiwaLimit + 1);
        for (int I = 2; I <= IsSosuuArr.Count - 1; I++) {
            IsSosuuArr[I] = true;
        }
        for (int I = 2; I * I <= IsSosuuArr.Count - 1; I++) {
            if (IsSosuuArr[I] == false) continue;
            for (int J = I * 2; J <= IsSosuuArr.Count - 1; J += I) {
                IsSosuuArr[J] = false;
            }
        }
        var SosuuList = new List<int>();
        for (int I = 2; I <= IsSosuuArr.Count - 1; I++) {
            if (IsSosuuArr[I]) SosuuList.Add(I);
        }
        SosuuArr = SosuuList.ToArray();
    }

    static void Main()
    {
        var sw = System.Diagnostics.Stopwatch.StartNew();

        Console.WriteLine("定和の上限を、{0}として検証します", TeiwaLimit);
        Eratosthenes();

        //2は使用できない
        SosuuArr = Array.FindAll(SosuuArr, X => X != 2);

        //深さ優先探索で魔方陣を列挙
        List<JyoutaiDef> LeafNodeList = ExecDFS();

        //回転した魔方陣を除外
        Console.WriteLine("回転した魔方陣の除外前の件数={0}", LeafNodeList.Count);
        RemoveKaitenKai(LeafNodeList);
        Console.WriteLine("回転した魔方陣の除外後の件数={0}", LeafNodeList.Count);
        Console.WriteLine();

        //ふたつの素数魔方陣が作れるもののなかで、最小の定和
        int AnswerTeiwa = LeafNodeList.GroupBy(A => A.Teiwa).Where(X => X.Count() > 1).Min(A => A.Key);

        Console.WriteLine("ふたつの素数魔方陣が作れるもののなかで、最小の定和={0}", AnswerTeiwa);
        Console.WriteLine("魔方陣は");
        foreach (JyoutaiDef EachLeafNode in LeafNodeList.Where(A => A.Teiwa == AnswerTeiwa)) {
            PrintBan(EachLeafNode.BanArr);
        }
        Console.WriteLine("経過時間={0}", sw.Elapsed);
    }

    struct JyoutaiDef
    {
        internal int CurrX;
        internal int CurrY;
        internal int[,] BanArr;
        internal int Teiwa;
    }

    static int UB = 3 - 1;

    //深さ優先探索で魔方陣を列挙
    static List<JyoutaiDef> ExecDFS()
    {
        var WillReturn = new List<JyoutaiDef>();

        var stk = new Stack<JyoutaiDef>();
        JyoutaiDef WillPush;
        WillPush.CurrX = WillPush.CurrY = 0;
        WillPush.BanArr = new int[UB + 1, UB + 1];
        WillPush.Teiwa = 0;
        stk.Push(WillPush);

        while (stk.Count > 0) {
            JyoutaiDef Popped = stk.Pop();

            //X座標の繰上げ処理
            if (Popped.CurrX > UB) {
                Popped.CurrX = 0;
                Popped.CurrY++;
            }

            //最終行を超えた場合
            if (Popped.CurrY > UB) {
                WillReturn.Add(Popped);
                Console.WriteLine("定和={0}の葉ノードを発見。現在の葉ノード数={1}",
                    Popped.Teiwa, WillReturn.Count);
                PrintBan(Popped.BanArr);
                continue;
            }

            //Push処理
            foreach (int EachInt in SosuuArr.Except(Popped.BanArr.Cast<int>())) {
                //対象解の除外で、左上の数値 < 右上の数値とする
                if (Popped.CurrX == 2 && Popped.CurrY == 0) {
                    if (Popped.BanArr[0, 0] > EachInt) continue;
                }

                WillPush.CurrX = Popped.CurrX + 1;
                WillPush.CurrY = Popped.CurrY;
                WillPush.BanArr = (int[,])Popped.BanArr.Clone();
                WillPush.BanArr[Popped.CurrX, Popped.CurrY] = EachInt;

                if (Popped.CurrX == 2 && Popped.CurrY == 0) {
                    WillPush.Teiwa = DeriveSum(WillPush.BanArr, 0, 0, 1, 0, 2, 0);
                    //定和の上限で枝切り
                    if (WillPush.Teiwa > TeiwaLimit) break;
                    //定和が3の倍数でなかったら枝切り
                    if (WillPush.Teiwa % 3 > 0) continue;
                }
                else WillPush.Teiwa = Popped.Teiwa;

                if (WillEdakiri(WillPush.BanArr, Popped.Teiwa, Popped.CurrX, Popped.CurrY))
                    continue;

                stk.Push(WillPush);
            }
        }
        return WillReturn;
    }

    //座標3つを引数として、和を求める
    static int DeriveSum(int[,] pBanArr, int pX1, int pY1, int pX2, int pY2, int pX3, int pY3)
    {
        return pBanArr[pX1, pY1] + pBanArr[pX2, pY2] + pBanArr[pX3, pY3];
    }

    //枝切りするかを判定
    static bool WillEdakiri(int[,] pBanArr, int pTeiwa, int pX, int pY)
    {
        if (pX == 0 && pY == 0) {
            if (pBanArr[0, 0] * 2 > TeiwaLimit) return true;
        }
        if (pX == 1 && pY == 0) {
            if (pBanArr[0, 0] * 2 + pBanArr[1, 0] > TeiwaLimit) return true;
        }
        if (pX == 0 && pY == 1) {
            int wkSum = pBanArr[0, 0] + pBanArr[0, 1];
            if (wkSum != pTeiwa / 3 + pBanArr[2, 0]) return true;
            if (Array.BinarySearch(SosuuArr, pTeiwa - wkSum) < 0) return true;
        }
        if (pX == 1 && pY == 1) {
            int wkSum1 = pBanArr[1, 0] + pBanArr[1, 1];
            if (wkSum1 > pTeiwa) return true;
            if (Array.BinarySearch(SosuuArr, pTeiwa - wkSum1) < 0) return true;

            int wkSum2 = pBanArr[0, 0] + pBanArr[0, 1];
            int wkSum3 = pBanArr[2, 0] + pBanArr[1, 1];
            if (wkSum2 != wkSum3) return true;
        }
        if (pX == 2 && pY == 1) {
            int wkSum1 = pBanArr[0, 0] + pBanArr[1, 1];
            int wkSum2 = pBanArr[2, 0] + pBanArr[2, 1];
            if (wkSum1 != wkSum2) return true;
            if (Array.BinarySearch(SosuuArr, pTeiwa - wkSum1) < 0) return true;

            int wkTeiwa = DeriveSum(pBanArr, 0, 1, 1, 1, 2, 1);
            if (pTeiwa != wkTeiwa) return true;
        }
        if (pX == 0 && pY == 2) {
            int wkTeiwa = DeriveSum(pBanArr, 0, 0, 0, 1, 0, 2);
            if (pTeiwa != wkTeiwa) return true;
        }
        if (pX == 1 && pY == 2) {
            int wkTeiwa = DeriveSum(pBanArr, 1, 0, 1, 1, 1, 2);
            if (pTeiwa != wkTeiwa) return true;
        }
        if (pX == 2 && pY == 2) {
            int wkTeiwa = DeriveSum(pBanArr, 2, 0, 2, 1, 2, 2);
            if (pTeiwa != wkTeiwa) return true;
        }
        return false;
    }

    //回転した解を除外
    static void RemoveKaitenKai(List<JyoutaiDef> pLeafNodeList)
    {
        Predicate<int> IsExist = (pCurrInd) =>
        {
            for (int I = 0; I <= pCurrInd - 1; I++) {
                //定和が違ってたら、違う魔方陣
                if (pLeafNodeList[pCurrInd].Teiwa != pLeafNodeList[I].Teiwa)
                    continue;

                bool IsOK1 = false, IsOK2 = false, IsOK3 = false, IsOK4 = false;
                bool IsOK5 = false, IsOK6 = false, IsOK7 = false; //回転1から7

                for (int X = 0; X <= UB; X++) {
                    for (int Y = 0; Y <= UB; Y++) {
                        int CurrVal = pLeafNodeList[pCurrInd].BanArr[X, Y];
                        int[,] wkP = pLeafNodeList[I].BanArr;
                        if (wkP[UB - X, Y] != CurrVal) IsOK1 = true;
                        if (wkP[UB - X, UB - Y] != CurrVal) IsOK2 = true;
                        if (wkP[X, UB - Y] != CurrVal) IsOK3 = true;
                        if (wkP[Y, X] != CurrVal) IsOK4 = true;
                        if (wkP[UB - Y, X] != CurrVal) IsOK5 = true;
                        if (wkP[UB - Y, UB - X] != CurrVal) IsOK6 = true;
                        if (wkP[Y, UB - X] != CurrVal) IsOK7 = true;
                    }
                }
                if (IsOK1 == false || IsOK2 == false || IsOK3 == false || IsOK4 == false
                 || IsOK5 == false || IsOK6 == false || IsOK7 == false)
                    return true;
            }
            return false;
        };

        for (int I = pLeafNodeList.Count - 1; I >= 0; I--) {
            if (IsExist(I)) pLeafNodeList.RemoveAt(I);
        }
    }

    //盤面を表示
    static void PrintBan(int[,] pBanArr)
    {
        var sb = new System.Text.StringBuilder();
        for (int Y = 0; Y <= UB; Y++) {
            for (int X = 0; X <= UB; X++) {
                sb.AppendFormat("{0,3},", pBanArr[X, Y]);
            }
            sb.AppendLine();
        }
        Console.WriteLine(sb.ToString());
    }
}

実行結果

省略
定和=177の葉ノードを発見。現在の葉ノード数=35
 17,113, 47,
 89, 59, 29,
 71,  5,101,

定和=177の葉ノードを発見。現在の葉ノード数=36
 17, 89, 71,
113, 59,  5,
 47, 29,101,

回転した魔方陣の除外前の件数=36
回転した魔方陣の除外後の件数=9

ふたつの素数魔方陣が作れるもののなかで、最小の定和=381
魔方陣は
157, 43,181,
151,127,103,
 73,211, 97,

157, 13,211,
181,127, 73,
 43,241, 97,

経過時間=00:00:03.7402711

解説

定和の上限を決めて、深さ優先探索で魔方陣を列挙し、
定和の上限未満の定和が解として存在するかを検証してます。

魔方陣を
ABC
DEF
GHI
とすると
 A+B+C +D+E+F +G+H+I +A+D+G +B+E+H +C+F+I +A+E+I +C+E+G = 定和*8
Eを含まない加算を定和に変形して、
 定和  +D+E+F +定和  +定和  +B+E+H +定和  +A+E+I +C+E+G = 定和*8
移行して、
D+E+F +B+E+H +A+E+I + C+E+G = 定和*4
B+E+Hも定和として、
D+E+F +A+E+I + C+E+G = 定和*3
ADGに着目して
  E+F   + E+I + C+E  = 定和*2
CFIに着目して
E+E+E = 定和
よって、定和は、3の倍数と分かりますので、枝切り条件に使ってます。

また、魔方陣に2を含むと、そのラインを含む定和が偶数となり、
2を含まないラインは、必ず奇数ですので、
2を候補から除外してます。