NumericのVALUE埋め込み/ポインタ参照の比較とRationalを用いた高精度な小数演算

今回ちょっと内容が多いからとりあえず目次というか概要というか。
主にこんな話をまとめてる。
VALUEとはいったい何なのか。
VALUE埋め込み表現出来るものとポインタ参照になるもの。
Ruby上でVALUE埋め込みなのかポインタなのかの調べ方。
VALUE埋め込み型とポインタ型の速度比較。
・Rationalで小数演算の精度を高める。
・HpSqrtの精度が高くなった。

VALUEとは

Rubyのオブジェクトを表すVALUEはunsigned longで、そのまま実体として値を持ったりヒープに確保されたメモリのポインタを表したり、持つ値によって意味合いが変わる。

VALUEがそのまま実体として値を持つものはruby_special_constsとして定義されている。
Qfalse,Qtrue,Qnil,Qundefはunsigned long値の一致、FIXNUM,FLONUM,SYMBOLは末尾ビットの一致で表される。
allocが返すポインタは必ず4の倍数である事を利用して、立っているビットを見れば埋め込みかポインタかを判断出来る仕組みになっている。

include/ruby/ruby.h

enum ruby_special_consts {
#if USE_FLONUM
    RUBY_Qfalse = 0x00,     /* ...0000 0000 */
    RUBY_Qtrue  = 0x14,     /* ...0001 0100 */
    RUBY_Qnil   = 0x08,     /* ...0000 1000 */
    RUBY_Qundef = 0x34,     /* ...0011 0100 */

    RUBY_IMMEDIATE_MASK = 0x07,
    RUBY_FIXNUM_FLAG    = 0x01, /* ...xxxx xxx1 */
    RUBY_FLONUM_MASK    = 0x03,
    RUBY_FLONUM_FLAG    = 0x02, /* ...xxxx xx10 */
    RUBY_SYMBOL_FLAG    = 0x0c, /* ...0000 1100 */
#else
    RUBY_Qfalse = 0,        /* ...0000 0000 */
    RUBY_Qtrue  = 2,        /* ...0000 0010 */
    RUBY_Qnil   = 4,        /* ...0000 0100 */
    RUBY_Qundef = 6,        /* ...0000 0110 */

    RUBY_IMMEDIATE_MASK = 0x03,
    RUBY_FIXNUM_FLAG    = 0x01, /* ...xxxx xxx1 */
    RUBY_FLONUM_MASK    = 0x00, /* any values ANDed with FLONUM_MASK cannot be FLONUM_FLAG */
    RUBY_FLONUM_FLAG    = 0x02,
    RUBY_SYMBOL_FLAG    = 0x0e, /* ...0000 1110 */
#endif
    RUBY_SPECIAL_SHIFT  = 8
};

詳しい事はここに書いてあるので割愛。
Ruby 1.7の頃の解説だからFLONUMに関する記述がなかったり若干古いけど、基本的な作りは変わっていない。
http://i.loveruby.net/ja/rhg/book/object.html

flagで使われるbit分、VALUE埋め込みで表せる表現の幅が狭くなる

C言語のlongよりもVALUE埋め込みのIntegerの方が1bit分表現の幅が狭い。
C言語のdoubleよりもVALUE埋め込みのFloatの方が2bit分表現の幅が狭い。

Integerが桁あふれした場合、内部ではstruct RBignumがallocされ、値はbignumとして保持される。
そのためC言語のlongの最大値を超えてメモリの許す限りの大きな数値を扱う事が出来る仕組みになっている。

Floatが桁あふれした場合、内部ではstruct RFloatがallocされ、値はdoubleとして保持される。
(ちなみにバージョンは定かではないが、確かRuby1.9辺りまではFloat値がVALUEに埋め込まれる事はなくすべてヒープメモリに置かれていた。)

なのでRubyのInteger/Floatの演算精度がCよりも劣るという事はない。
しかしVALUEが桁あふれした場合は以下を実行するためにオーバーヘッドが増える。
・オブジェクト生成時にallocによるヒープメモリ確保。
・オブジェクト参照時にポインタ参照。
GCによるヒープメモリ開放。

具体的にどんな値がVALUE埋め込み/ポインタ参照になるのか

internal.hのrb_float_new_inline(double)を元に、doubleからVALUEへの変換するプログラムを作った。

#include <stdio.h>

struct bits {
    unsigned char b1 : 1;
    unsigned char b2 : 1;
    unsigned char b3 : 1;
    unsigned char b4 : 1;
    unsigned char b5 : 1;
    unsigned char b6 : 1;
    unsigned char b7 : 1;
    unsigned char b8 : 1;
};
 
