Introduction

ブログ内検索

  • このサイトの記事を検索 by Google

おすすめの一冊!

無料ブログはココログ

« 素数判定 | トップページ | RSA 暗号 »

2010-04-13

素数判定(その2)

素数判定について調べたのですが、難しさがいまいち実感できなかったので、
ちょろっとコードを書いてみました。

多倍長演算の実装がめんどくさかったので、Perl です。

#!/usr/bin/perl -w

use strict;
use warnings;

use Math::BigInt;
use Math::Random::MT qw(srand rand);

#=============================================================================

sub rand32 {
	return (int(rand(0x10000)) << 16) | (int(rand(0x10000)));
}

sub rand_inf {
	my $arg = shift;
	if (ref($arg)) {
		# object of Math::BigInt
		my $x = substr($arg->as_hex(), 2);  # 先頭の "0x" は削っておく
		my $ntimes = int((length($x) / 8));
		my $y = '0x';
		if (length($x) % 8 != 0) {
			my $limit = hex(substr($x, 0, length($x)-8*$ntimes));
			$y .= sprintf("%x", rand32() % $limit);
		}
		for (my $i =0; $i < $ntimes; $i++) {
			$y .= sprintf("%08x", rand32());
		}
		return Math::BigInt->new($y);
	} else {
		return int(rand($arg));
	}
}

sub is_prime {
	my $n = shift;
	die unless ($n > 1);
	return 1 if ($n == 2);
	return 0 if $n->is_even();

	my $d = $n->copy();
	$d->bdec();
	$d->brsft(1) while $d->is_even();

	my $n1 = $n->copy()->bdec();
	
	foreach my $k (1..20) {
		my $a = rand_inf($n-2) + 1;
		my $t = $d->copy();
		my $y = $a->bmodpow($t, $n);
		while (($t != $n1) and ($y != 1) and ($y != $n1)) {
			$y = $y->bmodpow(2, $n);
			$t->blsft(1);
		}
		if (($y != $n1) and $t->is_even()) {
			print sprintf("%3d", $k);
			return 0;
		}
	}
	return 1;
}

sub get_prime {
	my $bit_length = shift;
	die unless ($bit_length % 4 == 0);
	my $limit = Math::BigInt->new('0x1' . ('0' x ($bit_length / 4)));
	my $count = 0;
	while (1) {
		my $x = rand_inf($limit);
		if (is_prime($x)) {
			print "\n";
			print "OK " . $bit_length . ' ' . $count . ' ' . sprintf("%258s", $x->as_hex()) . "\n";
			return $x;
		}
		$count++;
	}
}

{
	for (my $i=32; $i<=1024; $i*=2) {
		foreach (1..16) {
			print '### ' . time() . ' ' . `/bin/date` . "\n";
			my $value = get_prime($i);
		}
	}
}
実行すると、32bit から 1024bit までの素数を 16個ずつ見つけてくれます。 うちの PC で実行すると、 512bit の素数を1個みつけるのに 3分くらいかかっていました(打ち切ったので 1024bit は不明)。 ちなみに PC のスペックは以下の通り:
    AMD Phenom X4 9750 / 2.4GHz / WindowsXP SP3
    Perl v5.10.0 built for cygwin-thread-multi-64int
ところで、CPAN で "BigInt" で検索すると、 Math::BigInt::GMP とか Math::BigInt::Pari とか なにやら意味ありげなモジュールがみつかります。 GMP というのは GNU で開発している多倍長演算ライブラリらしい。 "GNU Multi-Precision Library" とのこと。 ⇒ GMP -- Wikipedia cygwin なので使えるかわからなかったのですが、 /usr/lib/libgmp.dll.a とかがすでに置いてあったので、 Let's Try! インストールは cpan コマンドで "install Math::BigInt::GMP" で一発。 使うには上記コードの6行目を変えるだけ。超カンタン♪
< use Math::BigInt;
--
> use Math::BigInt lib => 'GMP';
実行してみると、動作速度にびっくり。 512bit の素数が1秒未満で見つかります。1024bit でも数秒。 正直、ここまで速くなるとは思っていませんでした。GMP すごい! さて、これからちょっと高速化にチャレンジしてみようかと思います。 たぶん GMP は超えられないでしょうけど・・・

« 素数判定 | トップページ | RSA 暗号 »