整数除算を浮動小数点演算でエミュレートできるか

古典的な（BigInt以前の）JavaScriptは数値型が倍精度浮動小数点数のみでした。5.2までのLuaも同様です。

このような言語で整数除算を行いたい場合は、浮動小数点演算を経由して行うことになります。例えばJavaScriptで32ビット整数の除算をやるならこんな感じです：

// x, yは整数値とする
function Int_div(x, y) {
    return Math.floor(x / y);
}
function Int_quot(x, y) {
    return Math.trunc(x / y);
}
// 0除算やオーバーフローはここでは無視する

さて、浮動小数点演算と言えば誤差です。この「浮動小数点演算による整数除算のエミュレート」は誤差の影響を受けないのでしょうか？

目次

定式化

「浮動小数点演算による整数除算のエミュレート」が誤差というか丸めの影響を受けるのは「\(n=\lfloor x/y\rfloor\)とおいた時に\(x/y\)の真の値が\(n+1\)に近すぎて最近接丸めの結果\(R_{\mathrm{nearest}}(x/y)\)が\(n+1\)になってしまう」場合、と言えるでしょう。

ここでは浮動小数点数の基数\(b\ge 2\), 精度\(p\ge 1\), 指数部の範囲を\(\mathit{emin}\)以上\(\mathit{emax}\ge 1\)以下とします。JavaScriptで使われる倍精度浮動小数点数の場合は\(b=2\), \(p=53\), \(\mathit{emin}=-1022\), \(\mathit{emax}=1023\)です。

指数部の範囲についてはIEEEでは\(\mathit{emin}=1-\mathit{emax}\)と定めているので\(\mathit{emin}\le 0\)です。まあ常識的には\(\mathit{emax}\ge p\)でしょう。

\(x/y>0\)の場合に\(R_{\mathrm{nearest}}(x/y)=n+1\)を不等式で書くと、\(i\)を\(b^i<n+1\le b^{i+1}\)とした時に\[n+1-b^{i-p+1}/2\le x/y<n+1\]が成り立つこと、となります。もっと抽象的にいうと「ある整数\(m\)について\((m+\mathrm{nextDown}(m))/2\le x/y<m\)となることはあるか」です。

\(x/y<0\)の場合は絶対値を考えて、「ある整数\(m\)について\(m<\lvert x/y\rvert\le(m+\mathrm{nextUp}(m))/2\)となることはあるか」となるでしょう。

整数の範囲についてですが、\(b=2\), \(p=53\)に対して\(x\), \(y\)が32ビットくらいならうまくいきそうな気がします。一方、精度の桁数と同じくらい、53ビットくらいになるとどうでしょうか？浮動小数点形式で正確に表現できる整数の範囲は\([-b^p,b^p]\)です（「safe integer」と呼ばれるものはこの端点を除くというのは「2^53+1 は素数か」に書いた通りです）。この範囲の整数の除算は正確に行えるのでしょうか？

以後、この問題では\(0\le x\le b^p\), \(0<y\le b^p\)とし、「ある整数\(m\ge 0\)について\[\frac{m+\mathrm{nextDown}(m)}{2}\le \frac{x}{y}<m\]または\[m<\left| \frac{x}{y}\right|\le\frac{m+\mathrm{nextUp}(m)}{2}\]となることはあるか」を考えます。

反例

実は\(b>2\)の場合は反例があります。

例えば\(b=10\), \(p=3\)とします。十進3桁です。999と2はそれぞれこの形式で正確に表現できる整数です。一方、999/2=499.5は正確に表現できません。これを丸めると、方式によっては（偶数丸めや四捨五入など）500が結果として返ってきて、\(\lfloor 999/2\rfloor=499\)と一致しません。

なので、「浮動小数点演算による整数演算のエミュレート」は無条件でうまくいくわけではなく、一般には\(b^p\)よりも小さいところに上限がありそうだということがわかります。

補題

次の補題を用意します。

補題．整数\(m\), \(k\), \(x\), \(y\)（ただし\(k>0\), \(y>0\)）に対して\[m-\frac{1}{k}\le\frac{x}{y}<m\]または\[m<\frac{x}{y}\le m+\frac{1}{k}\]が成り立つならば、\(k\le y\)である。

直観的に明らかですが、一応証明しておきます。

\(m-\frac{1}{k}\le\frac{x}{y}<m\)とします。両辺に\(k\)や\(y\)をかけて\[y(km-1)\le kx,\quad x<ym\]を得ます。これらは整数についての不等式なので\[y(km-1)\le kx,\quad x\le ym-1\]が成り立ちます。\(x\)を消去すると\[y(km-1)\le k(ym-1)\]となり、\(k\le y\)が得られます。

\(m<\frac{x}{y}\le m+\frac{1}{k}\)の場合も同様です。【証明終】

場合分け

以下の3つに場合分けします：

  1. \(m=0\), \(0<x/y\le\mathrm{nextUp}(0)/2\)の場合
  2. \(1\le m\), \(\frac{m+\mathrm{nextDown}(m)}{2}\le x/y<m\)の場合
  3. \(1\le m\), \(m<x/y\le\frac{m+\mathrm{nextUp}(m)}{2}\)の場合

場合1

\(m=0\), \(0<x/y\le\mathrm{nextUp}(0)/2\)の場合です。

直観的には、浮動小数点数は0の近くでは目盛りが細かいので、\(x/y\)は正確に表現しやすそうです。

