理想の浮動小数点数である「浜田のURR」を理解する
はじめに
コンピュータでの標準的な「実数」の表現のしかたに「浮動小数点数」がある。浮動小数点数の符号化の様式は、IEEE 754で標準化されている。ここに定義されている符号化のやりかたは、実用的な意味では、よくできている。ただ、数学的な意味で、より洗練された符号化のやりかたがある。それが「浜田のURR(Universal Representation of Real numbers)」だ。この符号化を、Haskellのコードを書きながら理解していこう。
注記
ここで書くコードは、URRのおもに数学的な定義を理解することを目的としたものなので、必ずしも効率的なコードではない。「わかりやすさ」と「効率」とでは前者を優先する。
読むかたへ
とりあえず、書ききりました。ただ、最後のほうは、だいぶつかれていますので、やや説明が雑かもしれません。まちがいなどありましたら「やさしく」ご指摘ください。わからないところなどありましたら、お気軽にご質問ください。
コードはGitHubに置いてあります。
GitHub: test_haskell/tribial/qiita/urr
参考
- 元論文「二重指数分割に基づくデータ長独立実数値表現法II」(PDF)
- Wikipedia「浮動小数点数 - 浜田のURR」
浮動小数点数とは
指数と仮数(と符号)による実数の表現が浮動小数点数だ。たとえば10進法で、つぎのような数を考える。
123.456
これは指数と仮数を使って、つぎのように表現できる。
1.23456 * 10 ^ 2
この値の指数は2であり、仮数は123456である。べつの例をみてみよう。
0.00987654
9.87654 * 10 ^ (- 3)
指数: - 3
仮数: 987654
このように浮動小数点数では、大きな数も小さな数も(相対的な意味で)おなじくらいの正確さで表現できる。実際には10進法ではなく2進法が使われているが、基本的な考えかたは、おなじだ。
Haskellには、浮動小数点数としてFloatとDoubleがある。
> 1.23 :: Float
1.23
> 1.23 :: Double
1.23
> pi :: Float
3.1415927
> pi :: Double
3.141592653589793
Floatは32ビットのデータであり、つぎのようになっている。
符号部: 1ビット
指数部: 8ビット
仮数部: 23ビット
Doubleは64ビットのデータであり、つぎのようになる。
符号部: 1ビット
指数部: 11ビット
仮数部: 52ビット
ここで、指数部と仮数部へのビットのわりふりかたは、「人間の都合」であり、数学的な根拠はない(と思う)。また、128ビットや256ビットなど、さらに大きな浮動小数点数を定義しようとしたら、新たな規格を作成する必要がある(ちなみに128ビットの規格はすでにある)。
浜田のURR
浜田のURR(Universal Representation of Real number)では、指数部と仮数部のそれぞれのビット幅は固定ではない。大きな数や小さな数では指数部が長くなり、+1や-1に近い数では指数部が短くなる。指数部の幅は恣意的なものではなく、符号化のなかで自ずから決まってくるものだ。また、数を表すデータのビット長をのばしたとしても、新たな規格は必要なく、自然な延長によって「細かさ」と「最小と最大の値」との両方が拡張される。
また、しばしば使われる「ふつうの大きさの数」のあたりでは、正確な値を表現することができ、頻度が低いだろう「すごく小さい数」や「すごく大きい数」については「おおざっぱではあっても、表現できる」といった特徴もある。
数をビットで表現するということ(符号なし8ビット整数の例で)
たとえば、8ビットで非負整数を表現することを考える。しばしば使われる符号化を考えると、つぎのようになるだろう。
ビットの並び 十進数
00001011 11
10010010 146
最下位ビットから0, 1, 2, 3, 4, 5, 6, 7と番号をふる。それぞれのビットについて、その番号をnとしたとき、0ならば0を、1ならば2^nを、その値としてあたえて、それらの総和がそのビット列が表す値になっている。
この符号化の様式を、べつの角度からみてみよう。表せる全体の数の範囲(ここでは0以上256未満)を、つぎつぎに2分していっていると考えることができる。上記の00001011の例でみてみよう(m以上n未満を[m, n)と表現する)。
[0, 256)
0: [0, 128)
0: [0, 64)
0: [0, 32)
0: [0, 16)
1: [8, 16)
0: [8, 12)
1: [10, 12)
1: [11, 12)
はじめの区間は[0, 256)である。この区間を半分にする。つまり[0, 128)と[128, 256)にわける。最上位ビットが0なので、小さいほうの区間[0, 128)を選ぶ。さらに、つぎのビットも0なので、半分にした小さいほうの区間を選び[0, 64)となり、[0, 32)、[0, 16)と続く。この[0, 16)を半分にした[0, 8)と[8, 16)からは、対応するビットが1なので、[8, 16)が選ばれる。同様に[8, 12)と[12, 16)からは[8, 12)が選ばれ、さらに[10, 12)、[11, 12)のように範囲がせばめられる。ここで、[11, 12)の範囲内の整数は11のみなので、値は11に決定する。10010010の例も挙げておく。
[0, 256)
1: [128, 256)
0: [128, 192)
0: [128, 160)
1: [144, 160)
0: [144, 152)
0: [144, 148)
1: [146, 148)
0: [146, 147)
コードにする
このアルゴリズムをコードにしてみよう。
{-# OPTIONS_GHC -Wall -fno-warn-tabs #-}
import Data.Word
data Bit = O | I deriving Show
fromString :: String -> [Bit]
fromString = map fc
where fc '0' = 0; fc '1' = I; fc _ = error "Oops!"
toWord8 :: Int -> Int -> [Bit] -> Word8
toWord8 mn _ [] = fromIntegral mn
toWord8 mn mx (O : bs) = toWord8 mn ((mn + mx) `div` 2) bs
toWord8 mn mx (I : bs) = toWord8 ((mn + mx) `div` 2) mx bs
試してみる。
> :load word8.hs
> toWord8 0 256 $ fromString "00001011"
11
> toWord8 0 256 $ fromString "10010010"
146
上述のアルゴリズムをそのまま実装しただけだ。逆向きの変換も実装しておこう。
fromWord8 :: Int -> Int -> Word8 -> [Bit]
fromWord8 mn mx _ | mn + 1 >= mx = []
fromWord8 mn mx w
| fromIntegral w < sp = O : fromWord8 mn sp w
| otherwise = I : fromWord8 sp mx w
where sp = (mn + mx) `div` 2
はじめに、区間の幅が1以下(mn + 1 >= mx)のときの、値が決定し、それ以上のビットが不要である場合を定義する。そのうえでmnとmxの中間であるspと、表現したい値とを比較して、適切なビット(OまたはI)を追加したうえで、新しい範囲について、演算をくりかえしている。
コードを抽象化する
8ビット整数とビット列との相互変換のコードは、より抽象化することができる。「範囲をせばめていくことで、ビット列を値にする/値をビット列にする」というやりかたそのものをコードにすることができる。
ビット列から値を導出する
「ビット列から値を導出する関数」には、つぎのようなパラメータを指定する。
- 「範囲」から「分割する点」を導出する関数
- 最小値を値に変換する関数
コードは、つぎのようになる。
toValue :: (r -> r -> r) -> (r -> a) -> r -> r -> [Bit] -> a
toValue _ gv mn _ [] = gv mn
toValue sp gv mn mx (O : bs) = toValue sp gv mn (sp mn mx) bs
toValue sp gv mn mx (I : bs) = toValue sp gv (sp mn mx) mx bs
これを利用して、関数toWord8を書き換える。
toWord8 :: Int -> Int -> [Bit] -> Word8
toWord8 = toValue (\a b -> (a + b) `div` 2) fromIntegral
値からビット列を導出する
「値からビット列を導出する関数」には、つぎのようなパラメータを指定する。
- ビット列の長さ: 「範囲」からビット列の終端を決められない場合に必要
- 「範囲」から「終端かどうか」を判断する関数
- 「範囲」から「分割する点」を導出する関数
- 値が「分割する点」よりも小さいかどうかを判断する関数
コードは、つぎのようになる。
fromValue ::
Maybe Word -> (r -> r -> Bool) -> (r -> r -> r) -> (a -> r -> Bool) ->
r -> r -> a -> [Bit]
fromValue ln fn _ _ mn mx + | Just 0 <- ln = [] | fn mn mx = []
fromValue ln fn sp lt mn mx x
| x `lt` s = O : fromValue (subtract 1 <$> ln) fn sp lt mn s x
| otherwise = I : fromValue (subtract 1 <$> ln) fn sp lt s mx x
where s = sp mn mx
これを使って、関数fromWord8を書き換える。
fromWord8 :: Int -> Int -> Word8 -> [Bit]
fromWord8 = fromValue Nothing
(\a b -> a + 1 >= b)
(\a b -> (a + b) `div` 2)
(\w s -> fromIntegral w < s)
いろいろな符号化
1から無限大までを符号化してみる
「1から無限大までを符号化する」やりかたを考えてみよう。[1, 無限大)という範囲を何で分割すればいいだろうか。1と無限大のあいだを分割するんだから、すごく大きい数でわけるのがよさそうだ。と思うかと思う。けど、ここでは2で分割することにする。そして[2, 無限大)は4で分割し、[4, 無限大)は8で分割する。2倍2倍にしていくと、すごいいきおいで大きくなっていくので、これでかまわないだろう。例として1110101をみてみよう。
[1, 無限大)
1: [2, 無限大)
1: [4, 無限大)
1: [8, 無限大)
0: [8, 16)
1: [12, 16)
0: [12, 14)
1: [13, 14)
13になった。上限が無限大ではない場合には平均(相加平均)で分割している。このような符号化には、つぎのような特徴がある。
- はじめの0をはさんで、前後のビット数が等しい
- はじめの1の数をnとして、後半のビットパターンを2進数として解釈した値をmとすると表現される値は
+ 2 ^ n + m
この符号化方式での、デコード、エンコードのコードは、つぎのようになる。
まず、共通して使われる「分割点を導出する関数」を定義する。
positiveIntegerSplitter :: Maybe Integer -> Maybe Integer -> Maybe Integer
positiveIntegerSplitter (Just mn) Nothing = Just $ mn * 2
positiveIntegerSplitter (Just mn) (Just mx) = Just $ (mn + mx) `div` 2
positiveIntegerSplitter _ _ = error "Oops!"
「無限大」をNothingで表現した。これを使って、それぞれの関数をつくる。
toPositiveInteger :: [Bit] -> Integer
toPositiveInteger = toValue
positiveIntegerSplitter
(\case Just mn -> mn; _ -> error "Oops!")
(Just 1) Nothing
fromPositiveInteger :: Integer -> [Bit]
fromPositiveInteger = fromValue Nothing
(\a b -> case (a, b) of
(Just mn, Just mx) -> mn + 1 >= mx
_ -> False)
positiveIntegerSplitter
(\w -> maybe True (w <))
(Just 1) Nothing
LambdaCase拡張をつかったので、先頭につぎのプラグマを追加する。
{-# LANGUAGE LambdaCase #-}
ビット列を文字列に変換する関数も追加しておく。
toString :: [Bit] -> String
toString = map $ \case O -> '0'; I -> '1'
試してみる。
> :load urr.hs
> toPositiveInteger $ fromString "1110101"
13
> toString $ fromPositiveInteger 10
"1110010"
> toString $ fromPositiveIntger 300
"11111111000101100"
負の無限大から正の無限大までの整数
1以上の整数ではなく負の値や0も含めた整数全体を符号化してみよう。まず0で負の区間と正の区間を分割し、それらを、それぞれ-1と1で分割するという操作が追加で必要になる。例をみてみよう。
[-Inf, Inf)
1: [0, Inf)
1: [1, Inf)
1: [2, Inf)
1: [4, Inf)
1: [8, Inf)
0: [8, 16)
1: [12, 16)
0: [12, 14)
1: [13, 14)
111110101は13になる。
[-Inf, Inf)
0: [-Inf, 0)
0: [-Inf, -1)
0: [-Inf, -2)
0: [-Inf, -4)
0: [-Inf, -8)
1: [-16, -8)
0: [-16, -12)
1: [-14, -12)
1: [-13, -12)
000001011は-13になる。デコーダ、エンコーダを書いてみよう。
integerSplitter :: Maybe Integer -> Maybe Integer -> Maybe Integer
integerSplitter Nothing Nothing = Just 0
integerSplitter Nothing (Just 0) = Just $ - 1
integerSplitter (Just 0) Nothing = Just 1
integerSplitter Noting (Just mx) = Just $ mx * 2
integerSplitter (Just mn) Nothing = Just $ mx * 2
integerSplitter (Just mn) (Just mx) = Just $ (mn + mx) `div` 2
toWholeInteger :: [Bit] -> Maybe Integer
toWholeInteger = toValue integerSplitter id Nothing Nothing
fromWholeInteger :: Integer -> [Bit]
fromWholeInteger = fromValue Nothing
(\a b -> case (a, b) of
(Just mn, Just mx) -> mn + 1 >= mx
_ -> False)
integerSplitter
(\n -> maybe (error "Oops!") (n <))
Nothing Nothing
試してみよう。
> :load urr.hs
> toWholeInteger $ fromString "111110101"
Just 13
> toWholeInteger $ fromString "000001011"
Just (-13)
> toString $ fromWholeInteger 18
"11111100010"
> toString $ fromWholeInteger (- 6)
"0000110"
浜田のURR
さて、ようやく本題だ。浜田のURRでは、つぎのように分割していく。
区間の分割
まずは前処理的な分割。これは前述の例での0や1での分割のように、前処理的な分割だ。
- -Infから+Infまでを0で分割
- -Infから0を-1で分割
- 0からInfを1で分割
- ここまでで4区間になる[-Inf, -1), [-1, 0), [0, 1), [1, +Inf)
- これらを、それぞれ-2, -0.5, 0.5, 2で分割する
- ここまでで8区間になる
- [-Inf, -2)
- [-2, -1)
- [-1, -0.5)
- [-0.5, 0)
- [0, 0.5)
- [0.5, 1)
- [1, 2)
- [2, +Inf)
この8区間について本番の分割をしていく。つぎのようなルールで分割する。
- かたほうが0または無限大のとき、もう一方の値をxとして
- xの値の絶対値を2 ^ mとしたときの2 ^ (2 * m)を絶対値とする
- xの値と符号がおなじ
- そのような値で分割する
- 上限(2^m)と下限(2^n)の、一方がもう一方のk倍(k > 2)であるとき
- 上限と下限の相乗平均で分割する
- つまり2^((m + n) / 2)で分割する
- 上限と下限の、一方がもう一方の2倍のとき
- 上限と下限の平均(相加平均)値で分割する
浜田のURRのデコーダ、エンコーダを作る
まずは、範囲の下限、上限をあらわすデータ型を作成する。
data NP = N | P deriving Show
fromNP :: Num a => NP -> a
fromNP N = - 1
fromNP P = 1
data UrrRange = NInf | Zero | PInf | Exp NP Int | Raw Rational deriving Show
getRational :: UrrRange -> Rational
getRational Zero = 0
getRational (Exp np n) = fromNP np * 2 ^^ n
getRational (Raw r) = r
getRational _ = error "Oops!"
lessThan :: (Ord a, Fractional a) => a -> UrrRange -> Bool
lessThan _ NInf = False
lessThan x Zero = x < 0
lessThan _ PInf = True
lessThan x (Exp np n) = x < fromNP np * 2 ^^ n
lessThan x (Raw r) = x < fromRational r
範囲をせばめていく過程で、途中までは指数で表現したほうが、あつかいやすいので、指数表現と有理数表現の両方で表現できるようにした。分割点を導出する関数を定義する。
urrSplitter :: UrrRange -> UrrRange -> UrrRange
urrSplitter NInf PInf = Zero
urrSplitter NInf Zero = Exp N 0
urrSplitter Zero PInf = Exp P 0
urrSplitter NInf (Exp N 0) = Exp N 1
urrSplitter (Exp N 0) Zero = Exp N (- 1)
urrSplitter Zero (Exp P 0) = Exp P (- 1)
urrSplitter (Exp P 0) PInf = Exp P 1
urrSplitter NInf (Exp N mx) = Exp N (2 * mx)
urrSplitter (Exp N mn) Zero = Exp N (2 * mn)
urrSplitter Zero (Exp P mx) = Exp P (2 * mx)
urrSplitter (Exp P mn) PInf = Exp P (2 * mn)
urrSplitter (Exp np mn) (Exp _np mx)
| mn + 1 < mx = Exp np ((mn + mx) `div` 2)
urrSplitter mn mx = Raw $ (getRational mn + getRational mx) / 2
うえから、つぎのようになっている。
- 全体を負と正に分割
- 負の区域を-1(Exp N 0)で分割
- 正の区域を1(Exp P 0)で分割
- -1より小さい区域を-2(Exp N 1)で分割
- [-1, 0)を-1/2(Exp N (- 1))で分割
- [0, 1)を1/2(Exp P (- 1))で分割
- 1以上の区域を2(Exp P 1)で分割
- [-Inf, -2^mx)の区域を-2^(2*mx)で分割
- [-2^mn, 0)の区域を-2^(2*mn)で分割
- [0, 2^mx)の区域を2^(2*mx)で分割
- [2^mn, Inf)の区域を2^(2*mn)で分割
- [2^mn, 2^mx)の区域を
- mnとmxの差が1より大のときのみ、2^((mn+mx)/2)で分割
- それ以外のときは[mn, mx)について(mn + mx)/2で分割
デコーダ、エンコーダは、つぎのようになる。
fromUrr :: Fractional a => [Bit] -> a
fromUrr = toValue urrSplitter (fromRational . getRational) NInf PInf
toUrr :: (Ord a, Fractional a) => Word -> a -> [Bit]
toUrr ln = fromValue (Just ln) (\_ _ -> False) urrSplitter lessThan NInf PInf
試してみる。
> :load urr.hs
> fromUrr $ fromString "10101000"
0.625
> fromUrr $ fromString "01011000"
-0.625
> toString $ toUrr 8 0.625
"10101000"
> toString $ toUrr 8 (- 0.625)
"01011000"
> toUrr 64 pi
[I,I,I,O,I,O,O,...,O,O,O,O]
> fromUrr it
3.141592653589793
> toUrr 64 6.0221409e23
[I,I,I,I,I,I,I,I,I,I,O,O,O,I,I,I,O,...O,O,O,I,O,I]
> fromUrr it
6.022140899999995e23
構造としては、つぎのように考えることができる。
33.5
1111100100001100
1: 符号ビット
1: 指数部が0以上
1: 指数部が1以上
110: 指数部の桁数が2(連続する1の数が2なので)
01: 指数部が5(101)。はじめの1は省略される
00001100: 仮数部。2進法で1.00001100。はじめの1は省略される
ここで、符号ビットは負のときに1であってほしいので、つぎのようにする。
flipBit :: Bit -> Bit
flipBit O = I
flipBit I = O
heading :: (a -> a) -> [a] -> [a]
heading f (x : xs) = f x : xs
heading _ _ = error "Oops!"
fromUrr :: Fractional a => [Bit] -> a
fromUrr = toValue urrSplitter (fromRational . getRational) NInf PInf
. heading flipBit
toUrr :: (Ord a, Fractional a) => Word -> a -> [Bit]
toUrr ln = heading flipBit
. fromValue (Just ln) (\_ _ -> False) urrSplitter lessThan NInf PInf
まとめ
思ったより、記事を書くのに手間がかかった。浜田のURRについて、数学的な定義をHaskellで記述できた。領域を2分したうえで、かたほうを選ぶという作業が、それぞれのビットに対応しているのが、おもしろい。正負をわける作業が、符号ビットに対応し、0や無限大との分割が指数部のビット数を決めるビット列に対応。相乗平均によって分割していく作業が「指数部」に、相加平均によって分割していく作業が「仮数部」に対応している。
浜田のURRでは、これだけではなく、負の無限大を符号なし無限大と考え、あらわせるなかで最大の絶対値をもつものを正負の無限大と解釈し、最小の絶対値をもつものを正負の0(+0と-0)として解釈するといった、追加の定義があるけれど、そこまで「きちんと」説明するまえに、僕は力つきることにします。