有効数字で浮動小数点数を丸める

浮動小数点数の時系列データを圧縮するのに xor encoding という手法があるそうだ。

• Gorilla の論文 http://www.vldb.org/pvldb/vol8/p1816-teller.pdf
• 時系列データベースに関する基礎知識と時系列データの符号化方式について - クックパッド開発者ブログ

値が時系列に 15.5, 14.0625, 3.25, 8.625 というように並んでいたとする。

それぞれを倍精度浮動小数点で表すと、ビット列は下のようになる。

decimail	double precision
15.5	01000000 00101111 00000000 00000000 00000000 00000000 00000000 00000000
14.0625	01000000 00101100 00100000 00000000 00000000 00000000 00000000 00000000
3.25	01000000 00001010 00000000 00000000 00000000 00000000 00000000 00000000
8.625	01000000 00100001 01000000 00000000 00000000 00000000 00000000 00000000

時系列データというのは、同じ値が続いたり、前後の値が近いという特徴がある。そのとき、近い数字はビット列の並びも近いとすると、前の値ビット列と XOR を取ると、ビットが立っているところが集中するようになる。

00000000 00000011 00100000 00000000 00000000 00000000 00000000 00000000
00000000 00100110 00100000 00000000 00000000 00000000 00000000 00000000
00000000 00101011 01000000 00000000 00000000 00000000 00000000 00000000

ビットが立っているところのデータだけ記録すれば、保存に必要なディスク容量が少なくなるだろうというアイデアです。ビットが立っているところだけ記録する方法は Gorilla の論文を参照してください。

しかし、この手法は本当にいいのだろうか？

浮動小数点数のバイナリデータをダラダラと吐き出すシステムもないわけではないのだろうが、ほとんどの場合10進数のテキストデータが渡されるのではなかろうか。10進数を浮動小数点数に厳密に変換することはできません。例えば、 3.14 を倍精度浮動小数点にしたときのビット列は

01000000 00001001 00011110 10111000 01010001 11101011 10000101 00011111

となります。これを10進数に戻すと 3.14000000000000012434… となります。

3.13 を倍精度浮動小数点にしたら、

01000000 00001001 00001010 00111101 01110000 10100011 11010111 00001010

となるのですが、3.14 と XOR をとったら

00000000 00000000 00010100 10000101 00100001 01000000 01010010 00010101

です。2進数でキリが悪い数字となっているので、値は近いにも関わらず、ビット列はそこまで一致していません。

そもそも、3.14, 3.13 は有効数字は３桁であり、倍精度浮動小数点は10進数で約16桁の精度なので、10進数をそのまま浮動小数点数に直すと無効な桁も表現してしまっていることになります。有効数字だけ残して、無効な桁を表現するビットを丸めてしまえば、もっと効果的にデータを圧縮できそうです。

3.14 の場合は、下のように有効数字を示す右側のビットを落としてしまってよいでしょう。

丸めなし 01000000 00001001 00011110 10111000 01010001 11101011 10000101 00011111
丸めあり 01000000 00001001 00100000 00000000 00000000 00000000 00000000 00000000

有効数字３桁で丸めた浮動小数点を10進数に戻すと、 3.140625 となります。これが２進数ではキリのよい数ということです。記事の最後にこのような丸めを実現するプログラム例を載せておきます。

テキストの 10進数 → 有効数字で丸めた浮動小数点 → xor encoding

という風に変換していけば有効数字を維持したまま時系列データを圧縮できるのではないかと考えました。

実際にやってみました。Gorilla の論文に載っている符号化は実装が面倒すぎたので、代わりに前の値との xor 値を並べたバイト列を圧縮してどの程度効果があるのか試した。

sar でとった 1秒おきの PC の CPU 使用率のデータで試してみたら、以下のような結果になりました。

	サイズ [bytes]
もとのテキストデータ	70062
倍精度浮動小数点数に変換後 XOR 値をとって gzip 圧縮	29357
有効数字で丸めた浮動小数点数に変換後 XOR 値をとって gzip 圧縮	26887

一応、効果はあるようです。しかし、普通にテキストファイルを gzip 圧縮したら 18431 bytes となって、最も有効だった。PC の CPU の使用率は変動が激しすぎて、そもそも xor encoding が有効でないのだろう。滑らかな変動をするようなデータで試した方が良かったのかもしれない。

---

import argparse
import math
import re
import struct
import sys


