結果

問題 No.8030 ミラー・ラビン素数判定法のテスト
ユーザー iwkjosec
提出日時 2026-08-11 13:34:35
言語 C#
(.NET 10.0.201)
コンパイル:
dotnet_c
実行:
/usr/bin/dotnet_wrap
結果
AC  
実行時間 188 ms / 9,973 ms
+ 237µs
コード長 7,552 bytes
記録
記録タグの例:
初AC ショートコード 純ショートコード 純主流ショートコード 最速実行時間
コンパイル時間 14,486 ms
コンパイル使用メモリ 172,276 KB
実行使用メモリ 196,356 KB
最終ジャッジ日時 2026-08-11 13:35:19
合計ジャッジ時間 14,168 ms
ジャッジサーバーID
（参考情報） 		judge2_0 / judge1_0
このコードへのチャレンジ
（要ログイン）
ファイルパターン 結果
other AC * 10
権限があれば一括ダウンロードができます
コンパイルメッセージ
  復元対象のプロジェクトを決定しています...
  /home/judge/data/code/main.csproj を復元しました (82 ミリ秒)。
/home/judge/data/code/Main.cs(244,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/

ソースコード

diff #
raw source code 

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} {(MillerRabin(x) ? '1' : '0')}");
        }
        Console.Out.Flush();
    }
    public static ulong MontgomeryReduction
    (ulong xlo, ulong xhi, ulong m, ulong mInv)
    {
        ulong result = xhi - Math.BigMul(xlo * mInv, m, out _);
        return result > xhi ? result + m : result;
    }
    public static ulong BigMul(ulong a, ulong b, out ulong lo)
    {
        ulong alo = a & 0xffffffff, ahi = a >> 32;
        ulong blo = b & 0xffffffff, bhi = b >> 32;

        ulong lolo = alo * blo;
        ulong lohi = alo * bhi;
        ulong hilo = ahi * blo;
        ulong hihi = ahi * bhi;

        lo = lolo + ((lohi + hilo) << 32);
        ulong carry = ((lolo >> 32) + (lohi & 0xffffffff) + (hilo & 0xffffffff)) >> 32;
        ulong hi = hihi + (lohi >> 32) + (hilo >> 32) + carry;

        return hi;
    }
    public static ulong MultiplicativeInverse(ulong a)
    {
        Debug.Assert((a & 1) != 0);

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

        ulong x1 = x0 * (1 + y);
        y *= y;

        ulong x2 = x1 * (1 + y);
        y *= y;

        ulong x3 = x2 * (1 + y);
        y *= y;

        ulong x4 = x3 * (1 + y);

        return x4;
    }
    public static ulong MontgomeryMultiplication(ulong a, ulong b, ulong m, ulong mInv)
    {
        ulong xhi = Math.BigMul(a, b, out ulong xlo);
        return MontgomeryReduction(xlo, xhi, m, mInv);
    }
    public static ulong MontgomeryAddition(ulong a, ulong b, ulong m)
    {
        ulong add = a + b;
        return add < a || add >= m ? add - m : add;
    }

    public static ulong MontgomerySubtraction(ulong a, ulong b, ulong m)
    {
        ulong sub = a - b;
        return a < b ? sub + m : sub;
    }

    public static (ulong modinv, ulong rmod, ulong r2mod)
    MontgomeryConstant(ulong mod)
    {
        ulong modinv = MultiplicativeInverse(mod);
        ulong rmod = (0ul - mod) % mod;

        ulong r2mod = rmod;
        r2mod = MontgomeryAddition(r2mod, r2mod, mod);
        r2mod = MontgomeryAddition(r2mod, r2mod, mod);
        r2mod = MontgomeryMultiplication(r2mod, r2mod, mod, modinv);
        r2mod = MontgomeryMultiplication(r2mod, r2mod, mod, modinv);
        r2mod = MontgomeryMultiplication(r2mod, r2mod, mod, modinv);
        r2mod = MontgomeryMultiplication(r2mod, r2mod, mod, modinv);
        r2mod = MontgomeryMultiplication(r2mod, r2mod, mod, modinv);

        return (modinv, rmod, r2mod);
    }
    public static ulong ModPow(ulong value, ulong exponent, ulong mod)
    {
        var (modinv, rmod, r2mod) = MontgomeryConstant(mod);

        ulong power = MontgomeryMultiplication(value, r2mod, mod, modinv);
        ulong result = 1;

        while (exponent > 0)
        {
            if ((exponent & 1) != 0)
            {
                result = MontgomeryMultiplication(result, power, mod, modinv);
            }

            power = MontgomeryMultiplication(power, power, mod, modinv);
            exponent >>= 1;
        }

        return result;
    }
    private static readonly byte[] TrailingZeroLookup = [
    0, 1, 59, 2, 60, 40, 54, 3, 61, 32, 49, 41, 55, 19, 35, 4, 62, 52, 30, 33, 50, 12, 14, 42, 56, 16, 27, 20, 36, 23, 44, 5, 63,
        58, 39, 53, 31, 48, 18, 34, 51, 29, 11, 13, 15, 26, 22, 43, 57, 38, 47, 17, 28, 10, 25, 21, 37, 46, 9, 24, 45, 8, 7, 6
];

    public static int TrailingZeroCount(ulong value)
    {
        if (value == 0)
        {
            return 64;
        }

        return TrailingZeroLookup[((value & (0 - value)) * 0x03F566ED27179461) >> 58];
    }

    public static ulong ModMul(ulong a, ulong b, ulong mod)
    {
        // if mod is even then  
        if ((mod & 1) == 0)
        {
            int evenBits = TrailingZeroCount(mod);
            ulong oddMod = mod >> evenBits;
            var (modinv, rmod, r2mod) = MontgomeryConstant(oddMod);
            ulong mask = (1ul << evenBits) - 1;


            // MR(a * b)  
            ulong thi = BigMul(a, b, out ulong tlo);
            ulong tOdd = thi - BigMul(tlo * modinv, oddMod, out _);
            if (tOdd > thi)
            {
                tOdd += oddMod;
            }

            // == (a * b) & mask  
            ulong tEven = tlo & mask;

            // MR(ab * R^2)  
            thi = BigMul(tOdd, r2mod, out tlo);
            tOdd = thi - BigMul(tlo * modinv, oddMod, out _);
            if (tOdd > thi)
            {
                tOdd += oddMod;
            }

            // Garner's algorithm  
            ulong t = tOdd + (((tEven - tOdd) * modinv) & mask) * oddMod;
            return t;
        }
        else
        {
            // if mod is odd then  
            // do normal Montgomery Multiplication  
            var (modinv, rmod, r2mod) = MontgomeryConstant(mod);

            var mont = MontgomeryMultiplication(MontgomeryMultiplication(
                a, b, mod, modinv), r2mod, mod, modinv);
            return mont;
        }
    }

    // returns true if value is prime  
    public static bool MillerRabin(ulong value)
    {
        // 2 以下・偶数の場合の前処理  
        if (value <= 2)
        {
            return value == 2;
        }
        if ((value & 1) == 0)
        {
            return false;
        }

        ulong n1 = value - 1;
        int s = TrailingZeroCount(n1);
        ulong d = n1 >> s;

        // 「注意深く選ばれた」定数  
        // 2^64 までならこの 7 つで対応可能  
        ReadOnlySpan<ulong> MillerRabinConstants =
            [2, 325, 9375, 28178, 450775, 9780504, 1795265022];

        foreach (var a in MillerRabinConstants)
        {
            // 割り切れた場合は続行  
            if (a % value == 0)
            {
                continue;
            }

            // a^d mod value == 1 or value-1 なら続行  
            ulong t = ModPow(a, d, value);
            if (t == 1 || t == n1)
            {
                continue;
            }

            // t^(2^i) mod value について、  
            // == value-1 の要素があれば続行、なければ合成数  
            int i;
            for (i = 0; i < s; i++)
            {
                t = ModMul(t, t, value);
                if (t == n1)
                {
                    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;
    }
}
0