union dbbit {
    double d;
    unsigned long l;
    struct bits b[8];
};

#define RUBY_BIT_ROTL(v, n) (((v) << (n)) | ((v) >> ((sizeof(v) * 8) - n)))

void bits2chars(char *c, struct bits *ba) {
    struct bits *b;
    for (int i=7; 0<=i; i--) {
        b = &ba[i];
        sprintf(c, "%d%d%d%d%d%d%d%d ", b->b8, b->b7, b->b6, b->b5, b->b4, b->b3, b->b2, b->b1);
        c += 9;
    }
    //for (int i=0; i<=7; i++) {
    //    b = &ba[i];
    //    sprintf(c, "%d%d%d%d%d%d%d%d ", b->b1, b->b2, b->b3, b->b4, b->b5, b->b6, b->b7, b->b8);
    //    c += 9;
    //}
    c--;
    *c = '\0';
}

int main(void) {

    union dbbit u;
    struct bits* b;
    char bitschars[80];

    // heap
    u.l = 0;
    u.b[7].b7 = 1;
    u.b[7].b5 = 1;

    // VALUE
    //u.d = 123.0;

    int bits = (int)((u.l >> 60) & 0x7);
    /* bits contains 3 bits of b62..b60. */
    /* bits - 3 = */
    /*   b011 -> b000 */
    /*   b100 -> b001 */

    bits2chars(bitschars, u.b);
    printf("before\n");
    printf("  double: %lf\n", u.d);
    printf("  ulong : %ld\n", u.l);
    printf("  bits  : %s\n", bitschars);
    printf("  flags : %d\n\n", bits);

    printf("after\n");

    if (u.l != 0x3000000000000000 && !((bits-3) & ~0x01)) {
        u.l = (RUBY_BIT_ROTL(u.l, 3) & ~(unsigned long)0x01) | 0x02;

        bits2chars(bitschars, u.b);
        printf("  double: %lf\n", u.d);
        printf("  ulong : %ld\n", u.l);
        printf("  bits  : %s\n", bitschars);
    } else if (u.l==0) {
        printf("  0x8000000000000002\n");
    } else {
        printf("  heap\n");
    }

    return 0;
}

倍精度浮動小数点数のビットパターンはWikipediaを参照。
https://ja.wikipedia.org/wiki/%E5%80%8D%E7%B2%BE%E5%BA%A6%E6%B5%AE%E5%8B%95%E5%B0%8F%E6%95%B0%E7%82%B9%E6%95%B0

doubleをVALUEで表現出来る値

C言語のdouble値123.0をRubyVALUE表現へ変換した結果、RUBY_FLONUM_FLAG(末尾2つのbit 10)が立っている事が確認出来る。

before
  double: 123.000000
  ulong : 4638355772470722560
  bits  : 01000000 01011110 11000000 00000000 00000000 00000000 00000000 00000000
  flags : 4

after
  double: 0.000000
  ulong : 213358032346677250
  bits  : 00000010 11110110 00000000 00000000 00000000 00000000 00000000 00000010

irbで123.0のobject_idを調べると上のプログラムの変換後のulongと一致している事がわかる。
VALUE埋め込みの値なので何度やってもobject_idが変わる事はない。

2.7.0-preview1 :001 > 123.000000.object_id
 => 213358032346677250 
2.7.0-preview1 :002 > 123.000000.object_id
 => 213358032346677250 
doubleをVALUEで表現出来ない値

bit表記でx101、x110、x111から始まるものは桁あふれしてVALUEに埋め込めないのでヒープに置かれる。

before
  double: 231584178474632390847141970017375815706539969331281128078915168015826259279872.000000
  ulong : 5764607523034234880
  bits  : 01010000 00000000 00000000 00000000 00000000 00000000 00000000 00000000
  flags : 5

after
  heap

上の結果でVALUE埋め込み表現出来なかったdouble値をRubyで確認する。
irbでobject_idを調べると毎回変わるのでヒープに置かれているオブジェクトである事がわかる。

2.7.0-preview1 :001 > 231584178474632390847141970017375815706539969331281128078915168015826259279872.000000.object_id
 => 70246480427780 
2.7.0-preview1 :002 > 231584178474632390847141970017375815706539969331281128078915168015826259279872.000000.object_id
 => 70246464940200 
Ruby上でVALUE埋め込みなのかポインタなのか調べる

