tomabouの日記

Haskellなどで勉強したことなどを書いていきます

セルオートマトンで遊ぶ

はじめに

皆さんパソコンの背景は何にしていますか?なんか背景を変えたいなと思い、色彩センスはゼロですがいい感じの規則的な模様を勝手に生成させて背景にすることができたので紹介します。

一次元セルオートマトン

一次元セルオートマトンをご存知でしょうか
セル・オートマトン - Wikipedia
二次元の場合はライフゲームと言われたりします。簡単に言うと一列に0と1が並んでいて、1ステップごとに0と1が周りの0,1の配置に応じて変化するというものです。単純なものではi番目の数字は1ステップ前のi-1, i, i+1番目の数字(8通り考えられる)のみに依存して変わるもので、ルールを次のような表にまとめることが出来ます。

1ステップ前の数字 0になるか1になるか
0,0,0 0
0,0,1 1
0,1,0 1
0,1,1 1
1,0,0 1
1,0,1 0
1,1,0 0
1,1,1 0

このようなルールは256通り考えることが出来て、上の表を二進数だと思うことでルール0からルール255と名前がついています。上の場合はルール30(00011110)です。

興味深い性質

ほとんど意味のないルールもありますが(例えばルール0は1ステップですべて0になって変化しません)、ルール18,22,30,54,60,90,110などは面白い挙動を示すらしいです。ちなみにルール110はマリオメーカーと同様にチューリング完全であるらしいです。チューリング完全オタクの力は怖いですね。*1

画像生成

一次元の数字の列が時系列に沿って変化していくので、それを並べると二次元の0,1の表になります。せっかくなのでルールを入力して画像を出力できるようにしてみました。
pythonで書かれていて、実行時引数としてルールと初期条件を与えることが出来ます。

simpleのときは真ん中に1が一つで残り0。randomのときはその次に与えた引数の確率で1になります。確率が0に近いときと0.5に近いときで異なる見た目のものが生成されるのでいろいろ試してみてください。最後に自分が背景に使用している画像を張っておきます。
ライブラリとしてはPILを使用しています。

実行結果

入力

python automaton.py 30 simple

f:id:tomabou:20171008173553p:plain

python automaton.py 90 random 0.02

f:id:tomabou:20171008173116p:plain

デスクトップ背景

コード

from PIL import Image
import sys
import numpy as np
import random

