Cマガ電脳クラブ(第053回) 平方数の逆襲

問題

異なる平方数の逆数の和で、2分の1を作る。
平方数とは自然数を2乗したもので、1、4、9、16・・・のことである。式で表すと

   n
 シグマ( 1/ (a(I)の2乗) ) = 0.5
  I=1

となる。ただし、a(I)は自然数で、
1 <= a(1) < a(2) < ・・・ a(n)
である。

これだけだと何通りもの表し方があるが、このうちでa(n)の値が最小となるものを求めてほしい。
ヒント：加える個数は11個以下。余裕のある人はヒントを利用せずに求めてください。

ソース

using System;
using System.Collections.Generic;

class Program
{
    const int NJyougen = 35; //Nの上限

    //解の下限と上限
    const decimal AnswerKagen = 0.4999999999999999999M;
    const decimal AnswerJyougen = 0.5000000000000000001M;

    struct JyoutaiDef
    {
        internal int CurrP;
        internal decimal SumVal;
        internal List<int> BunboList;
    }

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

        Console.WriteLine("{0}の2乗をa(n)の上限として検証します", NJyougen);

        int[] HeihouSuuArr = new int[NJyougen + 1];
        for (int I = 2; I <= NJyougen; I++)
            HeihouSuuArr[I] = I * I;

        decimal[] GyakuSuuArr = new decimal[NJyougen + 1];
        for (int I = 2; I <= NJyougen; I++)
            GyakuSuuArr[I] = 1M / (decimal)HeihouSuuArr[I];

        DeriveRunSum(GyakuSuuArr);

        var Stk = new Stack<JyoutaiDef>();
        JyoutaiDef WillPush;

        for (int I = 2; I <= NJyougen; I++) {
            //添字以降の逆数の総和が、解の下限未満ならBreak
            if (mRunSumArr[I] < AnswerKagen) break;

            WillPush.CurrP = I;
            WillPush.SumVal = GyakuSuuArr[I];

            WillPush.BunboList = new List<int>() { HeihouSuuArr[I] };
            Stk.Push(WillPush);
        }

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

            //解の下限以上、かつ、解の上限以下であることが必要条件
            if (AnswerKagen <= Popped.SumVal && Popped.SumVal <= AnswerJyougen) {
                //真の値が0.5であることが十分条件
                if (IsAnswer(Popped.BunboList)) {
                    Console.WriteLine("解候補を発見");
                    PrintAnswer(Popped.BunboList);
                }
            }

            //添字以降の逆数の総和を足しても、解の下限未満なら枝切り
            int ValidMinInd = Popped.CurrP + 1;
            for (int I = Popped.CurrP + 1; I <= NJyougen; I++) {
                if (Popped.SumVal + mRunSumArr[I] < AnswerKagen) break;
                ValidMinInd = I;
            }

            for (int I = NJyougen; I >= Popped.CurrP + 1; I--) {
                if (ValidMinInd < I) continue;

                WillPush.CurrP = I;
                WillPush.SumVal = Popped.SumVal + GyakuSuuArr[I];
                WillPush.BunboList = new List<int>(Popped.BunboList);
                WillPush.BunboList.Add(HeihouSuuArr[I]);

                //解の上限を超えたら枝切り
                if (WillPush.SumVal > AnswerJyougen) break;

                Stk.Push(WillPush);
            }
        }
        Console.WriteLine("経過時間={0}", sw.Elapsed);
    }

    //単位分数の分子のListを引数として、真の値が0.5かを判定
    static bool IsAnswer(List<int> pBunboList)
    {
        long CurrBunsi = 0;
        long CurrBunbo = 1;

        foreach (int EachTaniBunsuuBunbo in pBunboList) {
            CurrBunsi = CurrBunsi * EachTaniBunsuuBunbo + CurrBunbo;
            CurrBunbo *= EachTaniBunsuuBunbo;
        }
        if (CurrBunbo % CurrBunsi > 0) return false;
        return CurrBunbo / CurrBunsi == 2;
    }

    //添字以降の逆数の総和を求める
    static decimal[] mRunSumArr = new decimal[NJyougen + 1];
    static void DeriveRunSum(decimal[] pGyakuSuuArr)
    {
        decimal RunSum = 0M;
        for (int I = NJyougen; 2 <= I; I--) {
            RunSum += pGyakuSuuArr[I];
            mRunSumArr[I] = RunSum;
        }
    }

    //平方根を求める
    static int GetHeihouKon(int pVal)
    {
        return (int)Math.Sqrt(pVal);
    }

    //解を出力
    static void PrintAnswer(List<int> pIntList)
    {
        var sb = new System.Text.StringBuilder();
        sb.Append("1/2 = ");
        for (int I = 0; I <= pIntList.Count - 1; I++) {
            sb.AppendFormat("1/({0}の2乗)", GetHeihouKon(pIntList[I]));
            if (I < pIntList.Count - 1) {
                if (I % 5 == 4) {
                    sb.AppendLine();
                    sb.Append("   ");
                }
                sb.Append(" + ");
            }
        }
        Console.WriteLine(sb.ToString());
    }
}

実行結果

35の2乗をa(n)の上限として検証します
解候補を発見
1/2 = 1/(2の2乗) + 1/(3の2乗) + 1/(4の2乗) + 1/(5の2乗) + 1/(7の2乗)
    + 1/(12の2乗) + 1/(15の2乗) + 1/(20の2乗) + 1/(28の2乗) + 1/(35の2乗)
経過時間=00:00:53.1619881

解説

a(n)の上限を決めてから深さ優先探索してます。

Decimal型の有効桁は、28桁ですので
Decimal型の単位分数40個の加算で、真の値が0.5ならば、
0.5- (10の-28乗)*40 <= 計算結果 >= 0.5 + (10の-28乗)*40
を満たします。(加算において、累計誤差は加算されるため)

上記の不等式を満たせば、
0.4999999999999999999 <= 計算結果 >= 0.5000000000000000001
も満たすので、

真の値が0.5ならば、
0.4999999999999999999 <= 計算結果 >= 0.5000000000000000001
を導出でき、この命題の対偶は、真の値が0.5である必要条件となります。

そして、必要条件を満たした分母の組み合わせの、
真の値が0.5かを調べて、十分性を検証してます。