以下Rubyプログラムでobject_idからVALUE埋め込み型なのかヒープに置かれているのか確認できる。
今回はSymbolについては割愛。

class Object
  RUBY_Qfalse = 0x00      # ...0000 0000
  RUBY_Qtrue  = 0x14      # ...0001 0100
  RUBY_Qnil   = 0x08      # ...0000 1000
  RUBY_Qundef = 0x34      # ...0011 0100

  RUBY_IMMEDIATE_MASK = 0x07
  RUBY_FIXNUM_FLAG    = 0x01  # ...xxxx xxx1
  RUBY_FLONUM_MASK    = 0x03
  RUBY_FLONUM_FLAG    = 0x02  # ...xxxx xx10
  RUBY_SYMBOL_FLAG    = 0x0c  # ...0000 1100
  RUBY_SPECIAL_SHIFT  = 8

  def memory_type
    case object_id
    when RUBY_Qfalse
      "Qfalse"
    when RUBY_Qtrue
      "Qtrue"
    when RUBY_Qnil
      "Qnil"
    when RUBY_Qundef
      "Qundef" # Rubyからは確認出来ない
    else
      if (object_id & RUBY_FIXNUM_FLAG)==RUBY_FIXNUM_FLAG
        "stack integer"
      elsif (object_id & RUBY_FLONUM_MASK)==RUBY_FLONUM_FLAG
        "stack flonum"
      else
        "heap or symbol"
      end
    end
  end
end

p false.memory_type
p true.memory_type
p nil.memory_type
p 123.memory_type
p 0xFFFFFFFFFFFFFFFF.memory_type
p 12.3.memory_type
p 231584178474632390847141970017375815706539969331281128078915168015826259279872.000000.memory_type
p "hello".memory_type

実行結果。

"Qfalse"
"Qtrue"
"Qnil"
"stack integer"
"heap or symbol"
"stack flonum"
"heap or symbol"
"heap or symbol"

実行速度の差

IntegerのVALUE埋め込みとヒープメモリでの速度比較。

require 'benchmark'

repeat = 1000
n = 1000

stack_time = Benchmark.realtime do
  n.times {
    repeat.times.map{|i| 1}.inject(:+)
  }
end

heap_time = Benchmark.realtime do
  n.times {
    repeat.times.map{|i| 0xFFFFFFFFFFFFFFFF}.inject(:+)
  }
end

puts "Stack time : #{stack_time}"
puts "Heap time  : #{heap_time}"
puts "Ratio      : #{heap_time / stack_time}"

速度比較結果。

Stack time : 0.07842599996365607
Heap time  : 0.13512099999934435
Ratio      : 1.7229107701777688

FloatのVALUE埋め込みとヒープメモリでの速度比較。

require 'benchmark'

repeat = 1000
n = 1000

stack_time = Benchmark.realtime do
  n.times {
    repeat.times.map{|i| 123.0}.inject(:+)
  }
end

heap_time = Benchmark.realtime do
  n.times {
    repeat.times.map{|i| 231584178474632390847141970017375815706539969331281128078915168015826259279872.000000}.inject(:+)
  }
end

puts "Stack time : #{stack_time}"
puts "Heap time  : #{heap_time}"
puts "Ratio      : #{heap_time / stack_time}"

速度比較結果。
bignumと違い固定長で中身doubleのため、そこまで極端に遅くはならない。

Stack time : 0.11943600000813603
Heap time  : 0.13702500006183982
Ratio      : 1.1472671560710808

小数演算の精度をもっと高めたい

Float同士の計算は丸め誤差がよく生まれる。
以下の2の平方根の1000倍を求めるプログラムでは、1000を掛けた場合は誤差がなく、1000回足した場合は誤差が生じている。

2.7.0-preview1 :001 > Math.sqrt(2)
 => 1.4142135623730951 
2.7.0-preview1 :002 > Math.sqrt(2)*1000
 => 1414.213562373095 
2.7.0-preview1 :003 > 1000.times.map{|i| Math.sqrt(2)}.inject(:+)
 => 1414.213562373105 

一般的に言えばBigDecimalを使えって話になるが、BigDecimal有理数虚数に対応していない。
平方根を計算する場合は虚数も考慮しなければいけない。
まずはRationalとして計算してみる。
以下の結果から、Rationalに変換してから計算するとBigDecimalや掛け算と同じ結果が得られる事がわかる。

2.7.0-preview1 :016 > f_add = 1000.times.map{|i| Math.sqrt(2).to_r}.inject(:+).to_f
 => 1414.213562373095 
