結果
| 問題 |
No.8030 ミラー・ラビン素数判定法のテスト
|
| ユーザー |
 MotSak
|
| 提出日時 | 2023-05-24 18:43:46 |
| 言語 | C# (.NET 8.0.404) |
| 結果 |
AC
|
| 実行時間 | 192 ms / 9,973 ms |
| コード長 | 6,169 bytes |
| コンパイル時間 | 10,118 ms |
| コンパイル使用メモリ | 166,644 KB |
| 実行使用メモリ | 187,720 KB |
| 最終ジャッジ日時 | 2024-12-23 19:22:36 |
| 合計ジャッジ時間 | 12,576 ms |
|
ジャッジサーバーID （参考情報） |
judge1 / judge3 |
（要ログイン）
| ファイルパターン | 結果 |
|---|---|
| other | AC * 10 |
コンパイルメッセージ
復元対象のプロジェクトを決定しています... /home/judge/data/code/main.csproj を復元しました (100 ms)。 MSBuild のバージョン 17.9.6+a4ecab324 (.NET) main -> /home/judge/data/code/bin/Release/net8.0/main.dll main -> /home/judge/data/code/bin/Release/net8.0/publish/
ソースコード
using System;
using System.IO;
using System.Numerics;
namespace PrimalityTest
{
class Program
{
static void Main()
{
Console.SetOut(new StreamWriter(Console.OpenStandardOutput()) { AutoFlush = false });
var n = FastIStream.Long();
for (int i = 0; i < n; i++)
{
var x = (ulong)FastIStream.Long();
Console.WriteLine(x.ToString() + " " + (x.IsPrime() ? "1" : "0"));
}
Console.Out.Flush();
}
}
static class FastIStream
{
static readonly Stream str = Console.OpenStandardInput();
const int size = 1024;
static readonly byte[] buffer = new byte[size];
static int ptr = size;
static byte Read()
{
if (ptr == size)
{
str.Read(buffer, 0, size);
ptr = 0;
}
return buffer[ptr++];
}
public static long Long()
{
var c = Read();
while (c < 0x21)
{
c = Read();
}
var n = false;
if (c == '-')
{
n = true;
c = Read();
}
var ret = 0L;
while (c > 0x20)
{
ret = ret * 10 + c - '0';
c = Read();
}
return n ? -ret : ret;
}
}
static class PrimalityTestExtension
{
public static bool IsPrime(this ulong n)
{
if (n == 2) { return true; }
if (n < 2 || (n & 1) == 0) { return false; }
var mont = new U64Mont(n);
return mont.BailliePSWTest();
}
private readonly struct U64Mont
{
private ulong N { get; }
private ulong Ninv { get; }
private ulong Nhalf { get; }
private ulong R { get; }
private ulong Rneg { get; }
public U64Mont(ulong n)
{
N = n;
Ninv = n;
for (int i = 0; i < 5; ++i)
{
Ninv *= 2 - n * Ninv;
}
Nhalf = (n >> 1) + 1;
R = (0 - n) % n;
Rneg = n - R;
}
public bool BailliePSWTest()
{
return MillerRabinTestBase2() && LucasTest();
}
private bool MillerRabinTestBase2()
{
int s = BitOperations.TrailingZeroCount(N - 1);
ulong d = N - 1 >> s;
ulong tR = PowR(Add(R, R), d);
if (tR == R || tR == Rneg) { return true; }
for (int j = 1; j < s; ++j)
{
tR = MR(tR, tR);
if (tR == Rneg) { return true; }
}
return false;
}
private bool LucasTest()
{
long d = 5;
for (int i = 0; i < 64; ++i)
{
if (Jacobi(d, N) == -1) { break; }
if (i == 32 && IntSqrtU64Remain(N) == 0) { return false; }
if ((i & 1) == 0) { d = -(d + 2); } else { d = 2 - d; }
}
ulong dm = AR(d);
ulong qm = AR((1 - d) / 4);
ulong qn = qm;
ulong um = R;
ulong vm = R;
for (ulong k = (N + 1) << (BitOperations.LeadingZeroCount(N + 1) + 1); k > 0; k <<= 1)
{
um = MR(um, vm);
vm = Sub(MR(vm, vm), Add(qn, qn));
qn = MR(qn, qn);
if ((k >> 63) != 0)
{
ulong uu = Div2(Add(um, vm));
vm = Div2(Add(MR(dm, um), vm));
um = uu;
qn = MR(qn, qm);
}
}
if (um == 0 || vm == 0) { return true; }
for (ulong x = ((N + 1) & (~N)) >> 1; x > 0; x >>= 1)
{
um = MR(um, vm);
vm = Sub(MR(vm, vm), Add(qn, qn));
if (vm == 0) { return true; }
qn = MR(qn, qn);
}
return false;
}
// a + b (mod n)
private ulong Add(ulong a, ulong b)
{
return AddOverflow(a, b, out ulong t) || t >= N ? t - N : t;
}
// a - b (mod n)
private ulong Sub(ulong a, ulong b)
{
return a < b ? N + a - b : a - b;
}
// a / 2 (mod n)
private ulong Div2(ulong a)
{
return (a & 1) == 0 ? a >> 1 : (a >> 1) + Nhalf;
}
// (ar * br) / r (mod n)
private ulong MR(ulong ar, ulong br)
{
return Sub(Math.BigMul(ar, br, out ulong t), Math.BigMul(t * Ninv, N, out _));
}
// a * r (mod n)
private ulong AR(ulong a)
{
ulong t = R;
ulong ar = (a & 1) != 0 ? R : 0;
a >>= 1;
while (a != 0)
{
t = Add(t, t);
if ((a & 1) != 0) { ar = Add(ar, t); }
a >>= 1;
}
return ar;
}
// a * r (mod n)
private ulong AR(long a)
{
return a < 0 ? N - AR((ulong)-a) : AR((ulong)a);
}
// ((ar / r) ** k) * r (mod n)
private ulong PowR(ulong ar, ulong k)
{
ulong tr = (k & 1) != 0 ? ar : R;
k >>= 1;
while (k != 0)
{
ar = MR(ar, ar);
if ((k & 1) != 0) { tr = MR(tr, ar); }
k >>= 1;
}
return tr;
}
private static bool AddOverflow(ulong lhs, ulong rhs, out ulong result)
{
return lhs > (result = lhs + rhs);
}
private static readonly BigInteger _maskLead = (((BigInteger)0xfffffffffffffffeUL) << 64) | 0x0000000000000000UL;
private static readonly BigInteger _maskCarry = (((BigInteger)0x0000000000000001UL) << 64) | 0x0000000000000000UL;
private static readonly BigInteger _maskTail = (((BigInteger)0x0000000000000000UL) << 64) | 0xffffffffffffffffUL;
private static ulong IntSqrtU64(ulong x, out ulong remain)
{
const int BITS = 64;
BigInteger tmpRemain = x;
BigInteger tmpAnswer = 0x4000000000000000UL;
for (uint i = 0; i < (BITS >> 1); ++i) {
BigInteger next = ((tmpAnswer << 1) & _maskLead) | (tmpAnswer & _maskTail);
if (tmpRemain >= tmpAnswer)
{
tmpRemain -= tmpAnswer;
tmpAnswer = next + _maskCarry;
}
else
{
tmpAnswer = next;
}
tmpRemain <<= 2;
}
remain = (ulong)(tmpRemain >> BITS);
return (ulong)(tmpAnswer >> BITS);
}
private static ulong IntSqrtU64Remain(ulong x)
{
_ = IntSqrtU64(x, out ulong remain);
return remain;
}
private static int Jacobi(long sa, ulong n)
{
ulong a; int t;
if (sa >= 0) { a = (ulong)sa; t = 1; } else { a = (ulong)-sa; t = ((n & 3) == 3) ? -1 : 1; }
while (a > 0)
{
int ba = BitOperations.TrailingZeroCount(a);
a >>= ba;
if (((n & 7) == 3 || (n & 7) == 5) && (ba & 1) != 0) { t = -t; }
if ((a & n & 3) == 3) { t = -t; }
ulong r = n; n = a; a = r; a %= n;
if (a > (n >> 1)) { a = n - a; if ((n & 3) == 3) { t = -t; } }
}
return (n == 1) ? t : 0;
}
}
}
}