砂場で遊ぼう

c++/python/mathematica などの練習帳

Project Euler 209 / 循環論法

問題概略

 k 入力の 2 進真理値表は  k 個の入力ビット(2 進数で 0 は「偽」,1 は「真」)から 1 個の出力ビットへの写像である。
たとえば論理和(AND)と排他的論理和(XOR)の 2 入力真理値表は次のようになる。

f:id:variee:20211219123938p:plain

6 ビットの入力( a,  b,  c,  d,  e,  f)に対して以下の式をみたす 6 入力の 2 進真理値表  \tau はいくつあるか。

 \tau(a,\, b,\, c,\, d,\, e,\, f)\, \mathtt{AND}\, \tau(b,\, c,\, d,\, e,\, f,\, a\,\mathtt{XOR}\, (b\,\mathtt{AND}\, c)) = 0

https://projecteuler.net/problem=209

解説の pdf も作りました。きれいなレイアウトで読みたい方はこちらをどうぞ。

drive.google.com

何を求めればいい?

左右の  \tau の引数を  x,  p(x) で表して AND と XOR をそれぞれ  \land,  \oplus で表すと,与式は次のようになります。

\begin{align*}
&\tau(x)\land \tau\left(p(x)\right)=0\\
&p(a,\, b,\, c,\, d,\, e,\, f)=(b,\, c,\, d,\, e,\, f,\, a\oplus (b\land c))
\end{align*}

以下,真偽を  \bigcirc\times で表します。

与式が成立するのは  x p(x) が「両方とも  \times」または「片方が  \bigcirc,もう片方が  \times」のときです。
この条件をみたすように  p(p(x)),  p(p(p(x)))……を決めるときのルールは次のようになります。

  •  x p(x) が両方とも  \times」のとき  p(p(x)) \bigcirc でも  \times でもいい
  •  x \times p(x) \bigcirc」のとき  p(p(x)) \times
  •  x \bigcirc p(x) \times」のとき  p(p(x)) \bigcirc でも  \times でもいい
  •  x p(x) が両方とも  \bigcirc」になる  x は存在しない

入試数学でよくある「 \times は連続してもいいが, \bigcirc は連続しない順列」の問題を連想しますね。

 x などがとる値は全部で  2^6=64 個なので, x\to p(x)\to p(p(x))\cdots の連鎖は循環しているはずです。この循環を円で表すと,本問は

「円周上に  \bigcirc が連続しないようにうまく  \bigcirc\times を置く方法が何通りあるか」

を求める問題です。

mathematica で実験

64 個の  x すべてが 1 つの円周上に並ぶとは限らないので,まずこれを mathematica で調べました。

 x p(x) を dsu の同じ集合に入れて,部分集合の要素数を数えました。
小さい方から順に 1, 2, 3, 6, 6, 46 個でした。

    In[]:= AbsoluteTiming[
        p[{a_, b_, c_, d_, e_, f_}] := {b, c, d, e, f, 
          BitXor[a, BitAnd[b, c]]};
        ds = CreateDataStructure["DisjointSet"];
        For[i = 0, i < 2^6, i++,
         x = IntegerDigits[i, 2, 6]; y = p@x;
         ds["Insert", x]; ds["Insert", y];
         ds["Unify", x, y]];
        Sort[Length /@ ds["Subsets"]]]
       
    Out[]= {0.0005681, {1, 2, 3, 6, 6, 46}}

 x p(x) の間に有向辺をはってグラフにすると視覚的にわかりやすいです。

    p[{a_, b_, c_, d_, e_, f_}] := {b, c, d, e, f, BitXor[a, BitAnd[b, c]]};
    lst = Tuples[{0, 1}, 6];
    f[lst_] := lst -> p@lst;
    Graph[f /@ lst]

f:id:variee:20211219124158p:plain

実はリュカ数列

これまでの考察をもとに,条件をみたす真理値表の埋め方が何通りあるか計算します。

