ハノイの塔

Pythonハノイの塔を解くプログラムを書きました。
GUIの表示にはmatplotlibの棒グラフを使いました。

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

def plot():
    plt.cla()
    plt.xlim(-0.5,2.5)
    plt.ylim(0,N)
    left = np.array([0,1,2])
    numl = [0,0,0]
    plt.bar(left,lis[N-1],color="0.0",linewidth = 0.5)
    numl = lis[N-1]
    for i in (reversed(range(N-1))):
        plt.bar(left,lis[i],bottom=numl,color=str(0.95-0.95*(i/N)),linewidth = 1)
        numl = lisum(numl,lis[i])
    plt.pause(0.3)

def hanoi(num,fro,des,wor):
    if num == 1:
        move(1,des)
        plot()
    else:
        hanoi(num-1,fro,wor,des)
        move(num,des)
        plot()
        hanoi(num-1,wor,des,fro)
def move(n,tar):
    s = [0,0,0]
    s[tar] = 1
    lis[n-1] = s

def lisum(a,b):
    s = []
    for i in range(len(a)):
        s.append(a[i]+b[i])
    return s
fig = plt.figure()
N = int(input())
lis = []
for i in range(N):
    lis.append([1,0,0])
plot()
hanoi(N,0,2,1)

実際にハノイの塔を解くアルゴリズムはこの部分

def hanoi(num,fro,des,wor):
    if num == 1:
        move(1,des)
        plot()
    else:
        hanoi(num-1,fro,wor,des)
        move(num,des)
        plot()
        hanoi(num-1,wor,des,fro)

一見めんどくさそうなゲームに見えるハノイの塔がこれだけのコードで解けるのは驚きです。
def hanoi(int,int,int,int)は、
num...解くハノイの塔のディスクの数
fro...棒の番号を左から0,1,2と振ったときのディスクがnum個ある開始位置の番号
des...num個のディスクを移動させたい棒の番号
wor...操作において利用可能な棒の番号
fro != des != worが成り立ちます。
move(n,des)関数は、n番目のディスク、(小さい方から1..nとしています)をdesに移動させる操作です。
手続きを見ていきます。
基底: num = 1のとき
ディスクが一つのみなので、1のディスクをdesに移動させて終わりです。
帰納: num = nのとき
まずnum-1までのディスクを操作において利用可能な棒、worにすべて移動させます。これは、hanoi(num-1,fro,wor,des)で可能です。
次に、動かせるようになったnum番目のディスクをdesに移動させます。これはmoveの一手で可能です。
最後に、desに移動したnum番目のディスクの上にworに乗っているnum-1個のディスクをのせて終わりです。これはhanoi(num-1,wor,des,fro)で可能です。

emacsの模様替え

emacsは文字や背景の色を設定できるプリセットテーマを搭載しています。今まで初期設定でやってきたので、emacsのテーマとついでにフォントも変更してみました。

before
f:id:xxxasdfghjk999:20180207003438p:plain

after
f:id:xxxasdfghjk999:20180207003433p:plain

フォントの変更、テーマのロードはemacs.d/init.elに

(load-theme 'wombat t)

;; Basic font
(custom-set-faces
 '(default ((t (
                :foundry "Hack" :family "Hack")))))

