連分数展開でネイピア数を100桁まで求める

　連分数展開を使ってネイピア数を100桁まで求めてみました。

ネイピア数を連分数展開すると、以下のように規則的な列になります。

e = [2; 1, 2, 1, 1, 4, 1, 1, 6, 1, 1, 8, 1, .... ]

上式の連分数を分数に変換し、さらに分数を小数に変換するという手順で計算します。

以下の機能が標準でサポートされているPythonで書きました。

• 多倍長整数の乗算
• ジェネレータを利用した無限リストの生成
• 任意精度の小数演算

import itertools
from decimal import Decimal, getcontext

def pi_expansion():
    yield 2
    for i in itertools.count(2, 2):
        yield 1
        yield i
        yield 1

getcontext().prec = 100

cf = pi_expansion()

p0, q0 = 0, 1
p1, q1 = 1, 0

for _ in range(100):
    a = cf.next()
    p1, p0 = a * p1 + p0, p1
    q1, q0 = a * q1 + q0, q1
    
    print Decimal(p1) / Decimal(q1)

100番目のconvergentまで計算して、
e ~ 2.718281828459045235360287471352662497757247093699959574966967627724076630353547594571382178525166427となりました。90番目あたりのconvergentからこの値に収束していました。
ググってみると、正しいみたいです。

単語が辞書に存在するかどうかを調べる

　暗号化された英語の文章を復号するという問題に挑戦しています。
復号化された文章が正しい英単語から構成されていることを確認するために、単語が辞書に存在するかどうかをチェックする機能を実装しました。Ubuntuの場合は、/usr/share/dict/american-englishにアメリカ英語の単語リストが存在するので、それを使いました。 　

#include <iostream>
#include <fstream>
#include <set>

using namespace std;

class dictionary {
    set<string> container;

public:
    void load(const string &path) {
        ifstream fs(path);

        string word;
        while (fs >> word)
            container.insert(word);
    }

    bool contains(const string &word) {
        return container.count(word);
    }
};


int main(int argc, char **argv) {

    dictionary dict;
    dict.load("/usr/share/dict/american-english");

    for (string s; cin >> s; ) {
        if (dict.contains(s))
            cout << s << " is in the dictionary." << endl;
        else
            cout << s << " is not in the dictionary." << endl;
    }

    return 0;
}

以下実行結果です。固有名詞も辞書に含まれているみたいです。

hello
hello is in the dictionary.
world
world is in the dictionary.
soccer
soccer is in the dictionary.
Linux
Linux is in the dictionary.
Beatles
Beatles is in the dictionary.
Fibonacci
Fibonacci is in the dictionary.
totient
totient is not in the dictionary.
Yankees
Yankees is in the dictionary.

Python Idioms (3) 文字列をreverseする

　リストにはreverseメソッドがありますが、文字列にはreverseメソッドがありません。
文字列をreverseしたい場合は、以下のようにスライスを使うのがPythonicなやり方らしいです。

>>> s = "hello, world"
>>> s[::-1]
'dlrow ,olleh'

のようにs[::-1]でsをreverseした文字列を取得出来ます。

これを使って整数nをreverseする関数を書いてみました。

>>> def revInt(n):
...   return int(str(n)[::-1])
... 
>>> 
>>> revInt(100)
1
>>> revInt(1234)
4321

エンジニア日記（3）1000倍速くなった

　某システムのとある機能が大変に遅かった。
ユーザーからも「遅い。遅い。」と苦情が多発していたらしい。

　中身を見てみた。たしかにこれは遅い。HashSetを使えば高速に処理出来るのだが、ArrayListを使っている。

　書き直してところ、ボトルネックとなっていた処理が1000倍速くなった。

　こういうのって結構あるんじゃないだろうか。僕はred coderでもないし、数学オリンピックの代表選手でもない。しかし、そんな僕にも改良できることが、世の中にはたくさんあるんじゃないだろうか。

エラトステネスの篩の計算量

　蟻本には、エラトステネスの篩の計算量はO(n log log n) と書いてあります。
log log nってなんぞや？どこから出てきたの？とずっと疑問でしたが、何故こうなるか分かったのでメモしておきます。

エラトステネスの篩のコード

C++で書く場合は、概ね以下のようなコードが一般的かと思います。

bool isPrime[N];

int main() {
    for (int i = 2; i < N; i++)
        isPrime[i] = true;

    for (int i = 2; i < N; i++) {
        if (!isPrime[i]) continue;
        for (int j = 2 * i; j < N; j += i)
            isPrime[j] = false;
    }
}

計算量の見積もり

まず簡単のため、上のソースの”素数でないなら飛ばす”という処理（※）が無かった場合を考えてみます。
（※）if (!isPrime[i] ... ) の行のこと。

二重ループのところが一番重い処理で、計算量は、

∑_{ i < N} N / i = N * ∑_{i < N} 1 / i = N ln (N)

となります。ここまでは簡単です。

次に、”素数でないなら飛ばす”がある場合を考えてみます。
計算量は、

∑_{ p < N, p is prime} N / p = N * ∑_{p < N, p is prime} 1/p

となりますが、∑_{p < N, p is prime} 1/pがもしやln ln(N)になるのでは？と予想をつけてみたのです。実験してみたところ、以下のようになりました。

N	∑_{p < N, p is prime} 1/p	ln ln N
1e6	2.887	2.626
1e7	3.041	2.780
1e8	3.175	2.913

