Детерминированная реализация Миллера-Рабина

Я пытаюсь реализовать функцию проверки простоты с помощью детерминированного алгоритма Миллера-Рабина, но результаты не всегда правильные: при проверке первых 1 000 000 чисел он находит только 78 495 вместо 78 498.

Это получается с использованием [2, 7, 61] в качестве базы, которая, согласно Википедии, всегда должна быть правильной для значений до 4 759 123 141. 2, 7 и 61).

Почему это происходит? Код, который я использую, следующий:

T modular_power(T base, T exponent, T modulo) {
    base %= modulo;
    T result = 1;

    while (exponent > 0) {
        if (exponent % 2 == 1)
            result = (result * base) % modulo;
        base = (base * base) % modulo;
        exponent /= 2;
    }

    return result;
}

bool miller_rabin(const T& n, const vector<T>& witnesses) {
    unsigned int s = 0;
    T d = n - 1;
    while (d % 2 == 0) {
        s++;
        d /= 2;
    }

    for (const auto& a : witnesses) {
        if (modular_power<T>(a, d, n) == 1)
            continue;

        bool composite = true;
        for (unsigned int r = 0; r < s; r++) {
            if (modular_power<T>(a, (T) pow(2, r) * d, n) == n - 1) {
                composite = false;
                break;
            }
        }

        if (composite)
            return false;
    }

    return true;
}

bool is_prime(const T& n) {
    if (n < 4759123141)
        return miller_rabin(n, {2, 7, 61});
    return false; // will use different base
}

person Becks    schedule 13.01.2018    source источник
comment
Далее - М-Р утверждает: {n in odd | n >= 3}. Таким образом, одно простое отклонение в is_prime: if ((n & 0x1) == 0) return (n == 2); ... отбрасывает все четные как составные, кроме (2). Детерминистический тест явно проходит случай: if ((a %= n) == 0) (поскольку этот тест ничего нам не говорит) и далее гарантирует, что (0 < a < n) с этой точки.   -  person Brett Hale    schedule 14.01.2018
comment
Мой детерминированный тест для всех беззнаковых 64-битных значений. Вы можете вернуться к uint64_t, если тестируете только 32-битные значения.   -  person Brett Hale    schedule 14.01.2018
comment
@BrettHale Я использую это для задач Project Euler, и для них требуются довольно большие числа, поэтому uint64_t может быть недостаточно. Ваши советы очень хороши, я уже реализовал четную проверку и, вероятно, воспользуюсь еще парой ваших оптимизаций, спасибо за помощь!   -  person Becks    schedule 14.01.2018


Ответы (1)


Миллер-Рабин действительно не работает, когда база и вход одинаковы. Что происходит в этом случае, так это то, что ad mod n равен нулю (поскольку mod n равен нулю, так что на самом деле это возведение нуля в нерелевантную степень), а остальная часть алгоритма не может " escape" с нуля и делает вывод, что вы имеете дело с композитом.

Как частный случай этого, Миллер-Рабин никогда не работает со входом 2, потому что нет базы, которую можно выбрать. 2 сам по себе бесполезен, как и 1, который ничего не оставляет.

person harold    schedule 13.01.2018