Тест Миллера-Рабина: как проверить простоту числа до 10^18 (обзор и рабочий код)?

Миллера-Рабина — вероятностный тест простоты, но для 64-битных чисел его можно сделать детерминированным, проверив фиксированный набор оснований. Идея: для простого p (нечётного) запишем p - 1 = d·2^s. Тогда для любого основания a (не кратного p) выполняется: либо a^d ≡ 1, либо a^{d·2^r} ≡ -1 для какого-то 0 ≤ r < s (следствие малой теоремы Ферма и того, что ±1 — единственные корни из 1 по простому модулю). Если для основания a это не так — n составное (a — «свидетель»).

Для n < 3.3·10^24 достаточно набора оснований {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37} — это покрывает все 64-битные числа детерминированно.

using u64 = unsigned long long; using u128 = __uint128_t;
u64 mulmod(u64 a, u64 b, u64 m){ return (u128)a*b%m; }
u64 powmod(u64 a, u64 b, u64 m){ u64 r=1; a%=m; while(b){ if(b&1) r=mulmod(r,a,m); a=mulmod(a,a,m); b>>=1; } return r; }
bool check(u64 n, u64 a, u64 d, int s){
    u64 x = powmod(a, d, n);
    if (x == 1 || x == n-1) return true;
    for (int i = 1; i < s; i++){ x = mulmod(x, x, n); if (x == n-1) return true; }
    return false;
}
bool is_prime(u64 n){
    if (n < 2) return false;
    for (u64 p : {2,3,5,7,11,13,17,19,23,29,31,37}) { if (n % p == 0) return n == p; }
    u64 d = n-1; int s = 0; while ((d & 1) == 0){ d >>= 1; s++; }
    for (u64 a : {2,3,5,7,11,13,17,19,23,29,31,37}) if (!check(n, a, d, s)) return false;
    return true;
}

Сложность — O(k·log n) умножений, где k — число оснований (12). Обязательно используй __uint128_t в mulmod, иначе a*b при n ~ 10^18 переполнит u64.