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

• 問題: No.8030 ミラー・ラビン素数判定法のテスト
• ユーザー: iwkjosec
• 提出日時: 2026-08-11 18:49:16
• 言語: C# (.NET 10.0.201)
• 結果: AC
• 実行時間: 187 ms / 9,973 ms
• コード長: 9,402 bytes
• コンパイル時間: 5,008 ms
• コンパイル使用メモリ: 171,084 KB
• 実行使用メモリ: 193,148 KB
• 最終ジャッジ日時: 2026-08-11 18:49:31
• 合計ジャッジ時間: 7,615 ms
• ジャッジサーバーID （参考情報）: judge2_0 / judge1_0

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

コンパイルメッセージ

  復元対象のプロジェクトを決定しています...
  /home/judge/data/code/main.csproj を復元しました (80 ミリ秒)。
  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;
using System.Runtime.CompilerServices;
using System.Text;

class Program
{
    static void Main()
    {
        Console.SetOut(new StreamWriter(Console.OpenStandardOutput()) { AutoFlush = false });

        var sb = new StringBuilder();
        var n = Scanner.Long();
        for (int i = 0; i < n; i++)
        {
            var x = Scanner.ULong();
            sb.Append(x).Append(' ').Append(PrimalityTest.MillerRabin(x) ? '1' : '0').AppendLine();
        }
        Console.Write(sb.ToString());
        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;
    }
}

static class Scanner
{
    private static readonly byte[] Buffer;
    private static int Index;

    static Scanner()
    {
        using var input = Console.OpenStandardInput();

        int capacity = 1 << 20;
        Buffer = new byte[capacity];

        int length = 0;

        while (true)
        {
            if (length == capacity)
            {
                capacity *= 2;
                Array.Resize(ref Buffer, capacity);
            }

            int n = input.Read(Buffer, length, capacity - length);

            if (n == 0) break;

            length += n;
        }

        if (length < Buffer.Length) Array.Resize(ref Buffer, length + 1);

        Buffer[length] = 0;
    }

    [MethodImpl(MethodImplOptions.AggressiveInlining)]
    public static int Int()
    {
        ref byte p = ref Buffer[Index];

        while (p <= 32)
        {
            p = ref Buffer[++Index];
        }

        bool negative = false;

        if (p == '-')
        {
            negative = true;
            p = ref Buffer[++Index];
        }

        int value = 0;

        while ((uint)(p - '0') < 10)
        {
            value = value * 10 + p - '0';
            p = ref Buffer[++Index];
        }

        return negative ? -value : value;
    }

    [MethodImpl(MethodImplOptions.AggressiveInlining)]
    public static long Long()
    {
        ref byte p = ref Buffer[Index];

        while (p <= 32)
        {
            p = ref Buffer[++Index];
        }

        bool negative = false;

        if (p == '-')
        {
            negative = true;
            p = ref Buffer[++Index];
        }

        long value = 0;

        while ((uint)(p - '0') < 10)
        {
            value = value * 10 + p - '0';
            p = ref Buffer[++Index];
        }

        return negative ? -value : value;
    }

    [MethodImpl(MethodImplOptions.AggressiveInlining)]
    public static uint UInt()
    {
        ref byte p = ref Buffer[Index];

        while (p <= 32)
        {
            p = ref Buffer[++Index];
        }

        uint value = 0;

        while ((uint)(p - '0') < 10)
        {
            value = value * 10u + (uint)(p - '0');
            p = ref Buffer[++Index];
        }

        return value;
    }

    [MethodImpl(MethodImplOptions.AggressiveInlining)]
    public static ulong ULong()
    {
        ref byte p = ref Buffer[Index];

        while (p <= 32)
        {
            p = ref Buffer[++Index];
        }

        ulong value = 0;

        while ((uint)(p - '0') < 10)
        {
            value = value * 10UL + (ulong)(p - '0');
            p = ref Buffer[++Index];
        }

        return value;
    }

    public static string String()
    {
        while (Buffer[Index] <= 32) ++Index;

        int start = Index;

        while (Buffer[Index] > 32) ++Index;

        return System.Text.Encoding.ASCII.GetString(Buffer, start, Index - start);
    }
}