hyperion/80-math/miller-rabin.md at incoming

2025-12-28 19:00:03 +03:00

19 KiB

Raw Permalink Blame History

пример кода, близкого к эталонному:

fn main() {
    let n: u64 = std::env::args().nth(1)
        .and_then(|s| s.parse().ok())
        .unwrap_or(35);

    println!("{} -> {}", n, if is_prime(n) { "prime" } else { "composite" });
}

pub fn is_prime(n: u64) -> bool {
    if n < 2 { return false; }
    if n == 2 || n == 3 { return true; }
    if n % 2 == 0 { return false; }

    let s = (n - 1).trailing_zeros();
    let d = (n - 1) >> s;

    let witnesses = get_witnesses(n);

    for &a in witnesses {
        if !miller_rabin_test(a, d, n, s) {
            return false;
        }
    }
    true
}

// O(1) lookup based on input size
fn get_witnesses(n: u64) -> &'static [u64] {
    const WITNESSES: &[(u64, &[u64])] = &[
        (2_046, &[2]),
        (1_373_652, &[2, 3]),
        (9_080_190, &[31, 73]),
        (25_326_000, &[2, 3, 5]),
        (4_759_123_140, &[2, 7, 61]),
        (1_112_004_669_632, &[2, 13, 23, 1662803]),
        (2_152_302_898_746, &[2, 3, 5, 7, 11]),
        (3_474_749_660_382, &[2, 3, 5, 7, 11, 13]),
        (341_550_071_728_320, &[2, 3, 5, 7, 11, 13, 17]),
        (3_825_123_056_546_413_050, &[2, 3, 5, 7, 11, 13, 17, 19, 23]),
        (u64::MAX, &[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37]),
    ];

    WITNESSES.iter()
        .find(|(limit, _)| *limit >= n)
        .map(|(_, w)| *w)
        .unwrap()
}

fn miller_rabin_test(a: u64, d: u64, n: u64, s: u32) -> bool {
    let mut x = mod_pow(a, d, n);
    if x == 1 || x == n - 1 {
        return true;
    }

    for _ in 1..s {
        x = mod_mul(x, x, n);
        if x == n - 1 {
            return true;
        }
    }
    false
}

// === Math Core ===

#[inline]
fn mod_mul(a: u64, b: u64, m: u64) -> u64 {
    ((a as u128 * b as u128) % m as u128) as u64
}

fn mod_pow(mut base: u64, mut exp: u64, m: u64) -> u64 {
    let mut res = 1;
    base %= m;
    while exp > 0 {
        if exp % 2 == 1 {
            res = mod_mul(res, base, m);
        }
        base = mod_mul(base, base, m);
        exp /= 2;
    }
    res
}

Что такое Miller-Rabin (с нуля, метафора "Судья и Лжецы")

Представь простое число как честного судью, а составное — как лжеца. У тебя есть несколько свидетелей (маленькие числа 2, 3, 7...), которые задают судье сложные вопросы.

Правило игры: Если судья честный (число простое), он всегда отвечает правильно. Если лжец (составное), то с высокой вероятностью хотя бы один свидетель поймает его на лжи.

Miller-Rabin — это математический допрос: мы заставляем число доказать свою честность, вычисляя сложные степени по модулю.

Шаг 1: nth(n) — Простой перебор с MR-тестом

pub fn nth(n: u32) -> u32 {
    let mut num = 1;  // начинаем с 2 (1-е простое)
    for _ in 0..=n {  // ищем n+1 простое число
        loop {
            num += 1;
            if miller_rabin(num as u64) { break; }  // MR говорит "простое!"
        }
    }
    num
}

Логика: "Ходи по числам, пока не найдешь n простых". Медленно (O(n log n)), но корректно. Для Exercism (n до 10k) работает мгновенно.

Шаг 2: Подготовка MR-теста (d и s)

let mut d = n - 1;
let mut s = 0;
while d % 2 == 0 { d /= 2; s += 1 }

