## 2014年12月25日木曜日

### 知ってると便利なSTL（12） reverse_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;
}
```

## 2014年12月21日日曜日

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

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 = -2;   // 1 is not prime

for (int i = 1; i * i < N; i += 2) {
continue;
for (int j = i * i; j < N; j += 2 * i)
}

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;
}
```

## 2014年12月7日日曜日

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

import句にasを付けるとモジュールに別名をつけることが出来る。

これを使って、

```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)
```

## 2014年11月30日日曜日

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

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
```

## 2014年11月29日土曜日

### 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;
}
```

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

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

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

## 2014年11月27日木曜日

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

商社は何をやっているのか？どうやって儲けているのか？

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

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

そこで登場するのが商社。商社は持ち前の営業力でモノを上手に宣伝し、顧客に売り込む。もちろん現代においては、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;
}
```

## 2014年10月31日金曜日

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

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

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

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

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

これを使って、0 <= x < 100にてsin(x)をテイラー展開近似するのに必要な項数を見積もってみる。

sin(x) = 1 - x3/3! + x5/5! - x7/7! + x9/9! - ....

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

## 2014年10月28日火曜日

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

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

S = 1 + 2 + 3 + 4 + ......
S1 = 1 - 1 + 1 - 1 + 1 ....
S2 = 1 - 2 + 3 - 4 + 5 ....

と定義する。

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

S2 = 1 - 2 + 3 - 4 + 5 ....
S2 =      1 - 2 + 3 - 4 + 5 ....

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

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

S - S2 = 4 + 8 + 12 + 16 + ... = 4Sより、3S = - S2となる。
よって、S = -1/12である。

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

## 2014年10月5日日曜日

### ペル方程式を解く

Pythonでペル方程式を解くプログラムを書きました。
ペル方程式とはディオファントス方程式の一つで以下のような形をしているものです。 ペル方程式を変形すると、 となります。よって、Dの平方根の有理数近似をx/yとすると、（x, y）が解の候補となりそうです。 とすると、[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
```

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

## 2014年9月30日火曜日

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

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

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

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

• 多倍長整数の乗算
• ジェネレータを利用した無限リストの生成
• 任意精度の小数演算
```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からこの値に収束していました。
ググってみると、正しいみたいです。

## 2014年9月23日火曜日

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

暗号化された英語の文章を復号するという問題に挑戦しています。

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

using namespace std;

class dictionary {
set<string> container;

public:
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;

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メソッドがありません。

```>>> 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
```

## 2014年9月20日土曜日

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

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

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

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

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

## 2014年9月14日日曜日

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

蟻本には、エラトステネスの篩の計算量は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

∑_{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

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

## 2014年9月13日土曜日

### 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']
>>>
```

## 2014年9月12日金曜日

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

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

これをPythonで。
```>>> x = 123456789
>>> sum(map(int, str(x)))
45
```

## 2014年9月4日木曜日

### Euler's totient functionとMobius functionの関係 ぱっと見よく分からなかったので、ちょっと考えてみた。

これは、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となる。

φ(36)
= 36 - (36以下の2の倍数の個数) - (36以下の3の倍数の個数) + (36以下の6の倍数の個数)
= 36 - 36/2 - 36/3 + 36/6
= 12

となる。

36 = 22 * 32と素因数分解できるので、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

となる。

(追記)
メビウスの反転公式を使うと、上式は となる。同じく最近見かけて何故成り立つか分からなかった式だけど、ここでこう繋がったか。。

## 2014年9月1日月曜日

### ピタゴラス数の列挙

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)
を満たす。

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

ただし、すべての辺の長さは整数である。

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;
}
```

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