結果

問題 No.3030 ミラー・ラビン素数判定法のテスト
ユーザー MotSakMotSak
提出日時 2023-05-20 11:00:46
言語 C#
(.NET 8.0.203)
結果
AC  
実行時間 247 ms / 9,973 ms
コード長 6,397 bytes
コンパイル時間 15,684 ms
コンパイル使用メモリ 167,084 KB
実行使用メモリ 187,628 KB
最終ジャッジ日時 2024-06-01 06:01:50
合計ジャッジ時間 17,077 ms
ジャッジサーバーID
（参考情報） 		judge1 / judge5
このコードへのチャレンジ
（要ログイン）

テストケース

テストケース表示
入力 結果 実行時間
実行使用メモリ
testcase_00 AC 44 ms
23,040 KB
testcase_01 AC 43 ms
22,912 KB
testcase_02 AC 42 ms
22,784 KB
testcase_03 AC 42 ms
22,912 KB
testcase_04 AC 170 ms
26,984 KB
testcase_05 AC 186 ms
26,880 KB
testcase_06 AC 151 ms
27,132 KB
testcase_07 AC 149 ms
26,624 KB
testcase_08 AC 149 ms
26,752 KB
testcase_09 AC 247 ms
187,628 KB
権限があれば一括ダウンロードができます
コンパイルメッセージ
  復元対象のプロジェクトを決定しています...
  /home/judge/data/code/main.csproj を復元しました (95 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 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;
			}
		}
	}
}
0