知ってると便利なSTL（12） reverse_iterator

基本

逆方向に進むiterator。
rbeginはend - 1と同じ場所を指す。rendはbegin - 1と同じ場所を指す。

上の説明を理解するための簡単なサンプル。

#include <iostream>
#include <vector>

using namespace std;

int main(int argc, char **argv) {
    
    vector<int> v{1,2,3,4,5,6,7,8};
    for (auto itr = v.rbegin(); itr != v.rend(); itr++) {
        cout << *itr << endl;
    }

    return 0;
}

実行すると、
8
7
6
5
4
3
2
1
と表示される。v.rbegin()はv.end()-1の場所を指しており、v.rend()はv.begin()-1の場所を指していることが分かる。また、reverse_iteratorは名前の通り逆向きに進むことが分かる。

これだけだと何が嬉しいか分からないので、以下に便利な使い方を示す。

逆順ソート

vectorを逆順にソートしたいときに、reverse_iteratorを使うと短く書ける。

#include <iostream>
#include <vector>
#include <algorithm>

using namespace std;

int main(int argc, char **argv) {
    
    vector<int> v{1,3,5,7,2,4,6,8};

    // functorを使って書くこともできるが、
    sort(v.begin(), v.end(), greater<int>());
    
    // reverse_iteratorを使うと、短く書ける
    sort(v.rbegin(), v.rend());

    return 0;
}

setの最小値、最大値

beginで最小値、rbeginで最大値を取得出来る。

#include <iostream>
#include <set>

using namespace std;

int main(int argc, char **argv) {
    
    set<int> s{3,2,1,4,8,7,6,5};

    cout << "min = " << *s.begin() << endl;
    cout << "max = " << *s.rbegin() << endl;

    return 0;
}

エラトステネスの定数倍高速化術

　Project Eulerのとある問題でエラトステネスの篩を高速化している人がいた。
なかなか興味深かったので、コードを読みといて自分なりにアレンジしてみた。

ローカル環境で約2倍の高速化に成功した。1e8までのふるいが0.7s程度で実行出来た。

普通のふるい

#include <iostream>
#include <algorithm>

using namespace std;

const int N = 1e8;
bool pr[N];

int main(int argc, char **argv) {
    
    fill(pr + 2, pr + N, true);

    for (int i = 0; i * i < N; i++) {
        if (!pr[i])
            continue;

        for (int j = i * i; j < N; j += i)
            pr[j] = false;
    }

    int x = 0;
    int y = 0;
    for (int i = 0; i < N; i++)
        if (pr[i])
            x ^= i, y = 23 * y + i;
    
    cout << x << endl;
    cout << y << endl;

    return 0;
}

定数倍高速版ふるい

アイディアとしては、

1. 2以外の偶数は素数ではないため、ふるいの配列は奇数に対してのみ用意する。
2. 素数か合成数かのフラグをboolではなくintのビットで管理する。

というもの。

1.で約2倍高速になった。これは直感的にも理解しやすい。一番内側の処理の実行回数が半分程度になるため2倍程度速くなりそうである。

1.の後に2.を施すことで、さらに1.4倍ほど高速になった。32ビット毎に何らかの処理を纏めたというわけではなく、単純にフラグの管理をビットで行うようにしただけであるが意外と速くなった。使用するメモリがコンパクトになったためキャッシュに乗りやすくなったということだろうか。

#include <iostream>
#include <algorithm>

using namespace std;

const int N = 1e8;
const int BITS = 32;

#define BACKET(n) (n / 2 / BITS)
#define MASK(n) (1 << n / 2 % BITS)

int pr[BACKET(N) + 1];    // sieve arrray containing only odd numbers  [1,3,5,7,,9,....]

int main(int argc, char **argv) {
    
    fill(pr, pr + BACKET(N) + 1, -1);
    pr[0] = -2;   // 1 is not prime

    for (int i = 1; i * i < N; i += 2) {
        if (!(pr[BACKET(i)] & MASK(i)))
            continue;
        for (int j = i * i; j < N; j += 2 * i)
            pr[BACKET(j)] &= ~MASK(j);
    }

    int x = 0;
    int y = 0;
    for (int i = 0; i < N; i++)
        if (i == 2 || i % 2 && (pr[BACKET(i)] & MASK(i)))
            x ^= i, y = 23 * y + i;

    cout << x << endl;
    cout << y << endl;

    return 0;
}

Python Idioms (4) モジュールに別名をつける

import句にasを付けるとモジュールに別名をつけることが出来る。
同様に、from import句にasを付けるとモジュール内の関数に別名を付けることが出来る。

