結果

問題 No.3030 ミラー・ラビン素数判定法のテスト
ユーザー MotSakMotSak
提出日時 2023-05-24 18:43:46
言語 C#
(.NET 8.0.203)
結果
AC  
実行時間 188 ms / 9,973 ms
コード長 6,169 bytes
コンパイル時間 13,937 ms
コンパイル使用メモリ 165,884 KB
実行使用メモリ 187,580 KB
最終ジャッジ日時 2024-06-06 06:50:43
合計ジャッジ時間 15,920 ms
ジャッジサーバーID
（参考情報） 		judge3 / judge2
このコードへのチャレンジ
（要ログイン）

テストケース

テストケース表示
入力 結果 実行時間
実行使用メモリ
testcase_00 AC 39 ms
23,296 KB
testcase_01 AC 38 ms
23,168 KB
testcase_02 AC 39 ms
23,268 KB
testcase_03 AC 40 ms
23,296 KB
testcase_04 AC 144 ms
26,876 KB
testcase_05 AC 150 ms
27,008 KB
testcase_06 AC 124 ms
26,492 KB
testcase_07 AC 126 ms
26,752 KB
testcase_08 AC 126 ms
26,624 KB
testcase_09 AC 188 ms
187,580 KB
権限があれば一括ダウンロードができます
コンパイルメッセージ
  復元対象のプロジェクトを決定しています...
  /home/judge/data/code/main.csproj を復元しました (104 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/

ソースコード

diff # 

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;
			}
		}
	}
}
0