Математика: Каждое простое удовлетворяет малой теореме Ферма: a^{n-1} \equiv 1 \pmod n.

Но вместо прямого вычисления a^{n-1} мы раскладываем показатель:

 n-1 = d \cdot 2^s

где d — нечетное. Почему? Чтобы проверять промежуточные квадраты.

Метафора: "Не спрашивай сразу 'ты честный?', а сначала проверь базовое свойство, потом возводи в квадрат шаг за шагом".

Шаг 3: Witness Table (функциональный стиль)

let witnesses = WITNESSES.iter()
    .find(|&&(hi, _)| hi >= n)  // ищем нужный диапазон
    .map(|&(_, wtnss)| wtnss)   // извлекаем список свидетелей
    .unwrap();

Гениальность: Для каждого диапазона известен минимальный набор свидетелей, гарантирующий 100% точность.

Диапазон	Свидетели	Тестов
< 2047		1
< 4.7e9	[1]	3
< u64::MAX	[7][1][2][3][4][5][6]	12

Функциональный стиль: iter().find().map() — чистый, без мутабельности, легко тестировать.

Шаг 4: Основной цикл MR (сердце алгоритма)

for &a in witnesses.iter() {  // для каждого свидетеля
    let mut power = mod_exp(a, d, n);  // x0 = a^d mod n
    
    if power == 1 || power == n-1 { continue }  // тест пройден
    
    for _r in 0..s {  // s раз возводим в квадрат
        power = mod_sqr(power, n);            // xr+1 = xr^2 mod n
        if power == n-1 { continue 'next_witness }  // тест пройден!
        if power == 1 { return false }        // составное!
    }
    return false;  // все свидетели поймали на лжи
}

Логика теста:

Вычисляем x_0 = a^d \mod n
Если x_0 = 1 или x_0 = n-1 → свидетель доволен
Иначе делаем s квадратов: x_{r+1} = x_r^2 \mod n
Если хоть раз получили n-1 → свидетель доволен
Если дошли до конца → лжец пойман!

Шаг 5: Быстрое возведение в степень (Binary Exponentiation)

fn mod_exp(mut x: u64, mut d: u64, n: u64) -> u64 {
    let mut ret: u64 = 1;
    while d != 0 {
        if d % 2 == 1 {  // если бит степени = 1
            ret = mod_mul(ret, x, n);  // ret *= x
        }
        d /= 2;                // сдвиг степени
        x = mod_sqr(x, n);     // x = x^2 (для следующего бита)
    }
    ret
}

Метафора: "Не умножай 2^1000 = 2×2×2... 1000 раз. Используй (2^2)^2 = 2^4, ((2^4)^2)^2 = 2^16 и т.д."

O(log d) умножений вместо O(d).

Шаг 6: Безопасное умножение по модулю (u64 → u128)

fn mod_mul(a: u64, b: u64, m: u64) -> u64 {
    match a.checked_mul(b) {
        Some(r) => if r >= m { r % m } else { r },
        None => mod_mul_(a, b, m),  // fallback на u128
    }
}

Проблема: 2^{64}-1 \times 2^{64}-1 = 2^{128}, переполнение! Решение: checked_mul + fallback на u128::from(a) * u128::from(b) % u128::from(m).

Оптимизация mod_sqr: Для a < 2^{32} обычное a*a быстрее u128.

Итог: Почему код Production-Ready

✅ Детерминированный (не вероятностный)
✅ O(√n) на свидетель (с log n на mod_exp)
✅ Безопасен от overflow (u128 fallback)
✅ Функциональный стиль (иммутабельные итераторы)
✅ Таблица свидетелей покрывает весь u64
✅ Микрооптимизации (mod_sqr для малых чисел)

Когда использовать: Проверка конкретных больших чисел (10^18 + 7?).
Когда НЕ использовать: Генерация последовательности простых — тут sieve в 1000x быстрее.

Этот код — золотой стандарт для криптографии и системного программирования.[4]

1 2 3 4 5 6 7

Малая теорема Ферма (шаг 2 MR)

Формулировка:
Если p — простое число и a не делится на p, то:

 a^{p-1} \equiv 1 \pmod{p}

или эквивалентно:
[ a^p \equiv a \pmod{p} ] [1]

Метафора "Кольцо часов"

Представь p = 5 простое число как циферблат с 5 делениями (0,1,2,3,4). Умножение на a = 3 — это поворот стрелки:

Стрелка на 1 → ×3 = 3
Стрелка на 2 → ×3 = 6 ≡ 1 (mod 5)  
Стрелка на 3 → ×3 = 9 ≡ 4 (mod 5)
Стрелка на 4 → ×3 = 12 ≡ 2 (mod 5)

Важно: Поворот на 1,2,3,4 дал все возможные позиции (1,3,4,2) — полный цикл! Умножение на 3 проходится по всему кольцу.

Теорема говорит: Если p простое, то любое a (не кратное p) за p-1 повортов вернется в исходную точку.

Пример вычисления

// p = 13 (простое), a = 5
5^12 mod 13 = ?  // малая теорема: должно быть 1

5^1 ≡ 5
5^2 ≡ 25 ≡ 12 ≡ -1  
5^4 ≡ (-1)^2 ≡ 1
5^8 ≡ 1^2 ≡ 1
5^12 = 5^8 × 5^4 ≡ 1 × 1 ≡ 1 ✓

Почему это работает (доказательство через перестановки)

Ключевой инсайт: Остатки 1, 2, ..., p-1 при умножении на a дают перестановку тех же чисел:

p = 7, a = 3
1×3 = 3
2×3 = 6  
3×3 = 9 ≡ 2
4×3 = 12 ≡ 5
5×3 = 15 ≡ 1
6×3 = 18 ≡ 4

Получили: [3,6,2,5,1,4] — все числа 1..6 в другом порядке!

Произведение всех элементов одинаково:

 1 \times 2 \times \dots \times (p-1) = a \times 2a \times \dots \times (p-1)a \pmod{p}

 (p-1)! \equiv a^{p-1} \times (p-1)! \pmod{p}

Сокращаем (p-1)! (оно не 0 mod p):

 1 \equiv a^{p-1} \pmod{p}

Почему MR использует именно эту теорему

Прямая проверка a^{p-1} \mod p неэффективна — показатель степени огромный (до 2^{64}).

Хитрость MR: Вместо a^{p-1} проверяем промежуточные квадраты:

 n-1 = d \times 2^s

 a^{n-1} = (a^d)^{2^s}

Если a^d \equiv 1 или хотя бы один квадрат дает n-1, то по свойствам степеней теорема выполняется.

Метафора: Вместо "сделай 1 трлн повортов стрелки" говорим "сделай 40 повортов, но проверяй промежуточные позиции".

Проверка на практике

fn fermat_test(n: u64, a: u64) -> bool {
    mod_exp(a, n-1, n) == 1  // прямая малая теорема
}

// Но! 561 — псевдопростое по основанию 2:
// 2^560 ≡ 1 mod 561, хотя 561 = 3×11×17

Вывод: Малая теорема — необходимое условие простоты, но недостаточное (есть "лжецы"). MR усиливает тест структурированными промежуточными проверками + несколькими основаниями.

Теперь понятен шаг 2? Идем к witness table? 🚀

1 2 3 4 5 6 7 8 9 10

Ты нащупал самую суть! Это называется "Binary Exponentiation" (Бинарное возведение в степень).

Твоя догадка: "мы ищем минимальную степень... домножаем... переиспользуем" — абсолютно верна.

Давай разложим это на механику "Сборщика" и "Генератора".

1. Генератор: Как мы переиспользуем `base`?

Посмотри на строку:

rust

base = mod_mul(base, base, m);

В каждой итерации цикла эта строка делает одно и то же: возводит текущую базу в квадрат.

Представь, что на входе base — это число X.

Итерация 1: base становится X^2.
Итерация 2: Мы берем этот X^2 и умножаем сам на себя: (X^2) \cdot (X^2) = X^4.
Итерация 3: Берем X^4 и умножаем сам на себя: (X^4) \cdot (X^4) = X^8.
Итерация 4: ... X^{16}.

Вот оно, переиспользование! Вместо того чтобы умножать X \cdot X \cdot X \dots (линейно), мы прыгаем по степеням двойки: 1, 2, 4, 8, 16, 32...
Мы генерируем эти "кирпичики" (X^1, X^2, X^4, X^8) на лету.

2. Сборщик: Как мы собираем результат `res`?

Теперь посмотри на условие:

rust

if exp % 2 == 1 { res = mod_mul(res, base, m); }

Здесь мы решаем: нужен нам этот кирпичик или нет?

Допустим, нам нужно вычислить 3^{13}.
Двоичное представление 13 — это 1101.
Это значит:
13=1⋅8+1⋅4+0⋅2+1⋅113 = 1 \cdot 8 + 1 \cdot 4 + 0 \cdot 2 + 1 \cdot 113=1⋅8+1⋅4+0⋅2+1⋅1
То есть:
313=38⋅34⋅(пропускаем 32)⋅313^{13} = 3^8 \cdot 3^4 \cdot (\text{пропускаем } 3^2) \cdot 3^1313=38⋅34⋅(пропускаем 32)⋅31

Мы просто берем только те степени, где в двоичном коде стоит 1.

Визуализация процесса (Трассировка)

Давай пройдем этот цикл для 3^{13} \pmod{100}.
Изначально: res = 1, base = 3. exp = 13 (1101 bin).

Итерация	exp (dec)	exp (bin)	Бит (четность)	Действие с `res` (Сборщик)	Состояние `res`	Действие с `base` (Генератор)	Новое `base`
1	13	`...1101`	1 (нечет)	Нужен! Умножаем на `3^1`	`1 \cdot 3 = \mathbf{3}`	`3 \cdot 3` (квадратируем)	`3^2 = 9`
2	6	`...0110`	0 (чет)	Не нужен. Пропускаем.	`\mathbf{3}`	`9 \cdot 9` (квадратируем)	`3^4 = 81`
3	3	`...0011`	1 (нечет)	Нужен! Умножаем на `3^4`	`3 \cdot 81 = 243 \to \mathbf{43}`	`81 \cdot 81 \to 61`	`3^8 = 61`
4	1	`...0001`	1 (нечет)	Нужен! Умножаем на `3^8`	`43 \cdot 61 \to \mathbf{23}`	`61 \cdot 61 \dots`	`3^{16}` (не важно)

Итог: Мы получили 23.
Проверка: 3^{13} = 1,594,323.
1,594,323 \pmod{100} = 23. Сошлось!

Доказательство формулы остатков

Ты спрашиваешь: "Как мы доказываем, что взятие остатка на каждом шаге не ломает математику?"

Это фундаментальное свойство модульной арифметики:
Остаток произведения равен произведению остатков.

(A⋅B)(modM)=((A(modM))⋅(B(modM)))(modM)(A \cdot B) \pmod M = ((A \pmod M) \cdot (B \pmod M)) \pmod M(A⋅B)(modM)=((A(modM))⋅(B(modM)))(modM)

Пример на пальцах:
Пусть M=10 (мы смотрим только на последнюю цифру).
Мы хотим узнать последнюю цифру 12 \cdot 13.

По-честному: 12 \cdot 13 = 156. Последняя цифра 6.
По-модульному:
- От 12 берем остаток \to 2.
- От 13 берем остаток \to 3.
- Умножаем остатки: 2 \cdot 3 = 6.
- Результат тот же: 6.

Именно это правило разрешает нам делать base %= m и mod_mul на каждом шаге. Мы просто отбрасываем "лишние круги", которые никак не влияют на конечный остаток. Без этого правила вся современная криптография (RSA, Diffie-Hellman) просто не работала бы на компьютерах.

19 KiB Raw Permalink Blame History Unescape Escape