浮動小数点数を分数に変換するアルゴリズムについて

はじめに

浮動小数点数を分数に変換したいと思った時，すぐによい資料に辿りつけなかったので，
自分なりにまとめ直してみる．
簡潔な答えは他の資料に譲るとして，結論ラストのぐたぐだ長々と書いていく．

問題としては全然難しくはないが，単純なアルゴリズムであっても何も考えずにプログラムを書くと，
誤差のせいで思ったような動作にならないことがある．
なので，そういう所にも少しだけ触れたい．

ただ，厳密なことを言い始めると中々面倒なので，ヒントや考える参考程度にして欲しい．

また，C言語のコードが載っているが，入力値に対する精査がなく，例えば負の数を渡すと予期した動作にならないので注意が必要である．
実際にはもっときちんとガードした方がよい．

これらのコードは以下の環境でコンパイルした．

ArchLinux
A10-7850K(Kaveri) AMD64
gcc==10.2.0

ちなみに，math.hのmodfを使っているが私の環境だと-lmオプションがなくてもコンパイルできた．

定義

始めに問題や文字の定義をする．
ある浮動小数点$f>0$があったとする．(符号は後でどうとでもできるため)
この時，$f>1$であってもよい．
また，2進数かつ指数表現で保存されているものとする．
単純にはIEEE 754のbinary64をイメージしてもらいたい．

そして，$n$,$m$を正の整数とした時，$f$は$n$,$m$を使って次のような分数に分解したい．

f = \frac{n}{m}

また，$f$の整数部を$I$，少数部を$F$として表すことにする．

f = I + F

10倍アルゴリズム

まず簡単に思いつくのが，人がよく使う手法だ．
ある数を10倍すれば小数点が一つ右にズレる．
したがって，小数点より右辺が0になるまで10倍すればよい．

ただし，プログラムで書く場合はちょっと注意が必要である．
単純にf=modf(f, &i)という風に更新していき，while(f!=0)とかすると，
誤差部分まで引っ張り出してしまい上手くいかない．

その辺を回避すれば一応求められるのだが，10倍というのもあまりに人間寄りだ．
これは最後の砦としたい．

C言語で書くとこんな感じになる．

#include <stdio.h>
#include <stdint.h>
#include <math.h>

uint64_t gcd_inner(uint64_t m, uint64_t n) {
  if (n == 0) {
    return m;
  }
  return gcd_inner(n, m % n);
}

uint64_t gcd(uint64_t m, uint64_t n) {
  /*  search a such that m = a*x and n = a*y
   * m = n * q + r = a * x
   * thus  q = a* z and r = a * w            */
  if (m < n) {
    uint64_t tmp = n;
    n = m;
    m = tmp;
  }
  return gcd_inner(m, n);
}

void frac_reduction(uint64_t *m, uint64_t *n) {
  uint64_t g = gcd(*m, *n);
  *m /= g;
  *n /= g;
}

void float2frac_10algorithm(double f, uint64_t *ret_n, uint64_t *ret_m) {
  double ans = f;
  uint64_t n = 0;
  uint64_t m = 1;
  while ((double)n / m != ans) {
    double i;
    f = modf(f, &i);
    n += i;
    n *= 10;
    m *= 10;
    f *= 10;
    printf("%f : %lu/%lu = %f\n", ans, n, m, (double)n / m);
  }
  frac_reduction(&n, &m);
  *ret_n = n;
  *ret_m = m;
}

int main(void) {

  double f = 0.0234;
  uint64_t n = 0;
  uint64_t m = 0;

  puts("10algorithm.");
  float2frac_10algorithm(f, &n, &m);
  printf("%f : %lu/%lu = %f\n", f, n, m, (double)n / m);
  return 0;
}

2倍アルゴリズム

ところで固定小数点数だと2倍すれば全ビットを上の桁に一つだけシフトすることと同義になる．
だから，10倍アルゴリズムではなく2倍アルゴリズムにすれば，もっと効率が良くなるんじゃないかと考えた．

ところが，浮動小数点数の場合指数と仮数に分かれており，
普通は指数の底は2である．
なので，2倍した所で指数の数字がインクリメントされるだけで，仮数は変化しない．

