[ruby-dev:31312] faster Bignum#to_s

遠藤と申します。

Karatsuba の基数変換アルゴリズムで Bignum#to_s を高速化するパッチを
書きました。何千桁もあるような大きな Bignum の表示が速くなります。
(そういうことをしたいことがどのくらいあるのかはわかりませんが)

数値演算のアルゴリズムは私も苦手なのでバグがあるかもしれませんが、
たたき台にでもなればと思います。よろしければ採用をご検討ください。



$ ruby -e '100.times { puts "puts((#{ rand(50000) }**#{ rand(50000)
}).to_s(#{ rand(35) + 2 }))" }' > test_bignum.rb

$ head test_bignum.rb
puts((29421**21521).to_s(10))
puts((14497**43375).to_s(7))
puts((1991**14216).to_s(23))
puts((7915**41353).to_s(27))
puts((22736**33284).to_s(6))
puts((24354**25848).to_s(16))
puts((28541**7140).to_s(29))
puts((2618**23001).to_s(14))
puts((4934**23263).to_s(26))
puts((29565**10694).to_s(29))

$ time ./ruby.org test_bignum.rb > result.org.txt # パッチ前

real    16m58.310s
user    0m19.590s
sys     16m35.400s

$ time ./ruby.new test_bignum.rb > result.new.txt # パッチ後

real    2m47.719s
user    0m5.820s
sys     2m41.090s



Karatsuba の基数変換アルゴリズムは、私の理解では以下のようなものです。

普通のアルゴリズムでは、10 で割りまくって各除算での余りを文字列にします
(今の Ruby の Bignum#to_s の実装もそうなっています) 。このとき 1 回の
除算の計算量は割られる数の桁数に対して線形の時間がかかるため、全体では
O(n^2) の時間がかかります。

Karatsuba の基数変換アルゴリズムでは、10 の累乗数で、桁数をおよそ半分に
するような値で割りまくります。12345678 の例で考えると、

   12345678 (を 10000 で割る)
-> 1234 と 5678 (をそれぞれ 100 で割る)
-> 12 と 34 と 56 と 78 (をそれぞれ 10 で割る)
-> 1 と 2 と 3 と 4 と 5 と 6 と 7 と 8

となり、各数字が文字列となります。除算回数は普通のアルゴリズムと同じ
ですが、各除算での割られる数の桁数が減っているため、全体としては
O(n log n) の時間ですむようです。



Index: bignum.c
===================================================================
--- bignum.c    (revision 12854)
+++ bignum.c    (working copy)
@@ -595,11 +595,97 @@
 }

 const char ruby_digitmap[] = "0123456789abcdefghijklmnopqrstuvwxyz";
+static inline int
+big2str_normal(VALUE x, long j, int base, int hbase, char *s, int trim)
+{
+    long i = RBIGNUM(x)->len;
+    BDIGIT *ds = BDIGITS(x);
+
+    while (i && j > 1) {
+       long k = i;
+       BDIGIT_DBL num = 0;
+
+       while (k--) {
+           num = BIGUP(num) + ds[k];
+           ds[k] = (BDIGIT)(num / hbase);
+           num %= hbase;
+       }
+       if (trim && ds[i-1] == 0) i--;
+       k = SIZEOF_BDIGITS;
+       while (k--) {
+           s[--j] = ruby_digitmap[num % base];
+           num /= base;
+           if (!trim && j < 1) break;
+           if (trim && i == 0 && num == 0) break;
+       }
+    }
+    return j;
+}
+
+#define KARATSUBA_DIGITS 1024
+static VALUE big2str_table[37][16];
+
+static VALUE bigsqr(VALUE x);
+static void bigdivmod(VALUE x, VALUE y, VALUE *divp, VALUE *modp);
+
+static inline int
+big2str_karatsuba(VALUE x, int n, int base, int hbase, char *s, int trim)
+{
+    long i, j, k, bs_len, sign = RBIGNUM(x)->sign;
+    volatile VALUE t = big2str_table[base][0], t2;
+    VALUE *as, *bs, q, r;
+
+    j = RBIGNUM(t)->len;
+    for (i=0; j <= RBIGNUM(x)->len; i++) j *= 2;
+
+    as = ALLOCA_N(VALUE, i);
+
+    for (i=0,j=1; ; i++,j*=2) {
+       as[i] = t;
+       if(big2str_table[base][i + 1]) {
+           t2 = big2str_table[base][i + 1];
+       }
+       else {
+           t2 = bigsqr(t);
+           if(i + 1 < 16) {
+               big2str_table[base][i + 1] = t2;
+               rb_global_variable(&big2str_table[base][i + 1]);
+           }
+       }
+       if(RBIGNUM(x)->len < RBIGNUM(t2)->len) break;
+       t = t2;
+    }
+
+    bs_len = j * 2;
+    bs = ALLOCA_N(VALUE, bs_len);
+    bs[0] = x;
+    RBIGNUM(x)->sign = 1;
+    for (; j>0; i--, j/=2) {
+       for (k=0; k<bs_len; k+=j*2) {
+           bigdivmod(bs[k], as[i], &q, &r);
+           bs[k] = q;
+           bs[k + j] = r;
+       }
+    }
+    RBIGNUM(x)->sign = sign;
+
+    j = 0;
+    while(RBIGNUM(bs[j])->len == 1 && BDIGITS(bs[j])[0] == 0) j++;
+    for (i=bs_len-1; i>j; i--) {
+       k = big2str_normal(
+           bs[i], KARATSUBA_DIGITS, base, hbase,
+           s + n - KARATSUBA_DIGITS, Qfalse);
+       n -= KARATSUBA_DIGITS - k;
+    }
+    n = big2str_normal(bs[j], n, base, hbase, s, trim);
+
+    return n;
+}
+
 VALUE
 rb_big2str0(VALUE x, int base, int trim)
 {
     volatile VALUE t;
-    BDIGIT *ds;
     long i, j, hbase;
     VALUE ss;
     char *s;
@@ -646,29 +732,16 @@
 #endif

     t = rb_big_clone(x);
-    ds = BDIGITS(t);
     ss = rb_str_new(0, j+1);
     s = RSTRING_PTR(ss);

     s[0] = RBIGNUM(x)->sign ? '+' : '-';
-    while (i && j > 1) {
-       long k = i;
-       BDIGIT_DBL num = 0;
-
-       while (k--) {
-           num = BIGUP(num) + ds[k];
-           ds[k] = (BDIGIT)(num / hbase);
-           num %= hbase;
-       }
-       if (trim && ds[i-1] == 0) i--;
-       k = SIZEOF_BDIGITS;
-       while (k--) {
-           s[--j] = ruby_digitmap[num % base];
-           num /= base;
-           if (!trim && j < 1) break;
-           if (trim && i == 0 && num == 0) break;
-       }
+    if (RBIGNUM(x)->len > RBIGNUM(big2str_table[base][0])->len * 4) {
+       j = big2str_karatsuba(t, j, base, hbase, s, trim);
     }
+    else {
+       j = big2str_normal(t, j, base, hbase, s, trim);
+    }
     if (trim) {while (s[j] == '0') j++;}
     i = RSTRING_LEN(ss) - j;
     if (RBIGNUM(x)->sign) {
@@ -683,6 +756,22 @@
     return ss;
 }

+static void
+init_big2str_table(void)
+{
+    int i, j;
+    VALUE v;
+
+    for (i=0; i<37; i++) {
+       v = rb_big_pow(rb_int2big(i), INT2FIX(KARATSUBA_DIGITS));
+       big2str_table[i][0] = v;
+       rb_global_variable(&big2str_table[i][0]);
+       for (j=1; j<16; j++) {
+           big2str_table[i][j] = Qfalse;
+       }
+    }
+}
+
 VALUE
 rb_big2str(VALUE x, int base)
 {
@@ -2215,4 +2304,6 @@
     rb_define_method(rb_cBignum, "to_f", rb_big_to_f, 0);
     rb_define_method(rb_cBignum, "abs", rb_big_abs, 0);
     rb_define_method(rb_cBignum, "size", rb_big_size, 0);
+
+    init_big2str_table();
 }

-- 
Yusuke ENDOH <[e-mail]>