#1176252 (C#) No.8030 ミラー・ラビン素数判定法のテスト

• 問題: No.8030 ミラー・ラビン素数判定法のテスト
• ユーザー: iwkjosec
• 提出日時: 2026-08-11 17:47:57
• 言語: C# (.NET 10.0.201)
• 結果: AC
• 実行時間: 189 ms / 9,973 ms
• コード長: 7,102 bytes
• コンパイル時間: 9,501 ms
• コンパイル使用メモリ: 173,016 KB
• 実行使用メモリ: 192,576 KB
• 最終ジャッジ日時: 2026-08-11 17:48:33
• 合計ジャッジ時間: 8,960 ms
• ジャッジサーバーID （参考情報）: judge2_0 / judge1_0

ファイルパターン	結果
other	AC * 10

コンパイルメッセージ

  復元対象のプロジェクトを決定しています...
  /home/judge/data/code/main.csproj を復元しました (96 ミリ秒)。
/home/judge/data/code/Main.cs(223,13): warning CA2022: 'System.IO.Stream.Read(byte[], int, int)' による不正確な読み取りを避ける (https://learn.microsoft.com/dotnet/fundamentals/code-analysis/quality-rules/ca2022) [/home/judge/data/code/main.csproj]
  main -> /home/judge/data/code/bin/Release/net10.0/main.dll
  main -> /home/judge/data/code/bin/Release/net10.0/publish/

ソースコード

using System.Diagnostics;
using System.Numerics;

class Program
{
    static void Main()
    {
        Console.SetOut(new StreamWriter(Console.OpenStandardOutput()) { AutoFlush = false });
        var n = FastIO.Long();
        for (int i = 0; i < n; i++)
        {
            var x = (ulong)FastIO.Long();
            Console.WriteLine($"{x} {(PrimalityTest.MillerRabin(x) ? '1' : '0')}");
        }
        Console.Out.Flush();
    }
}

public static class PrimalityTest
{
    // m は奇数。
    // x * 2^(-64) mod m を計算する。
    private static ulong MontgomeryReduction(ulong xlo, ulong xhi, ulong m, ulong mInv)
    {
        ulong q = xlo * mInv;
        ulong result = xhi - Math.BigMul(q, m, out _);
        return result > xhi ? result + m : result;
    }

    // 奇数 a の 2^64 を法とする逆元を Newton 法で求める。
    private static ulong MultiplicativeInverse(ulong a)
    {
        Debug.Assert((a & 1) != 0);

        ulong x = (3 * a) ^ 2;
        ulong y = 1 - a * x;

        x *= 1 + y;
        y *= y;

        x *= 1 + y;
        y *= y;

        x *= 1 + y;
        y *= y;

        x *= 1 + y;

        return x;
    }

    // Montgomery 表現同士を乗算する。
    // MontMul(a, b) = a * b * R^(-1) mod m, R = 2^64。
    private static ulong MontgomeryMultiplication(ulong a, ulong b, ulong m, ulong mInv)
    {
        ulong hi = Math.BigMul(a, b, out ulong lo);
        return MontgomeryReduction(lo, hi, m, mInv);
    }

    // a + b mod m。
    // a,b < m を仮定する。
    private static ulong MontgomeryAddition(ulong a, ulong b, ulong m)
    {
        ulong sum = a + b;
        return sum < a || sum >= m ? sum - m : sum;
    }

    // Montgomery 演算に必要な定数。
    private readonly struct MontgomeryContext
    {
        public readonly ulong Mod;
        public readonly ulong ModInv;
        public readonly ulong R2Mod;

        public MontgomeryContext(ulong mod)
        {
            Debug.Assert(mod > 1 && (mod & 1) != 0);

            Mod = mod;
            ModInv = MultiplicativeInverse(mod);

            ulong rMod = (0ul - mod) % mod;

            // 4R から始め、Montgomery 平方を5回行う。
            // xR -> x^2 R なので、
            // 4^32 R = 2^64 R = R^2 となる。
            ulong r2 = MontgomeryAddition(rMod, rMod, mod);
            r2 = MontgomeryAddition(r2, r2, mod);

            r2 = MontgomeryMultiplication(r2, r2, mod, ModInv);
            r2 = MontgomeryMultiplication(r2, r2, mod, ModInv);
            r2 = MontgomeryMultiplication(r2, r2, mod, ModInv);
            r2 = MontgomeryMultiplication(r2, r2, mod, ModInv);
            r2 = MontgomeryMultiplication(r2, r2, mod, ModInv);

            R2Mod = r2;
        }

        // 通常表現 -> Montgomery 表現。
        public ulong ToMontgomery(ulong x) => MontgomeryMultiplication(x, R2Mod, Mod, ModInv);

        // Montgomery 表現 -> 通常表現。
        public ulong FromMontgomery(ulong x) => MontgomeryMultiplication(x, 1, Mod, ModInv);

        // Montgomery 表現同士の乗算。
        public ulong Multiply(ulong a, ulong b) => MontgomeryMultiplication(a, b, Mod, ModInv);

        // 通常表現 a,b の積。
        //
        // まず abR^-1 を求め、それに R^2 を掛けることで
        // 最終的に通常表現の ab mod m を得る。
        public ulong MultiplyNormal(ulong a, ulong b)
        {
            ulong x = MontgomeryMultiplication(a, b, Mod, ModInv);
            return MontgomeryMultiplication(x, R2Mod, Mod, ModInv);
        }
    }

    // a^e mod m。
    // power は Montgomery 表現、result は通常表現のまま処理する。
    private static ulong ModPow(ulong a, ulong e, in MontgomeryContext ctx)
    {
        ulong power = ctx.ToMontgomery(a);
        ulong result = 1;

        while (e != 0)
        {
            if ((e & 1) != 0) result = ctx.Multiply(result, power);
            e >>= 1;
            if (e != 0) power = ctx.Multiply(power, power);
        }

        return result;
    }

    // a * b mod m。
    // 奇数 modulus 用の高速経路。
    private static ulong ModMulOdd(ulong a, ulong b, in MontgomeryContext ctx) => ctx.MultiplyNormal(a, b);

    // a * b mod m。
    // m が偶数の場合は m = 2^k * oddMod と分解し、
    // Garner 法で結果を再構成する。
    private static ulong ModMulEven(ulong a, ulong b, ulong mod)
    {
        int k = BitOperations.TrailingZeroCount(mod);
        ulong oddMod = mod >> k;
        ulong mask = (1ul << k) - 1;

        var ctx = new MontgomeryContext(oddMod);

        ulong hi = Math.BigMul(a, b, out ulong lo);

        // oddMod における積。
        ulong tOdd = MontgomeryReduction(lo, hi, oddMod, ctx.ModInv);
        tOdd = ctx.Multiply(tOdd, ctx.R2Mod);

        // 2^k における積。
        ulong tEven = lo & mask;

        // Garner 法による CRT 再構成。
        ulong c = ((tEven - tOdd) * ctx.ModInv) & mask;
        return tOdd + c * oddMod;
    }

    // a * b mod mod。
    public static ulong ModMul(ulong a, ulong b, ulong mod)
    {
        if ((mod & 1) != 0) return ModMulOdd(a, b, new MontgomeryContext(mod));

        return ModMulEven(a, b, mod);
    }

    // 64bit unsigned 整数に対する決定的 Miller-Rabin。
    public static bool MillerRabin(ulong value)
    {
        if (value <= 2) return value == 2;

        if ((value & 1) == 0) return false;

        // ここから先は必ず奇数なので、Montgomery 定数を一度だけ作る。
        var ctx = new MontgomeryContext(value);

        ulong nMinusOne = value - 1;
        int s = BitOperations.TrailingZeroCount(nMinusOne);
        ulong d = nMinusOne >> s;

        // 2^64 未満の整数に対して決定的な7つの底数。
        ReadOnlySpan<ulong> bases = [2, 325, 9375, 28178, 450775, 9780504, 1795265022];

        foreach (ulong a in bases)
        {
            if (a % value == 0) continue;

            ulong x = ModPow(a, d, ctx);

            if (x == 1 || x == nMinusOne) continue;

            int i;
            for (i = 1; i < s; i++)
            {
                x = ctx.MultiplyNormal(x, x);
                if (x == nMinusOne) break;
            }

            if (i == s) return false;
        }

        return true;
    }
}

public static class FastIO
{
    static Stream str = Console.OpenStandardInput();
    const int size = 1024;
    static 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;
    }
}