Haskell/GHCでの浮動小数点数の扱い

Haskell/GHCでの浮動小数点数の扱いを見ていきます。

浮動小数点数特有の落とし穴や、Haskell/GHC特有の落とし穴（バグ）にも若干触れます。

型

Haskell標準で規定されている浮動小数点数型は

• Float
• Double
• CFloat
• CDouble

の4つです。

大抵の環境では、

• Float はIEEE binary32（単精度）
• Double はIEEE binary64（倍精度）
• CFloat は Float のnewtype
• CDouble は Double のnewtype

と思って良いでしょう。というか、GHCがそうでない環境をサポートしているか疑問です。

リテラル

Haskellにももちろん小数リテラルがあります。

Haskellの小数リテラルはオーバーロードされており、 Fractional a => a という型がつきます。つまり、0.1 :: Float にも 0.1 :: Double にも 0.1 :: Rational にもなれます。

内部的には有理数 n % m と fromRational 関数を使って

0.1 → fromRational (1 % 10)
3.14 → fromRational (157 % 50)

という風に脱糖されます。fromRational に渡される値は正確なので、型次第で好きに調理できます。

GHC 8.4では HexFloatDecimals 拡張として十六進小数リテラルも導入されました。また、GHC 8.6で導入された NumericUnderscores 拡張を使うと数値リテラルに区切り文字 _ を入れることができます。

例：

