結果
| 問題 |
No.8030 ミラー・ラビン素数判定法のテスト
|
| ユーザー |
 MotSak
|
| 提出日時 | 2023-05-20 11:00:46 |
| 言語 | C# (.NET 8.0.404) |
| 結果 |
AC
|
| 実行時間 | 257 ms / 9,973 ms |
| コード長 | 6,397 bytes |
| コンパイル時間 | 16,663 ms |
| コンパイル使用メモリ | 167,468 KB |
| 実行使用メモリ | 188,128 KB |
| 最終ジャッジ日時 | 2024-12-21 09:35:09 |
| 合計ジャッジ時間 | 19,233 ms |
|
ジャッジサーバーID （参考情報） |
judge5 / judge1 |
（要ログイン）
| ファイルパターン | 結果 |
|---|---|
| other | AC * 10 |
コンパイルメッセージ
復元対象のプロジェクトを決定しています... /home/judge/data/code/main.csproj を復元しました (105 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 class U64Mont
{
private ulong N { get; }
private ulong Ninv { get; }
private ulong Nhalf { get; }
private ulong R { get; }
private ulong Rneg { get; }
private ulong D { get; }
private int S { 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;
S = BitOperations.TrailingZeroCount(n - 1);
D = n - 1 >> S;
}
public bool BailliePSWTest()
{
return MillerRabinTestBase2() && LucasTest();
}
public bool MillerRabinTestBase2()
{
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;
}
public bool LucasTest()
{
ulong n = N;
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 qm = AR((d < 0) ? ((((ulong)(1 - d)) >> 2) % n) : (n - (((ulong)(d - 1)) >> 2) % n));
ulong k = (n + 1) << BitOperations.LeadingZeroCount(n + 1);
ulong um = R;
ulong vm = R;
ulong qn = qm;
ulong dm = AR((d < 0) ? (n - (ulong)-d % n) : ((ulong)d % n));
for (k <<= 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; }
ulong x = (n + 1) & (~n);
for (x >>= 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(BigMul(ar, br, out ulong t), 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;
}
// ((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 ulong BigMul(ulong a, ulong b, out ulong low)
{
return Math.BigMul(a, b, out low);
}
private static ulong IntSqrtU64(ulong x, out ulong remain)
{
const int BITS = 64;
BigInteger tmpRemain = x;
BigInteger tmpAnswer = (((BigInteger)0x0000000000000000UL) << 64) | 0x4000000000000000UL;
BigInteger maskLead = (((BigInteger)0xfffffffffffffffeUL) << 64) | 0x0000000000000000UL;
BigInteger maskCarry = (((BigInteger)0x0000000000000001UL) << 64) | 0x0000000000000000UL;
BigInteger maskTail = (((BigInteger)0x0000000000000000UL) << 64) | 0xffffffffffffffffUL;
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 j;
if (sa >= 0) { a = (ulong)sa; j = 1; } else { a = (ulong)-sa; j = ((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) { j = -j; }
if ((a & n & 3) == 3) { j = -j; }
ulong t = n; n = a; a = t; a %= n;
if (a > (n >> 1))
{
a = n - a;
if ((n & 3) == 3) { j = -j; }
}
}
return (n == 1) ? j : 0;
}
}
}
}