無理数を有理数で近似

無理数xを、有理数N/Dで近似するというのは、よくあるテーマ。

現実的には、適当に有理数の分子、分母に上限を付ける。例えば、1 <= N, D <= Lのような範囲のN, Dに対して、|x - N/D|が最小となるようなN, Dを求める問題を考えることになる。

何も考えないで実装すると、N, Dを全探索するはめになり、O(L^2)の擬多項式時間計算量が必要となる。アルゴリズムをやったことがある人なら、Dを全部まわして、Nを二分探索するという方法を思いつくだろう。これだと、O(L * log L)で L = 1,000,000程度でも一瞬で答えが出る。

実は二分探索よりも、もっと面白い方法がある。

N/D > x なら 探索点を N/(D+1)に移動、そうでなければ探索点を(N+1)/Dに移動。これを逐次繰り返していき、最良点を保持すれば答えが出る。

以下、サンプルソース。無理数じゃないけど、x = 100007/100009として解いてみる。

#define MAX_VAL 1000000

int main() {
   int N, D;
   double x = 0.99998000179983804;

   double s, t;
   s = t = 1.0;

   while (s <= MAX_VAL && t <= MAX_VAL) {
       if (fabs(x - s/t) < fabs(x - 1.*N/D)) {
           N = s;
           D = t;
       }
       if (s / t > x)
           ++t;
       else
           ++s;
   }

   printf("%.18lf = %d/%d\n", x, N, D);
}

見事、元の分数の分子、分母を再現できる。

厳密な証明ではないが、数学的帰納法的なイメージで、上記の方法で適切なN, Dが求まることが分かる。（ざっくりした証明的なものです。。数学得意な人がいたら、厳密な証明方法教えてください。）

|x - N/D|が最小となるN, DをそれぞれN*, D*とおく。

i) (N* - 1, D*)、または、(N*, D* -1)が探索されれば、最適解にたどりつく。

ii) (N* - 1, D*) を (N~, D~) と改めて置きなおす。(N~-1, D~)、または、(N~, D~-1)が探索されれば、(N~, D~)に辿り着く。(N*, D*-1)に対しても同様。

iii)N, Dは1以上の整数なので、点列(D, D-1, D-2, ...)、および、(N, N-1, N-2, ..)は有限回の試行で、それぞれ1に収束する。

i)、ii)、iii)、より(N, D) = (1, 1)からはじめて、上記の方法を繰り返せば、最適解(N*, D*)に辿りつく。□

辿り着いた後も、N, D < Lの間は探索は続くので、最良解はメモリ上に保持しておかないといけないので注意。