def make_image(rule,mode = "simple",p=0.5):
    color1 = (0x66,0x99,0x99) #0の色
    color2 = (0xff,0x99,0x99) #1の色
    x = 1920
    y = 1080 #出力画像サイズ 今はフルHDの画像を出力します
    dotsize = 2 #一つのセルがドットいくつ分か
    filename = "rule"+str(rule)+"_"+str(x)+"x"+str(y)+"_"+str(dotsize)+"_"+mode+".png"
    img = Image.new('RGB',(x,y),color1) 

    rule_list = list()
    for i in range(8):
        rule_list.append(rule %2)
        rule = rule // 2

    print (rule_list)
    
    nx = x//dotsize + 1
    ny = y//dotsize + 1
    cell_map =[[0 for i in range(ny)] for j in range(nx+2) ] 
    if mode=="simple":
        cell_map[nx//(2)][0] = 1
    elif mode =="random":
        for i in range(nx):
            cell_map[i][0]=1 if random.random()<p else 0

    cell_map[nx][0]=cell_map[0][0]
    cell_map[nx+1][0]=cell_map[1][0]

    for j in range(ny-1):
        for i in range(nx):
            num = 4*cell_map[i][j] + 2* cell_map[i+1][j]+cell_map[i+2][j]
            cell_map[i+1][j+1] = rule_list[num]

        cell_map[0][j+1]=cell_map[nx][j+1]
        cell_map[nx+1][j+1]=cell_map[1][j+1]
        

    for i in range(x):
        for j in range(y):
            if cell_map[i//dotsize][j//dotsize]==1:
                img.putpixel((i,j),color2)
    
    img.save(filename)
    return img

if __name__ == "__main__":
    argv = sys.argv
    if len(argv)==1:
        rule = 30 
    else:
         rule = int(argv[1])
    if len(argv)<=2:
        mode = "simple"
    else:
        mode = argv[2]
    
    if len(argv)>3:
        p = float(argv[3])
    else:
        p = 0.5

    img = make_image(rule,mode,p)
    img.show()

*1:ちなみにマリオメーカーチューリング完全であることを示すのに使用されたcyclic tag systemチューリング完全性はこのルール110のチューリング完全性を示すために用いられたものなので深い関係が実はあります

AIの「語彙力」を制御する

「よくある単語」判定

りんなというMicrosoftの作ったチャットAIがあります。りんなにしりとりを挑むと「そんな単語知るかよ」という単語を連発してきて人間はとてもではないですが勝てません。しりとりのような言葉を用いた遊びをするAIを作成するには適切にAIの語彙力を制限する必要があります。今回は日本語コーパスを利用することで「簡単な単語」かどうかを簡単に判定できる、ということを紹介したいと思います。

日本語コーパス

コーパスというのは言語の研究などをするために大量の文章を単語ごとに分けて、単語の働きで分類してタグ付けしたものです。みるからに死ぬほど手間がかかりますが、国立国語研究所が作成しており、単語リストは無料で公開されています。
その単語リストを利用すると各単語がどれほどの頻度で出現するかを知ることが出来るので、簡単に「よくある単語」と「難しい単語」の識別が出来ます。

具体例

試しに次のようなアプリを作ってみました。ある三文字の単語から別の三文字の単語まで、一文字づつ変えることで到達するルート(例えば、課題→指定→視点→資金→付近→布団)を教えてくれます。
以下のコードをダウンロードしたBCCWJ_frequencylist_suw_ver1_0.tsvと同じディレクトリにおいて実行すると、できるだけ簡単な単語を使ってできるだけ短いルートでスタートからゴールまで変化させてくれます。
なお出力される二番目の数字は一億単語中何回現れるかの頻度を表しています。

出力例
('カダイ', 10, 'スタート')
('カテイ', 13306.0, '家庭')
('シテイ', 13030.0, '指定')
('シテン', 3687.0, '視点')
('シキン', 11181.0, '資金')
('フキン', 3552.0, '付近')
('フトン', 10, 'ゴール')
コード

3000という値はパスの長さと単語の簡単さどちらを重視するかのパラメータです。小さくすると短いものを重視します。最短パスはダイクストラで求めています。

import re
import heapq

start_word = "カダイ"
goal_word = "フトン"
length_or_frequency = 3000

PATTERN = re.compile("普通名詞")
print("hello")
KATAKANA = [chr(i) for i in range(12449, 12532+1)]
FILE = open('BCCWJ_frequencylist_suw_ver1_0.tsv', 'r', errors='replace', encoding='utf-8')
index_profile = list()
yomi_to_index = dict() 

def make_node(index, yomi, frequency, kanji):
    """make node for Dijkstra"""
    index_profile.append((yomi, frequency, kanji))
    if not yomi in yomi_to_index:
        yomi_to_index[yomi] = [index]
    else:
        yomi_to_index[yomi].append(index)


i = 0
for line in FILE:
    l = line.split('\t')
    mOB = PATTERN.search(l[3])
    if mOB:
        if len(l[1]) == 3:
            make_node(i, l[1], float(l[6]), l[2])
            i = i+1


START = i
GOAL = i+1

def search_path():
    """do Dijkstra"""
    make_node(START, goal_word, 10, "ゴール")
    make_node(GOAL, start_word, 10, 'スタート')

    que = []
    heapq.heappush(que, (0, START))

    costs = dict()
    preindex = dict()

    while que != []:
        (cost, v) = heapq.heappop(que)
        (oyomi, ofreq, ka) = index_profile[v]

        for nyomi in [oyomi[0]+oyomi[1]+c for c in KATAKANA] \
            + [oyomi[0]+c+oyomi[2] for c in KATAKANA] \
            + [c+oyomi[1]+oyomi[2] for c in KATAKANA]:

            if oyomi != nyomi and nyomi in yomi_to_index:
                for nindex in yomi_to_index[nyomi]:
                    ncost = cost + 1 + length_or_frequency/index_profile[nindex][1]
                    if not nindex in costs:
                        costs[nindex] = ncost
                        preindex[nindex] = v
                        heapq.heappush(que, (ncost, nindex))
                    elif ncost < costs[nindex]:
                        heapq.heappush(que, (ncost, nindex))
                        costs[nindex] = ncost
                        preindex[nindex] = v

    if GOAL not in costs:
        print("cannot reach")
    else:
        ind = GOAL
        while True:
            print(index_profile[ind])
            if ind==START:
                break
            ind = preindex[ind]
        
search_path()

忘れてしまうsegment木

競技プログラミングは忘れたころにやりたくなります
セグメント木は便利なので忘れたころに使いたくなります
しかしライブラリを作ろうにも使いたくなるころには忘れているので紛失してしまいます。
なので備忘録変わりにセグメント木をあげておきます

普通のセグメント木

プログラミングコンテストチャレンジブックに書いてあるセグメント木は以下のようなものです。
数列と、結合法則が成り立ち零元が存在する(モノイドという)演算に対して次のことが出来ます。(例えば最小値をとる演算とします)

  • ある区間の最小値をとる
  • ある値を更新する
#include<algorithm>
#include<functional>
#include<vector>

template<typename T>
class segtree {
public:
	segtree(int n) {//数列の長さを指定してセグメント木を作る
		for (size = 1; size < n; size *= 2);
		vec.resize(2 * size -1, zero);
	}
	segtree(std::vector<T> v) {//vectorで数列を与えてセグメント木を作る
		int n = v.size();
		for (size = 1; size < n; size *= 2);
		vec.resize(2 * size -1, zero);
		for (int i = 0; i < n; i++) {
			vec[i + size-1] = v[i];
		}
		for (int i = size - 2; i >= 0; i--) {
			vec[i] = func(vec[i * 2 + 1], vec[i * 2 + 2]);
		}
	}
	void set(int n, T x) {//値の更新
		n += size - 1;
		while (n >= 0) {
			vec[n] = func(vec[n], x);
			n = (n + 1) / 2 - 1;
		}
	}
	int get(int a, int b) {//値の取得([a,b)で指定)
		return get_func(0, 0, size, a, b);
	}
private:
	T zero =2147483647;//モノイドの零元 ここを下の関数と一緒に適切なものに変えればよい
	std::function<T(T, T)> func = [](T a, T b) {return std::min(a, b); };//演算
	std::vector<T> vec;
	int size;
	
	int get_func(int n, int p, int q, int x, int y) {
		if (q - p == 1)return vec[n];
		if (x == p && q == y) return vec[n];
		int mean = (p + q) / 2;
		T l = zero;
		T r = zero;
		if (x < mean)l = get_func(n * 2 + 1, p, mean, x, std::min(y,mean));
		if (mean < y)r = get_func(n * 2 + 2, mean, q, std::max(x,mean), y);
		return func(l, r);
	}
};

普通のセグメント木の双対

数列と可換なモノイドに対して次のことが出来ます(例えば和とるとします)

  • ある区間全体に同じ数字を足す
  • ある値を取得する

これは普通のセグメント木のデータの流れを逆向きにすると実装できます(このコードはその双対が分かりやすいと思います)

template<typename T>
class segtree2 {
public:
	segtree2(int n) {
		for (size = 1; size < n; size *= 2);
		vec.resize(2 * size - 1, zero);
	}
	segtree2(std::vector<T> v) {
		int n = v.size();
		for (size = 1; size < n; size *= 2);
		vec.resize(2 * size - 1, zero);
		for (int i = 0; i < n; i++) {
			vec[i + size - 1] = v[i];
		}
		for (int i = size - 2; i >= 0; i--) {
			vec[i] = func(vec[i * 2 + 1], vec[i * 2 + 2]);
		}
	}
	T get(int n) {
		T ans = zero;
		n += size - 1;
		while (n >= 0) {
			ans = func(vec[n], ans);
			n = (n + 1) / 2 - 1;
		}
		return ans;
	}
	void set(int a, int b, T x) {
		return set_func(0, 0, size, a, b,x);
	}
private:
	T zero = 0;
	std::function<T(T, T)> func = [](T a, T b) {return a+ b; };
	std::vector<T> vec;
	int size;

	void set_func(int n, int p, int q, int x, int y, T val) {
		if (q - p == 1) {
			vec[n] = func(vec[n],val);
			return;
		}
		if (x == p && q == y) {
			vec[n] = func(vec[n], val);
			return ;
		}
		int mean = (p + q) / 2;
		if (x < mean)set_func(n * 2 + 1, p, mean, x, std::min(y, mean),val);
		if (mean < y)set_func(n * 2 + 2, mean, q, std::max(x, mean), y,val);
		return;
	}
};

c++コンパイラの頭の良さの比較

はじめに

windows10は、bash on ubuntu on windows なるものが登場し簡単にgccが導入できます。
visual studio で使用しているマイクロソフトコンパイラ(以下vc++と呼称)とどちらが高速なコードを生成できるのか調べてみたくなりました。
簡単なコードを比較し、どちらが高速か比較してみましょう

計算させるもの

次のような計算をさせたいと思います。(大体一秒程度で終わるように調整しました)

    auto start = chrono::system_clock::now();

    /*unsigned */long long x=0;
    for(/*unsigned */long long i=0;i< 1+(1ll<<28);i++){
        x+=i;
        x%=13;
    }
    auto end = chrono::system_clock::now();
    
    cout<<x<<endl;
コンパイラ かかった時間(10回平均)(signed) かかった時間(10回平均)(unsigned)
g++ on bash 914 ms 811 ms
vc++ 1089 ms 896 ms

符号付きの場合2割近く符号なしの場合一割ほどvc++が遅いです。また、符号なしの場合が符号付きの場合に比べて早くなっています。なぜこのような差が生まれたのかアセンブリを出力させて調べてみます。

vc++とg++の違い

g++では -S オプションで、vc++では /FA オプションでアセンブリを出力させることができます。
アセンブリにはAT&T構文とintel構文があるようですが、個人的にはintel構文の方が読みやすいです。g++では-masm=intelintel構文のアセンブリを出力させることができます。
上記の計算部分に対応しているところだけを抜き出します。

vc++(unsigned)のアセンブリ


ループの内部のアセンブリを抜き出すと以下のようになります(レジスタを0に初期化したりしているところは省略)

	mov	r8, 5675921253449092805	 ; 謎の数字
$LL64@main:
; Line 11
	add	rsi, rcx     ; rcxはiに対応 rsiはxに対応 この行は x+=i
; Line 12
	mov	rax, r8      ; r8に入っている謎の数字をraxに格納
	mul	rsi          ; mulはraxとの積を上位64bitをrdx、下位64bitをraxに格納
	inc	rcx          ; iをインクリメント
	shr	rdx, 2       ; 上位64bitを右に2ビットシフトし、掛け算の結果の上位66bitを求める
	imul	rax, rdx, 13 ; 掛け算の結果の上位66bitを13倍
	sub	rsi, rax     ; xから上で計算した数字を引いている
	cmp	rcx, 268435457    ; for文の条件文
	jb	SHORT $LL64@main

謎の数字 5675921253449092805があります。これは何なのでしょうか。
以下のような関係式が成り立ちます
5675921253449092805 * 13 = 2^{66} + 1
つまり
\displaystyle \frac{5675921253449092805}{2^{66}}  \approx \frac{1}{13}
2のべき乗での割り算はビットシフトで行うことが出来るため、時間のかかるdiv命令を使わなくて済むという高速化です。

上のアセンブリは以上のテクニックを生かして13で割った商を高速に求め、
その商を13倍して元の数から引くという方法であまりを計算していることがわかります。
頭いい…

g++(unsigned)のアセンブリ

ではそのvc++のコードより一割早いg++の生成したコードは何が違うのでしょうか

	movabs	rdi, 5675921253449092805 ; 同じ
.L2:
	lea	rsi, [r12+rcx]   ; r12+rcxをrsiに格納 x+=i
	add	rcx, 1           ;インクリメント
	mov	rax, rsi         ; xをraxに移して
	mul	rdi              ; 謎の数をかける
	shr	rdx, 2           ; 上位66bitをとる
	lea	rax, [rdx+rdx*2] ; raxにrdxの3倍を格納
	lea	rax, [rdx+rax*4] ; raxにrdx+rax*4を格納 つまり13倍している
	sub	rsi, rax         ; xから商の13倍を引いて
	cmp	rcx, 268435457   ; 1+(1<<28)と比較する
	mov	r12, rsi
	jne	.L2

13倍するコードを素直にmul命令を使わずlea命令に置き換えています。
lea命令はもともとアドレス計算用の命令ですが、算術演算にも使用でき、
フラグを立てず、使用する演算ユニットが違うため、*1もろもろ良いみたいです。
1,2,4倍にしかできないため
13 = 1+ 4*(1+2*1)
といううまい変形をすることでlea命令で計算できるようにしています。
めっちゃ頭いい…
なおjumpする命令が違ったりしますが、それはcmp命令が建てたフラグのうちどれを使っているかが違うだけです。

結論

どちらも頭が良いが、gccの方が頭がいい

unsignedとsignedの違い

なぜ符号付きにすると大幅に遅くなってしまうのでしょうか。
符号付き整数を表現するために2の補数表現を用いているということ、ビットシフトには論理シフトと算術シフトがあることを思い出します。

vc++の出力アセンブリ
	mov	r8, -5675921253449092805		
$LL64@main:
; Line 11
	add	rsi, rcx
; Line 12
	mov	rax, r8
	imul	rsi
	inc	rcx
	sar	rdx, 2
	mov	rax, rdx       ;rdxをrdaにコピー
	shr	rax, 63    ;63ビットシフトするという謎の処理が挟み込まれている 
	add	rdx, rax       ;63ずらしたものを足している
	imul	rax, rdx, 13
	add	rsi, rax
	cmp	rcx, 268435457	
	jl	SHORT $LL64@main

63右に論理シフトすると、もともとの数の最上位のビットが1つまり負であれば1、正であれば0です。
商が負であるときに1足すという処理をしています。

g++の出力したアセンブリ
	movabs	rdi, 5675921253449092805

.L2:
	lea	rsi, [r12+rcx]
	add	rcx, 1
	mov	rax, rsi
	imul	rdi
	mov	rax, rsi
	sar	rax, 63         ;同じく63右にビットシフトしている
	sar	rdx, 2
	sub	rdx, rax
	lea	rax, [rdx+rdx*2]
	lea	rax, [rdx+rax*4]
	sub	rsi, rax
	cmp	rcx, 268435457
	mov	r12, rsi
	jne	.L2

sarは算術シフトです。(shift arithmetic right)最上位ビットが1のとき、63ビットシフトすると全部1のビットが立った数字になり、2の補数表現では-1です。
つまり割られる数xが負ならば-1を引くという処理です

これらの計算はc++の余りの計算の仕様が数学的に不自然なため生じた処理です。c++で商が負のときは0に近い方に切り捨て、余りは負になると定められています(c++03までは処理系依存でした)これは余りが正になったり負になったりして不自然です。
割り算を掛け算に変更するテクニックを用いると、商が負のときも数学的に自然な商を超えない最大の整数が求まり、c++の求める答えと1ずれます。
そこのずれを修正するために無駄な計算をしているようです。

結論

人間の定めた仕様は頭悪い

Appendix

環境
CPU core i7-4510U
memory 8GB
OS windows 10
コード全文
#include<iostream>
#include<chrono>

using namespace std;

int main(){
    auto start = chrono::system_clock::now();

    long long x=0;
    for(long long i=0;i<1 + (1ll<<28);i++){
        x+=i;
        x%=13;
    }
    auto end = chrono::system_clock::now();
    
    cout<<x+2<<endl;

    auto diff = end - start;
    std::cout << "elapsed time = "
              <<chrono::duration_cast<std::chrono::milliseconds>(diff).count()
              << " msec."
              << endl;

    return 0;
}

Haskellで10を作るゲームを解く

はじめに

切符に書かれた四つの数字を使用して10を作る有名な遊びがあります*1
目的の数字を与えられた数字と四則演算を使用して作るプログラムを書いてみました。言語は最近ハマっているHaskellを使用してみました。

構想

式を定義
→与えられた数字から式を列挙
→式を計算
→目的の数字が出たら勝利
このような流れで書きます。工夫の余地がある点は式を列挙するパートで、できるだけ重複しないで列挙したいです。例えば(3+2)*4と4*(2+3)どちらも試すのはアホです。なので次のように正規形を定義して正規形のみ列挙することにします。(与えられた数字は以下狭義単調増加に並んでいるとします)

正規形::任意の演算子について左側にある数字の最小値<右側にある数字の最小値

例えば(1*3+2)*4は正規形ですが、(2+3*1)*4は正規形ではないです。任意の式を正規形にするためには3-2や3/2を2 [演算子] 3と書く必要があります。a - b = b $- a , a / b = b $/ a として新しい演算子を定めることで任意の式をくるくるひっくり返して正規形にすることができます。

実装

式の定義

data Equation = Number Integer | Operator String Equation Equation 

日本語で書けば式とは数字単体、もしくは一つの演算子と二つの式の組である。非常に明快です。二分木の変形ですね。

式の表示

instance Show Equation where 
    show (Number x) = show x
    show (Operator ('$':xs) a b) = "(" ++ show b ++ " " ++ xs ++" "++show a ++")"
    show (Operator xs a b) ="("++ (show a) ++" "++ xs++" " ++(show b)++")"

数字はそのまま表示し、演算子があったら括弧でくくってから左右の式の間にその演算子を書きます。$という記号がついた演算子は左右をひっくり返して表示します。

式の計算

calculate :: Equation -> Maybe Rational
calculate (Number x) = Just(x % 1)
calculate (Operator "+" a b) = (+) <$> calculate a <*> calculate b
calculate (Operator "*" a b) = (*) <$> calculate a <*> calculate b
calculate (Operator "-" a b) = (-) <$> calculate a <*> calculate b
calculate (Operator "$-" a b) = (-) <$> calculate b <*> calculate a
calculate (Operator "/" a b) 
    |calculate b /= Just 0 = (/) <$> calculate a <*> calculate b
    |otherwise = Nothing
calculate (Operator "$/" a b) 
    |calculate a /= Just 0 = (/) <$> calculate b <*> calculate a
    |otherwise = Nothing

Rationalは割り算を表す型です。数字だったらそのままにして演算子があったら左右にその演算子を適用し、$があったらひっくり返しています。割り算の失敗を表すためにMaybeモナドを用います。アプリカティブファンクターとして書いているところが個人的に綺麗で気に入っています。

式の列挙

makeEq :: [Integer] ->[Equation]
makeEq [] = []
makeEq [x] = [Number x]
makeEq (x:xs) = do
    let n = 1 + length xs
    i <-[0..(n-2)]
    ope <- ["+","*","-","$-","/","$/"]
    (ps,qs)<- combi i xs
    a <- makeEq (x:ps)
    b <- makeEq (qs)
    return (Operator ope a b)

リストモナドを用いて複数通りを列挙しています。リストの内包表記と本質的に同じものなのでリストモナドが良くわからなくても内包表記から類推すればなんとなくわかると思います。与えられた数字の一番最初は正規形の制約より必ず左側に渡します。任意に左側の数字の個数を選び、演算子を選び、先ほど決めた個数数字をチョイスして、左右に配置しています。チョイスする関数の実装は次のようなものです。

combi :: Int -> [a] -> [([a],[a])]
combi 0 xs = [([],xs)]
combi n (x:xs) 
    |length xs + 1 == n = [(x:xs,[])]
    |otherwise = [(x:ps,qs)|(ps,qs)<-combi (n-1) xs] 
        ++ [(ps,x:qs)|(ps,qs)<- combi n xs]

個数とリストを受け取って、その個数分取り出したものとその残りを全通り列挙しています。素直に再帰的に定義できていると思います。このとき順序を保っているのは重要です。常に使用する数字は最初の入力と同じ順序で並んでいるため先頭の数字を左に配置すれば正規形になるのです。

インターフェース

もう論理部分は完成しています。あとは入出力だけです。

main = forever $ do
    putStrLn "What Integer do you want to make?"
    i <- readLn ::IO Integer
    putStrLn "Input numbers like [1,2,3,4]"
    numbers <- readLn :: IO [Integer]
    let ans = headMaybe . (filter (\x -> calculate x == Just (i%1))) $ makeEq numbers
    putStrLn (show i ++" = " ++ showMaybe ans)
    putStrLn ""

実際に探索しているのはletの行だけです。headMaybeは長さゼロのリストに対してリストの先頭を取る関数を対応させたもの、showMaybeはMaybeモナドに文字列に変換する関数を対応させたものです。

headMaybe::[a] -> Maybe a
headMaybe [] = Nothing
headMaybe xs = Just (head xs) 

showMaybe ::(Show a)=> Maybe a -> String
showMaybe (Just a) = show a
showMaybe Nothing = "No Answer"

これで完成です!

テスト

難問の[1,1,5,8]を解かせてみましょう(答えを知らない人はネタバレに注意)

What Integer do you want to make?
10
Write numbers you should use like [1,2,3,4]
[1,1,5,8]
10 = (8 / (1 - (1 / 5)))

やりました!
速度としては遅延評価により答えが見つかればすぐに表示されるため、目的の数字が小さめですぐ答えが見つかるものは使用する数字が多くても大丈夫です。もしも作れない場合、使用する数字が6個を超えてくると全通り試してNoAnswerを表示するまでに非常に時間がかかってしまいます。組み合わせは爆発するものなので仕方ないですね。

感想、改善点

非常に直感的に書くことができたので非常に満足感がありました。計算式が本質的に二分木であることがここまで明示的に書けるとは想像以上に気持ち良いです。
改善点としてはまだ無駄な列挙をしている点です。入力に同じ数字がある場合は非常に無駄な探索をしてしまいます。また、結合則を無視しているため、(1+2)+3と1+(2+3)どちらも調べています。まあ速さを求めるならC++などで動的計画法を活用して書くべきである気がするのでそこまで気にしないでおいておきます。
コードの改善点がある場合ぜひ教えてください。
ソースとexeファイルを上げておきます。
findEquation - Google ドライブ

*1:[3,4,7,8] [1,1,5,8]などが難しいことで有名です

Haskellでフィボナッチ

去年の六月ごろに「すごいHaskell楽しく学ぼう」(以下すごいH)を買って読んでみようと思ったのですが、当時はほとんどプログラミングをしたことがなく、良くわからないまま終わりました。8か月ほど経って経験を積んだのでようやく本棚でほこりをかぶっていた本を読み始めてみました。
とりあえず再帰フィボナッチ数列を求めてみるコードを書いてみます

fibo1 ::(Num a)=> Int -> a
fibo1 0 = 1
fibo1 1 = 1
fibo1 n = fibo1 (n-1) + fibo1 (n-2)

いやぁ気持ちが良いですね
なんか漸化式をそのまま書くと動くので嬉しいです

このままだとn番目のフィボナッチ数を求めるオーダーが指数オーダーなので動的計画法っぽいことをして、O(n)で計算させてみます。

applyN :: Int ->(a->a) ->a ->a
applyN 0 _ a = a
applyN n f a = f $ applyN (n-1) f a

fibo2 :: Int -> Integer
fibo2 n = fst $ applyN n (\(a,b) -> (a+b, a)) (1,0)

ラムダ計算で書くとなんかそれっぽくてカッコいいですね
それにバグっていないという安心感がありますね
これが副作用がないということなんでしょうか
折角なので行列を繰り返し二乗法で計算してO(log n)で計算するやつをやってみました

data Mat a = Mat a a a a deriving (Show)

mmul:: (Num a)=>Mat a -> Mat a -> Mat a 
Mat x y z w `mmul` Mat i j k l 
    = Mat (x*i+y*k) (x*j+y*l) (z*i+w*k) (z*j+z*l)

getnum :: Mat a -> a
getnum (Mat x _ _ _) = x

fiboMat = Mat 1 1 1 0

appself ::(a-> a->a) -> a -> a
appself f x = f x x

foldn:: Int -> (a -> a-> a)-> a-> a
foldn 1 _ a = a
foldn n f a 
    |even n = appself f (foldn (n `div` 2) f a )
    |otherwise = f a $ appself f (foldn (div n 2) f a )

fibo3 :: Int -> Integer
fibo3 n = getnum (foldn n mmul fiboMat)

関数の名前を付けるセンスがないですねぇ。
いや非常にテンションが上がりますね。行列を定義するのが非常に自然です。
これめちゃくちゃ自然に書けますねぇ
appselfという関数が実は重要で、これを使わないと再帰展開が分岐してしまってO(n)になってしまいます(しまいました)
すごいH本を読み進めていきましょう