結果

問題 No.8030 ミラー・ラビン素数判定法のテスト
ユーザー yuppe19 😺
提出日時 2019-02-15 15:41:42
言語 C++14
(gcc 13.3.0 + boost 1.87.0)
結果
AC  
実行時間 9,496 ms / 9,973 ms
コード長 2,333 bytes
コンパイル時間 792 ms
コンパイル使用メモリ 67,200 KB
実行使用メモリ 5,248 KB
最終ジャッジ日時 2024-11-16 23:12:43
合計ジャッジ時間 24,256 ms
ジャッジサーバーID
（参考情報） 		judge1 / judge4
このコードへのチャレンジ
（要ログイン）
ファイルパターン 結果
other AC * 10
権限があれば一括ダウンロードができます

ソースコード

diff # 

#include <iostream>
#include <vector>
using namespace std;
using U64 = uint64_t;
using S64 = int64_t;

const S64 Primes[] = { 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37 };

U64 ModMul(U64 a, U64 b, U64 modulus)
{
  /*
	{
		S64 a2 = (S64)a;
		if (a2 >= 0)
		{
			S64 b2 = (S64)b;
			if (b2 >= 0)
			{
				if (!MulAsm(&a2, b2))
					return (U64)a2 % modulus;
			}
		}
	}
  */
	U64 result = 0;
	while (a != 0)
	{
		if ((a & 1) != 0)
			result = (result + b) % modulus;
		a >>= 1;
		b = (b << 1) % modulus;
	}
	return result;
}

U64 ModPow(U64 value, U64 exponent, U64 modulus)
{
	U64 w = 1;
	while (exponent > 0)
	{
		if ((exponent & 1) != 0)
			w = ModMul(w, value, modulus);
		value = ModMul(value, value, modulus);
		exponent >>= 1;
	}
	return w;
}


bool is_prime(S64 n)
{
	if (n <= 1)
		return false;
	if ((n & 1) == 0)
		return n == 2;
	if (n <= 1920000)
	{
    for(U64 p : {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37}) {
      if(n == p)     { return true; }
      if(n % p == 0) { return false; }
    }
    /*
		if (n == 3)
			return true;
		if (n % 6 != 1 && n % 6 != 5)
			return false;
		S64 m = n / 6 * 2 + (n % 6 == 1 ? 0 : 1);
		size_t size;
		const U8* primes_bin = GetPrimesBin(&size);
		return (primes_bin[m / 8] & (1 << (m % 8))) != 0;
    */
	}

	// Miller-Rabin primality test.
	U64 enough;
	if (n < 2047)
		enough = 1;
	else if (n < 1373653)
		enough = 2;
	else if (n < 25326001)
		enough = 3;
	else if (n < 3215031751)
		enough = 4;
	else if (n < 2152302898747)
		enough = 5;
	else if (n < 3474749660383)
		enough = 6;
	else if (n < 341550071728321)
		enough = 7;
	else if (n < 3825123056546413051)
		enough = 9;
	else
	{
		// n < 2^64 < 318665857834031151167461
		enough = 12;
	}
	{
		U64 d = (U64)n - 1;
		U64 s = 0;
		while ((d & 1) == 0)
		{
			s++;
			d >>= 1;
		}
		for (U64 i = 0; i < enough; i++)
		{
			U64 x = ModPow(Primes[i], d, (U64)n);
			U64 j;
			if (x == 1 || x == (U64)n - 1)
				continue;
			bool probablyPrime = false;
			for (j = 0; j < s; j++)
			{
				x = ModPow(x, 2, (U64)n);
				if (x == (U64)n - 1)
				{
					probablyPrime = true;
					break;
				}
			}
			if (!probablyPrime)
				return false;
		}
		return true;
	}
}

int main(void) {
  int n; scanf("%d", &n);
  for(int loop=0; loop<n; ++loop) {
    S64 x; scanf("%ld", &x);
    printf("%ld %d\n", x, is_prime(x));
  }
  return 0;
}
0