結果

問題 No.8030 ミラー・ラビン素数判定法のテスト
ユーザー Jashinchan
提出日時 2022-08-29 12:45:24
言語 C
(gcc 13.3.0)
結果
WA  
実行時間 -
コード長 4,534 bytes
コンパイル時間 1,017 ms
コンパイル使用メモリ 33,280 KB
実行使用メモリ 6,820 KB
最終ジャッジ日時 2024-11-06 10:51:23
合計ジャッジ時間 1,322 ms
ジャッジサーバーID
（参考情報） 		judge3 / judge2
このコードへのチャレンジ
（要ログイン）
ファイルパターン 結果
other AC * 9 WA * 1
権限があれば一括ダウンロードができます
コンパイルメッセージ
main.c: In function 'in':
main.c:72:16: warning: implicit declaration of function 'getchar_unlocked' [-Wimplicit-function-declaration]
   72 |     while (c = getchar_unlocked(), c < 48 || c > 57);
      |                ^~~~~~~~~~~~~~~~
main.c: In function 'out':
main.c:83:5: warning: implicit declaration of function 'putchar_unlocked' [-Wimplicit-function-declaration]
   83 |     putchar_unlocked(x - ((((__uint128_t)x * 14757395258967641293ull) >> 3) >> 64) * 10 + 48);
      |     ^~~~~~~~~~~~~~~~

ソースコード

diff # 

#include <stdint.h>
#include <stdio.h>

uint64_t gcd64(uint64_t a, uint64_t b)
{
    if (!a || !b)
        return a | b;
    uint64_t t;
    uint64_t sh = __builtin_ctzll(a | b);
    a >>= __builtin_ctzll(a);
    do {
        b >>= __builtin_ctzll(b);
        if (a > b) t = a, a = b, b = t;
        b -= a;
    } while (b);
    return a << sh;
}

uint32_t mr32(uint64_t A, uint32_t n, uint32_t ninv)
{
    uint32_t y = (uint32_t)(A >> 32) - (uint32_t)(((uint64_t)((uint32_t)A * ninv) * n) >> 32);
    return (int32_t)y < 0 ? y + n : y;
}

uint32_t mul32(uint32_t a, uint32_t b, uint32_t n, uint32_t ninv)
{
    return mr32((uint64_t)a * b, n, ninv);
}

uint32_t pow32(uint32_t a, uint32_t k, uint32_t n, uint32_t ninv, uint32_t r)
{
    uint32_t ret = r;
    while (k > 0)
    {
        if (k & 1)
            ret = mul32(ret, a, n, ninv);
        a = mul32(a, a, n, ninv);
        k >>= 1;
    }
    return ret;
}


uint64_t mr64(__uint128_t A, uint64_t n, uint64_t ninv)
{
    uint64_t y = (uint64_t)(A >> 64) - (uint64_t)(((__uint128_t)((uint64_t)A * ninv) * n) >> 64);
    return (int64_t)y < 0 ? y + n : y;
}

uint64_t mul64(uint64_t a, uint64_t b, uint64_t n, uint64_t ninv)
{
    return mr64((__uint128_t)a * b, n, ninv);
}

uint64_t pow64(uint64_t a, uint64_t k, uint64_t n, uint64_t ninv, uint64_t r)
{
    uint64_t ret = r;
    while (k > 0)
    {
        if (k & 1)
            ret = mul64(ret, a, n, ninv);
        a = mul64(a, a, n, ninv);
        k >>= 1;
    }
    return ret;
}


uint64_t in(void)
{
    uint64_t c, x = 0;
    while (c = getchar_unlocked(), c < 48 || c > 57);
    while (47 < c && c < 58)
    {
        x = x * 10 + c - 48;
        c = getchar_unlocked();
    }
    return x;
}
void out(uint64_t x)
{
    if (x >= 10) out((((__uint128_t)x * 14757395258967641293ull) >> 3) >> 64);
    putchar_unlocked(x - ((((__uint128_t)x * 14757395258967641293ull) >> 3) >> 64) * 10 + 48);
}

int miller_rabin32(uint32_t n)
{
    static uint32_t base[] = { 2u, 7u, 61u };
    int s = __builtin_ctz(n - 1);
    uint32_t d = (n - 1) >> s;

    uint32_t r = (uint32_t)-1u % n + 1;
    uint32_t r2 = (uint64_t)(int64_t)-1 % n + 1;
    uint32_t ninv = n;
    for (int _ = 0; _ < 4; ++_) ninv *= 2 - ninv * n;
    
    for (int i = 0; i < 3; ++i)
    {
        if (base[i] >= n) break;
        uint32_t a = mr32((uint64_t)r2 * base[i], n, ninv);
        uint32_t c = pow32(a, d, n, ninv, r);
        if (c == r) continue;
        int f = 0;
        for (int q = 0; q < s && f == 0; q++)
        {
            f |= (c == (n - r));
            c = mul32(c, c, n, ninv);
        }
        if (f == 0) return 0;
    }

    return 1;
}
int miller_rabin64(uint64_t n)
{
    static uint64_t base[] = { 325, 9375, 28178, 450775, 9780504, 1795265022 };
    
    int s = __builtin_ctzll(n - 1);
    uint64_t d = (n - 1) >> s;

    uint64_t r = (uint64_t)-1ull % n + 1;
    uint64_t r2 = (__uint128_t)(__int128_t)-1 % n + 1;
    uint64_t ninv = n;
    for (int _ = 0; _ < 5; ++_) ninv *= 2 - ninv * n;

    {
        uint64_t l = (n - 1) << __builtin_clzll(n - 1);
        uint64_t t = r << 1;
        if (t >= n) t -= n;
        for (l <<= 1; l; l <<= 1)
        {
            t = mul64(t, t, n, ninv);
            if (l >> 63)
            {
                t <<= 1;
                if (t >= n) t -= n;
            }
        }
        if (t != r)
        {
            uint64_t x = (n - 1) & -(n - 1);
            uint64_t rev = n - r;
            for (x >>= 1; t != rev; x >>= 1)
            {
                if (x == 0)
                    return 0;
                t = mul64(t, t, n, ninv);
            }
        }
    }

    for (int i = 0; i < 6; ++i)
    {
        if (base[i] >= n) break;
        uint64_t a = mr64((__uint128_t)r2 * base[i], n, ninv);
        uint64_t c = pow64(a, d, n, ninv, r);
        if (c == r) continue;
        int f = 0;
        for (int q = 0; q < s && f == 0; q++)
        {
            f |= (c == (n - r));
            c = mul64(c, c, n, ninv);
        }
        if (f == 0) return 0;
    }
    return 1;
}
int is_prime(uint64_t n)
{
    if (n < 2ull) return 0;
    if (n < 4ull) return 1;
    if (n & 1 == 0) return 0;

    if (gcd64(n, 15) != 1)
        return 0;
    
    if (n < 4294967296ull)
        return miller_rabin32(n);
    
    else
        return miller_rabin64(n);
}

int main(void)
{
    uint64_t Q = in();
    while (Q--) {
        uint64_t x = in();
        out(x);
        putchar_unlocked(' ');
        out(is_prime(x));
        putchar_unlocked('\n');
    }
    return 0;
}
0