Problem329 素数ガエル

問題

Susanは素数ガエルを飼っている.
このカエルは, 1から500までの番号がつけられた500個の正方形の上を跳びまわる.
左右のいずれかの正方形に等確率で跳ぶことができ,
範囲 [1;500] の外に跳ぶことはできない.
(いずれかの端に着地したなら, 次は自動的に移動可能な正方形に跳ぶ.)

素数である数が書かれた正方形にいる場合, 次の正方形に跳ぶ直前に,
カエルは 2/3 の確率で 'P' (PRIME) と鳴き,
1/3 の確率で 'N' (NOT PRIME) と鳴く.

素数でない数が書かれた正方形にいる場合, 次の正方形に跳ぶ直前に,
カエルは 1/3 の確率で 'P' と鳴き,
2/3 の確率で 'N' と鳴く.

カエルの開始位置は全ての正方形に対し等確率でランダムであり,
また最初の15回の鳴き声を聞いたとすると,
PPPPNNPPPNPPNPNと聞こえる確率は何か.

約分済みの分数 p/q の形で答えを入力せよ.

ソース

using System;
using System.Collections.Generic;

class Program
{
    const int UB = 500;

    static int[] mSosuuArr;

    //エラトステネスの篩
    static void Eratosthenes()
    {
        var IsSosuuArr = new System.Collections.BitArray(UB + 1);

        for (int I = 2; I <= IsSosuuArr.Count - 1; I++) {
            IsSosuuArr[I] = true;
        }
        for (int I = 2; I * I <= IsSosuuArr.Count - 1; I++) {
            if (I != 2 && I % 2 == 0) continue;

            if (IsSosuuArr[I]) {
                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);

        mSosuuArr = SosuuList.ToArray();
    }

    //仮分数クラス
    internal class CKabunSuu
    {
        private decimal mBunshi; //分子
        private decimal mBunbo; //分母

        //コンストラクタ (値指定)
        internal CKabunSuu(decimal pBunshi, decimal pBunbo)
        {
            mBunshi = pBunshi;
            mBunbo = pBunbo;
        }

        //コンストラクタ (インスタンス指定)
        internal CKabunSuu(CKabunSuu pIns)
        {
            mBunshi = pIns.mBunshi;
            mBunbo = pIns.mBunbo;
        }

        //加算用メソッド (インスタンス指定)
        internal void Add(CKabunSuu pIns)
        {
            decimal SavedBunbo = mBunbo; //通分する前の分母

            //通分してから、加算
            mBunbo *= pIns.mBunbo;
            mBunshi = this.mBunshi * pIns.mBunbo + pIns.mBunshi * SavedBunbo;
            Yakubun();
        }

        //乗算用メソッド (インスタンス指定)
        internal void Prod(CKabunSuu pIns)
        {
            this.mBunshi *= pIns.mBunshi;
            this.mBunbo *= pIns.mBunbo;
            Yakubun();
        }

        //除算用メソッド (値指定)
        internal void Div(decimal pVal)
        {
            mBunbo *= pVal;
            Yakubun();
        }

        //約分
        private void Yakubun()
        {
            if (mBunshi == 0) return;

            decimal GCD; //最大公約数
            GCD = DeriveGCD(mBunshi, mBunbo);
            if (GCD == 1) return;
            mBunshi /= GCD;
            mBunbo /= GCD;
        }

        //ユークリッドの互除法で2数の最大公約数を求める
        private decimal DeriveGCD(decimal pVal1, decimal pVal2)
        {
            decimal WarareruKazu = Math.Max(pVal1, pVal2);
            decimal WaruKazu = Math.Min(pVal1, pVal2);

            while (true) {
                decimal Amari = WarareruKazu % WaruKazu;
                if (Amari == 0) return WaruKazu;
                WarareruKazu = WaruKazu;
                WaruKazu = Amari;
            }
        }

        internal bool IsZero()
        {
            return mBunshi == 0;
        }

        //分子と分母を返す
        internal string ToStr()
        {
            return string.Format("{0}/{1}", mBunshi, mBunbo);
        }
    }

    static void Main()
    {
        Eratosthenes();

        //正しい鳴き順である確率[現在位置]なDP表
        var PrevDP = new CKabunSuu[UB + 1];

        for (int I = 1; I <= UB; I++) {
            PrevDP[I] = new CKabunSuu(1, UB);
        }

        for (int I = 1; I <= 15; I++) {
            var CurrDP = new CKabunSuu[UB + 1];
            for (int J = 1; J <= UB; J++) {
                CurrDP[J] = new CKabunSuu(0, 1);
            }

            for (int J = 1; J <= UB; J++) {
                if (PrevDP[J].IsZero()) continue;

                //移動可能な方向の数
                int CanMoveHoukou = 2;
                if (J == 1 || J == UB) CanMoveHoukou = 1;

                Action<int> MergeAct = (pNewJ) =>
                {
                    if (pNewJ < 1 || UB < pNewJ) return;

                    var wkIns = new CKabunSuu(PrevDP[J]);
                    var CorrectKakuritu = DeriveCorrectKakuritu(I, J);
                    wkIns.Prod(CorrectKakuritu);

                    wkIns.Div(CanMoveHoukou);
                    CurrDP[pNewJ].Add(wkIns);
                };

                MergeAct(J - 1);
                MergeAct(J + 1);
            }
            PrevDP = CurrDP;
        }

        var AnswerIns = new CKabunSuu(0, 1);
        for (int I = 1; I <= UB; I++) {
            AnswerIns.Add(PrevDP[I]);
        }
        Console.WriteLine(AnswerIns.ToStr());
    }

    //ジャンプ回数と、現在位置を引数として、正しく鳴く確率を返す
    static CKabunSuu DeriveCorrectKakuritu(int pJumpCnt, int pCurrPos)
    {
        const string CorrectStr = "PPPPNNPPPNPPNPN";
        char CurrChar = CorrectStr[pJumpCnt - 1];

        bool IsPrime = (Array.BinarySearch(mSosuuArr, pCurrPos) >= 0);

        bool IsMatch = false;
        if (CurrChar == 'P' && IsPrime) IsMatch = true;
        if (CurrChar == 'N' && IsPrime == false) IsMatch = true;

        return new CKabunSuu(IsMatch ? 2 : 1, 3);
    }
}

実行結果

199740353/29386561536000

解説

正しい鳴き順である確率[現在位置]を、
確率の加法定理と、確率の乗法定理で更新していく、
動的計画法で解いてます。