2.7.0-preview1 :017 > 1000.times.map{|i| BigDecimal(Math.sqrt(2).to_s)}.inject(:+).to_f
 => 1414.213562373095 

Rationalは分母VALUE、分子VALUEを持つ。
つまり分母と分子にBignum値を持つ事が出来る。
そして虚数にも対応するためにComplexとする。
Complex(Rational(Bignum, Bignum), Rational(Bignum, Bignum))とすれば、実数/虚数ともに小数値を可変長精度で表せる事になる。
分子Bignum分母Bignumの有理数を実数と虚数に持つ複素数同士の足し算をする事で、有理数にも虚数にも対応した高精度な小数演算を実現する。

2.7.0-preview1 :005 > 1000.times.map{|i| Math.sqrt(Complex(1, 1))}.inject(:+)
 => (1098.6841134678039+455.08986056222534i) 
2.7.0-preview1 :006 > c = 1000.times.map{|i| c = Math.sqrt(Complex(1, 1)); Complex(c.real.to_r, c.imag.to_r)}.inject(:+)
 => ((618504170501439125/562949953421312)+(1024771263224069125/2251799813685248)*i) 
2.7.0-preview1 :007 > Complex(c.real.to_f, c.imag.to_f)
 => (1098.68411346781+455.0898605622274i) 
Rationalを内包するComplexの速度

Floatの足し算とRationalを内包するComplexの足し算の速度比較。

require 'benchmark'

repeat = 1000
n = 1000

f_mul = Math.sqrt(2) * repeat
r_mul = Math.sqrt(2).to_r * repeat

f_add = nil
r_add = nil

f_time = Benchmark.realtime do
  n.times {
    f_add = repeat.times.map{|i| Math.sqrt(2)}.inject(:+)
  }
end

r_time = Benchmark.realtime do
  n.times {
    r_add = repeat.times.map{|i| Math.sqrt(2).to_r}.inject(:+)
  }
end

puts "Float root 2        : #{Math.sqrt(2)}"
puts

puts "Float mul           : #{f_mul}"
puts "Rational mul        : #{r_mul.to_f}"
puts "Diff                : #{f_mul - r_mul}"
puts

puts "Float repeat add    : #{f_add}"
puts "Rational repeat add : #{r_add.to_f}"
puts "Diff                : #{f_add - r_add}"
puts

puts "Float time          : #{f_time}"
puts "Rational time       : #{r_time}"
puts "Ratio               : #{r_time / f_time}"

結果。
約6.35倍かかるようになったけど誤差はなくなった。

Float root 2        : 1.4142135623730951

Float mul           : 1414.213562373095
Rational mul        : 1414.213562373095
Diff                : 0.0

Float repeat add    : 1414.213562373105
Rational repeat add : 1414.213562373095
Diff                : 1.000444171950221e-11

Float time          : 0.1697259999345988
Rational time       : 1.0777420001104474
Ratio               : 6.349893360626763

HpSqrtでRationalを内包するComplexを採用

高精度平方根演算ライブラリHpSqrtでRationalを内包するComplexを採用した。
上記のベンチマークは一番差が出る極端な例で、この変更前後でのTravisCIのユニットテスト実行時間(assertions/s)を見ると1.22倍程度となっている。
バージョンによってassertion数も増えているので一概には比較出来ないが、テストで書いたような使い方をするなら処理速度に問題はないかなと。

1.8.0
Finished in 0.002686s, 2978.0063 runs/s, 18612.5394 assertions/s.

1.10.0
Finished in 0.004189s, 2148.2426 runs/s, 15276.3917 assertions/s.

計算速度は遅くなったけど、その分精度は上がった。
変更前も掛け算の精度は高かったけど、割り算・足し算・引き算の精度は高くなかった。(と言ってもRuby標準と同じ精度の意)
今回の更新で四則演算すべてでRuby標準のものよりも高精度に演算出来るようになった。

遅くなった理由と精度は確かだという事を言いたいだけで書いた技術記事なんで、これ読んだ人は高精度なHpSqrtを使ったらいいと思う。

GitHub: https://github.com/yoshida-eth0/ruby-sqrt
RubyGems: https://rubygems.org/gems/hpsqrt

環境

$ ruby -v
ruby 2.7.0preview1 (2019-05-31 trunk c55db6aa271df4a689dc8eb0039c929bf6ed43ff) [x86_64-darwin18]

$ gem list | grep hpsqrt
hpsqrt (1.10.0)