ruby-dev

むらけんです．

まつもとさんの日記で話題になっていた Bignum#to_s の
Karatsuba 基数変換による高速化を C で実装してみました．
ベキ数のキャッシュアルゴリズムが少し変な気がしますが，パッチを添付します．

一応，添付した test_karatsuba.rb でテストもして，きちんと動いているようです．

処理速度的には，べき乗の演算時間が含まれていますが 1/4 以下になりました．
$ time ./ruby -e 'p((65536**65536).to_s_orig(10))' > r1.txt

real    4m6.110s
user    1m1.368s
sys     0m0.016s
$ time ./ruby -e 'p((65536**65536).to_s(10))' > r2.txt

real    0m41.312s
user    0m11.421s
sys     0m0.012s

-- 
Kenta Murata
OpenPGP FP = FA26 35D7 4F98 3498 0810 E0D5 F213 966F E9EB 0BCC

Attachments (2)

bignum_to_s_karatsuba.patch (8.25 KB, text/x-diff)

Index: bignum.c
===================================================================
--- bignum.c	(revision 12859)
+++ bignum.c	(working copy)
@@ -19,6 +19,8 @@
 #include <ieeefp.h>
 #endif
 
+#include <st.h>
+
 VALUE rb_cBignum;
 
 #if defined __MINGW32__
@@ -727,12 +729,283 @@
     return ss;
 }
 