{-# LANGUAGE HexFloatLiterals #-}
{-# LANGUAGE NumericUnderscores #-}

main = do
  print 0x1.0000_0000_0000_1 -- 末尾の指数部分は省略可能
  print 0x1p-1024

実行結果：

1.0000000000000002
5.562684646268003e-309

型クラスと関数

各種型クラスと関数を確認します。

Eqクラス

instance Eq a where
  (==), (/=) :: Double -> Double -> Bool

等価性を比較します。浮動小数点数に関して要注意なのは

• x == x が成り立つとは限らない（x がNaNの可能性があるので）
• x == y だからといって f x == f y が成り立つとは限らない（0 vs -0）

です。一方で、浮動小数点数だからといって何でも成り立たないというわけではなく、

• x == y と y == x は同値
• not (x == y) と x /= y は同値

です。

Ordクラス

instance Eq a => Ord a where
  compare :: Double -> Double -> Ordering
  (<), (<=), (>), (>=) :: Double -> Double -> Bool
  max, min :: Double -> Double -> Double

浮動小数点数には順序づけられない値、NaNがあります。なので、

• x < y と not (x >= y) は同値ではない
• x <= y と not (x > y) は同値ではない
• x <= x は成り立つとは限らない

です。一方で

• x < y と y > x は同値
• x <= y と y >= x は同値
• x <= y && y <= x と x == y は同値

です。

GHCでは Float, Double に対する compare は

compare x y | x < y     = LT
            | x == y    = EQ
            | otherwise = GT

と定義されているようです。つまり、一方がNaNなら常にGTを返します。

min, max はデフォルト定義

max x y | x <= y    = y
        | otherwise = x
min x y | x <= y    = x
        | otherwise = y

が採用されます。つまり、

• (min x y, max x y) は常に (x, y) または (y, x) のいずれかを返す
• min 0.0 (-0.0) は 0.0 を返し、 max 0.0 (-0.0) は -0.0 を返す

となります。

Numクラス

instance Num a where
  (+), (-), (*) :: a -> a -> a
  negate, abs :: a -> a
  signum :: a -> a
  fromInteger :: Integer -> a

(+), (-), (*), abs はおなじみのソレです。

signum は、

• 引数が正の場合は 1.0 を
• 引数が負の場合は -1.0 を
• それ以外（±0.0, NaN）の場合は引数をそのまま

返します。

Float/Double に関する fromInteger は要注意で、現在のGHCでは丸め方法が一貫していません。引数が小さい場合は最近接偶数丸めを、大きい時は0方向への丸めを行います。

• conversion from Integer to Double rounds differently depending on the value (#15926) · Issues · Glasgow Haskell Compiler / GHC · GitLab

実行例：

Prelude Numeric> showHFloat (fromInteger 0xFFFFFFFFFFFFFC0 :: Double) "" -- 最近接偶数丸めが使用される
"0x1p60"
Prelude Numeric> showHFloat (fromInteger 0xFFFFFFFFFFFFFC00 :: Double) "" -- 切り捨てが行われる
"0x1.fffffffffffffp63"

Fractionalクラス

instance Num a => Fractional a where
  (/) :: a -> a -> a
  recip :: a -> a
  fromRational :: Rational -> Double

(/) はおなじみの除算です。recip x は 1.0 / x として定義されています。

注意点は

• x / y == x * recip y は一般には成立しない
  • 例：x = 0.01 :: Double, y = 0.1 :: Double とすると x / y == 9.999999999999999e-2, x * recip y == 0.1 となります。
  • どの組み合わせが反例となるかは精度に依存します。例えば、x = 0.01 :: Float, y = 0.1 :: Float の場合は x / y == x * recip y となります。

です。

fromRational は fromInteger と異なり、正しく最近接偶数丸めを行うようになっています。

Prelude Numeric> showHFloat (fromRational 0xFFFFFFFFFFFFFC00 :: Double) "" -- 最近接偶数丸めが使用される
"0x1p64"

よって、注意点としては

• fromRational (toRational x) == fromInteger x は（fromInteger の問題により）常に成り立つとは限らない

となります。

Floatingクラス

instance Fractional a => Floating a where
  pi :: a
  exp, log, sqrt :: a -> a
  (**), logBase :: a -> a -> a
  sin, cos, tan :: a -> a
  asin, acos, atan :: a -> a
  sinh, cosh, tanh :: a -> a
  asinh, acosh, atanh :: a -> a

  -- GHC.Float (since base-4.9.0.0):
  log1p, expm1 :: a -> a
  log1pexp, log1mexp :: a -> a

おなじみの初等関数です。sqrt の精度は特に規定はありませんが、大抵の環境では真の値に最も近い浮動小数点数を返すと思います。

GHC.Float モジュールからは見慣れない関数がexportされています。

log1p x は $\log(1+x)$ を、 expm1 x は $\exp x-1$ を計算します。これらは、x が 0 に近いときに log (1 + x) や exp x - 1 よりも正確な値を返すことが期待されます。

log1pexp x は $\log (1 + \exp x)$ を、 log1mexp x は $\log (1 - \exp x)$ を計算します。

Realクラス

instance (Num a, Ord a) => Real a where
  toRational :: a -> Rational

引数を（多倍長）有理数に変換します。整数型、有理数型、浮動小数点数型などがこのクラスのインスタンスになれます。

「有理数に正確に変換できる型」なので、実数の部分集合であっても $\mathbf{Q}(\sqrt{2})$ のように無理数を正確に表現できる型はこのクラスのインスタンスになれません。まあこの記事の本題は浮動小数点数なので、この件には深く突っ込まないでおきましょう。

浮動小数点数についての注意点としては、

• 引数が無限大やNaNの場合はめちゃくちゃな値を返す

ことです。toRational する際は引数が有限の浮動小数点数であることを確認しましょう。

RealFracクラス

instance (Real a, Fractional a) => RealFrac a where
  properFraction :: Integral b => a -> (b, a)
  truncate, round, ceiling, floor :: Integral b => a -> b

浮動小数点数や有理数を整数に変換します。引数が整数でなかった場合は、引数の両隣の整数のうち、

• truncate は絶対値が小さい方を
• round は最近接偶数丸めで（C言語の round 関数とは異なります）
• ceiling は大きい方を
• floor は小さい方を

返します。

properFraction x は $x=n+f$ となるような整数 $n$ と実数 $f$ を返します。$n$ は truncate x と等しく、$f$ の符号は $x$ と同一です。

浮動小数点数に関しての注意点は

• 引数が無限大やNaNの場合はめちゃくちゃな値が返る（toRational と同様）
• 変換先が固定長整数型の場合に範囲外の値を渡したら何が返ってくるかわからない

ことです。これらの関数を使う際は引数が有限の浮動小数点数であって変換先の型で表現できることを確認しましょう。

RealFloatクラス

instance (RealFrac a, Floating a) => RealFloat a where
  floatRadix :: a -> Integer
  floatDigits :: a -> Int
  floatRange :: a -> (Int, Int)
  decodeFloat :: a -> (Integer, Int)
  encodeFloat :: Integer -> Int -> a
  exponent :: a -> Int
  significand :: a -> a
  scaleFloat :: Int -> a -> a
  isNaN, isInfinite :: a -> Bool
  isDenormalized, isNegativeZero :: a -> Bool
  isIEEE :: a -> Bool
  atan2 :: a -> a -> a

浮動小数点数の中身を操作できるクラスです。

floatRadix, floatDigits, floatRange は浮動小数点数形式についての情報を与えます。引数型がアレですが、引数を評価しない（undefined を渡しても問題ない）ことを期待できます。

IEEE 754準拠の浮動小数点数型なら floatRadix は 2 または 10 を返すはずです。以後この値を $b$ とおきます。

floatDigits は仮数部の桁数 $d$ です。Float なら24、Double なら53です。ケチ表現は関係ないので注意してください。

floatRange は正規化数の仮数部の絶対値を $[1/b,1)$ に正規化した時の指数部の範囲です。Float なら (-125,128), Double なら (-1021,1024) です。いずれも両端を含みます。

decodeFloat は浮動小数点数 $x$ を仮数部 $m$ と指数部 $e$ に分解します。$x = mb^e$ という関係が成り立ちます。$x=0$ ならば $(m,e)=(0,0)$ で、そうでなければ $b^{d-1}\le \lvert m\rvert\lt b^d$ です。

encodeFloat は仮数部と指数部から浮動小数点数を構築します。$m$, $e$ が与えられたときに $mb^e$ を返します。丸め方法の規定はなさそうなので、正確な値を返して欲しければ $m$ の絶対値は（正規化数の場合） $b^d$ 以下とするのが無難でしょう。

exponent, significand は仮数部の絶対値を $[1/b,1)$ に正規化した時の指数部と仮数部を返します。引数が正規化数の場合は exponent は floatRange の範囲の値を返しますが、非正規化数の場合はその範囲をはみ出ます。

scaleFloat は浮動小数点数に基数の冪乗を乗算します。正規化数どうしで変換する場合は正確な値を返しますが、変換後の値が非正規化数の場合（アンダーフロー）は丸めが発生します。

isﾅﾝﾄｶ 系のやつは特殊な浮動小数点数についての性質を確認する関数です。

atan2 y x は (sin t, cos t) の逆関数みたいなやつです。

冪乗

Haskellには3種類の冪乗演算子があります。

(^) :: (Num a, Integral b) => a -> b -> a
(^^) :: (Fractional a, Integral b) => a -> b -> a
(**) :: Floating a => a -> a

• (^) は自然数乗で、任意の Num について適用できます。繰り返し乗算で実装されます。
• (^^) は整数乗で、任意の Fractional について適用できます。
• (**) は実数の実数乗あるいは、複素数の複素数乗です。

絶対値が小さい数を (^^) で構築する際は注意が必要です。計算途中でオーバーフロー・アンダーフローが起こる可能性があります。非正規化数の記述にはリテラルや encodeFloat を使うべきでしょう。

Prelude> :set -XHexFloatLiterals
Prelude> encodeFloat 1 (-1074) :: Double
5.0e-324
Prelude> 0x1p-1074 :: Double
5.0e-324
Prelude> (1/2)^^(1074 :: Int) :: Double
5.0e-324
Prelude> 2^^(-1074 :: Int) :: Double -- !!
0.0

• (^^) is not correct for Double and Float (#4413) · Issues · Glasgow Haskell Compiler / GHC · GitLab

型変換

整数と浮動小数点数の変換、浮動小数点数同士の変換のための関数がいくつか用意されています。

-- 整数からの変換
fromInteger  :: Num b =>         Integer -> b
toInteger    :: Integral a =>          a -> Integer
fromIntegral :: (Integral a, Num b) => a -> b

-- 有理数・浮動小数点数から整数への変換
truncate :: (RealFrac a, Integral b) => a -> b
ceiling  :: (RealFrac a, Integral b) => a -> b
floor    :: (RealFrac a, Integral b) => a -> b
round    :: (RealFrac a, Integral b) => a -> b

-- 有理数を介した変換
fromRational :: Fractional b =>   Rational -> b
toRational   :: Real a =>                a -> Rational
realToFrac  :: (Real a, Fractional b) => a -> b

-- GHC.Float
float2Double :: Float -> Double
double2Float :: Double -> Float

fromIntegral は fromInteger . toInteger として定義されていますが、色々と書き換え規則が定義されており、Integer を介さない変換が可能なことも多いです。

（それでも、前述の Float / Double の fromInteger の問題のために書き換え規則で動作が変わるコーナーケースが存在します。）

要注意なのは realToFrac です。この関数は fromRational . toRational として定義されていますが、色々と書き換え規則が定義されており、Rational を介さない変換が可能な場合もあります。浮動小数点数同士の場合は、無限大やNaNを変換先の型のそれへ変換される可能性もあります。

ただ、最適化が効かない場合は toRational が使用されます。 toRational は前述の通り、入力が無限大やNaNの場合に無意味な値（何らかの有限値）を返します。それを fromRational に与えると、その有限値を結果の型で解釈した値を返します。

要するに、 realToFrac に無限大やNaNを与えるとヤバイ、ということです。次のコードを最適化ありとなしでそれぞれコンパイル・実行してみましょう：

import Numeric

nanF :: Float
nanF = 0 / 0

nanD :: Double
nanD = 0 / 0

infinityF :: Float
infinityF = 1 / 0

main = do
  -- 最適化の有無で出力が変わる：
  putStrLn $ showHFloat (realToFrac nanF :: Double) ""
  putStrLn $ showHFloat (realToFrac infinityF :: Double) ""
  putStrLn $ showHFloat (realToFrac nanD :: Float) ""

最適化なし：

$ stack ghc -- RealToFrac.hs
[1 of 1] Compiling Main             ( RealToFrac.hs, RealToFrac.o )
Linking RealToFrac ...
$ ./RealToFrac
-0x1.8p128
0x1p128
-Infinity

最適化あり：

$ stack ghc -- -O2 RealToFrac.hs
[1 of 1] Compiling Main             ( RealToFrac.hs, RealToFrac.o ) [Optimisation flags changed]
Linking RealToFrac ...
$ ./RealToFrac
NaN
Infinity
NaN

最適化なしの場合は無限大やNaNが化けてしまっていることがわかります。

対策は、

• realToFrac の入力が常に有限値であることを確認する（無限大やNaNの場合は自前で処理する）
• GHC.Float モジュールの float2Double / double2Float を直接呼び出す（この場合、無限大やNaNは正しく処理される）

となるでしょう。

関連issue：

• realToFrac doesn't sanely convert between floating types (#3676) · Issues · Glasgow Haskell Compiler / GHC · GitLab
• realToFrac between Float and Double gives different results depending on compile flags (#16519) · Issues · Glasgow Haskell Compiler / GHC · GitLab

文字列との変換

Prelude の関数でも文字列との相互変換ができますが、フォーマット等を指定したい場合は Numeric モジュールの関数が使用できます。

面倒なので詳しくは触れません。

拙作ライブラリー

Haskell標準ではIEEE準拠の浮動小数点数を扱う道具がそれなりに揃っていますが、IEEE 754-2019で規定された演算と比較すると足りないものも色々あります。そこで、筆者が作っているのがfp-ieeeおよびrounded-hwというライブラリーです。

• minoki/haskell-floating-point: Haskell libraries for floating point numbers
  • fp-ieee: IEEE準拠の各種演算（未リリース）
    • 丸め方法を指定できる fromInteger / fromRational 等
    • fused multiply-add
    • 四捨五入を行う roundAway（C言語の round 関数に相当）
    • truncate, ceiling, floor, round, roundAway の、浮動小数点数を返すバージョン
    • 「まともな」（最適化の有無で挙動が変わらない）浮動小数点数型同士の変換：realFloatToFrac
    • NaNのペイロードの扱い
  • rounded-hw: 丸め方法を指定した各種演算
    • 区間演算ができるようになる

fp-ieeeの方針としては、独自の型クラスはなるべく提供せずに、ジェネリックな型を持つ関数に対して書き換え規則を定義することによって個々の型に最適化された実装を用意します。また、基本的な機能はHaskellだけで完結するようにして、C FFIは最適化のためだけに利用します。

まだ足りない機能（sqrt の引数の型と返り値の型が異なるやつ）がありますが、Apple Silicon Mac上で動作する見込みが立ったらそろそろリリースしても良いかなと思っています。