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

• 問題: No.8030 ミラー・ラビン素数判定法のテスト
• ユーザー: iwkjosec
• 提出日時: 2026-08-11 17:41:40
• 言語: C# (.NET 10.0.201)
• 結果: AC
• 実行時間: 201 ms / 9,973 ms
• コード長: 8,797 bytes
• コンパイル時間: 11,524 ms
• コンパイル使用メモリ: 169,856 KB
• 実行使用メモリ: 196,616 KB
• 最終ジャッジ日時: 2026-08-11 17:41:56
• 合計ジャッジ時間: 13,771 ms
• ジャッジサーバーID （参考情報）: judge1_0 / judge2_0

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

コンパイルメッセージ

  復元対象のプロジェクトを決定しています...
  /home/judge/data/code/main.csproj を復元しました (91 ミリ秒)。
/home/judge/data/code/Main.cs(303,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
{
    /// <summary>
    /// Montgomery reduction.
    ///
    /// Returns x * R^-1 (mod m), where R = 2^64.
    /// m must be odd.
    /// </summary>
    private static ulong MontgomeryReduction(ulong xlo, ulong xhi, ulong m, ulong mInv)
    {
        ulong q = xlo * mInv;
        ulong result = xhi - Math.BigMul(q, m, out _);

        // Add m when the subtraction underflowed.
        return result > xhi ? result + m : result;
    }

    /// <summary>
    /// Computes a^-1 (mod 2^64) for odd a.
    ///
    /// Newton iteration:
    ///     x' = x * (2 - a*x)
    ///
    /// If y = 1 - a*x, the new error is y^2.
    /// Thus the number of correct bits approximately doubles
    /// at every iteration.
    /// </summary>
    private static ulong MultiplicativeInverse(ulong a)
    {
        Debug.Assert((a & 1) != 0);

        // Initial approximation.
        ulong x = (3 * a) ^ 2;

        // y = 1 - a*x
        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;
    }

    /// <summary>
    /// Computes a * b * R^-1 (mod m), where R = 2^64.
    /// m must be odd.
    /// </summary>
    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);
    }

    /// <summary>
    /// Computes (a + b) mod m.
    /// Requires a,b < m.
    /// </summary>
    private static ulong MontgomeryAddition(ulong a, ulong b, ulong m)
    {
        ulong sum = a + b;

        // Either the addition overflowed, or sum >= m.
        return sum < a || sum >= m ? sum - m : sum;
    }

    /// <summary>
    /// Returns constants required for Montgomery arithmetic:
    ///
    /// modInv = m^-1 (mod 2^64)
    /// rMod   = R mod m
    /// r2Mod  = R^2 mod m
    ///
    /// m must be odd.
    /// </summary>
    private static (ulong modInv, ulong rMod, ulong r2Mod) MontgomeryConstants(ulong m)
    {
        ulong modInv = MultiplicativeInverse(m);

        // R mod m, where R = 2^64.
        ulong rMod = (0ul - m) % m;

        // Compute R^2 mod m.
        //
        // Starting from R mod m:
        //   x <- 4R mod m
        //   x <- x^2 repeatedly
        //
        // Five squarings give:
        //   (4R)^(2^5) = R^(34)
        //
        // But because the Montgomery multiplication itself
        // contributes R^-1 at each step, the resulting value
        // is R^2 mod m.
        ulong r2Mod = MontgomeryAddition(rMod, rMod, m);
        r2Mod = MontgomeryAddition(r2Mod, r2Mod, m);

        for (int i = 0; i < 5; i++)r2Mod = MontgomeryMultiplication(r2Mod, r2Mod, m, modInv);

        return (modInv, rMod, r2Mod);
    }

    /// <summary>
    /// Computes value^exponent mod m using Montgomery multiplication.
    ///
    /// 'power' is kept in Montgomery representation:
    ///     value * R mod m
    ///
    /// 'result' deliberately remains in ordinary representation.
    /// This works because MontgomeryMultiplication(normal, mont)
    /// returns the ordinary product.
    /// </summary>
    private static ulong ModPow(ulong value, ulong exponent, ulong m)
    {
        var (modInv, rMod, r2Mod) = MontgomeryConstants(m);

        // value -> value * R (mod m)
        ulong power = MontgomeryMultiplication(value, r2Mod, m, modInv);

        // Deliberately ordinary representation.
        ulong result = 1;

        while (exponent != 0)
        {
            if ((exponent & 1) != 0)
                // ordinary * Montgomery -> ordinary
                result = MontgomeryMultiplication(result, power, m, modInv);

            // Montgomery * Montgomery -> Montgomery
            power = MontgomeryMultiplication(power, power, m, modInv);

            exponent >>= 1;
        }

        return result;
    }

    /// <summary>
    /// Computes a * b mod m.
    ///
    /// Handles both odd and even moduli.
    /// </summary>
    private static ulong ModMul(ulong a, ulong b, ulong m)
    {
        if ((m & 1) != 0) return ModMulOdd(a, b, m);
        return ModMulEven(a, b, m);
    }

    /// <summary>
    /// Computes a * b mod m for odd m.
    ///
    /// Uses:
    ///     MontMul(a,b) = abR^-1
    ///
    /// followed by multiplication by R^2, giving abR,
    /// and one final Montgomery reduction implicitly through
    /// the second multiplication.
    /// </summary>
    private static ulong ModMulOdd(ulong a, ulong b, ulong m)
    {
        var (modInv, _, r2Mod) = MontgomeryConstants(m);

        ulong result = MontgomeryMultiplication(a, b, m, modInv);

        result = MontgomeryMultiplication(result, r2Mod, m, modInv);

        return result;
    }

    /// <summary>
    /// Computes a * b mod m for even m.
    ///
    /// Write:
    ///     m = 2^k * oddMod
    ///
    /// Compute the product modulo both factors, then reconstruct
    /// the result using Garner's algorithm.
    /// </summary>
    private static ulong ModMulEven(ulong a, ulong b, ulong m)
    {
        int k = BitOperations.TrailingZeroCount(m);
        ulong oddMod = m >> k;
        ulong mask = (1ul << k) - 1;

        var (modInv, _, r2Mod) = MontgomeryConstants(oddMod);

        // --------------------------------------------------------
        // Product modulo oddMod
        // --------------------------------------------------------

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

        // First Montgomery reduction:
        //     abR^-1 mod oddMod
        ulong tOdd = MontgomeryReduction(lo, hi, oddMod, modInv);

        // Convert back:
        //     (abR^-1) * R^2 * R^-1 = ab
        tOdd = MontgomeryMultiplication(tOdd, r2Mod, oddMod, modInv);

        // --------------------------------------------------------
        // Product modulo 2^k
        // --------------------------------------------------------

        ulong tEven = lo & mask;

        // --------------------------------------------------------
        // Garner reconstruction
        //
        // t = tOdd + c * oddMod
        //
        // Choose c so that:
        //     t == tEven (mod 2^k)
        //
        // c = (tEven - tOdd) * oddMod^-1 (mod 2^k)
        // --------------------------------------------------------

        ulong c = ((tEven - tOdd) * modInv) & mask;

        return tOdd + c * oddMod;
    }

    /// <summary>
    /// Deterministic Miller-Rabin test for all 64-bit unsigned integers.
    /// </summary>
    public static bool MillerRabin(ulong value)
    {
        // 0, 1 and 2.
        if (value <= 2) return value == 2;

        // Even numbers other than 2 are composite.
        if ((value & 1) == 0) return false;

        ulong nMinusOne = value - 1;

        // n - 1 = d * 2^s, d odd.
        int s = BitOperations.TrailingZeroCount(nMinusOne);
        ulong d = nMinusOne >> s;

        // Deterministic for every 64-bit unsigned integer.
        ReadOnlySpan<ulong> bases = [2, 325, 9375, 28178, 450775, 9780504, 1795265022];

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

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

            // a^d == 1 or -1 (mod n)
            if (x == 1 || x == nMinusOne) continue;

            int i;
            for (i = 1; i < s; i++)
            {
                x = ModMul(x, x, value);
                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;
    }
}