これを使って、
重複なし組み合わせ（nCk）と重複組み合わせ（nHk）を計算する関数を書いてみた。

from math import factorial as f

def c(n, k):
    return f(n) / f(k) / f(n - k)

def h(n, k):
    return c(n - 1 + k, k)

おいら、レベル4に上がったっぺ

大きな数の最初の10桁を求める

大きな数の最初の10桁を効率的に計算しないといけないような問題に遭遇した。

下n桁を求めるのは剰余で出来るけど、上n桁ってどうやるのだろう。
logを使うとうまく出来そうだ。ということでやってみた。
とりあえず2のべき乗の上10桁を求めてみた。

プログラム

#include <iostream>
#include <cmath>

using namespace std;

const double EPS = 1e-6;

// return the first 10 digits of 2^k
long long calc(int k) {
    double x = k * log(2);
    int d = x / log(10);
    x -= d * log(10);
    x += 9 * log(10);
    return exp(x) + EPS;
}


int main(int argc, char **argv) {
    
    cout << calc(1000) << endl;
    cout << calc(10000) << endl;
    cout << calc(100000) << endl;

    return 0;
}

結果確認

プログラムの実行結果。

1071508607
1995063116
9990020930

Pythonで2^kを実際に計算して検算してみた。正しく計算出来ていた。

$ python -c "print str(2**1000)[:10]"
1071508607
$ python -c "print str(2**10000)[:10]"
1995063116
$ python -c "print str(2**100000)[:10]"
9990020930

Graphvizで二分探索木を描画する

Graphvizという便利なツールがあることを知った。
AT&T研究所が開発したオープンソースのグラフ描画ツールである。とりあえず動かして遊んでみようということで、二分探索木を描画してみた。

インストール

Ubuntuマシンにインストール。

$ sudo apt-get install graphviz

二分木を作るプログラム

まず、20個の乱数を二分探索木に格納する。
その後、Graphvizに読み込ませるDOT言語のスクリプトをdump関数で出力する。

#include <iostream>
#include <cstdlib>

using namespace std;

struct node {
    int x;
    node *lch, *rch;

    node (int x) : x(x), lch(NULL), rch(NULL) {}
};

node *insert(node *v, int x) {
    if (v == NULL)
        return new node(x);
    if (x < v->x) 
        v->lch = insert(v->lch, x);
    if (x > v->x) 
        v->rch = insert(v->rch, x);
    return v;
}

void dump(node *v, node *par = NULL) {
    if (v == NULL)
        return;

    if (par == NULL) {
        cout << "digraph G{" << endl;
        cout << "  graph [ordering=\"out\"];" << endl;
    }
    else {
        cout << "  " << par->x << " -> " << v->x;
        if (v->x > par->x)
            cout << " [style = dotted]";
        cout << ";" << endl;
    }

    dump(v->lch, v);
    dump(v->rch, v);

    if (par == NULL)
        cout << "}" << endl;
}

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

    srand(time(NULL));
    node *root = NULL;

    for (int i = 0; i < 20; i++)
        root = insert(root, rand() % 100);

    dump(root);

    return 0;
}

上のソースはbst.cppという名前で保存しているとする。

以下を実行。

$ g++ bst.cpp -o bst
$ ./bst >sample.dot
$ dot -Tpng sample.dot -o sample.png
$ xdg-open sample.png 

すると、以下の画像が表示される。実戦が左の子への辺、破線が右の子への辺である。

これで木の回転をして遊んだり、平衡二分木を自分で実装したり、などなどするときに動作確認がしやすくなる。

商社は何をやっているのか

　商社は何をやっているのか？どうやって儲けているのか？
半年くらい働いてみて分かったことをまとめておく。
商社の規模や扱っている商品によって違いはあるかもしれないが、基本的な部分はだいたい共通していると思う。

販売額 - 仕入れ額のマージンが儲け

商社の利益は、販売額と仕入れ額の差から生まれる。メーカーから商品を仕入れて、それに値段を上乗せして、消費者や小売業者に販売する。このとき上乗せした値段が利益になる。

　ここまでは前から知っていた。いまいち分かっていなかったのは商社の存在意義。メーカーが直接売ればいいじゃん。なんか商社って楽して儲けてるだけで何の価値も創造していないような・・と思っていた。

いいモノを作れば売れるというわけではない

商社の存在意義を考えるうえでポイントとなるのが、「いいモノを作れば必ず売れる」というのは間違いということ。いいモノを作っても宣伝しなければ売れない。当たり前だ。

