Project Euler 198 題

Project Euler 的 198 題。簡單的說，要算出分母在 10^8 以內，大小在 0~0.01 之間的「混淆數」總共有幾個 。而所謂的混淆數 x，是指如果對於某個 m，分母不超過 m 的「最近似 x 的分數」不只一個時，就稱 x 為混淆數。比方對於 x=9/40 來說，m=6 時， 1/5 和 1/4 都是分母不超過 6 最近似 x 的分數。
這個題目其實如果沒有任何基礎知識的話，還不太容易，所以解出的人不多。但站在一些基礎知識之上， 10^8 這種數量級是很容易用程式來秒殺的。據說 Knuth 的 The Art of Computer Programming 裡面就有提供相關的演算法。
不過使用 Haskell，可以很直接的用遞迴解出：

f a b l m | a*b>l || (a==1 && b > m)  = 0
         | otherwise = (f a c l m)+1+(f c b l m) where c=a+b
main = putStrLn $ show result
      where result= (f 1 100 l2 m) + 49+l2-m
      l = 10^8
      l2 = l `div` 2
      m = floor (sqrt (fromIntegral l2)

七行就簡單解決了。如果你知道背後的數學，上面的程式碼其實是很直接的。但用 Python 解就麻煩了，因為上面這個遞迴解超出 Python 的遞迴深度。可以用迴圈來取代遞迴，如下面的演算法：

def p198(M=10**8, z=100, cnt=0):
 M2=M/2
 a=m=int(M2**0.5)
 stack=range(z,m)
 while stack:
     b=stack[-1]
     if a*b>M2:
         a=stack.pop()
     else:
         cnt+=1
         stack.append(a+b)
 return cnt+M2-z/2

(程式碼我從 Project Euler 直接貼過來，排版有點亂掉）
這個演算法是我參考別人的解答後寫出的，配上 psyco 後，其實還挺巧妙也挺快的，但不太容易讀懂，要跟蹤一下堆疊是怎麼跑的。
但有些簡單的工作， Haskell 寫起來又很累贅。比方有一個整數 n，想找出它的根號的整數部份，如後印出來，在 python 來說再簡單不過的

print int(x**0.5).

Haskell 的作法要先將整數轉成浮點數，開根號，然後再用 floor 轉成整數，最後再用 show 轉成字串才能印出。所以是

putStrLn.show $ (floor . sqrt) (fromIntegral n)

如果有許多浮點數與整數，那更複雜一點的運算就要轉來轉去。雖說這樣可以強迫正確的設計風格，但有時只是簡單的小計算，這樣也未免太累了。
果然 Haskell 能將複雜工作簡單化，簡單工作複雜化。

Read More

zhpy 的另類用法

配合上 zhpy ，我們可以寫出下面這樣的 python 程式:

# encoding: zhpy_utf8
from math import *
from operator import *
print Σ(range(100))
print sin(π/4) ≠ √(2)/2
print ㏒(100), ㏑(e)
print 5 × 30 ÷ 2 ≦ tan(√(5)*π/4)
Π=λ f:reduce(mul, f, 1)
print Π(range(1,6))
print set([1, 2, 3, 4, 5]) ∩ set([1, 3, 5, 7, 9])

(第一行是 MagicCodec 語法，普通的 zhpy 請去掉第一行)
不管實不實用，至少增加了一點可讀性。既然是λ，為什麼要寫成 lambda?
以上的程式碼只需要下面的 zhpy .ini 檔案（可以叫做 math.ini，放在執行的目錄下)

[math]
λ=lambda
π=pi
Σ=sum
≠=!=
÷=/
×=*
㏒=log10
㏑=log
√=sqrt
≧=>=
≦=<= ∩=& ∪=l

for 和 in 可以換成相對應的符號, set 的 <, <= 也可以換成 set 符號。
Haskell 的 ide Lesksah 就能做到這件事情。 如果能夠把 lyx 修改成 python, haskell 或者其他程式語言的編輯器，應該也挺有趣的。
相關的東西是 TeXmacs 的 python plugin tmPython，

Read More

根號 2 到兩百萬位的挑戰

Sphere Online Judge 有一個題目是在 20 秒內計算根號 2，越多位越好（上限是兩百萬位）。有不少人能達到這個目標，但都是使用 C 或者 C++。我想試試看 python 這種以速度慢聞名的語言，能達到什麼程度。試驗的結果是十五萬位。雖然不太多，但是也擠進了排行榜裡面前二十名。而目前 python 語言的第二名還不到五萬位。

我的程式沒有用什麼特殊的技巧。事實上， python 的乘法不算太慢，現在 python 已經使用了 Karatsuba 演算法來做乘法，在一萬位以內的效能都還不錯，算到十幾萬位也都還能接受。但再上去就有點吃力了。而最主要的瓶頸在於 long 到 str 的轉換速度。光是將一個兩百萬位的數字轉成字串，就會花上非常久的時間。

改用 gmpy 會快很多，應該能夠算到百萬位以上，雖然我沒試過，但 spoj 應該不能使用 gmpy。

所以，也許應該直接寫個簡單的 FFT 乘法模組。事實上，我找到了一個叫做 DecInt 的 Pure python 套件，裡面實現了幾種快速的乘法及除法演算法，而且因為基本上資料是用十進位表示，數字轉字串非常快。他的乘法也許不是最快的，但是如果要我自己寫，大概也會是類似的東西，所以相當適合先拿來測試一下效能。

果然，大致上可以讓算十五萬位的速度快兩三倍。如果自己再強化改進，也許能再快個兩倍，但估計最多只能勉強上到百萬位，而且需要相當長時間的測試和實驗，所以我放棄。

於是我改用 Haskell 來寫。我沒有寫過 Haskell 程式（現在也還不會），而且連 Tutorial 都沒有讀完，就抓了 ghc 用最笨的方式把程式寫出來了。我只需要知道怎麼用 = 來定義函數能寫了。結果算到了一百六十萬位，足足比 python 多了十倍(Haskell 現在第二名約一百二十萬位)。如果再多加微調，或者我如果對 Haskell 再多熟悉一點，也許還能多個不少位，甚至到兩百萬位也不無可能。

其實我只要算到五十萬位，然後再用個 4-way Toom Cook 就行了，而算到四十萬位不用三秒。不過問題是，就像前面說得，其實我還不會寫 Haskell。

演算法基本上和 python 版本的相同，不過改以 recursive 的寫法。 Python 雖然完全支援 recursive 函數，但是 Python 的函數呼叫很慢，又限制層數（內定1000，最多 5000)，所以，其實有點降低我把東西寫成遞迴的誘因。反之，Haskell 是 functional language，所以其實是獎勵這種寫法。

而 Haskell 之所以算得這麼快，其實有點勝之不武，因為 ghc 內建用 gmp 來做計算。

讓人疑惑的反而是，為什麼 python 不用 gmp 來做計算？

當演算法夠好，時間夠長， python 的速度劣勢就能被彌補，但 online judge 的時間限制都不太長，二十秒已經算是超級長的了，但這二十秒差不多只等於我電腦上的五、六秒，典型的時間限制是兩秒。用 pure python 實現很好的乘法演算法，也要到上萬位才能擊敗 python 內建的長整數乘法。所以，現在想用純 python（而且程式碼在4kb以內）上兩百萬位，還挺困難的，不過，也許還是有人能辦到也說不定。

(Update 2008-06-11) Haskell 達成兩百萬位。當然我 Haskell 的能力比兩天前進步一點，但主要還是用笨方法寫 recursive function。程式碼總共 11 行，而且其中有 3 行是用來定義一個函數，來達到 python 的 zfill 功能。

Read More