したがって，固定小数点と思って同様の処理を書くと，収束しないことがある．
例えばもっとも小さい桁が10進数の4だったりすると，
2倍を繰り返しても一番小さい桁は$4,8,6,2,4$と循環して10倍アルゴリズムのように桁が上がらない．

これでは使えない．

とあるページで見てしまったアイディア

色々探しているとよいアディアが見つかった．
ちょっとだけ問題を変形して$I=0$の時を考える．

f=F=\frac{1}{m}

この時，$n=1$になり$m$を求めればよいことになる．
ここで，$m$が求まれば問題は解決するのだが，これは単純に以下である．

m=\frac{1}{F}

ここで$m$が厳密な整数になる保証はないが，もし$\frac{n}{m}=f$になるのであれば，
正しい$n$,$m$が求まったといってよいだろう．

また，もし違った場合は$m$の小数部分のの誤差が無視できなかったということなので，
さらにその少数部分を同様の手法で求めればよい．
つまり，

m = I' + F'

という風に解釈して，$F'$を同様にして求めればよい．
ちなみに，$F\lt 1$であるため，$m \gt 1$になる．
したがって，$I' \ge 1$である．

アルゴリズム

前置きが長くなってしまったが，ここからが本題である．
少しややこしいいが，順を追って説明していく．

まず，fを次のように変形していく．

f=\text{floor}(f)+\left(f-\text{floor}(f)\right)=I+F

ここまではプログラミング言語風に書き換えただけで，やっていることは同じだ．
ここで，$F$は次のようになっていれば$f$は簡単に分数にできる．

F=\frac{1}{m}

したがって，$m$が正の整数という条件を満たすには，次のようにすればよい．

m = \text{floor}\left(\frac{1}{F}\right)

これで$m$が求まったので，次は$n$を求めよう．

I+F=I+\frac{1}{m}=\frac{Im+1}{m}

だから，

n = Im+1

とできる．
こうして求めた$\frac{n}{m}$が$f$と一致すれば万事解決だ．

しかし，実際には$m$を求める際にfloorを使っている．
したがって，この誤差が無視できない場合上手くいかない．

次にこの誤差について考えてみる．
$m$の導出の際にfloorを使わなければ，$\frac{1}{F}$は次のようになるはずである．

m \sim \frac{1}{F} = f' = \text{floor}(f') + (f' - \text{floor}(f')) = I' + F'

なぜなら，$F\lt 1$だからだ．

誤差が小さいということは$F'\sim 0$であればよいということだ．
この時，$m=f'=I'$となり一発で分数が求まる．
しかし，$F'$の値が大きいと$f'$と$I'$が乖離し誤差が大きくなってしまう．

ところで，この$f'$は実はそれ自身も分数で表わせるかもしれない．
となると，$m$自身も分数になってしまう．

しかし，思い出して欲しい．
今，

f=\frac{Im+1}{m}

であり，もし$m$が$\frac{a}{b}$という分数だったとする．
すると，

f=\frac{Im+1}{m}=\frac{I\frac{a}{b}+1}{\frac{a}{b}}=\frac{Ia+b}{a}

こんな感じで式変形できるから，$m$が分数であろうが，
$n$と$m$を新たに置き換えてやれば問題はない．

つまり，$f=\frac{n}{m}$を求めるには，$m$を求める必要がある．
しかし，$n$と$m$を求めるにはさらに別の分数を求める必要がある．そしてそのためにまた別の分数を求める必要がある．...，という入れ子構造になっているのだ．

$f$の分数を求めるには$f'$の分数表現が必要で，...，と続き，最終的にこれが止まるのは，

m_i = \frac{1}{F_i}

の$m_i$がちょうど整数になる．あるいは，十分な精度で近似できた時だ．

これで基本的な求め方は決まった．

次の問題は，この再帰的な分数を求めるアルゴリズムをどう実装するかだ．
ここがまたややこしい．

さて，一度問題を整理する．

f = I+F = \frac{n}{m_m}

I = \text{floor}(f)

F = f - I

そして，

f'=\frac{1}{F}

とすると，次のように変形できる．