+#define MAX_POWER_CACHE_ENTRIES 64
+
+struct power_cache_data
+{
+    VALUE power;
+    long  magic;
+};
+
+struct power_cache_entry
+{
+    st_table* cache;
+    long      len;
+};
+
+struct power_cache
+{
+    struct power_cache_entry entries[35];
+};
+
+static int
+mark_cache_entry(long n1, struct power_cache_data* data, void* _)
+{
+    rb_gc_mark(data->power);
+    return ST_CONTINUE;
+}
+
+static void
+power_cache_mark(struct power_cache* pc)
+{
+    int i;
+    for (i = 0; i < 35; ++i)
+        st_foreach(pc->entries[i].cache, mark_cache_entry, 0);
+}
+
+static int
+free_cache_entry(long n1, struct power_cache_data* data, void* _)
+{
+    if (data != 0)
+        free(data);
+    return ST_DELETE;
+}
+
+static void
+power_cache_free(struct power_cache* pc)
+{
+    int i;
+    for (i = 0; i < 35; ++i) {
+        st_foreach(pc->entries[i].cache, free_cache_entry, 0);
+        st_free_table(pc->entries[i].cache);
+    }
+    free(pc);
+}
+
+static VALUE big2str_power_cache = 0;
+
+static void
+power_cache_init(void)
+{
+    struct power_cache* pc;
+    int i;
+
+    big2str_power_cache = Data_Make_Struct(
+      rb_cData, struct power_cache, power_cache_mark, power_cache_free, pc);
+    rb_global_variable(&big2str_power_cache);
+    for (i = 0; i < 35; i++) {
+        pc->entries[i].cache = st_init_numtable();
+        pc->entries[i].len = 0;
+    }
+}
+
+struct power_cache_find_to_remove_data
+{
+    long n1;
+    long max_magic;
+};
+
+static int
+power_cache_find_to_remove(long n1, struct power_cache_data* pcd,
+                           struct power_cache_find_to_remove_data* data)
+{
+    if (pcd->magic > data->max_magic) {
+        data->n1 = n1;
+        data->max_magic = pcd->magic++;
+    }
+    return ST_CONTINUE;
+}
+
+static VALUE
+power_cache_lookup(base, n1)
+    int base;
+    long n1;
+{
+    struct power_cache* pc;
+    struct power_cache_data* pcd = 0;
+
+    if (big2str_power_cache == 0) power_cache_init();
+    Data_Get_Struct(big2str_power_cache, struct power_cache, pc);
+    if (pc->entries[base - 2].len == 0 ||
+        !st_lookup(pc->entries[base - 2].cache, n1, (st_data_t*)&pcd))
+        return Qnil;
+    return pcd->power;
+}
+
+static void
+power_cache_insert(base, n1, power)
+    int base;
+    long n1;
+    VALUE power;
+{
+    struct power_cache* pc;
+    struct power_cache_data* pcd;
+
+    Data_Get_Struct(big2str_power_cache, struct power_cache, pc);
+    if (pc->entries[base - 2].len == MAX_POWER_CACHE_ENTRIES) {
+        struct power_cache_find_to_remove_data data = { 0, 0 };
+        st_foreach(pc->entries[base - 2].cache,
+                   power_cache_find_to_remove, (st_data_t)&data);
+        st_delete(pc->entries[base - 2].cache, &data.n1, 0);
+        pc->entries[base - 2].len--;
+    }
+
+    pcd = ALLOC(struct power_cache_data);
+    pcd->power = power;
+    pcd->magic = 0;
+    st_insert(pc->entries[base - 2].cache, n1, (st_data_t)pcd);
+    pc->entries[base - 2].len++;
+}
+
+static VALUE
+power_cache_get_power(base, n1)
+    int base;
+    long n1;
+{
+    VALUE power = power_cache_lookup(base, n1);
+    if (NIL_P(power)) {
+        power = rb_big_pow(rb_int2big(base), LONG2NUM(n1));
+        power_cache_insert(base, n1, power);
+    }
+    return power;
+}
+
+/* big2str_karatsuba_find_n1
+ *
+ * Let a natural number x is given by:
+ * x = 2^0 * x_0 + 2^1 * x_1 + ... + 2^(B*n_0 - 1) * x_{B*n_0 - 1},
+ * where B is BITSPERDIG (i.e. BDIGITS*CHAR_BIT) and n_0 is
+ * RBIGNUM(x)->len.
+ *
+ * Now, we assume n_1 = min_n \{ n | 2^(B*n_0/2) <= b_1^(n_1) \}, so
+ * it is realized that 2^(B*n_0) <= {b_1}^{2*n_1}, where b_1 is a
+ * given radix number. And then, we have n_1 <= (B*n_0) /
+ * (2*log_2(b_1)), therefore n_1 is given by ceil((B*n_0) /
+ * (2*log_2(b_1))).
+ */
+static long
+big2str_karatsuba_find_n1(VALUE x, int base)
+{
+    static const double log_2[] = {
+        1.0,              1.58496250072116, 2.0,
+        2.32192809488736, 2.584962500721,   2.58496250072116,
+        2.8073549220576,  3.0,              3.16992500144231,
+        3.32192809488736, 3.4594316186373,  3.58496250072116,
+        3.70043971814109, 3.8073549220576,  3.90689059560852,
+        4.0,              4.08746284125034, 4.16992500144231,
+        4.24792751344359, 4.32192809488736, 4.39231742277876,
+        4.4594316186373,  4.52356195605701, 4.58496250072116,
+        4.64385618977472, 4.70043971814109, 4.75488750216347,
+        4.8073549220576,  4.85798099512757, 4.90689059560852,
+        4.95419631038688, 5.0,              5.04439411935845,
+        5.08746284125034, 5.12928301694497, 5.16992500144231
+    };
+    long bits;
+
+    if (base < 2 && 36 < base)
+      rb_bug("illegal radix %d", base);
+
+    if (FIXNUM_P(x)) {
+        bits = (SIZEOF_LONG*CHAR_BIT - 1)/2 + 1;
+    }
+    else if (BIGZEROP(x)) {
+        return 0;
+    }
+    else {
+        bits = BITSPERDIG*RBIGNUM(x)->len;
+    }
+
+    return (long)(bits/(2*log_2[base - 2]) + 1);
+}
+
+static long
+big2str_karatsuba_without_sign(x, base, ptr, n1)
+    VALUE x;
+    int base;
+    char* ptr;
+    long n1;
+{
+    long len = 2*n1;
+    VALUE b, qr;
+
+    if (FIXNUM_P(x)) {
+        VALUE str = rb_fix2str(x, base);
+        memset(ptr, '0', len - RSTRING(str)->len);
+        MEMCPY(ptr + len - RSTRING(str)->len, RSTRING(str)->ptr,
+               char, RSTRING(str)->len);
+        return;
+    }
+    if (BIGZEROP(x)) {
+        memset(ptr, '0', len);
+        return;
+    }
+
+    if (n1 % 2 == 1) ++n1;
+    b = power_cache_get_power(base, n1);
+    qr = rb_big_divmod(x, b);
+    big2str_karatsuba_without_sign(
+        RARRAY(qr)->ptr[0], base, ptr,              (len - n1)/2);
+    big2str_karatsuba_without_sign(
+        RARRAY(qr)->ptr[1], base, ptr + (len - n1), n1/2);
+}
+
 VALUE
