砂場で遊ぼう

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

トップ > mathematica > Project Euler 186 / ネットワークの接続性

Project Euler 186 / ネットワークの接続性

mathematica project euler python union-find (dsu)

問題概略

100 万人のユーザをもつ電話システムの通信記録がある。

RecNr Caller Called
1 200007 100053
2 600183 500439
3 600863 701497
... ... ...

n 番目の記録のかけ手と受け手の電話番号は次のラグつきフィボナッチ生成器で決定される。


S_k= \mathrm{Mod}(100003 - 200003k + 300007k^3,\, 10^6)\quad (1\leqq k\leqq 55)

S_k=\mathrm{Mod}(S_{k-24} + S_{k-55},\, 10^6)\quad (k\geqq 56)


\mathrm{Caller}(n)=S_{2n-1}, \mathrm{Called}(n) = S_{2n} である。

\mathrm{Caller}(n) = \mathrm{Called}(n) のとき，ユーザは間違って電話をかけたとされて通信は切断される。

電話をかけた人と受けた人は友達であるとする。 X が Y の友達であり Y が Z の友達であるとき X は Z の友達である。

首相の電話番号は 524287 である。何回電話をかけると 99% のユーザ（首相自身も含む）が首相の友達になるだろうか？ なお，切断された通信は数えない。

https://projecteuler.net/problem=186

これは union-find (dsu) を使う問題です。
電話をかけた人と受けた人を辺で結んでいきます。
辺で結ばれた人の集合が友達のグループで，首相を含むグループの大きさが 99 万になるのに電話は何回必要かを求めます。

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

drive.google.com

mathematica で解く

まず mathematica で解きました。mathematica にも dsu は一応組み込まれていますが，なぜか size（ある要素が含まれる集合の大きさを返す関数）は実装されていません。

自分で書いてもよかったのですが正直面倒でしたし，仮に書いても mathematica では大した速度は出ないと思ったので，あるものをそのまま使いました。

首相を含む部分集合の大きさを計算

最初に考えたのは，mathematica に組み込まれている dsu を使ってまともに部分集合を求めてその大きさを計算する解法です。

電話を 1 回かけるたびにグループの大きさを計算するとものすごく時間がかかるので，二分探索しました。計算時間は約 6 分です。*1

無効なものも含めて電話を 一応 n 回かけたときの首相を含むグループの大きさを len(n) とします。

  1. low, high の初期値として \mathrm{len(low)< target\leqq len(high)} をみたす値をセットする。差は仮に 100 万にしました
  2. 二分探索で \mathrm{len(low)< target\leqq len(high)} かつ \mathrm{high}=\mathrm{low}+1 をみたす low, high を探す
  3. \mathrm{len(high)} から無効な電話の回数を引いたものが答え
    In[]:= AbsoluteTiming[
        percentage = 99;
        target = 10^6*percentage/100;
        pm = 524287;
        s[k_] := s[k] = Piecewise[
           {{Mod[100003 - 200003 k + 300007 k^3, 10^6], 1 <= k <= 55},
            {Mod[s[k - 24] + s[k - 55], 10^6], k >= 56}}];
        len[n_] := Module[{ds = CreateDataStructure["DisjointSet"], i},
          For[i = 1, i <= n, i++,
           ds["Insert", s[2 i - 1]];
           ds["Insert", s[2 i]];
           ds["Unify", s[2 i - 1], s[2 i]]];
          Length@SelectFirst[ds["Subsets"], MemberQ[#, pm] &]];
        high = NestWhile[# + 10^6 &, 0, len@# < target &];
        low = high - 10^6;
        While[high > low + 1,
         mid = Floor[(low + high)/2];
         If[len@mid < target, low = mid, high = mid]];
        ans = high - Length@Select[Range@high, s[2 # - 1] == s[2 # ] &];]
       
    Out[]= {358.637, Null}

グラフの連結成分を求める

次は mathematica に組み込まれている dsu を使わないコードを書きました。

  1. 下のコード中の calls で作った辺を集めてグラフにする
  2. グラフの連結成分を求める組み込み関数 ConnectedComponents を使って，首相を含む友達グループの大きさを求める

前のコードと同じ二分探索で約 4 分です。組み込み関数の dsu を使うメリットはなかったわけですね。

    In[]:= AbsoluteTiming[
        percentage = 99;
        target = 10^6*percentage/100;
        pm = 524287;
        s[k_] := s[k] = Piecewise[
           {{Mod[100003 - 200003 k + 300007 k^3, 10^6], 1 <= k <= 55},
            {Mod[s[k - 24] + s[k - 55], 10^6], k >= 56}}];
        calls[n_] := s[2 n - 1] <-> s[2 n];
        len[n_] := Length@First[ConnectedComponents[calls /@ Range@n, pm], {}];
        high = NestWhile[# + 10^6 &, 0, len@# < target &];
        low = high - 10^6;
        While[high > low + 1,
         mid = Floor[(low + high)/2];
         If[len@mid < target, low = mid, high = mid]];
        ans = high - Length@Select[Range@high, s[2 # - 1] == s[2 # ] &];]
       
    Out[]= {248.302, Null}

python で解く

本命の python で解きます。

ラグつきフィボナッチ生成器はネットでみつけたコードをほぼそのまま使いました。
S_n の値は高々 3 回しか使わないので，長さ 55 の \verb|deque| に入れて使い終わったものはどんどん捨てています。

dsu は ACL(AtCoder Library) の python 移植版を使いました。

github.com

caller と called を merge して size を調べるだけで，二分探索するまでもなく約 5 秒で終わりました。

    import time
    from collections import deque
    from atcoder.dsu import DSU
    
    MOD = 10 ** 6
    PERCENTAGE = 99
    N = 10 ** 6
    TARGET = N * PERCENTAGE / 100
    PM = 524287
    
    def lfg():
        dq = deque(maxlen=55)
        for k in range(1, 56):
            t = (100003 - 200003 * k + 300007 * pow(k, 3, MOD)) % MOD
            yield t
            dq.appendleft(t)
        while True:
            t = (dq[23] + dq[54]) % MOD
            yield t
            dq.appendleft(t)
    
    def solve():
        s = lfg()
        g = DSU(N)
        res = 0
        while g.size(PM) < TARGET:
            caller, called = next(s), next(s)
            if caller == called:
                continue
            res += 1
            g.merge(caller, called)
        return res
    
    if __name__ == '__main__':
        start = time.time()
        print(solve())
        elapsed_time = time.time() - start
        print("elapsed_time: {0} [sec]".format(elapsed_time))

通話記録の中身はどうなってるか

mathematica を使って通話記録をいろいろ調べました。

  • 無効な電話は 3 回だけ
  • 首相は一度も電話をかけない。かかってくるのは 2 回
  • 首相が一人目の友達を得るまでに 100 万回上の電話が必要。この 1 回の電話で首相は一気に 90 万人以上の友人を得る

首相はハブられまくっている内気な上級国民です。

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