Hatena::ブログ(Diary)

まめめも このページをアンテナに追加 RSSフィード

2008-01-16

[][][] 多倍長整数演算の速度比較

rubyperlpython の多倍長整数演算の速度を適当に比較してみました。フィボナッチ数の計算速度比較以上にどーでもいい比較です。

前置き

使用した処理系のバージョンはそれぞれ以下のとおりです。

perl は多倍長整数を扱う方法が何個かあるようなので、標準装備らしい bigint と、別途 CPANインストールしないといけない Math::BitInt::GMP (GMP の binding) の 2 種類を試しました。python の psyco はなんか動かなくて面倒だったので試してません。この比較にはあまり影響しないような気がしてます。僕は PerlPython は素人なので、もっと速くて美しい書き方・やり方があったら教えてください。

実験 1: 5 の累乗の計算

5 の 30,000 乗、100,000 乗、300,000 乗、1,000,000 乗の計算にかかる時間を調べました。単位は秒。

300001000003000001000000
ruby 0.09 0.13 0.78 8.21
python 0.14 0.12 0.31 1.37
perl(bigint) 3.67 37.05 N/A N/A
perl(GMP) 0.06 0.06 0.09 0.21

GMP は当然として python も結構速いです。python は乗算に Karatsuba 法を実装しているようです。perl の bigint は圧倒的に遅いです。100 秒超えたら N/A です。

以下ソース。n を数字に置き換えてください。

# ruby
5**n
# python
5**n
# perl (bigint)
use bigint;
5**n;
# perl (GMP)
use Math::BigInt lib => "GMP";
Math::BigInt->new(5)->bpow(n);

実験 2: 5 の累乗の表示

実験 1 の数字を 10 進表示するのにかかる時間を調べました。多倍長整数を文字列化するということは基数変換するということであり、これが結構時間のかかる処理なのです。

300001000003000001000000
ruby 0.01 0.05 4.05 42.88
python 0.24 3.16 28.92 N/A
perl(bigint) 0.01 1.74 N/A N/A
perl(GMP) 0.02 0.06 0.21 1.20

やはり GMP が圧倒的ですが、Ruby もかろうじて動きます。Ruby は基数変換に Karatsuba 法を用いているためです (ruby-dev:31323) 。ちょっと手前味噌

以下ソース。上の表の数値は下のソースの実行時間から実験 1 の数値を引いてあります (つまり累乗の計算時間は含みません) 。

# ruby
p(5**n)
# python
print(5**n)
# perl (bigint)
use bigint;
print(5**n);
print "\n";
# perl (GMP)
use Math::BigInt lib => "GMP";
print Math::BigInt->new(5)->bpow(n);
print "\n";

実験 3: 円周率の計算

最後に、多倍長整数を使って円周率を 100 桁、1,000 桁、10,000 桁計算する時間を調べました。

100100010000
ruby 0.03 0.25 5.77
python 0.09 0.13 6.11
perl(bigint) 0.25 3.44 N/A
perl(GMP) 0.11 0.59 6.56

あれほど圧倒的だった GMP がなぜか普通です。もっと大きな桁で威力を発揮するのかな。rubypython はほぼ同じ。bigint はやっぱり全然ダメです。

以下ソース。python わからないのでなんだかかっこわるい。GMP の演算は破壊的なので copy() とか必要でうっとうしすぎます。

# ruby
a = b = 10 ** n
(n * 8 + 1).step(3, -2) do |i|
  a = (i / 2) * (a + b * 2) / i
end
puts "3.#{ a - b }"
# python
a = b = 10 ** n
for i in range(n * 8 + 1, 1, -2):
  a = divmod(divmod(i, 2)[0] * (a + b * 2), i)[0]
print("3." + str(a - b))
# perl (bigint)
use strict;
use bigint;

my $a;
my $b;
my $i;

$a = $b = 10 ** n;
for($i = n * 8 + 1; $i >= 3; $i -= 2) {
  $a = int(int($i / 2) * ($a + $b * 2) / $i);
}
$a -= $b;
print("3.$a\n");
# perl (GMP)
use strict;
use Math::BigInt lib => "GMP";

