ProjectEuler 361 Subsequence of Thue-M - 拼圖

寫一下我解這題的心路歷程好了 XD

引文後有雷

※ 引述《babufong (嗶嗶)》之銘言：
: 361. Subsequence of Thue - Morse sequence
: http://projecteuler.net/problem=361
: Thue - Morse sequence {Tn} 是一個二進位制的數列，他符合以下條件：
: T0 = 0
: T2n = Tn
: T2n+1 = 1 - Tn
: 我們可以知道 {Tn} 的前幾項是這樣的：
: 01101001100101101001011001101001....
: 定義 {An} 是一個整數數列，而他的每一項用二進位制表示都有出現在 {Tn} 之中
: 舉例來說，十進位制的 18 用二進位制表示就是 10010
: 10010 出現在 T8 到 T12 之間，所以 18 是 {An} 的其中一項
: 14 用二進位制表示為 1110，而 1110 從未出現在 {Tn} 中，故 14 不是 {An} 其中一項
: An 的前幾項我們也知道了，如下：
: n 0 1 2 3 4 5 6 7 8 9 10 11 12 ...
: An 0 1 2 3 4 5 6 9 10 11 12 13 18 ...
: 我們還能證明出 A100 = 3251 而 A1000 = 80852364498
: 18
: 請計算出 Σ A10^k 的末九位
: k=1

一開始我是用了另一個等價的 T_n 的定義在做:

從"0"開始 把目前有的字串 0 1 反轉之後接到後面 無限接下去就是 T_n

於是在想到底什麼型式是不會出現的

很容易找到了 111 000 11011 00100

但發現 1101011 0010100 這種也要時我果斷放棄

後來回到原來的定義才發現原來直接用"延伸"的方式去找就好

例如 10010 這一串 表示前面必然有個地方有 10010 => 100 這個字串在

用類似的道理就可以知道

所有夠長的字串都能夠由某個較短的字串做以下四個動作之一得到:

(1) 做 1 -> 10, 0 -> 01
(2) 做 1 -> 10, 0 -> 01 但去尾
(3) 做 1 -> 01, 0 -> 10 但截頭
(4) 做 1 -> 01, 0 -> 10 但截頭去尾

例如 100 做 (1) 得到 100101
做 (2) 得到 10010
做 (3) 得到 11010
做 (4) 得到 1101

仔細比較得到的這四個數 會發現永遠有 (4) < (2) < (3) < (1) < 長一點字串的 (4)

因此 我們會得到結論

長度 4 的字串前半段是長度 2 做 (1)
後半段是長度 3 做 (4)
長度 5 的字串前半段是長度 3 做 (2)
後半段是長度 3 做 (3)
長度 6 的字串前半段是長度 3 做 (1)
後半段是長度 4 做 (4)
依此類推

這樣一來 長度 n 的字串的個數 S(n) 就是如下的數列

1 2 3 5 6 8 10 11 12 14 16 18 20 21 22 23 24 ...

(這就是我在前篇推文提到的數列)

它的遞迴關係是 S(2n) = S(n)+S(n+1) n≧2
S(2n+1) = 2S(n+1) n≧2
S(1) = 1, S(2) = 2, S(3) = 3

那麼根據上面得到的這個規則 A_n 就只要去求 S(n) 的前幾項和到正好超過

再去計算它是那個長度的哪個位置 是前半還是後半

是由長度多少字串的第幾個延伸來的 就可以知道它是什麼字串

結果 A_{10^15} 因為用 Mathematica 直接生字串生到記憶體炸掉了...XD

這中間為了節省時間我做了非常多數學計算

甚至把上面這數列的前 N 項和寫了一個公式出來

(由於這數列等於 1.5n-1.5 加上一個很有規律升降的數列

前 N 項和可以由 Σ(1.5n-1.5) 再去調整 這樣就有公式了)

但無論如何總是快不起來 就算記憶體不炸掉要算出 A_{10^18} 也要半小時

回頭看到題目發現範例特別把 18 (10010) 出現在 T_8 ~ T_12 給說出來

這才打醒我不必要去記整個字串 去記它在字串的哪裡就好了...

而根據給定的 T_n 的定義 這個延伸的動作就變得像是乘以 2 一樣簡單

不過要從在哪裡的字串去求出那一位是什麼也讓我想了一會兒

然後我才想起一件很重要的事: T_n 其實還有另一個等價定義

T_n = {1, 若 n 的二進位表示法中 1 的個數為奇數
{0, 若 n 的二進位表示法中 1 的個數為偶數

也就是俗稱(?)的 Parity

利用這一點 給定範圍之後我能夠一位一位算過去求餘數

發現這一點之後我改用 C 寫

求 Parity 的部份我特別寫了一支小小的組合語言副程式嵌進去

利用了組合語言裡會根據結果的 Parity 設定旗標的功能直接讀出我要的 Parity

(這個副程式組譯之後只有 20 byte 算是滿自豪的 XD

所以就直接寫成字串後當成函式來呼叫了)

不過跑出來的結果還是不對 (這個時候 A_{10^18} 只要約一分鐘就能有結果)

一直到 utomaya 的推文 (A_{10^18} 約是 11 億位二進位) 之後我回頭檢查算式

才發現原來是計算上面那數列的前 N 項和的公式的中間結果溢位了...

把這個溢位改掉之後才終於得到正確答案 orz

總計從我開始解這題的星期日開始 邊忙要交的作業邊解

然後還跨過了這週四某科的期中考

一直到現在下一題快出來了才把這題幹掉...orz

--
來去把 360 也幹掉拿最新五題的成就 @_@

--
LPH [acronym]
= Let Program Heal us
-- New Uncyclopedian Dictionary, Minmei Publishing Co.

--