class DoublePrecisionFloat:
    def __init__(self):
        self.bits = None

    @staticmethod
    def from_float(f):
        instance = DoublePrecisionFloat()
        instance.bits = struct.unpack('>Q', struct.pack('>d', f))[0]
        return instance

    @staticmethod
    def from_string(string):
        f = float(string)
        return DoublePrecisionFloat.from_float(f)

    @staticmethod
    def from_bits(bits):
        instance = DoublePrecisionFloat()
        instance.bits = bits
        return instance

    @property
    def raw_float(self):
        return struct.unpack('>d', self.bits.to_bytes(8, byteorder='big'))[0]

    def to_binstr(self):
        bb = self.bits.to_bytes(8, byteorder='big')
        return ' '.join(['{:08b}'.format(b) for b in bb])

    def to_hexstr(self):
        bb = self.bits.to_bytes(8, byteorder='big')
        hex2bytes = ['{:02x}{:02x}'.format(bb[i], bb[i + 1])
                     for i in range(0, len(bb), 2)]
        return ' '.join(hex2bytes)

    @staticmethod
    def __significant_bits(significant_figures, base):
        return int(math.ceil(
            significant_figures * math.log2(base)))

    # 11111111 11110000 00000000 ... 00000000
    # |<----------><------  fraction  ------>
    # |  exponent
    # sign
    __FILTER = int.from_bytes(
        b'\xFF\xF0\x00\x00\x00\x00\x00\x00',
        byteorder='big',
        signed=True)

    def round_off(self, significant_figures, base=10):
        # 0         0.1         1
        # |<-------------------||
        #        round off
        sig_bits = self.__significant_bits(significant_figures, base)
        # left side is padded with 1
        filt = self.__FILTER >> sig_bits
        new_bits = self.bits & filt
        return DoublePrecisionFloat.from_bits(new_bits)

    def round_up(self, significant_figures, base=10):
        # 0         0.1         1
        # ||------------------->|
        #        round up

        sig_bits = self.__significant_bits(significant_figures, base)
        # left side is padded with 1
        filt = self.__FILTER >> sig_bits

        # TODO: the case of overflow
        new_bits = (self.bits + (~filt)) & filt
        return DoublePrecisionFloat.from_bits(new_bits)

    def round(self, significant_figures, base=10):
        # 0         0.1         1
        # |<--------||--------->|
        #  round off   round up

        sig_bits = self.__significant_bits(significant_figures, base)
        filt = self.__FILTER >> sig_bits
        round_up_bit = 0x0010000000000000 >> (sig_bits + 1)
        new_bits = (self.bits + round_up_bit) & filt
        return DoublePrecisionFloat.from_bits(new_bits)

    def xor(self, other):
        return DoublePrecisionFloat.from_bits(self.bits ^ other.bits)


class FloatString:
    def __init__(self, string):
        self.string = string

    @property
    def significant_figures(self):
        # TODO: This code is ugly

        string = self.string

        # 0.0
        f = float(string)
        if f == 0.0:
            return 0

        # 0.0123
        pattern1 = r'[+-]?0\.0*(\d+)'
        match1 = re.search(pattern1, string)
        if match1:
            return len(match1.group(1))

        # 123.4 or 123.40
        pattern2 = r'[+-]?0*([1-9]+\d*\.\d+)'
        match2 = re.search(pattern2, string)
        if match2:
            return len(match2.group(1)) - 1

        # +1E+10 or -1e-10
        pattern3 = r'[+-]?(\d+)[eE][+-]\d+'
        match3 = re.search(pattern3, string)
        if match3:
            return len(match3.group(1))

        # +1.23E+10 or -1.23e-10
        pattern4 = r'[+-]?0*(\d\.\d+)[eE][+-]\d+'
        match4 = re.search(pattern4, string)
        if match4:
            return len(match4.group(1)) - 1

        # 1234
        pattern5 = r'[+-]?0*(\d+)'
        match5 = re.search(pattern5, string)
        if match5:
            return len(match5.group(1))

        raise ValueError("could not get significant figures: '%s'" % string)

    def to_rounded_float(self):
        d = DoublePrecisionFloat.from_string(self.string)
        sigfig = self.significant_figures
        return d.round(sigfig)

    def to_nonrounded_float(self):
        d = DoublePrecisionFloat.from_string(self.string)
        return d


def main():
    parser = argparse.ArgumentParser()
    out_formats = parser.add_mutually_exclusive_group()
    out_formats.add_argument('-b', action='store_true',
                             help='print results with binary strings')
    out_formats.add_argument('-x', action='store_true',
                             help='print results with hex strings')
    out_formats.add_argument('-f', action='store_true',
                             help='print results with decimal number')
    out_formats.add_argument('-n', action='store_true',
                             help='not round results')
    args = parser.parse_args()

    if args.b is True:
        write_func = lambda f: print(f.to_binstr())
    elif args.x is True:
        write_func = lambda f: print(f.to_hexstr())
    elif args.f is True:
        write_func = lambda f: print(str(f.raw_float))
    else:
        write_func = lambda f: sys.stdout.buffer.write(
            struct.pack('>d', f.raw_float))

    def stdin():
        try:
            while True:
                yield input().strip()
        except EOFError:
                pass

    for line in stdin():
        if args.n is not True:
            rounded = FloatString(line).to_rounded_float()
        else:
            rounded = FloatString(line).to_nonrounded_float()
        write_func(rounded)

if __name__ == '__main__':
    main()