+rb_big2str0_karatsuba(x, base, trim)
+    VALUE x;
+    int base;
+    int trim;
+{
+    VALUE ss, xx;
+    long n1;
+
+    if (FIXNUM_P(x)) {
+        return rb_fix2str(x, base);
+    }
+    if (BIGZEROP(x)) {
+        return rb_str_new2("0");
+    }
+
+    if (base < 2 && 36 < base)
+        rb_raise(rb_eArgError, "illegal radix %d", base);
+
+    xx = rb_big_clone(x);
+    n1 = big2str_karatsuba_find_n1(xx, base);
+    ss = rb_str_new(0, 2*n1 + 1); /* plus one for sign */
+    RSTRING(ss)->ptr[0] = RBIGNUM(xx)->sign ? '+' : '-';
+    RBIGNUM(xx)->sign = 1;
+
+    big2str_karatsuba_without_sign(
+        xx, base, RSTRING(ss)->ptr + 1, n1);
+    if (trim) {
+        int i = 1;
+        while (RSTRING(ss)->ptr[i] == '0') ++i;
+        if (RBIGNUM(x)->sign) {
+            MEMMOVE(RSTRING(ss)->ptr, RSTRING(ss)->ptr + i,
+                    char, RSTRING(ss)->len - i); /* plus one for '\0' */
+        }
+        else {
+            MEMMOVE(RSTRING(ss)->ptr + 1, RSTRING(ss)->ptr + i,
+                    char, RSTRING(ss)->len - i); /* plus one for '\0' */
+        }
+        RSTRING(ss)->len -= i - !RBIGNUM(x)->sign;
+    }
+    RSTRING(ss)->ptr[RSTRING(ss)->len] = '\0';
+
+    return ss;
+}
+
+VALUE
 rb_big2str(VALUE x, int base)
 {
     return rb_big2str0(x, base, Qtrue);
 }
 
+VALUE
+rb_big2str_karatsuba(VALUE x, VALUE base)
+{
+    return rb_big2str0_karatsuba(x, FIX2INT(base), Qtrue);
+}
+
 /*
  *  call-seq:
  *     big.to_s(base=10)   =>  string
@@ -748,7 +1021,7 @@
  */
 
 static VALUE
-rb_big_to_s(argc, argv, x)
+rb_big_to_s_orig(argc, argv, x)
     int argc;
     VALUE *argv;
     VALUE x;
@@ -762,6 +1035,21 @@
     return rb_big2str(x, base);
 }
 
+static VALUE
+rb_big_to_s(argc, argv, x)
+    int argc;
+    VALUE *argv;
+    VALUE x;
+{
+    VALUE b;
+    int base;
+
+    rb_scan_args(argc, argv, "01", &b);
+    if (argc == 0) base = 10;
+    else base = NUM2INT(b);
+    return rb_big2str_karatsuba(x, base);
+}
+
 static unsigned long
 big2ulong(x, type)
     VALUE x;