;; Japanese font
(set-fontset-font t 'japanese-jisx0208 (font-spec :family "GenShinGothic"))

を書き込みました。
英数字のフォントはHack、日本語には源真フォントを導入しました。
両者とも太い字で見やすいです。

ポーカー遊び完結編

前回のコードに加え、リストをシャッフルする汎用関数(トランプデッキにも適用可能)などを実装しました。デッキを人数分に分けるdistributeも作りました。

randompick [] = return Nothing
randompick xs = do
                x <- randomRIO (0,length xs - 1) :: IO Int
                return (Just (xs !! x))
fromJust (Just a) = a
pickcard deck = do
                x <- randompick deck
                if x /= Nothing then return (x,delete (fromJust x) deck)
                                else return (Nothing,[])
shuffle deck = do
               (card,d) <- pickcard deck
               if card /= Nothing then
                do
                  f <- shuffle d
                  return ((fromJust card):f)
                  else return []
addone [] _ = []
addone [[]] (d:deck) =  [[d]]
addone (x:xs) [] = (x:xs)
addone (x:xs) (d:deck)  = (d:x):(addone xs deck)

distribute n deck = let s = replicate n []
                        distribi l [] = l
                        distribi l deck = let xs = addone l deck
                                          in
                                          distribi xs (drop n deck)
                    in
                    distribi s deck
draw deck 0 = []
draw (x:xs) n = x:(draw xs (n-1))
prepeat p n = do
     x <- shuffle deck
     let d = suitsort (draw x 5)
         (PR s _) = culcPrize d
     if s == p then do print p
                       print d
                       print n
               else prepeat p (n+1)

ポーカーをこれでちょっと遊べて満足したので、ひとまずはこれでおしまいにします。

Haskellでポーカー?

長いコードになりましたので続きから

続きを読む

Haskellの型クラスを利用した有理数型

タイトル通りです。

data Ration = Integer :/ Integer | AddID 
instance Eq Ration where
         (x :/ y) == (x':/y') = let s1 = gcd x y
                                    s2 = gcd x' y'
                     in
                 (x `div` s1 == x' `div` s2)
                                           && (y `div` s1 == y' `div` s2)
instance Show Ration where
         show (x :/ y) = let s = gcd x y
                         in
                        show (x`div`s) ++ "/" ++ show (y`div`s)

rationToFloat :: Ration -> Float
rationToFloat (x :/ y) = fromIntegral x / fromIntegral y
mul (x :/ y) (x' :/ y') = reduct ((x * x'):/(y * y'))
reduct (x :/ y) = let g = gcd x y
                  in
                  (x `div` g) :/ (y`div`g)
di (x :/ y) (x' :/  y') = mul (x :/ y) (y' :/ x')
add x AddID = x
add AddID y = y
add (x :/ y) (x' :/ y') = let l = lcm y y'
                              m1 = l `div` y
                              m2 = l `div` y'
                          in reduct ((m1*x + m2*x') :/ l)
sub (x :/ y) (x' :/ y') = add (x :/ y) ((-x') :/ (y'))

有理型は和も積も単位元が複数存在するのが厄介です。
さらに0除算のことなどいろいろな例外があります。
今回のコードは急ごしらえでそのようなことは考えていません。
和の単位元だけ試験的に実装しています。(AddID)
今回はこれを応用して、ライプニッツ級数
\[\sum^{\infty}_{n=0}\frac{(-1)^n}{2n+1} = \frac{\pi}{4} \]
より円周率を求めてみます。
利用する関数、定義

let lip = [ ((-1)^n):/(2*n+1) | n<-[0,1..]]
partsum b 0 _ = b
partsum b n (x:xs) = partsum (add b x) (n-1) xs

実行例

*Main> ((4:/1) `mul` (partsum AddID 100 lip))
8252079759413970386664454687621174435983101115012912631997769614579677862845786070667088/2635106162757236442495826303084698495565581115509040892412867358728390766099042109898375
*Main> ((4:/1) `mul` (partsum AddID 200 lip))
731070388038968202720510198283781775454721070454247644121830326235077353957426480846277182638522756913731374269936715193263560910895530282361264868908266204946805311497184/233077884665326398150824103062920077329319208418023678398907683319848459110106903491713919541172119871364032797568478630600311502003066966840474629472267552661106168111125
*Main> ((4:/1) `mul` (partsum AddID 300 lip))
5225645700732578660757222494473018434678943097092985751800195180838596610013780584278253709093305256797047031098916493425558192295042537198564257121082268188433032685860147637228857628560723477301587102016574531025314432495082724959309868725561540131084688/1665141453284292119501128176265135852020805040785019269013402590677638665112799180271757599577286560166540552724725936015136272490152610881122344030143206974706298054993350415948784680656377567820299938133978899260687047126681401325662561292696661412536625
(実数値で上から)
3.131592903558553
3.136592684838817
3.1382593295155905

収束遅い...

高速フィボナッチ

Haskellで書きました。

matrix2 (x1,y1,z1,w1) (x2,y2,z2,w2) = (x2*x1+y1*z2 , x1*y2+y1*w2,x2*z1+w1*z2,y2*z1 + w1*w2)
repeat' n f x y
       | n < 1 = f x y 
       | n == 1 = f x y
       | n > 1 = repeat' (n-1) f (f x y) y
third (x,y,z,w) = z
matrxpow m n
      | n == 0 = (1,0,0,1)
      | n `mod` 2 == 0 = matrxpow (matrix2 m m) (n `div` 2)
      | otherwise = matrix2 m (matrxpow (matrix2 m m) (n `div` 2))
fib x = third ( matrxpow (1,1,1,0) (x))

フィボナッチ数列、というか線形な漸化式は行列で表示することができます。
まず、

\[ \left(
\begin{array}{cc}
a_{n+2} \\
a_{n+1}
\end{array}
\right)
=
\left(
\begin{array}{cc}
1 & 1 \\
1 & 0
\end{array}
\right)
\left(
\begin{array}{cc}
a_{n+1} \\
a_{n}
\end{array}
\right)
\]
となります。
これを再帰的に用いると、


\[ \left(
\begin{array}{cc}
a_{n+2} \\
a_{n+1}
\end{array}
\right)
=
\left(
\begin{array}{cc}
1 & 1 \\
1 & 0
\end{array}
\right)^{n+1}
\left(
\begin{array}{cc}
a_{1} \\
a_{0}
\end{array}
\right)
\]
\[(a_0 = 0 , a_1 = 1)\]
となり、行列のべきと初期条件の積の形に変形できました。
今回のプログラムではこの行列のべきを計算し、2×2行列の左下の要素を取ってくることでフィボナッチ数列を計算しています。
行列のべきを求める関数 matrxpow m n (Int*Int*Int*Int -> Int -> Int*Int*Int*Int)には、行列の積の結合律を利用した高速化処理がなされています。
基底 : \[A^0 = I\]
行列の0乗は単位行列です。
帰納 : \[A^{2k+1} = A(AA)^k , k \in \mathbb{Z} , k≧0\]
\[A^{2k} = (AA)^k , k \in \mathbb{Z} , k≧0\]
べき数の偶奇によって場合分けします。

エラストテネスのふるい

あけましておめでとうございます。
エラストテネスのふるいを書きました(byHaskell)

erast t = let list = [2,3..t]
              erasti [] = []
              erasti (x:xs) = x : erasti (filter (\y -> (y `mod` x) /= 0) xs)
          in
          erasti list

ふるいに用いられる数字は素数であることが確定するので、再帰的に書くことができます。