非正規化数を考慮すると\(\mathrm{nextUp}(0)=b^{\mathit{emin}-p+1}\)なので、\[\frac{\mathrm{nextUp}(0)}{2}=\frac{1}{2b^{-\mathit{emin}+p-1}}\]です。\(\mathit{emin}\le 0\le p-1\)なので分母は整数です。よって補題より\(y\ge 2b^{-\mathit{emin}+p-1}\)となります。

\(y\)の範囲は\(2b^{-\mathit{emin}+p-1}\le y\le b^p\)となり、このような\(y\)が存在するならば\(2b^{-\mathit{emin}+p-1}\le b^p\)なので\(b^{\mathit{emax}}\le b^2/2\)です。よって\(\mathit{emax}=1\)でなくてはなりません。

逆に言うと、\(\mathit{emax}>1\)であれば場合1は起こらないということになります。

場合2

\(m\ge 1\), \(\frac{m+\mathrm{nextDown}(m)}{2}\le x/y<m\)の場合です。\(x/y\le b^p\)なので\(m\le b^p\)の場合を考えればよいです。

\(\mathrm{nextDown}(m)\)の具体的な値は、\(-1\le i\le p-1\)として\(b^i<m\le b^{i+1}\)の時\(\mathrm{nextDown}(m)=m-b^{i-p+1}\)となります。すると\[\frac{m+\mathrm{nextDown}(m)}{2}=m-\frac{b^{i-p+1}}{2}=m-\frac{1}{2b^{p-1-i}}\]です。\(i\le p-1\)なので\(2b^{p-1-i}\)は整数です。よって補題より、\(y\ge 2b^{p-1-i}\)です。

\(m-\frac{b^{i-p+1}}{2}\le\frac{x}{y}\)だったので\[y\left(m-\frac{b^{i-p+1}}{2}\right)\le x\]で、\(y\ge 2b^{p-1-i}\)を使うと\[2b^{p-1-i}\left(m-\frac{b^{i-p+1}}{2}\right)\le x\]が得られます。

\[2b^{p-1-i}\left(m-\frac{b^{i-p+1}}{2}\right)=2b^{p-1-i}m-1>2b^{p-1}-1\]なので、\(2b^{p-1}\le x\)が得られます。

等号成立条件について考えます。\(x=2b^{p-1}\)と仮定すると\[2b^{p-1-i}\left(m-\frac{b^{i-p+1}}{2}\right)\le 2b^{p-1}\]で、\(m\)について解くと\[m\le b^i\left(1+\frac{1}{2b^{p-1}}\right)=b^i+\frac{1}{2b^{p-1-i}}\]です。ここで\(i\le p-1\)と\(b^i<m\)を思い出すと\(b^i<m\le b^i+1/2\)を得て、\(m\)が整数という仮定に反します。よって、等号\(x=2b^{p-1}\)は成り立たないことがわかりました。

従って、場合2の結論は\(2b^{p-1}+1\le x\)です。

場合2の例は、\(b>2\)の場合に\(x=2b^{p-1}+1\), \(y=2\), \(m=b^{p-1}+1\)とすると\[\frac{m+\mathrm{nextDown}(m)}{2}=\frac{x}{y}\]が成り立ちます。

場合3

\(m\ge 1\), \(m<x/y\le\frac{m+\mathrm{nextUp}(m)}{2}\)の場合です。\(x/y\le b^p\)なので\(m<b^p\)の場合を考えればよいです。

\(\mathrm{nextUp}(m)\)の具体的な値は、\(0\le i\le p-1\)として\(b^i\le m<b^{i+1}\)の時\(\mathrm{nextUp}(m)=m+b^{i-p+1}\)となります。すると\[\frac{m+\mathrm{nextUp}(m)}{2}=m+\frac{b^{i-p+1}}{2}=m+\frac{1}{2b^{p-1-i}}\]です。\(i\le p-1\)なので\(2b^{p-1-i}\)は整数です。よって補題より、\(y\ge 2b^{p-1-i}\)です。

\(m<x/y\)だったので\(ym<x\)で、両辺が整数であることを使うと\[x\ge ym+1\ge 2b^{p-1-i}m+1\ge 2b^{p-1}+1\]を得ます。ただし\(y\ge 2b^{p-1-i}\)と\(b^i\le m\)を使いました。

よって場合3の結論は\(2b^{p-1}+1\le x\)です。

場合3の例は、\(b>2\)の場合に\(x=2b^{p-1}+1\), \(y=2\), \(m=b^{p-1}\)とすると\[\frac{m+\mathrm{nextUp}(m)}{2}=\frac{x}{y}\]が成り立ちます。

結論

\(\mathit{emax}\ge\max\{p,2\}\)とします。この時、絶対値が\(b^p\)以下の整数の除算を最近接丸めでエミュレートする時に問題になりうるのは以下の場合です：

  • \(b=2\): 問題なし
  • \(b>2\): \(x\ge 2b^{p-1}+1\)の場合。逆に、\(x\le 2b^{p-1}\)の場合は問題ない

特に、JavaScriptで絶対値\(2^{53}\)以下の整数の除算を行うために浮動小数点演算を利用しても問題ない、ということになります。

文献

このトピックに関しては

  • Lefèvre, Vincent. (2005). The Euclidean Division Implemented with a Floating-Point Division and a Floor.

というresearch report (?)があるようです。問題の定式化はこの記事とはちょっと違います。

このエントリーをはてなブックマークに追加