@@ -2314,6 +2602,7 @@
 {
     rb_cBignum = rb_define_class("Bignum", rb_cInteger);
 
+    rb_define_method(rb_cBignum, "to_s_orig", rb_big_to_s_orig, -1);
     rb_define_method(rb_cBignum, "to_s", rb_big_to_s, -1);
     rb_define_method(rb_cBignum, "coerce", rb_big_coerce, 1);
     rb_define_method(rb_cBignum, "-@", rb_big_uminus, 0);

test_karatsuba.rb (855 Bytes, text/x-ruby)

require 'test/unit'

class TestBignum < Test::Unit::TestCase
  def setup
    srand(Time.now.to_i)
  end

  def test_to_s
    100.times do
      x = rand(2<<1024 - 1)
      (2..36).to_a.each do |b|
        assert_equal(x.to_s_orig(b), x.to_s(b))
      end
    end
  end

  def test_to_s  # [ruby-core:10686]
    assert_equal("fvvvvvvvvvvvv" ,18446744073709551615.to_s(32))
    assert_equal("g000000000000" ,18446744073709551616.to_s(32))
    assert_equal("3w5e11264sgsf" ,18446744073709551615.to_s(36))
    assert_equal("3w5e11264sgsg" ,18446744073709551616.to_s(36))
    assert_equal("nd075ib45k86f" ,18446744073709551615.to_s(31))
    assert_equal("nd075ib45k86g" ,18446744073709551616.to_s(31))
    assert_equal("1777777777777777777777" ,18446744073709551615.to_s(8))
    assert_equal("-1777777777777777777777" ,-18446744073709551615.to_s(8))
  end
end

Thread

Prev Next

In This Thread

Prev Next

[#31320] Import RubyGems to Ruby 1.9 — "NAKAMURA, Hiroshi" <nakahiro@...>

[#31323] Bignum#to_s の Karatsuba 基数変換による高速化 — "Kenta Murata" <muraken@...>

[#31330] Rubyアソシエーション設立のご報告 — Shugo Maeda <shugo@...>

[#31333] Invalid error message by illegal regexp — KIMURA Koichi <kimura.koichi@...>

[#31344] IO.popen with an array in trunk — sheepman <sheepman@...>

[#31345] RDoc i18n化パッチ — KURODA Hiraku <hiraku@...>

[#31349] Array#[]=(pos, length, val) and Array#[]=(range, val) — sheepman <sheepman@...>

[#31351] set_trace_func NULL pointer given — eklerni <eklerni@...>

[#31356] nmake don't have $< macro — KIMURA Koichi <kimura.koichi@...>

[#31357] invalid string for Date.parse — Yukihiro Matsumoto <matz@...>

[#31360] infinite loop at peephole optimization — "Yusuke ENDOH" <mame@...>

[#31361] retry neglects signals — "Yusuke ENDOH" <mame@...>

[#31362] typo in enum.c — Kazuhiro NISHIYAMA <zn@...>

[#31371] simultaneous exceptions dump core — "Yusuke ENDOH" <mame@...>

[#31372] jump in rescue clause in eval dumps core — "Yusuke ENDOH" <mame@...>

[#31373] jump in ensure dumps core — "Yusuke ENDOH" <mame@...>

[#31374] ARGF.each_byte yields a string — Tanaka Akira <akr@...>

[#31376] Re: [ ruby-Bugs-9490 ] Date module, step method, infinite loop if +step+ is 0 should raise an exception? — Urabe Shyouhei <shyouhei@...>

[#31377] Re: [ ruby-Patches-11719 ] add a :passive option to open-uri's open method — Urabe Shyouhei <shyouhei@...>

[#31378] Re: [ ruby-Patches-11954 ] CGI::HTTP_STATUS is incomplete — Urabe Shyouhei <shyouhei@...>

[#31379] Re: [ ruby-Bugs-10646 ] Failure parsing dates returned by WIN32OLE — Urabe Shyouhei <shyouhei@...>

[#31380] Re: [ ruby-Bugs-1823 ] NUM2UINT and rb_num2ull do not rase RangeError for negative Fixnum — Urabe Shyouhei <shyouhei@...>

[#31397] File exists - /tmp/bootstraptest.tmpwd — Tanaka Akira <akr@...>

[#31403] Re: [ ruby-Bugs-3361 ] Ruby 1.8.4 parse errors — Urabe Shyouhei <shyouhei@...>

[#31404] negative allocation size in rb_f_sprintf — Tietew <tietew+ruby-dev@...>

[#31407] [BUG] Stack consistency error (sp: 11, bp: 12) — Tanaka Akira <akr@...>

[#31416] Dir.mktmpdir — Tanaka Akira <akr@...>

[#31417] Process.euid [BUG] Segmentation fault — Tanaka Akira <akr@...>

[#31420] method `method_missing' called on terminated object (0xb7cbd6dc) (NotImplementedError) — Tanaka Akira <akr@...>

[#31424] Re: [ruby-cvs:20177] Re: Ruby:r12931 (trunk): * io.c (argf_close): always close via method. — Nobuyoshi Nakada <nobu@...>

[#31434] LP64: -0x10000000000000001 << -0x80000001 — Tanaka Akira <akr@...>

[#31437] yield [0] to {|*,v|} — Tanaka Akira <akr@...>

[#31438] /bin/sh: BASERUBY@: command not found — Tanaka Akira <akr@...>

[#31439] typo in ChangeLog — Kazuhiro NISHIYAMA <zn@...>

[#31440] yield [1] to {|v,&b| p v } — Tanaka Akira <akr@...>

[#31448] Ruby's (new) Bizarre Operator(s) — Nobuyoshi Nakada <nobu@...>

[#31462] Dir.mktmpdir for 1.8 — Tanaka Akira <akr@...>

[#31464] lambda {|a| p a }.call(1,2) — Tanaka Akira <akr@...>

[#31470] nested fiber invocation — Yukihiro Matsumoto <matz@...>

[#31472] lambda {|&b|}.call(3) — Tanaka Akira <akr@...>

[#31473] setter of $! — SASADA Koichi <ko1@...>

[#31475] lambda {|(v0,v1),v2|}.call([1],2) — Tanaka Akira <akr@...>

[#31480] trunk/test/ruby/test_proc.rb:127 — SASADA Koichi <ko1@...>

[#31481] trunk/test/test_time.rb — SASADA Koichi <ko1@...>

[#31492] SEGV: def m() yield [] end; m(&lambda {|*v| p v}) — Tanaka Akira <akr@...>

[#31502] {|(a,a)|} — Tanaka Akira <akr@...>

[#31507] {|(a)|} assigns outer local variable — Tanaka Akira <akr@...>

[#31511] trunk/test/ruby/test_eval.rb — SASADA Koichi <ko1@...>

[#31522] a, a = 1, 2 — Tanaka Akira <akr@...>

[#31525] いくつかのバグ報告と提案(5点) — eklerni <eklerni@...>

[#31531] 多重代入の*とto_ary — eklerni <eklerni@...>

[#31539] strtod の精度 — Satoshi Nakagawa <snakagawa@...>

[#31547] duck typing — "NAKAMURA, Hiroshi" <nakahiro@...>

[#31561] ignore exception message — SASADA Koichi <ko1@...>

[#31566] test_set and test_hash — Tanaka Akira <akr@...>

[#31576] test/win32ole — SASADA Koichi <ko1@...>

[#31583] Fiber reviesed — SASADA Koichi <ko1@...>

[#31592] [BUG] iseq_compile_each: unknown node: NODE_METHOD — Tanaka Akira <akr@...>

[#31603] benchmark of loops — SASADA Koichi <ko1@...>

[#31609] catch without a tag argument — Tanaka Akira <akr@...>

[#31625] IO.sysdup2, IO.sysdup, IO.sysclose — Tanaka Akira <akr@...>

[#31643] Re: [ruby-list:43915] Re: Hashへの生成順は保障されないのか？ — SASADA Koichi <ko1@...>

[#31646] Re: [ruby-cvs:20498] Ruby:r13261 (trunk): * encoding.c: provide basic features for M17N. — Tanaka Akira <akr@...>

[#31650] ungetc — Tanaka Akira <akr@...>

[#31651] rb_enc_mbclen — Tanaka Akira <akr@...>

[#31652] test/ruby/ — SASADA Koichi <ko1@...>

[#31659] each_byte and break — Tanaka Akira <akr@...>

[#31661] String#inspectで無限ループ — akira yamada / やまだあきら <akira@...>

[#31663] return in rescue in while in eval dumps core — "Yusuke ENDOH" <mame@...>

[#31665] String#slice! — WATANABE Tetsuya <Tetsuya.WATANABE@...>

[#31666] IndexError in rindex — Nobuyoshi Nakada <nobu@...>

[#31669] Re: [ ruby-Bugs-11882 ] recursive invokation of garbage_collect() in dlptr_free() — Urabe Shyouhei <shyouhei@...>

[#31671] Array#slice! — Satoshi Nakagawa <snakagawa@...>

[#31672] -Kオプションとparser->encの関係 — akira yamada / やまだあきら <akira@...>

[#31674] flush_register_windowsの"ta 0x83" — akira yamada / やまだあきら <akira@...>

[#31678] optimize empty hash space — Tanaka Akira <akr@...>

[#31679] ruby_1_8_5 undefined reference to `_RSTRING_PTR' — Masaki Suketa <masaki.suketa@...>

[#31682] mkmf.rb:310:in `method_missing': undefined method `quote' for nil:NilClass — akira yamada / やまだあきら <akira@...>

[#31684] WindowsのUnicodeファイル名 — OJ <oj@24.am>

[#31687] instance variable as formal argument — "Yusuke ENDOH" <mame@...>

[#31351] set_trace_func　NULL pointer given — eklerni <eklerni@...>