　そこで登場するのが商社。商社は持ち前の営業力でモノを上手に宣伝し、顧客に売り込む。もちろん現代においては、ITを駆使した顧客分析、商品分析なども行う。今はこういう商品が人気だから、今度はこういう商品を開発しましょう！というように商品企画に関与することもある。

　つまり商社はメーカーが作ったいいモノを世の中により広く流通させるために営業力の面からサポートを行う仲介業者みたいなイメージ。これが分かったとき、商社は楽して儲けているだけというイメージは完全に無くなった。

モノは買い手にすぐ届くわけではない

モノを買いたい人が見つかったらそれで終わりというわけではない。
モノを顧客の元に届けなければならない。そのためには配送ルートが必要だ。それからモノを格納しておく倉庫も必要だ。何をいつどれだけ作ればいいかを考える生産計画も必要だ。

　こういう物流まわりも商社がやっている。メーカーが作ったモノを営業して売り込み、顧客のもとに届けるまで面倒をみる。これが商社。商品の総合プロデューサーという感じだろうか。

Optimal Strategy for Grundy's Game

 I learned about "Grundy's Game." This game is played by two players as follows:

1. The game starts with a single heap of objects.
2. The first player chooses a heap and splits the heap into two heaps of different sizes.
3. The second player does the same thing as 2.
4. Repeat 2. and 3. The player who cannot make a proper move loses.

For example, let's say the size of an initial heap is 9.

The first player splits a heap into (4, 5).
The second player splits the second heap into (2, 3), making the heap sizes (4, 2, 3).
The first player splits the first heap into (2, 2), making the heap sizes (2, 2, 2, 3).
The second player splits the last heap into (1, 2), which is the only allowed move here, making the heap sizes  (2, 2, 2, 1, 2).

Now there are no proper moves available.
You cannot split a heap of the size 1.
And You cannot split a heap of the size 2 into two heaps of different sizes.
Therefore the second players wins in the example.

Umm, sounds like a good way to kill time. You can play with some coins.

I wrote a program that returns the optimal strategy for this game. As you can guess from the name of this game, I used Grundy numbers. Just run the program, input the heap sizes, like, 10 or 5 10 20 or whatever. That's it.

#include <iostream>
#include <vector>
#include <set>
#include <sstream>

using namespace std;

const int MAX_V = 1000;
int g[MAX_V + 1];

void init() {
    for (int i = 2; i <= MAX_V; i++) {
        set<int> s;
        for (int j = 1; 2 * j < i; j++) {
            s.insert(g[j] ^ g[i-j]);
        }
        while (s.count(g[i]))
            ++g[i];
    }
}

void solve(vector<int> &a) {

    int sum = 0;
    for (int i = 0; i < a.size(); i++)
        sum ^= g[a[i]];

    for (int i = 0; i < a.size(); i++) {
        sum ^= g[a[i]];
        
        for (int j = 1; 2 * j < a[i]; j++) {
            if (!(sum ^ g[j] ^ g[a[i] - j])) {

                for (int k = 0; k < a.size(); k++) {
                    if (k == i)
                        cout << j << " " << (a[i] - j) << " ";
                    else
                        cout << a[k] << " ";
                }
                cout << endl;

                return;
            }
        }
        
        sum ^= g[a[i]];
    }
    
    cout << "Uncle." << endl;
}

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

    init();
    
    for (string line; getline(cin, line), line != "bye"; ) {

        istringstream iss(line);
        vector<int> a;

        int n;
        while (iss >> n)
            a.push_back(n);

        solve(a);

    }

    return 0;
}

スターリングの公式を使って、sinのテイラー展開近似の項数を見積る

　ζ(2) = π² / 6の証明を読んでいて、ふとsin(x)のテイラー展開について気になった。0の周りでテイラー展開するとして、例えば x = 100くらいでもそれなりの精度で近似するためにはどれくらいの項数がいるのだろうか。
sin(x)は周期関数なので、[0, 2π)に入るようにxを変換すれば、x = 2πのときにどれだけ項数が必要か考えれば十分だと思うけど、周期性は利用しないという前提で。

x < 1であれば、x^k → 0になるので高次の項は無視できるというのはすぐに分かるけど、x > 1のときって高次の項は無視できるのだろうか？

x^k/k! を考えると、kが大きくなるにつれ、分子の増加よりも分母の増加の方が速くなるため、x > 1のときも十分に大きなkを取れば、x^k/k! → 0 になりそうである。じゃあ十分に大きなkってどれくらいだろうか？スターリングの公式が使えそうな気がする。

k! ~ (k/e)^kなので、
x^k/k! ~ x^k/((k/e)^k) = (ex/k)^k
である。