f = I+F = I+\frac{1}{f'} = \frac{If'+1}{f'} =\frac{n}{m}

そして，今度は$f'$を同様にして分数表現に直していけば，いつかは十分な精度で収束するだろう．

しかし，ここで問題が発生する．
このプロセスをずっと繰り返し，その度に新たな$I_i$が求まったとすると，次のようになる．

f = I_0 + \frac{1}{I_1 + \frac{1}{ I_2 +\frac{1}{I_3+ \ldots} }}

ちなみに，このように分数の中に分数がある分数を連分数と呼ぶ．

この分数形式のままだと最後に$I_n$が求まった後，分母の一番下から順番に上へ上へと上がっていく必要がある．
つまり，$I_0,\ldots,I_n$まで全ての$I_i$をどこかに保存しておく必要がある．

これの何が問題かというと，先程while(f==0)は誤差を拾うため使えないといった．
だから，終了判定は最初のfをansとして保存していたとすると，while(ans == (double)n/m)を使う必要がある．
したがって，whileの終了判定の度にnとmを求める必要があり，それには全ての$I_i$を辿る必要がある．
これでは使えない．

というわけで，漸化式にして$n$と$m$が更新されるようにしたい．
ただ，この時点ではそもそも上手く漸化式に置き換えられるかどうかは分からない．
(本当をいうと連分数が漸化式にできるという話とおおよその形まで見てしまったている)

まずはそもそもどんな形にしたいのかから考える．
欲しいのは$f$の分数表現$\frac{n}{m}$である．
そして，ループを繰り返す毎に$n$と$m$が逐次更新されていき，その度に精度が上がっていって欲しい．

こんな風に．

n \leftarrow n,
m \leftarrow m

これを念頭に置いた上で漸化式を作るための法則を見出すために，順番に計算していく．

f = I_0 + F_0 = I_0 +  \frac{1}{I_1} = \frac{I_0 I_1 + 1}{I_1}

ここで，$I_1$は本当は$f_1$であり，

f_1 = I_1 + F_1

において$F_1 = 0$と置いた時という解釈である．
ところで，$F_1 = \frac{1}{I_2}$にできたとすると，
先程の式において$I_1$を次のように置き換えることと同じになる．

I_1 \rightarrow f_1 = I_1 +  F_1 = I_1 + \frac{1}{I_2}

したがって，一般化すれば更新式は次のようになる．

I_i \rightarrow I_i + \frac{1}{I_{i+1}} = \frac{I_i I_{i+1} + 1}{I_{i+1}}

それではもう一度最初から計算していく．

\frac{n_0}{m_0} =  I_0

\frac{n_1}{m_1} = I_0 +  \frac{1}{I_1} = \frac{I_0 I_1 + 1}{I_1}

\frac{n_2}{m_2}
=\frac{I_0 I_1 + 1}{I_1} = \frac{I_0 \left( I_1 + \frac{1}{I_2}\right) + 1}{I_1 + \frac{1}{I_2}}
= \frac{I_0 \left( I_1 I_2 + 1 \right) + I_2}{I_1 I_2 + 1}
= \frac{I_2 (I_0 I_1 + 1 )  +I_0}{I_1 I_2 + 1}
= \frac{I_2 n_1   + n_0 }{ I_2 m_1 + m_0}

\frac{n_3}{m_3}
= \frac{I_2 n_1   + n_0 }{ I_2 m_1 + m_0}
= \frac{\left(I_2 + \frac{1}{I_3}\right) n_1   + n_0 }{ \left(I_2 + \frac{1}{I_3}\right) m_1 + m_0}
= \frac{\left(I_3 I_2 + 1\right) n_1   + I_3 n_0 }{ \left(I_3 I_2 + 1 \right) m_1 + I_3 m_0}
= \frac{I_3 (I_2 n_1 + n_0 ) + n_1  }{ I_3 (I_2 m_1 + m_0 ) + m_1 }
= \frac{I_3 n_2 + n_1  }{ I_3 n_2 + m_1 }

何か法則が見えてきた気がする．
厳密な証明は後でするとして，おそらく次のようになっているはずだ．

\frac{n_i}{m_i}= \frac{I_i n_{i-1} + n_{i-2}  }{ I_i m_{i-1} + m_{i-2} }

そして，初期値は次のようになるはず．

n_0 = I_0, m_0 = 1

n_1 = I_0I_1 + 1, m_1 = I_1

ではこれを証明していこう．

ある正の整数$i \ge 2, i \in \mathbb{N}$において，以下の式が成り立つとする．

\frac{n_i}{m_i}= \frac{I_i n_{i-1} + n_{i-2}  }{ I_i n_{i-1} + m_{i-2} }

ただし，$i=0,1$においては以下である．

n_0 = I_0, m_0 = 1

n_1 = I_0I_1 + 1, m_1 = I_1

そして，ある数$i$までこの式が成立したとする．
そして，最初に定義の通り次のように置き換える．

I_i \rightarrow I_i + \frac{1}{I_{i+1}}

すると，$i+1$番目は次のようになるはずだ．

\frac{n_{i+1}}{m_{i+1}}
= \frac{\left(I_i + \frac{1}{I_{i+1}}\right) n_{i-1} + n_{i-2}  }{ \left(I_i + \frac{1}{I_{i+1}}\right) n_{i-1} + m_{i-2} }
= \frac{I_{i+1}\left(I_i + 1\right) n_{i-1} + I_{i+1}n_{i-2}  }{ I_{i+1}\left(I_i + 1\right) m_{i-1} + I_{i+1}m_{i-2} }
= \frac{I_{i+1}\left(I_i n_{i-1} + n_{i-2}\right) + n_{i-1} } { I_{i+1}\left(I_i m_{i-1} + m_{i-2}\right) + m_{i-1} }
= \frac{I_{i+1}n_i + n_{i-1} } { I_{i+1}m_i + m_{i-1} }

よって，任意の正の整数$i \in \mathbb{N}$の場合で次が成立する．

\frac{n_i}{m_i}= \frac{I_i n_{i-1} + n_{i-2}  }{ I_i m_{i-1} + m_{i-2} }

ただし，

n_0 = I_0, m_0 = 1

n_1 = I_0I_1 + 1, m_1 = I_1

最後にC言語で書くとこんな感じになる．

#include <stdio.h>
#include <stdint.h>
#include <math.h>

uint64_t gcd_inner(uint64_t m, uint64_t n) {
  if (n == 0) {
    return m;
  }
  return gcd_inner(n, m % n);
}

uint64_t gcd(uint64_t m, uint64_t n) {
  /*  search a such taht m = a*x and n = a*y
   * m = n * q + r = a * x
   * thus  q = a* z and r = a * w            */
  if (m < n) {
    uint64_t tmp = n;
    n = m;
    m = tmp;
  }
  return gcd_inner(m, n);
}

void frac_reduction(uint64_t *m, uint64_t *n) {
  uint64_t g = gcd(*m, *n);
  *m /= g;
  *n /= g;
}

void float2frac_algorithm(double f, uint64_t *ret_n, uint64_t *ret_m) {
  double ans = f;
  double i = 0;

  /* calculate I0 */
  f = 1.0 / modf(f, &i);
  uint64_t I[2] = {0};
  I[0] = i;
  if (ans == i) {
    *ret_n = I[0];
    *ret_m = 1;
    return;
  }

  /* calculate I1 */
  f = 1.0 / modf(f, &i);
  I[1] = i;

  uint64_t n[3] = {I[0], I[0] * I[1] + 1, 0};
  uint64_t m[3] = {1, I[1], 0};
  n[2] = n[1];
  m[2] = m[1];

  while ((double)n[2] / m[2] != ans) {
    printf("%f : %lu/%lu = %f\n", ans, n[2], m[2], 
	(double) n[2] / m[2]);
    f = 1.0 / modf(f, &i);
    n[2] = i * n[1] + n[0];
    n[0] = n[1];
    n[1] = n[2];
    m[2] = i * m[1] + m[0];
    m[0] = m[1];
    m[1] = m[2];
  }
  frac_reduction(&n[2], &m[2]);
  *ret_n = n[2];
  *ret_m = m[2];
  return;
}

int main(void) {

  double f = 0.0234;
  uint64_t n = 0;
  uint64_t m = 0;

  puts("algorithm.");
  float2frac_algorithm(f, &n, &m);
  printf("%f : %lu/%lu = %f\n", f, n, m, (double)n / m);
  return 0;
}