浮動小数点数」タグアーカイブ

TwoSumの証明に2週間かかった

前に「浮動小数点数小話」という同人誌を出しましたが、私としては浮動小数点数についてもっとちゃんとした本をいずれ出したいと思っています。執筆に着手するのは早い方が良いので、浮動小数点数についてのまとめノートみたいものを書き始めています。

続きを読む

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

古典的な(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除算やオーバーフローはここでは無視する

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

続きを読む

2^53+1 は素数か

ふと \(2^{53}+1\) が素数かどうか気になりました。

短い答え:No

Linuxだと factor コマンドで素因数分解できます。

$ factor $((2**53+1))
9007199254740993: 3 107 28059810762433

macOSだとHomebrewやMacPortsでcoreutilsを入れる必要があります。

筆者のMacには factor コマンドは入っていなかったので、最初はMaximaで確認しました:

$ maxima    
Maxima 5.45.1 https://maxima.sourceforge.io
using Lisp SBCL 2.2.9
Distributed under the GNU Public License. See the file COPYING.
Dedicated to the memory of William Schelter.
The function bug_report() provides bug reporting information.
(%i1) factor(2^53+1);
(%o1)                        3 107 28059810762433

というわけで、 \(2^{53}+1\) は\[2^{53}+1=3\times 107\times 28059810762433\]と素因数分解でき、素数ではありません。

暗算でもできるのでは?

今回調べたいのは \(2^{53}+1\) という形の数です。ちょっと頭を使えば、この数が3で割り切れることがわかります:\begin{align*}2^{53}+1&=2^{53}-(-1)^{53}\\&=(2-(-1))(2^{52}+\cdots+2^{52-i}\cdot(-1)^i+\cdots+(-1)^{52})\end{align*}と因数分解できます。

一般化すると、\(2^{2m+1}+1\)の形の整数はすべて3で割り切れることがわかります。つまり、\(m\neq 0\)で\(2^{2m+1}+1\)が3よりも大きくなる場合はこの形の数はすべて合成数です。

もっと一般化

もっと一般に、\(2^n+1\)の形の整数が素数となるのは\(n=0\)または\(n=2^k\ (k\geq 0)\)の場合に限られます。

\(n=2^k(2m+1)\)とおくと、\begin{align*}2^n+1&=2^{2^k(2m+1)}+1\\&=\left(2^{2^k}\right)^{2m+1}+1\\&=\left(2^{2^k}\right)^{2m+1}-(-1)^{2m+1}\\&=\left(2^{2^k}-(-1)\right)\left(\left(2^{2^k}\right)^{2m}+\cdots+\left(2^{2^k}\right)^{2m-i}\cdot(-1)^i+\cdots+(-1)^{2m}\right)\\
&=\left(2^{2^k}+1\right)\left(\left(2^{2^k}-1\right)\left(2^{2^k}\right)^{2m-1}+\cdots+\left(2^{2^k}-1\right)\left(2^{2^k}\right)^{2(m-i)+1}+\cdots+\left(2^{2^k}-1\right)\cdot 2^{2^k}+1\right)\\
&=\left(2^{2^k}+1\right)\left(\sum_{i=0}^{m-1} \left(2^{2^k}-1\right)\left(2^{2^k}\right)^{2i+1}+1\right)\end{align*}と因数分解できます。

\(k\geq 0\)なので、\(2^{2^k}+1\geq 3\)であり、最初の因数は非自明です。そして、2番目の因数に登場する\(2^{2^k}-1\)は1以上で、\(m\geq 1\)ならば2番目の因数も非自明です。

よって、\(2^{2^k(2m+1)}+1\)の形の数が素数となるのは\(m=0\)の場合に限ることがわかりました。

フェルマー数

\(2^{2^k}+1\)の形の数はフェルマー数と呼ばれています。\(k=0,1,2,3,4\)の場合はこれは素数となり、フェルマー素数と呼ばれています。

ところでなんで \(2^{53}+1\) を調べたかったのか

一部のプログラミング言語では、倍精度浮動小数点数型を唯一の数値型として提供しています。BigInt以前のJavaScriptや、5.2までのLuaなど。

こういう言語では、絶対値が\(2^{53}\)以下の整数は正確に表現できます。そして演算結果の絶対値がそれを超えると、丸めが発生して正しい答えが返ってきません。JavaScriptで試してみましょう:

$ node
Welcome to Node.js v17.9.1.
Type ".help" for more information.
> 2**53 - 1
9007199254740991
> 2**53
9007199254740992
> 2**53 + 1
9007199254740992
> 2**53 + 2
9007199254740994
> 2**53 + 3
9007199254740996

2**53 + 1 を表現できておらず、「表現可能な最も近い値のうち、仮数部の末尾が偶数な方」が表示されています(最近接偶数丸め)。

では、演算結果が「正確に表現できる範囲」を超えたかどうかはどうやったらわかるでしょうか?まず、演算結果の絶対値が\(2^{53}+2\)以上であれば間違いなく超えています。

では、演算結果の絶対値が\(2^{53}\)以下なら結果は正確と言えるでしょうか?否ですね。2**53 + 1 の演算結果は\(2^{53}\)以下であるにも関わらず、不正確です。

一方、演算結果の絶対値が\(2^{53}-1\)以下であれば結果は正確であることが保証されます。この範囲、\([-(2^{53}-1),2^{53}-1]\)はJavaScript界隈ではsafe integerと呼ばれています。

さて、筆者が作っているStandard ML処理系、LunarMLでもsafe integerに相当する整数型を提供できると便利です。提供するとするとビット数は符号ビットも含めて54ビット、 Int54 みたいなものになります。

ですが、Standard MLの固定長整数は2の補数表現を仮定しており、表現できる範囲は\([-2^n,2^n-1]\)の形である必要があります。これとsafe integerを比べると、\(-2^{53}\)がはみ出ます。

なので、LunarMLで Int54 を提供するには、演算結果が\(-(2^{53}+1)\)の場合を何らかの方法で検出し、オーバーフロー例外を起こさなければなりません。この時、\(2^{53}+1\)が素数だったなら乗算の場合に演算結果が\(-(2^{53}+1)\)になる可能性を考慮しなくてよかったのに、現実は残酷でした、という話です。

浮動小数点数の文字列化(基数変換)

動機付けと問題

計算機の内部では二進浮動小数点数が使われることが多い一方で、プログラムのソースコードやテキストベースのデータ形式(例:JSON)では十進小数が使われることが多い。

データのシリアライズ等で、内部的な二進浮動小数点数を十進小数に変換して、再度二進小数に戻すという操作が考えられる。この時、元々の二進小数の値が保持されることが望ましい。

有限桁の二進小数は原理的には有限桁の十進小数で表現できるが、指数部が大きかったり小さかったりすると仮数部の桁数も膨張するため実用的ではない。

そこで、基数変換の際に丸めが発生することを許容して、二進小数→十進小数→二進小数の変換が恒等写像となるようにしたい。ただし丸めの方法は最近接丸めであるとする。この時、

  • 途中の十進小数の仮数部は何桁あれば十分か?
  • なるべく短い桁数の十進小数を、正しい丸めで得るためのアルゴリズムはどのようなものか?

という問題が考えられる。

続きを読む

浮動小数点数の演算の丸めモードを制御する

以前、Haskellで浮動小数点数の丸めモードを制御したいという話を書いた。

以前の記事では、丸めモードを変える方法としてはC言語のfesetroundだけを扱ったが、高速化するなら各アーキテクチャ固有の方法を使えた方が有利である。また、動的な丸めモード指定ではない静的な丸めモードについてもこの記事では触れる。

続きを読む

Haskellでの浮動小数点数の方向付き丸めを考える

週刊 代数的実数を作る」の #5 で、区間演算と方向付き丸めの話を書いた。浮動小数点数の計算は不正確だと思われがちだが、方向付き丸め等をうまく使えばある種の「正しい結果」(この数は確実に1.0より大きい、等)を得ることができる、という話だ。

MPFRのようなソフトウェア実装の浮動小数点数だと引数で丸めモードを指定できる。しかし、ハードウェア組み込み(注)の浮動小数点数(floatやdouble)の丸めモードを指定する方法は、言語や実行環境に依存する。

(注:環境によってはCPUに浮動小数点演算器が組み込まれておらず、floatやdoubleの演算もソフトウェア実装だったりするが、我々が普段使うPCではfloatやdoubleはほぼ確実にハードウェア実装されているため、以下「ハードウェア実装」で通す)

CやFortranのような低レベルかつ数値計算のニーズがあるようなプログラミング言語だと、丸めモードを変更する手段が用意されているが、JavaScriptやHaskellなど、そういうニーズが薄い言語では丸めモードの変更には対応していないことが多い。せいぜい、「浮動小数点数演算には最近接丸めを使用する」と規定されているのが関の山である。

しかし、「Haskellではできない処理がある」というのはどうにも気にくわない。どうにかして、ハードウェア組み込みの浮動小数点数の方向付き丸めをHaskellで扱う方法を考えたい。

ちなみに、ハードウェア組み込みの浮動小数点数にこだわらないのであれば、MPFRの丸めモード指定を使うHaskellライブラリーがすでに存在する:

http://hackage.haskell.org/package/rounded

続きを読む

Lua で遊ぶ浮動小数点数

(倍精度)浮動小数点数でいろいろ遊ぶ際に、 Lua が便利なのではないかと思った。以下、特に断らない限り「浮動小数点数」と言ったら倍精度のものを指す。

なぜ Lua を使うのか?

Lua はC言語で実装されていて、数値の扱いについてはC言語と近い挙動を示すと考えられる。C言語と違って累乗の演算子 (x^y) があるのが地味に便利である。

Lua 5.2 以降では、浮動小数点数の16進表記をサポートするようになった。ソースコード中にリテラルで 0x1.fp2 と書けるし、文字列から数値に変換するときに tonumber("0x1.fp2") と書けるし、数値から文字列に変換するときは %a または %A を使って string.format("%a",7.75) と書ける。

あとは、標準ライブラリに足りない機能があったときに簡単にC言語で拡張ライブラリを書ける。

続きを読む