これを使って、0 <= x < 100にてsin(x)をテイラー展開近似するのに必要な項数を見積もってみる。
大事なことを書き忘れていた。sin関数のテイラー展開は以下のようになる。

sin(x) = 1 - x³/3! + x⁵/5! - x⁷/7! + x⁹/9! - ....

x = 100として、x^k/k! ~ (ex/k)^k < 10^-6となるようなkを概算してみる。
k = 300のとき、(ex/k)^k = 1.42 * 10^-13なので、300/2 = 150項あれば十分なのではないかと予想してみる。

以下に、[0, 100]でのsin関数とそのテイラー展開近似（150項まで）をグラフにしてみた。十分な近似ができているようだ。

1+2+3+4+5+ ..... = -1/12

最近摩訶不思議な式に出くわした。

1 + 2 + 3 + 4 + ..... = -1/12らしい。

詳しいことは分からないが、様々な物理の分野で使用される式らしい。

証明

S = 1 + 2 + 3 + 4 + ......
S₁ = 1 - 1 + 1 - 1 + 1 ....
S₂ = 1 - 2 + 3 - 4 + 5 ....

と定義する。

まず、
1 - S₁ = 1 - (1 - 1 + 1 - 1 ....) = S₁より、 2S₁ = 1となる。
よってS₁ = 1/2が成り立つ。

次に、
S₂ = 1 - 2 + 3 - 4 + 5 ....
S₂ =      1 - 2 + 3 - 4 + 5 ....

として（下段のS₂が1項右にずれているのがポイント）、両辺を足すと、
2S₂ = 1 - 1 + 1 - 1 + 1 .... = S₁より、
2S₂ = 1/2となる。よって、S₂ = 1/4が導かれる。

最後に、
S = 1 + 2 + 3 + 4 + .....
S₂ = 1 - 2 + 3 - 4 + .....
として上段から下段の式を引くと、

S - S₂ = 4 + 8 + 12 + 16 + ... = 4Sより、3S = - S_{2となる。
よって、S = -1/12である。}

参考URL

[1] ASTOUNDING: 1 + 2 + 3 + 4 + 5 + ... = -1/12

ペル方程式を解く

　Pythonでペル方程式を解くプログラムを書きました。
ペル方程式とはディオファントス方程式の一つで以下のような形をしているものです。

解き方

ペル方程式を変形すると、

となります。よって、Dの平方根の有理数近似をx/yとすると、（x, y）が解の候補となりそうです。
有理数近似には、連分数展開が使えそうです。
実際に、ペル方程式の解を(x, y)とすると、x/yは√Dのconvergentとなる[1]ことが証明できます。よってconvergentを列挙し、方程式を満たすかどうか調べれば不足なくペル方程式の解を求めることが出来ます。

上の計算手続きで方程式の解を求めることができますが、コンピュータを使って計算する場合は注意点があります。
浮動小数点演算を使って平方根の連分数展開を行うと誤差が大きくなってしまい、正しい結果が得られないということです。これを回避するために、整数演算のみで平方根の連分数展開を行う方法が知られています[2]。

とすると、[a0; a1, a2, a3, ...]が√Dの連分数展開となります。

プログラム

D=2の場合のペル方程式を解いています。解を10個列挙しています。

import itertools

# check if x^2 - D y^2 = k holds.
def check_pell_equation(D, k):
    def f(x, y):
        return x * x - D * y * y == k
    return f

# calculate the continued fraction of the square root of d.
def continued_fraction_of_sqrt(d):
    p, q = 0, 1
    x = int(d ** 0.5)

    while True:
        a = (x + p) / q
        yield a
        p = a * q - p
        q = (d - p * p) / q

# solve a Pell's equation x^2 - D y^2 = k, 
#  where D is a non-square number, k is +1 or -1.
def solve_pell_eqation(D, k):
    pell_eq = check_pell_equation(D, k)
    cf = continued_fraction_of_sqrt(D)
    
    p0, q0 = 0, 1
    p1, q1 = 1, 0

    while True:
        if pell_eq(p1, q1):
            yield p1, q1
        a = cf.next()
        p1, p0 = a * p1 + p0, p1
        q1, q0 = a * q1 + q0, q1

if __name__ == '__main__':
    ret = solve_pell_eqation(2, 1)
    
    for x, y in itertools.islice(ret, 10):
        print x, y

参考

[1] Solution of Pell's Equation is a Convergent
[2] c++ - Generating continued fractions for square roots - Stack Overflow
[3] MathForKids Pell's Equation (YouTube Video)

連分数展開でネイピア数を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;
}

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