要素数 1 の集合に属する  x p(x) で自分自身にうつる元なので  \times の 1 通りです。

要素数 2 の集合は「両方とも  \times」か「片方が  \bigcirc,もう片方が  \times」なので  1+2=3 通り。

要素数 3 の集合に含まれる  \bigcirc は高々 1 個。0 個のものが 1 通りで,1 個のものが 3 通りの計 4 通り。

要素数 6 の集合も同様の casework で求められますが,要素数 46 個の集合は無理なので一般化して要素数  n 個の集合を考えます。
一番若い数字の元が  \bigcirc \times かで場合分けします。

  •  \bigcirc のとき,その両隣は  \times で残り  n-3 個の元の  \bigcirc\times は任意
  •  \times のとき,残り  n-1 個の元の  \bigcirc\times は任意

 \times は連続してもいいが, \bigcirc は連続しない」をみたすように円周上に  n 個の  \bigcirc\times を並べる方法が  a_n 通りあり,直線上に並べる方法が  b_n 通りあるとします。上の場合分けより

 a_n=b_{n-1}+b_{n-3}\mbox{ ……(1)}

 b_n がフィボナッチ型の漸化式をみたすことは入試数学ではよく知られた事実です。
詳細は割愛しますが,最初の 1 個が  \bigcirc \times かで場合分けすると次の式が導けます。

 b_n=b_{n-1}+b_{n-2},\, b_1=2,\, b_2=3\mbox{ ……(2)}

(1)(2)から  a_n もフィボナッチ型の漸化式をみたすことがわかります。

\begin{align*}
a_n+a_{n+1} &=(b_{n-1}+b_{n-3})+(b_{n}+b_{n-2})\\
&=(b_{n-1}+b_{n})+(b_{n-3}+b_{n-2})\\
&=b_{n+1}+b_{n-1}\\
&=a_{n+2}
\end{align*}

要素数 1, 2 の集合のところで求めたように  a_1=1,  a_2=3 です。
これらを使って a_6,  a_{46} を求めれば解けますが,実はもうひと工夫できます。

 a_{n+2}=a_{n+1}+a_n,  a_1=1,  a_2=3 をみたす  \{a_n\} はリュカ数列  \{L_n\} です。

ja.wikipedia.org

実装

mathematica

組み込み関数を使って  L_1 L_2 L_3 L_6{}^2 L_{46} を求めると速いです。約0.0005秒でした。*1

    In[]:= AbsoluteTiming[
        p[{a_, b_, c_, d_, e_, f_}] := {b, c, d, e, f, 
          BitXor[a, BitAnd[b, c]]};
        ds = CreateDataStructure["DisjointSet"];
        For[i = 0, i < 2^6, i++,
         x = IntegerDigits[i, 2, 6]; y = p@x;
         ds["Insert", x]; ds["Insert", y];
         ds["Unify", x, y]];
        ans = Apply[Times, LucasL[Length@#] & /@ ds["Subsets"]];]
       
    Out[]= {0.0004947, NUll}

python

python でも解きました。計算時間は約 0.001 秒で,mathematica より遅いという意外な結果でした。
dsu は ACL(AtCoder Library) の python 移植版を使いました。
github.com

    import time
    from atcoder.dsu import DSU
    from sympy import lucas
    from math import prod
    
    def p(n):
        a, b, c, d, e, f = [int(x) for x in format(n, '06b')]
        lst = [b, c, d, e, f, a ^ (b & c)]
        return int(''.join(map(str, lst)), 2)
    
    def solve():
        g = DSU(2 ** 6)
        for x in range(0, 2 ** 6):
            g.merge(x, p(x))
        return prod(map(lambda x: lucas(len(x)), g.groups()))
    
    if __name__ == '__main__':
        start = time.time()
        print(solve())
        elapsed_time = time.time() - start
        print("elapsed_time: {0} [sec]".format(elapsed_time))

*1:101問以降なので答えの数値は省略。