my $a = Math::BigInt->new(10)->bpow(n);
my $b = $a->copy();
my $i;
my $t;

for($i = n * 8 + 1; $i >= 3; $i -= 2) {
  $t = $b->copy();
  $a->badd($t);
  $a->badd($t);
  $t = Math::BigInt->new(int($i / 2));
  $t->bmul($a);
  $t->bdiv($i);
  $a = $t;
}
$a->bsub($b);
print("3.$a\n");

いい加減なまとめ

現状では、気軽に多倍長演算を享受したい時は rubypython がおすすめです。perl の bigint は速度がどうでもいいときだけ使いましょう速度が重要なら GMP とあわせて使いましょう。ソースを書くのが面倒でもそれなりに速く多倍長演算したいときは GMP がよさそうです。

以下追記:

H.I. 2008/01/17 10:09

「use bigint lib => 'GMP'」と書くことができて、GMPで

operator/constant overload をしてくれるみたいです。

との情報をいただきました。Perl では記述性と速度を両取りする選択肢もあるようです。Perl はさすがですね。円周率 10000 桁で試したところ、10 秒でした。GMP を直接たたくより若干遅い *1 ようですが、記述性があがる方がいいですね。

あと、すべて実時間で計測してます。プロセス初期化時間やコンパイル時間が含まれますので、0.0x のような小さい数値はあてにしないでください。

Appendix

以上の実験を自動でやってくれるスクリプトです。清書してないので読みにくくてごめんなさい。

require "timeout"

RB="./ruby19/local/bin/ruby"
PY="./python30/local/bin/python"
PL="./perl5/local/bin/perl"
LIMIT = 100

puts `#{RB} -v`
puts `#{PY} --version`
puts `#{PL} -v`.split("\n")[1]
puts

def test(prog, src, offset = 0)
  open("test-code", "w") {|fh| fh.print src }
  pid = fork
  unless pid
    $stdout.reopen(open("/dev/null", "w"))
    Process.setrlimit(Process::RLIMIT_CPU, LIMIT + 10)
    exec(prog, "test-code")
  end
  t = Time.now
  Process.wait
  t = Time.now - t
  puts(t > LIMIT ? "n/a" : "%3.3f" % (t - offset))
  t
end

def test_rb(src, offset = 0)
  print "ruby: "
  test(RB, src, offset)
end

def test_py(src, offset = 0)
  print "python: "
  test(PY, src, offset)
end

def test_pl(src, offset = 0)
  print "perl(bigint): "
  test(PL, src, offset)
end

def test_plg(src, offset = 0)
  print "perl(Math::BigInt::GMP): "
  test(PL, src, offset)
end

[30000, 100000, 300000, 1000000].each do |n|
  puts "5^#{n}"
  trb = test_rb "(5**#{n})"
  tpy = test_py "(5**#{n})"
  tpl = test_pl <<SRC
use strict;
use bigint;
5**#{n};
SRC
  tplg = test_plg <<SRC
use strict;
use Math::BigInt lib => "GMP";

