7. 円周率の分数表現を探す¶
355/113が非常によい円周率の近似値であることはよく知られている。電卓を使って計算してやれば
\[355/113 = 3.14159292035\]
と驚くべき精度があることがわかる。しかもこの数字は分母から順に1,1,3,3,5,5と、最も小さな奇数から順に二個ずつ並んでいるのが憎らしい。これだけ綺麗な分数は見つかりそうもないが、他の分数表現を探してみようじゃないか。
7.1. 評価基準¶
単純に円周率に近い値が欲しければ例えば\(31415/10000\)のように\(10^n\)で割ってやればよいだけなのでつまらない。355/113が人を感動させるひとつの理由は分子、分母それぞれ3桁しかないのにそれを大きく超える精度で\(\pi\)と一致するからである。そこで、次の評価基準を採用する。
\[\text{err} = \frac{-\log|\pi-m/n|}{({\rm digit}(m)+{\rm digit}(n))/2}\]
ただし、digit関数は引数の桁を返す関数である。
7.2. 結果¶
次のプログラムを用いて100万未満を探索した。
#include <iostream>
#include <math.h>
#include <stdio.h>
#include <algorithm>
#include <set>
unsigned int digit(unsigned int n)
{
return (n==0 ? 1 : (int)log10(n)+1);
}
int gcd( int m, int n )
{
while(m!=n)
{
if (m > n) m = m - n;
else n = n - m;
}
return m;
}
class Frac
{
public:
Frac(unsigned int m, unsigned int n, double eval): m{m}, n{n}, eval{eval}
{
}
unsigned int m, n;
double eval;
bool operator>(const Frac &frac) const
{
return this->eval > frac.eval;
}
};
int main()
{
double C = M_PI;
int N = 1000000;
std::set<Frac, std::greater<Frac> > frac;
for(unsigned int n=1; n<N; n++)
{
unsigned int tmp_m = round(n * C);
unsigned int tmp_n = n;
//約分
int GCD = gcd(tmp_m,tmp_n);
tmp_m /= GCD;
tmp_n /= GCD;
double dig = (digit(tmp_n) + digit(tmp_m))/2.;
double eval = -log10(fabs((double)tmp_m/(double)tmp_n-C)) / dig;
frac.insert(Frac(tmp_m,tmp_n,eval));
}
for(auto &f: frac)
{
std::cout << f.m << " " << f.n << " " << f.eval << std::endl;
}
return 0;
}
結果発表!!
1位:355/113 評価値=2.19129
2位: 22/7 評価値=1.93206
3位: 312689/99532 評価値=1.91554
おおう!355/113すごいな!22/7も大健闘!円周率の連分数表現の最初の方が大健闘という結果。
このプログラムは円周率以外でも使える。例えばネイピア数。
49171/18089が最も評価値が高い。あと目立ったのは193/71。これもやはり連分数表現から出てくるらしい。