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。これもやはり連分数表現から出てくるらしい。