425 lines
19 KiB
Markdown
425 lines
19 KiB
Markdown
|
||
|
||
пример кода, близкого к эталонному:
|
||
```rust
|
||
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-тестом
|
||
|
||
```rust
|
||
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)
|
||
|
||
```rust
|
||
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 (функциональный стиль)
|
||
|
||
```rust
|
||
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 (сердце алгоритма)
|
||
|
||
```rust
|
||
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; // все свидетели поймали на лжи
|
||
}
|
||
```
|
||
|
||
**Логика теста:**
|
||
1. Вычисляем \(x_0 = a^d \mod n\)
|
||
2. Если \(x_0 = 1\) или \(x_0 = n-1\) → свидетель доволен
|
||
3. Иначе делаем \(s\) квадратов: \(x_{r+1} = x_r^2 \mod n\)
|
||
4. Если хоть раз получили \(n-1\) → свидетель доволен
|
||
5. Если дошли до конца → **лжец пойман!**
|
||
|
||
## Шаг 5: Быстрое возведение в степень (Binary Exponentiation)
|
||
|
||
```rust
|
||
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)
|
||
|
||
```rust
|
||
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](https://github.com/SergioBenitez/Rocket/issues/1522)
|
||
[2](https://gliderkite.github.io/posts/sr-cache/)
|
||
[3](https://stackoverflow.com/questions/16761898/redis-as-write-back-view-count-cache-for-mysql)
|
||
[4](https://github.com/jedisct1/rust-sieve-cache)
|
||
[5](https://dev.to/artslob/32-rust-crates-you-better-know-about-4aa7)
|
||
[6](https://lib.rs/crates/staticvec)
|
||
[7](https://www.perplexity.ai/search/444a1554-6943-4987-9a27-5923865aea84)
|
||
|
||
|
||
## Малая теорема Ферма (шаг 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 \) повортов **вернется в исходную точку**.
|
||
|
||
### Пример вычисления
|
||
|
||
```rust
|
||
// 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 \) дают **перестановку** тех же чисел:
|
||
|
||
```rust
|
||
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 повортов, но проверяй промежуточные позиции".
|
||
|
||
### Проверка на практике
|
||
|
||
```rust
|
||
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](https://ru.wikipedia.org/wiki/%D0%9C%D0%B0%D0%BB%D0%B0%D1%8F_%D1%82%D0%B5%D0%BE%D1%80%D0%B5%D0%BC%D0%B0_%D0%A4%D0%B5%D1%80%D0%BC%D0%B0)
|
||
[2](https://math.mosolymp.ru/upload/files/2022/other/1568/8-1/Malaya_teorema_Ferma.pdf)
|
||
[3](https://foxford.ru/wiki/matematika/malaya-teorema-ferma)
|
||
[4](https://www.lirmm.fr/~ashen/mathclass/09fermat-prob.pdf)
|
||
[5](https://wiki.algocode.ru/index.php?title=%D0%9C%D0%B0%D0%BB%D0%B0%D1%8F_%D1%82%D0%B5%D0%BE%D1%80%D0%B5%D0%BC%D0%B0_%D0%A4%D0%B5%D1%80%D0%BC%D0%B0)
|
||
[6](https://kvant.mccme.ru/pdf/2000/04/kv0400ferm.pdf)
|
||
[7](https://ru.wikipedia.org/wiki/%D0%A2%D0%B5%D0%BE%D1%80%D0%B5%D0%BC%D0%B0_%D0%A4%D0%B5%D1%80%D0%BC%D0%B0_%E2%80%94_%D0%AD%D0%B9%D0%BB%D0%B5%D1%80%D0%B0)
|
||
[8](https://journalpro.ru/articles/malaya-teorema-ferma-i-ee-primenenie-v-kriptosistemakh/)
|
||
[9](https://www.reddit.com/r/math/comments/2v6db7/intuition_for_fermats_little_theorem/)
|
||
[10](https://math.fandom.com/ru/wiki/%D0%9C%D0%B0%D0%BB%D0%B0%D1%8F_%D1%82%D0%B5%D0%BE%D1%80%D0%B5%D0%BC%D0%B0_%D0%A4%D0%B5%D1%80%D0%BC%D0%B0)
|
||
|
||
Ты нащупал самую суть! Это называется **"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$.
|
||
|
||
1. **По-честному:** $12 \cdot 13 = 156$. Последняя цифра **6**.
|
||
|
||
2. **По-модульному:**
|
||
|
||
- От 12 берем остаток $\to 2$.
|
||
|
||
- От 13 берем остаток $\to 3$.
|
||
|
||
- Умножаем остатки: $2 \cdot 3 = 6$.
|
||
|
||
- Результат тот же: **6**.
|
||
|
||
|
||
Именно это правило разрешает нам делать `base %= m` и `mod_mul` на каждом шаге. Мы просто отбрасываем "лишние круги", которые никак не влияют на конечный остаток. Без этого правила вся современная криптография (RSA, Diffie-Hellman) просто не работала бы на компьютерах.
|