Math::BigInt->new(5)->bpow(#{n});
SRC
  puts

  puts "5^#{n} with print"
  test_rb "p(5**#{n})", trb
  test_py "print(5**#{n})", tpy
  test_pl <<SRC, tpl
use strict;
use bigint;

print 5**#{n};
print "\\n";
SRC
  test_plg <<SRC, tplg
use strict;
use Math::BigInt lib => "GMP";

print Math::BigInt->new(5)->bpow(#{n});
print "\\n";
SRC
  puts
end

[100, 1000, 10000].each do |n|
  puts "#{n} digits of pi"
  test_rb <<SRC
a = b = 10 ** #{n}
#{n * 8 + 1}.step(3, -2) do |i|
  a = (i / 2) * (a + b * 2) / i
end
puts "3.\#{ a - b }"
SRC
  test_py <<SRC
a = b = 10 ** #{n}
for i in range(#{n * 8 + 1}, 1, -2):
  a = divmod(divmod(i, 2)[0] * (a + b * 2), i)[0]
print("3." + str(a - b))
SRC
  test_pl <<SRC
use strict;
use bigint;

my $a;
my $b;
my $i;

$a = $b = 10 ** #{n};
for($i = #{n * 8 + 1}; $i >= 3; $i -= 2) {
  $a = int(int($i / 2) * ($a + $b * 2) / $i);
}
$a -= $b;
print("3.$a\n");
SRC
  test_plg <<SRC
use strict;
use Math::BigInt lib => "GMP";

my $a = Math::BigInt->new(10)->bpow(#{n});
my $b = $a->copy();
my $i;
my $t;

for($i = #{n * 8 + 1}; $i >= 3; $i -= 2) {
  $t = $b->copy();
  $a->badd($t);
  $a->badd($t);
  $t = Math::BigInt->new(int($i / 2));
  $t->bmul($a);
  $t->bdiv($i);
  $a = $t;
}
$a->bsub($b);
print("3.$a\n");
SRC
  puts
end

*1:全然調べてませんが、Perlオブジェクトと GMP のオブジェクトを毎回変換してるとかかな?

実計算時間で計測?実計算時間で計測? 2008/01/17 09:44 Pythonあたりコンパイル時間が含まれているような

H.I.H.I. 2008/01/17 10:09 「use bigint lib => ’GMP’」と書くことができて、GMPで
operator/constant overload をしてくれるみたいです。

k.inabak.inaba 2008/01/17 15:45 divmod(a,b)[0] が整数の範囲で割り算したいということであれば、a//b でおkだと思いまする

ku-ma-meku-ma-me 2008/01/17 21:27 > 実計算時間で計測?
そうです。ひとこと書いておくべきでしたね。

> H.I. さん
おおお、情報ありがとうございます。まとめに追記しました。

> きーなばー
// ありがとうございまする

アルゴリズムが気になるアルゴリズムが気になる 2010/11/24 08:00 ここで使用されている円周率を計算するアルゴリズムはどういったものなのでしょうか?
名前など、宜しければ詳細を教えて頂けませんか?

succeedsucceed 2010/11/25 00:11 アルゴリズムは、4/π=1+1/(3+4/....)っていう公式を変形してπ=の形にしたものだと思います。
この式の連分数の分母の部分をaと置く、つまり4/π=1+1/aとすると、
π=1+4/(1+a)となりますが、後はこのaを展開していけばいいです。
http://ja.wikipedia.org/wiki/%E5%86%86%E5%91%A8%E7%8E%87

kikkik 2010/11/25 01:05 aもbも10^nで割っておけば、
a = (i/2)/i*(2+a) (i= 8n+3, 8n+1, ..., 3)
なんだから
pi-2 = 最後の a = (3/2)/3*(2+"ひとつ前のa")
= 2/3 + 2/3*"ひとつ前のa"
= 2/3 + 2/3*(5/2)/5*(2+"2つ前のa")
= 2/3 + (2/3)*(4/5) + (2/3)*(4/5)*"3つ前のa"
= 2^1*1/3 +2^2*1*2/3*5+2^3*1*2*3/3*5*7+2^4*1*2*3*4/3*5*7*9+...

だから 2nCn 系の公式のどれかだよ

ku-ma-meku-ma-me 2010/11/25 01:18 問題のコードを書いた人ですが、何と言うことか、どうやって書いたか覚えてないしわからなかったので、twitter で聞いてみたらいろんな人が教えてくれました。以下の公式っぽいです。

http://en.wikipedia.org/wiki/Numerical_approximations_of_%CF%80#Arctangent
http://twitter.com/#!/ainsophyao/status/7451872415715328
http://twitter.com/#!/ainsophyao/status/7463916833603584

みんなありがとうございます!

気になるって言った人気になるって言った人 2010/11/26 13:21 ありがとうございますm(_ _)m

スパム対策のためのダミーです。もし見えても何も入力しないでください
ゲスト


画像認証