両者の値はnearly equalであることが分かりました。ということで、
∑_{p < N, p is prime} 1/p = ln ln(N)
が漸近的に成り立つのではないか？そうであれば、エラトステネスの篩の計算量がO(N log log (N))というのも納得できるのだが。。
と思い調べていると、以下の情報を見つけました。

Euler then concluded
1/2 + 1/3 + 1/5 + 1/7 + 1/11 ..... = ln ln (+∞)

It is almost certain that Euler meant that the sum of the reciprocals of the primes less than n is asymptotic to ln(ln(n)) as n approaches infinity. It turns out this is indeed the case.

-- Divergence of the sum of the reciprocals of the primes - Wikipedia, the free encyclopedia

ということで、予想は正しかったみたいです。

Python Idioms (2) 副作用なしでリスト要素を順番に処理する

　リストの要素をソート順に処理する、もしくは、逆順に処理することを考える。

　リストにはsortメソッド、reverseメソッドがあるのでこれらを使えばいいのだが、これらのメソッドは元々のリスト自体を書き換えてしまう。

　元々のリストに副作用を与えずに処理したい場合はどうするか？そのようなときは、以下のようにsorted関数、reversed関数を使えばよい。

>>> employees = ['Charlie', 'David', 'Alice', 'Bob']
>>> 
>>> # リストの要素をソート順に処理する
... for emp in sorted(employees):
...     print emp
... 
Alice
Bob
Charlie
David
>>> 
>>> # もともとのリストがソートされたわけではない
... print employees
['Charlie', 'David', 'Alice', 'Bob']
>>> 
>>> # リストの要素を後ろから処理する
... for emp in reversed(employees):
...     print emp
... 
Bob
Alice
David
Charlie
>>> 
>>> # 同じくもともとのリストが逆順に並び替えられたわけではない
... print employees
['Charlie', 'David', 'Alice', 'Bob']
>>> 

Python Idioms (1) 数字の桁の合計を求める

123456789という数字の桁の合計は、
1 + 2 + 3 + 4 + 5 + 6 + 7 + 8 + 9 = 45.

これをPythonで。

>>> x = 123456789
>>> sum(map(int, str(x)))
45

Euler's totient functionとMobius functionの関係

最近見た式。

ぱっと見よく分からなかったので、ちょっと考えてみた。

例えばφ(10)について考えてみる。
これは、10と互いに素な10以下の正の整数の個数だ。

10以下の正の整数は10個ある。10 = 2 * 5なので、ここから2の倍数の個数10/2と、5の倍数の個数10/5を引く。

10 - 5 - 2 = 3となる。

ただし、10は2の倍数としても引かれているし、5の倍数としても引かれているので、1回余分に引かれていることになる。

よって
φ(10) = 10 - 5 - 2 + 1 = 4となる。

次に、n = 36を考えてみる。
φ(36)
= 36 - (36以下の2の倍数の個数) - (36以下の3の倍数の個数) + (36以下の6の倍数の個数)
= 36 - 36/2 - 36/3 + 36/6
= 12

となる。

36 = 2² * 3²と素因数分解できるので、4の倍数とか9の倍数とかについても考えたくなるが、これらは2の倍数、3の倍数ですでに数えられているので考えなくてよい。

つまり、平方数で割れる約数については考えなくてよい。
あとは、重複している部分を数えるために包除原理を使って、足したり引いたりすればよい。
この２つが感覚的に分かれば、最初に示した式はごく自然に見える。

一般に、
n = p1^k1 * p2^p2 * p3^k3 * ......
として

φ(n)
= n - n/p1  - n/p2 - n/p3 ..... + n / (p1 * p2) + n / (p1 * p3) + .... - n / (p1 * p2 * p3) .....
= ∑_{d | n} μ(d) n / d

となる。

(追記)
メビウスの反転公式を使うと、上式は


となる。同じく最近見かけて何故成り立つか分からなかった式だけど、ここでこう繋がったか。。

ピタゴラス数の列挙

a^2 + b^2 = c^2
を満たす正の整数(a, b, c)をPythagorean tripleと言う。

Pythagorean tripeを列挙するには、以下のユークリッドの公式を使うといい。

すべてのPythagorean tripleは正の整数n, m, kを用いて、

a = m^2 - n^2
b = 2 m n
c = m^2 + n^2

のように一意に表すことが出来る。ただし、m, nは
m > n,
(m, n) = 1,
m != n (mod 2)
を満たす。

証明はここのProof of Euclid's formulaが分かりやすい。

これを使って以下の問題を解いてみる。

三辺の長さの合計が1,000,000以下となるような直角三角系は何個存在するか？
ただし、すべての辺の長さは整数である。

a + b + c <= 1,000,000を満たすPythagorean triple (a, b, c)の個数を求めればよい。以下のようなプログラムを書いてみた。 808950個。あってるかな？

#include <iostream>

using namespace std;

const int PERIMETER = 1e6;

long long gcd(long long a, long long b) {
    if (b == 0)
        return a;
    return gcd(b, a%b);    
}

int main(int argc, char **argv) {
    
    long long ret = 0;
    
    for (int m = 1; m * m < PERIMETER; m++) {
        for (int n = 1; n < m && n * n + m * m < PERIMETER; n++) {
            if (m % 2 == n % 2) continue;
            if (gcd(m, n) > 1) continue;

            int a = m * m - n * n;
            int b = 2 * m * n;
            int c = m * m + n * n;

            ret += PERIMETER / (a + b + c);
        }
    }

    cout << ret << endl;

    return 0;
}

おまけ：ピタゴラスの定理の証明が面白かったので、リンクを張っておく。