Сложно ли генерировать 1024-битные простые числа?

Простые числа удивительны!

С одной стороны, их легко объяснить: это просто числа, которые делятся только на единицу и на себя; с другой стороны, они содержат в себе бесконечную сложность. Они встречаются во множестве разных сфер, от математических концепций и гипотез до любопытных визуализаций и криптографии, лежат в основе многих Интернет-стандартов и протоколов безопасности, которые мы используем ежедневно.

Несмотря на моё восхищение простыми числами, я никогда не исследовал их подробно, поэтому решил бросить себе вызов. А есть ли лучший способ изучить простые числа, чем воспользоваться моей любовью к кодингу для их генерации?

Вызов

Но какие простые числа мне генерировать? Нахождение миллиардного числа — это слишком просто, а попадание в таблицу рекордов самых больших известных простых чисел — слишком амбициозная задача для первой попытки. Соединив простые числа и интерес к криптографии, я придумал такую задачу:

Генерировать простые числа, способные генерировать ключи для алгоритма RSA

На момент написания этой статьи хорошей длиной ключей RSA считаются 2048 битов. Ключи RSA генерируются перемножением двух простых чисел, так что для получения 2048-битного ключа нам нужны два числа длиной примерно 1024 бита. Это ограничивает рамки задачи генерацией 1024-битных простых чисел. Теперь вы знаете, откуда взялось число из заголовка поста.

В дополнение к задаче я установил для себя некоторые правила:

• Код должен быть написан с нуля; в противном случае можно было бы просто написать openssl prime -generate -bits 1024! «С нуля» здесь означает отсутствие внешних зависимостей.

• Никакого мощного внешнего оборудования или облаков; то есть мы не можем решить задачу «грубой силой», вбросив в неё дополнительные вычислительные мощности. Я буду пользоваться своим ноутбуком с AMD Ryzen 7 CPU и 16 ГБ ОЗУ.

• Генерировать простые числа за «разумное» время: определение намеренно оставлено расплывчатым: так я буду я стремиться к тому, чтобы оптимизировать код, но не попаду в ловушку чрезмерной оптимизации.

В качестве языка программирования я выбрал Rust, в основном потому, что недавно его изучал, а моя задача показалась хорошей практикой. Я считаю, что Rust достаточно низкоуровневый, чтобы экспериментировать с глубокими концепциями, и в то же время достаточно высокоуровневый, чтобы код было относительно легко понять. Я не буду применять сложные фичи Rust, потому что пока незнаком с ними.

Итак, с условиями мы разобрались, давайте приступать!

16 битов — самое лёгкое!

Я планирую постепенно подбираться к 1024 битам, поэтому в качестве разогрева начал с 16 битов. Теоретически, процесс генерации любых N-битных простых чисел элементарен:

while true {
    let number = <<random N-bit integer>>
    if is_prime(number) {
        break
    }
}

Генерируем новые случайные N-битные числа, пока не найдём то, которое пройдёт тест на простоту. Ещё до того, как я приступил к созданию тестов, у меня возникла первая проблема: откуда вообще брать случайные числа? У Rust есть превосходный крейт (что-то типа библиотеки/пакета) под названием rand, который практически можно считать частью стандартной библиотеки. Но прежде чем нарушить своё правило «никаких зависимостей», я решил, что сначала стоит хотя бы попробовать сделать это самостоятельно.

Я вспомнил, что где-то слышал о /dev/urandom, и после изучения этого устройства оказалось, что оно идеально мне подходит. В ядре Linux есть встроенный Cryptographically Secure Pseudo Random Number Generator (CSPRNG), доступ к которому осуществляется чтением из файла псевдоустройства /dev/urandom. Оно собирает энтропию из пользовательского окружения и использует её для периодического создания порождающего значения детерминированного потокового шифра ChaCha20 (забавное название!), который способен генерировать «истинные» случайные биты. Поначалу я не решался использовать его, но одна статья убедила меня.

Вот придуманная мной реализация:

// rng.rs
use std::fs::File;
use std::io::Read;
fn insert_random_bytes(mut bytes: &mut[u8]) -> std::io::Result<()> {
    File::open("/dev/urandom")?.read_exact(&mut bytes)?;
    Ok(())
}
fn u16() -> u16 {
    let mut bytes = [0u8; 2];
    insert_random_bytes(&mut bytes).expect("Cannot access /dev/urandom");
    u16::from_le_bytes(bytes)
}

Примечание: выражения без точки с запятой, расположенные в конце функций Rust, используются в качестве возвращаемого значения функции.

insert_random_bytes() получает в качестве входных данных изменяемый массив байтов и заполняет его выводом из /dev/urandom. Функция u16() создаёт буфер из 2 байтов (16 битов), заполняет буфер случайными битами, а затем создаёт из этих битов integer u16 (в Rust u16 обозначает беззнаковое 16-битное число). Затем u16() используется следующим образом:

fn run() {
    println!("random no - {}", rng::u16() | 0b1000000000000001);
}

Возвращаемое случайное число подвергается OR с 0b1000000000000001, чтобы его первый и последний бит были равны 1. Последний бит, равный 1, делает его нечётным числом, а первая единица гарантирует, что число будет достаточно большим и охватывать весь нужный мне диапазон битов.

Вот несколько сгенерированных 16-битных случайных чисел:

random no - 36111
random no - 52205
random no - 45689
random no - 33631

Теперь, когда у меня есть свой собственный генератор случайных чисел, можно быстро разобраться с 16-битными простыми числами. Для начала создадим enum для хранения результатов:

enum PrimeResult {
    Prime,
    Composite,
}

Затем применим простой тест на простоту под названием «перебор делителей» (trial division). Он выполняет цикл от 3 до sqrt(num) и проверяет, окажется ли одно из чисел делителем num:

fn trial_division_simple(n: u16) -> PrimeResult {
    let root_n = (n as f64).sqrt() as u16;
    for x in 3..root_n {
        if n % x == 0 {
            return PrimeResult::Composite;
        }
    }
    PrimeResult::Prime
}

И, наконец, простой цикл:

fn run() {
    loop {
        let num = rng::u16() | 0b1000000000000001;
        if trial_division_simple(num) == PrimeResult::Prime {
            println!("Prime found: {num}");
            break;
        }
    }
}

➜ time cargo run --release
Prime found: 44809
cargo run --release  0.03s user 0.01s system 99% cpu 0.038 total

Всё хорошо работает, а для генерации 16-битных простых чисел в среднем требуется около 40 мс. Чтобы убедиться, что мои простые числа на самом деле простые, я использовал крутой онлайн-тестер, написанный на WebAssembly (хостится здесь) вместе с командой OpenSSL: openssl prime <number>.

Итак, мы освежили в памяти основные концепции и немного разогрелись. Пока делать следующий шаг!

64 бита, в 4 раза больше битов!

После 16 битов я перескочил сразу к 64-битным числам. Сегодня 64-битная архитектура используется в большинстве современных вычислительных устройств, и использовав 64 бита, мы заберёмся далеко на территорию 20-значных чисел (для контекста: 1 триллион — это 13 знаков). Способен ли будет простой алгоритм перебора делителей справиться с такими большими числами?

➜ time cargo run --release
Prime found: 14288847644715868907
cargo run --release  30.27s user 0.02s system 99% cpu 30.294 total

Ну, типа того. 30 секунд на генерацию 64-битного простого числа — не очень хороший результат, но мы ещё не раскрыли полностью потенциал перебора делителей. В этом разделе мы попробуем выжать из него максимум.

Для начала вот более оптимизированная версия перебора делителей:

fn trial_division(n: u64, start: u64) -> PrimeResult {
    // assumption: n is odd, n > 3 and start > 3
    if n % 3 == 0 {
        return PrimeResult::Composite;
    }
    let root_n = (n as f64).sqrt() as u64;
    for x in (start..(root_n + 1)).step_by(6) {
        if n % x == 0 || n % (x + 2) == 0 {
            return PrimeResult::Composite;
        }
    }
    PrimeResult::Prime
}

Изменения:

• Теперь алгоритм получает параметр start , позволяющий указать, с какого числа начинать цикл.

• Цикл движется с шагом 6, а не 1.

• Наряду с n % x == 0 выполняется проверка n % (x + 2) == 0.

По сути, здесь учитываются только делители от start до sqrt(n), имеющие вид 6k+1. Цитата из Википедии:

Это вызвано тем, что все целые числа можно выразить в виде (6k+i), где i = −1, 0, 1, 2, 3 или 4. Стоит отметить, что 2 является делителем (6k+0), (6k+2), и (6k+4), а 3 является делителем (6k+3). То есть будет эффективнее проверять, делится ли n на 2 или 3, а затем проверять все числа формата 6k±1 ≤ √n. Это в 3 раза быстрее, чем проверять все числа до √n.

Далее напишем функцию для генерации списка небольших простых чисел при помощи этого усовершенствованного перебора делителей:

fn generate_small_primes<const N: usize>() -> [u64; N] {
    let mut primes: [u64; N] = [0; N];
    primes[0] = 2;
    primes[1] = 3;
    let mut n: u64 = 3;
    let mut nth: u64 = 2;
    let mut i: usize = 2;
    let limit = N as u64;
    loop {
        n += 2;
        if trial_division(n, 5) == PrimeResult::Prime {
            primes[i] = n;
            i += 1;
            nth += 1;
            if nth == limit {
                return primes
            }
        }
    }
}

Добавим в наш результат ещё одно возможное состояние:

enum PrimeResult {
    Prime,
    Composite,
    Unknown,   // <----- new
}

А теперь используем список малых простых чисел для предварительной проверки на легко делимые числа, прежде чем возвращаться к перебору делителей:

fn primes_64bit() {
    const N: usize = 10000;
    let start = (N + 1) as u64;
    let primes = utils::generate_small_primes::<N>();
    loop {
        let num = rng::u64() | 0x8000000000000001u64;
        let mut result = PrimeResult::Unknown;
        for i in 0..N {
            if num % primes[i] == 0 {
                result = PrimeResult::Composite;
                break;
            }
        }
        if result == PrimeResult::Unknown {
            result = algos::trial_division(num, start)
        }
        if result == PrimeResult::Prime {
            println!("Prime found: {num}");
            break;
        }
    }
}

После всей этой работы результат оказался таким:

➜ time cargo run --release
Prime found: 12589778476493955313
cargo run --release  6.40s user 0.01s system 99% cpu 6.414 total

Значительно лучше, чем предыдущие 30 секунд. У нас ещё есть возможности для оптимизации, но если даже для 64-битных чисел выполнение занимает 6 секунд, то очевидно, что его нельзя масштабировать до 1024-битных.

Итак, нам придётся покинуть уютный мир детерминированных алгоритмов и войти во вселенную неопределённости с вероятностными алгоритмами!

128 битов, но с неожиданным поворотом!

Теперь всё становится интереснее. Поначалу концепция вероятностных тестов простоты показалась мне странной, и я попытался найти детерминированные алгоритмы, способные справляться с большими числами. У меня получилось найти два: APR-CL и ECPP. Оба они настолько математически сложны, что мне вообще не удалось разобраться в научных статьях о них, а в Интернете практически нет доступной информации для тех, кто плох в математике.

Почитав онлайн-обсуждения, изучив исходный код OpenSSL и рекомендации NIST, я понял, что почти везде, в том числе и в RSA, используются вероятностные алгоритмы. Тонкость в том, что при правильной реализации эти алгоритмы имеют чрезвычайно низкую частоту ошибок, которой можно пренебречь. Все приведённые ниже алгоритмы не «доказывают», что число простое, а сообщают, что это «возможно простое» с определённой долей точности. Первым из этих алгоритмов я исследовал малую теорему Ферма.

Малая теорема Ферма

Эта теорема Ферма гласит: если p простое, a — любое целое, не делящееся на p, то число a^p-1 делится на p. В модульной арифметике это можно выразить так:

a^p-1 = 1 (mod p)

Мы можем выбирать разные значения a при a < p, чтобы a по определению не могло делиться на p, и теоретически подстановка этих значений в соответствие покажет, является ли p простым.

Нам достаточно реализовать это соответствие в коде, начав с a^p-1 и a = 2:

fn run() {
    let num = rng::u128() | 0x80000000000000000000000000000001u128;
    let base = 2u128;
    println!("{}", base.pow(num - 1);
}

➜ cargo run
...
error[E0308]: mismatched types
--> src/lib.rs:71:29
|
71  |     println!("{}", base.pow(num - 1);
|                         --- ^^^^^^^ expected u32, found u128
|                         |
|                         arguments to this function are incorrect

Мне понадобилась минута, чтобы понять, что ошибка была не с моей стороны. Функция pow() намеренно получает u32, так как возвещение u128 в любую более высокую степень уже переполняет u128! К счастью, наше соответствие записано в модульной арифметике, то есть мы можем брать модуль каждого этапа, в итоге получив результат меньше u128.

То есть, по сути

a × b (mod m) = [ a (mod m) × b (mod m) ] (mod m)

а значит

a^p-1 (mod p) = ((((a × a (mod p)) × a (mod p)) × a (mod p)) × ...... p - 1 раз)

Алгоритм для реализации этого называется возведением в степень по модулю. Я реализовал его, напрямую повторив псевдокод из Википедии, использовав версию, реализующую «возведение в степень возведением в квадрат», чтобы получить более эффективный алгоритм.

Попытка номер 2:

fn run() {
    let num = rng::u128() | 0x80000000000000000000000000000001u128;
    let base = rng::u128_range(2, num - 1);
    println!("{}", mod_exp(base, num - 1, num));
}

➜ cargo run
...
thread 'main' panicked at 'attempt to multiply with overflow', src/utils.rs:38:16
note: run with RUST_BACKTRACE=1 environment variable to display a backtrace

Хм.

Я не сообразил, что даже результат перемножения двух u128 запросто может казаться слишком большим для хранения в u128. Потерпев поражение, я решил двигаться дальше, сохраняя внутри u128 только 64-битные числа. Грубо говоря, максимальное пространство, необходимое для перемножения двух N-битных чисел — это 2N, отсюда и решение хранить два 64-битных числа в u128. Эта идея распределения удвоенного количества битов от необходимого тоже встретится ниже. Любопытно, что в предыдущем этапе с 64 битами и перебором делителей не было умножений, поэтому я не столкнулся с этой проблемой.

Добавляем в enum ещё одно возможное состояние:

enum PrimeResult {
    Prime,
    Composite,
    Unknown,
    ProbablePrime,   // <----- new
}

А вот реализация теста Ферма. Она просто выполняет уравнение k раз со случайными основаниями:

fn fermat_test(num: u128, k: usize) -> PrimeResult {
    for _ in 0..k {
        let base = rng::u128_range(2, num - 1);
        if mod_exp(base, num - 1, num) != 1 {
            return PrimeResult::Composite;
        }
    }
    PrimeResult::ProbablePrime;
}

Здесь есть интересная подробность реализации: функция rng::u128_range() подразумевает, что равномерно выбирает случайное u128 от 2 до num - 1 , но мне показалось более практичным возвращать непосредственно случайное число, которое на несколько байтов короче num. Это сильно упрощает логику, в то же время обеспечивая нам по большей мере случайное достаточно большое число от 2 до num - 1. В дальнейшем этот трюк будет использоваться там, где понадобятся случайные значения из интервала.

А вот тест целиком!

fn primes_128bit() -> u128 {
    loop {
        let num = (rng::u64() | 0x8000000000000001u64) as u128;
        if fermat_test(num, 10) == PrimeResult::ProbablePrime {
            return num;
        }
    }
}

➜ time cargo run --release
Prime found: 9944209443870115157
cargo run --release  0.03s user 0.01s system 99% cpu 0.033 total

Это сильно быстрее, чем прогоны по 6 секунд, которые мы ранее получали для 64 битов, но то, что мы получаем «вероятно простой» результат, могло подсказать вам, что здесь есть тонкость. Изъяном малой теоремы Ферма оказываются «псевдопростые числа». Соответствие, определённое малой теоремой Ферма, справедливо для всех простых чисел, но также справедливо и для некоторых составных. Если RNG генерирует одно из таких особых составных чисел, то код скажет, что оно простое, даже если это не так. Эти составные числа, также называемые «псевдослучайными числами Ферма», встречаются редко, но их всё равно достаточно много, чтобы мы не могли полагаться на точность теста Ферма.

Тест простоты Миллера — Рабина

Тест Миллера — Рабина — это усовершенствованный вероятностный тест простоты, работающий на тех же принципах, что и тест Ферма; однако он гораздо надёжнее и практичнее благодаря некоторым ключевым отличиям. Например, в этом тесте ни одно составное число не будет сильным псевдопростым для всех оснований одновременно, в отличие от теста Ферма, для которого такие составные существуют (их называют числами Кармайкла). Кроме того, тест Миллера — Рабина обладает гораздо меньшей частотой ошибок, которую в большинстве случаев можно считать «несущественной». На самом деле, изучив то, чем пользуются другие люди, например, в исходном коде OpenSSL и в рекомендациях NIST, то можно понять. что многие источники рекомендуют или уже используют тест Миллера — Рабина!

Математика, лежащая в основе теста Миллера — Рабина, не так важна для реализации самого алгоритма, но я, тем не менее, постараюсь изложить то, что понял. Можете пропустить эту часть и перейти непосредственно к коду.

Изученное нами соответствие в тесте Ферма выглядело так:

a^n-1 = 1 (mod n),  если n простое        (1)

В более общем виде a^n-1 можно записать как a^{2^s × d}, где:

• d = нечётное число, оставшееся после вынесения за скобки всех степеней двойки из n.

• s = степень двойки как делитель n.

Из чего следует:

n - 1 = 2^s × d       (2)

Соединив (1) и (2), мы можем сказать, что n является сильно вероятным простым числом, если истинно одно из следующих условий:

• a^d = 1 (mod n).

• a^{2^r × d} = n - 1 (mod n) для некоторых 0 <= r < s.

Примечание: n - 1 (mod n) эквивалентно -1 (mod n). Выражение оставлено в развёрнутом виде, потому что я работаю с беззнаковыми int и не могу представить -1.

По сути, вместо выполнения одной проверки с a^n-1 код выполняет множество проверок, начиная с a^{2^0 × d}, то есть a^d, затем a^{2^1 × d}, a^{2^2 × d}, a^{2^3 × d} и так далее, пока не дойдёт до n в a^{2^s × d}.

Вот моя реализация, полученная сочетанием приведённой выше математики с псевдокодом из Википедии:

fn miller_rabin_test(n: u128, k: usize) -> PrimeResult {
    let mut s = 0;
    let mut d = n - 1;
    while d % 2 == 0 {
        d = d / 2;
        s += 1;
    }
    'main_loop: for _ in 0..k {
        let base = rng::u128_range(2, n - 2);
        let mut x = utils::mod_exp(base, d, n);
        if x == 1 || x == n - 1 { continue 'main_loop; }
        for _ in 0..(s - 1) {
            x = utils::mod_exp(x, 2, n);
            if x == n - 1 { continue 'main_loop; }
        }
        return PrimeResult::Composite;
    }
    PrimeResult::ProbablePrime
}

Первый цикл while выносит за скобки степени двойки, преобразуя n-1 в 2^s × d. Затем main_loop выполняет все описанные выше тесты, возводя в квадрат x (повышая степень на 2) и проводя тестирование, пока не достигнет 2^s-1.

А затем обычный цикл для нахождения простых чисел:

fn primes_128bit() -> u128 {
    loop {
        let num = (rng::u64() | 0x8000000000000001u64) as u128;
        if miller_rabin_test(num, 10) == PrimeResult::ProbablePrime {
            return num;
        }
    }
}

➜ time cargo run --release
Prime found: 15333511742700010117
cargo run --release  0.03s user 0.01s system 99% cpu 0.042 total

Скорость такая же, как у теста Ферма, но как насчёт «вероятно простых»? В наихудшем случае ошибка этого теста стремится к 4^-k, но при больших значениях n, в среднем ошибка гораздо меньше, порядка 8^-k. Какова вероятность того, что тест Миллера — Рабина вернёт составное число при k = 10?

➜ python3 -c "print(f'chance of error = {8 ** -10 :.15f}%')"
chance of error = 0.000000000931323%

Мне вполне достаточно. Для контекста: вероятность в точности равна вероятности тридцать раз подряд выбросить на монетке все орлы (2^-30). Однако при реальном применении в криптографии нужно быть аккуратнее при выборе оснований случайной генерации и предполагать наличие враждебных условий.

Итак, мы, наконец, получили способ генерации случайных чисел и тест простоты, который достаточно быстр и эффективен при работе с крупными числами. Единственное, чего нам не хватает — это сами большие числа, так что нам нужно идти глубже...

1024 бита — делаем небольшой крюк?

Теперь уже очевидно, что мы не можем двигаться дальше, чем на 64 бита, пользуясь только встроенными целочисленными типами данных Rust. Нам требуется bigint, то есть реализация арифметики с произвольной точностью, обычно называемой в большинстве языков bigint или bignum. Разумно было бы импортировать крейт bigint для rust и закончить на этом, но я поставил перед собой условие отсутствия внешних зависимостей и буду его соблюдать.

Итак, настало время реализовать BigInt.

Создание BigInt — это не ответ на вопрос «как генерировать большие простые числа?», но никаких больших простых чисел (и составных тоже) не было бы без BigInt, поэтому мы сделаем небольшой крюк и разберёмся, как работает BigInt.

Попытка номер 1 — BigInt как десятичные разряды

Вкратце прочитав статью Википедии об арифметике произвольной точности и проведя небольшое исследование, я выяснил, что эту задачу можно решить несколькими способами. Проще всего хранить все разряды большого числа в массиве, поэтому в качестве первой попытки я выбрал его.

Поначалу всё очень просто. Число представляется в виде списка цифр, поэтому для реализации сложения и умножения я просто скопировал в код обычные расчёты, которые мы делали в средней школе. Но когда добрался до деления, то понял, что всё будет не так просто, и после нескольких неудачных попыток сдался. Если вы хотите бросить себе вызов, то попробуйте придумать, как реализовать в коде алгоритм деления столбиком, который мы выполняли в тетради.

Попытка номер 2 — BigInt в двоичном виде

После неудачи с первой попыткой я подумал: почему бы не хранить числа в двоичном виде? В частности, почему бы не хранить число как список из нулей и единиц? С этого началась моя вторая попытка.

Таким был мой очень простой BigInt, всего лишь массивом булевых значений:

const N: usize = 2048;
struct BigInt {
    bits: [bool; N]
}

В качестве размера используется 2048, а не 1024, потому что, как мы уже видели, для перемножения двух N-битных чисел требуется максимум 2N битов пространства.

Дальше нам нужно реализовать арифметику. Я ещё немного помню университетский курс логики и микропроцессоров, поэтому для создания сложения и вычитания, по сути, реализовал в коде полный сумматор. Вот как выглядело сложение:

fn bigint_add(own: &[bool], other: &[bool]) -> [bool; N] {
    let mut bits = [false; N];
    let mut carry = false;
    for (i, (d1, d2)) in own.iter().zip(other.iter()).enumerate() {
        bits[i] = d1 ^ d2 ^ carry;
        carry = (d1 & d2) | (carry & (d1 ^ d2));
    }
    if carry { panic!("Attempt to add with overflow"); }
    bits
}
impl Add for BigInt {
    type Output = Self;
    fn add(self, other: Self) -> Self {
        Self { bits: bigint_add(&self.bits, &other.bits) }
    }
}
impl AddAssign for BigInt {
    fn add_assign(&mut self, other: Self) {
        self.bits = bigint_add(&self.bits, &other.bits);
    }
}

Примечание: здесь я использую трейты Add и AddAssign из Rust для переопределения операторов + и += моего типа BigInt и делаю то же самое для всех других операторов.

Далее я реализовал операторы сдвига влево (<<) и вправо (>>). Они просто сдвигают весь список битов влево или вправо на указанную величину, отбрасывая все переполнения.

fn bigint_shl(own: &[bool], amount: usize) -> [bool; N] {
    let mut bits = [false; N];
    if amount > N { return bits; }
    let mut i = amount;
    for bit in own.iter().take(N - amount) {
        bits[i] = *bit;
        i += 1;
    }
    bits
}
fn bigint_shr(own: &[bool], amount: usize) -> [bool; N] {
    let mut bits = [false; N];
    if amount > N { return bits; }
    let mut i = 0;
    for bit in own.iter().skip(amount) {
        bits[i] = *bit;
        i += 1;
    }
    bits
}

Перейдём к умножению. Удобство работы с двоичными значениями заключается в том, что умножение становится очень простым. Двоичные разряды могут принимать значения только 0 или 1, поэтому каким бы длинным ни было число, единственными двумя результатами его умножения на бит может быть или 0, или оно само. Это упрощает классический алгоритм умножения до гораздо более простого алгоритма сдвига и сложения, и теперь мы можем использовать только что реализованные сдвиг и сложение:

fn bigint_mul(own: &[bool], other: &[bool]) -> [bool; N] {
    let mut result = BigInt::zero();
    let n1 = BigInt::from(own);
    let mut current;
    for (shift, d2) in other.iter().enumerate() {
        if !(*d2) { continue; }
        for i in (N - shift)..N {
            if own[i] { panic!("Attempt to multiply with overflow"); }
        }
        current = n1 << shift;
        result += current;
    }
    result.bits
}

Теперь мы перейдём к делению. Единственная причина того, что я решил выбрать способ с двоичными числам, заключается в том, что двоичное деление в столбик тоже сильно упрощает задачу. В отличие от деления в столбик с цифрами 0-9, цифры частного могут иметь значения только 1 или 0, то есть в промежуточных делениях используется или делитель, или 0. Это похоже на представленный выше алгоритм сдвига и сложения и может быть реализовано с помощью одних лишь сдвига и вычитания:

fn bigint_div(own_bits: &[bool], other_bits: &[bool]) -> ([bool; N], [bool; N]) {
    let mut quotient = BigInt::zero();
    let mut dividend = BigInt::from(own_bits);
    let mut remainder = BigInt::zero();
    let divisor = BigInt::from(other_bits);
    if divisor == BigInt::zero() {
        panic!("Attempt to divide by zero");
    }
    let mut no_of_bits = N;
    while !dividend.bits[N - 1] {
        dividend <<= 1;
        no_of_bits -= 1;
    }
    for i in 0..no_of_bits {
        remainder <<= 1;
        remainder.bits[0] = dividend.bits[N - 1 - i];
        quotient <<= 1;
        if remainder >= divisor {
            remainder -= divisor;
            quotient.bits[0] = true;
        }
    }
    (quotient.bits, remainder.bits)
}

Немного потестировав код, чтобы убедиться, что арифметика работает, я заменил все  u128 в miller_rabin_test() и mod_exp() на свой BigInt, поменял RNG, чтобы он заполнял 1024 бита, и запустил код. Он не закончил работу через несколько минут, а время было позднее, поэтому я оставил его работать и ушёл спать. Проснувшись на следующий день, я увидел такую картину:

➜ time cargo run --release
Prime found: 1100001101000111101111110011001111000010011111010110100010101011111110100010001100001101110101100100100110011110101100001000110010100011100101001010100110010110001101110111001110100110000001010111100000110111100010010010100110101011101010100101000100111110000001100000000011101111010100100101001111111100010010000110101010101000101010000011011110111101100011000010111011010000110101100111001101011000000001000100011001011100100000011110011011000000101000010001010001010010001101111100110001000011100110111100000010100000011101011100101101110100010111110000110000010111110000110101001110100100110101011101100111000100010110101001101110100111100010000000110000001001011100100100101001110100101100110110001001101110010100011011011110111100010011111011101010100010111010110101101010011110010111100100110111010111100101111111110010100101010111100111010001010101011100001000100101001111110101100101011100100001101111000100000001111100001001001010100101101111010000101001111111010110111111011101010100111011100111100001010101011101
cargo run --release  1959.67s user 0.09s system 99% cpu 32:44.90 total

Если выразить по основанию 10:

137130462909417371581865483489043797725909059024661411704723085022816692663284008207826785132470756353352621332808019668785759110990576815741502628035997147255459016128105305451010585699069674494217365521467940783164171729442866016775055913991624626502191730619275815532321664270492537447637102633611801007453

Это первое найденное мной 1024-битное простое число! Я решил свою задачу!

Но оставалась одна небольшая проблема — счётчик времени исполнения показал, что для нахождения этого простого числа потребовалось около тридцати минут. Хотя, строго говоря, я решил задачу и понял, как это делается, по тридцать минут на одно простое число я бы не назвал «разумным временем», особенно учитывая, что OpenSSL для этого требуется всего 30 мс! Требование «разумного времени» было частью ограничений задачи, потому что я хотел и решить её, и сделать это эффективно.

Попытка номер 3 — BigInt как байты

Двоичная реализация, созданная мной при попытке номер 2, была невероятно ценной. Она не только дала мне первое 1024-битное простое число, но и стало правильно работающей реализацией, относительно которой можно тестировать все дальнейшие изменения. Это сильно помогло в существенном ускорении моих экспериментов и придало мне уверенности в попытках сделать что-то более сложное.

Когда я начал разбираться в причинах медленности двоичного решения, то первым делом обнаружил, что в массиве bool каждый bool занимал в памяти байт, а не один бит, как я считал. В ответе на Stackoverflow приведены причины этого. Это означало, что мой массив bool размера 2048 занимал в памяти не 2048 битов, а 2048 байтов! Это 2 КБ памяти на хранение лишь одного числа. Вероятно, моя двоичная реализация тратила всё своё время на ожидание чтения или записи чисел из ОЗУ из-за промахов кэша L1. В то время я не знал, как это протестировать, но подумал, что всё равно нужно попробовать более эффективно использующую память версию и проверить, улучшит ли она ситуацию.

Значит, вполне логично будет хранить биты как блоки размером в байт, а не как отдельные биты в списке. То есть мы сможем хранить все 2048 битов в массиве из 256 байтов. Как ни удивительно, сложение/вычитание и умножение работали с этим новым форматом без серьёзных изменений в алгоритмах. Вместо сложения бита за битом и использования дополнительного бита для переноса мы могли теперь складывать байт за байтом, используя для переноса дополнительный байт. Для умножения я заменил алгоритм сдвига и сложения снова на «тетрадный», но использовал байты вместо цифр 0-9, и он тоже заработал без изменений. Для деления я добавил несколько строк кода, чтобы обрабатывать список байтовых блоков как единый список битов, а остальное оставил без изменений.

После всех этих изменений я получил своё второе 1024-битное простое число за 4 минуты 43 секунды. Хорошее улучшение, но этого всё равно недостаточно.

Попытка номер 4 — BigInt как блоки по u64

Проведя дополнительные исследования арифметики с произвольной точностью в попытках найти новые способы оптимизации, я наткнулся на интересное открытие. Есть причина того, что мои арифметические алгоритмы в попытке номер 3 работали напрямую с байтами вместо битов. Невольно я реализовал BigInt на основе цифр, схожий с тем, что сделано в попытке номер 1, но с использованием цифр «высокого основания». В попытке номер 1 использовались цифры по основанию 10, от 0 до 9, которые нам привычны, но компьютер может работать с любым основанием, какое мы выберем. Вероятно, вы слышали об основании 16 (или шестнадцатеричном), в котором цифрами являются 0-9 и A-F, или даже об основании 64, цифры которого представляют все алфавитно-цифровые символы. В моём коде из попытки номер 3, по сути, использовалось основание 256, где каждый байт работал, как отдельная «цифра».

Например, вот одно и то же число, выраженное по разным основаниям:

base-10:    3,095,627                          (7 разрядов)
base-16:    2F3C4B                             (6 разрядов)
base-64:    LzxL                               (4 разряда)
base-256:   00101111 00111100 01001011         (3 разряда)

(Версия по основанию 256 — это, по сути, 3095627 в двоичном виде, разбитое на 3 байта.)

Любопытно то, что я читал об этом много раз, но каждый раз думал «это слишком сложная для меня концепция», но после её реализации в моей голове, наконец, что-то «щёлкнуло». Как только я разобрался, для меня стало логично, что нет никакой причины ограничиваться цифрами размером с байт и можно увеличить их размер до максимума.

Так что вот как выглядит мой последний BigInt:

const N: usize = 2048 / 64;
struct BigInt {
    chunks: [u64; N]
}

В нём используется массив из 32 блоков u64 для хранения до 2048 битов. Как обычно, он в два раза больше необходимого нам, чтобы оставалось достаточно места для хранения результата перемножения двух BigInt. Аналогично, каждая отдельная «цифра» может быть не больше u64, потому для хранения результата умножения двух «цифр» потребуется u128. Остальная часть кода по большей части не изменилась с попытки номер 3. Были внесены только несущественные изменения, например, переменная переноса теперь имеет размер u64, а не байт. На этом этапе BigInt использует основание (2⁶⁴-1), или 18446744073709551615, и для представления числа, которое в десятеричной системе занимает 309 цифры, ему нужно всего 16 «цифр»!

Теперь для генерации 1024-битных простых чисел требуется примерно 60-90 секунд, что существенно лучше, чем в случае с двоичными числами, но всё равно недостаточно быстро.

Попытка номер 5 — BigInt как блоки u64, но лучше

Дальше я решил выполнить несколько простых бенчмарков, чтобы попробовать разобраться, что же нас замедляет. Вот результаты:

	двоичные	блоки u64
a + b и a - b	5537,35 нс	123,57 нс
a * b	1292283,14 нс	842,32 нс
a / b и a % b	733446,76 нс	44440,12 нс
a << b и a >> b	276,85 нс	140,88 нс
a < b и a > b	2506,02 нс	58,91 нс

(Все результаты усреднены из 1000 прогонов, измеренных в наносекундах)

Это демонстрирует существенный прогресс, проделанный после двоичных чисел. Остальная часть попытки номер 5 будет списком того, что я сделал, чтобы ещё больше увеличить скорость.

Деление

Больше всего в бенчмарках бросается в глаза деление. Несмотря на то, что всё остальное сильно улучшилось, деление по-прежнему использует тот же алгоритм, что и в двоичных числах, продолжая выполнять деление в столбик по одному биту за раз. Я часто слышал, как люди жалуются на деление, и теперь знаю, почему. Деление — гораздо более сложная задача, чем сложение или умножение.

Многие найденные в поисках более совершенных алгоритмов статьи и источники указывали на книгу Handbook of Applied Cryptography. Я обнаружил её в Internet Archive, создал аккаунт и позаимствовал её на 14 дней. Вот, что я нашёл на странице 598:

С одной стороны, там говорится о «представлении по основанию b», и к этому моменту я уже понял достаточно, чтобы осознать, что это отсылка к тем же «цифрам с высоким основанием», которые я обнаружил ранее. С другой стороны, всё остальное мне было абсолютно непонятно, и книга не упрощала мне жизнь, потому что в ней не было никакого текста с объяснением алгоритма. Я три дня подряд перечитывал эту страницу и сделал множество неудачных попыток, в результате написав работающую реализацию. Понял я то, что алгоритм выполняет деление столбиком чисел по основанию N, используя первые три «цифры» делимого и первые две «цифры» делителя для подбора текущей «цифры» частного в цикле, пока не найдёт правильное значение; но пока я недостаточно освоился с лежащими в основе этого вычислениями, чтобы объяснить их. Ещё я слышал, что очень похожий алгоритм встречается в Искусстве программирования, но я никак не мог быстро свериться с ней, кроме как купив книгу, что я рано или поздно сделаю.

Впрочем, трата времени на это была не напрасной: это экономит примерно 40 тысяч нс (40 мкс) на каждое деление, а в каждом прогоне и в тысяче прогонов Миллера — Рабина довольно много операций деления и деления с остатком, пока не будет найдено простое число. Я внёс ещё одну дополнительную оптимизацию — проверку того, является ли делитель одной «цифрой» (одним блоком u64), после чего напрямую выполняю деление столбиком с использованием u128, чтобы перехватывать любые переполнения. Это позволяет полностью пропускать затратный алгоритм; в то же время это один из тех случаев, которые часто встречаются в тесте Миллера — Рабина.

fn bigint_div(mut dividend: BigInt, mut divisor: BigInt) -> (BigInt, BigInt) {
    if divisor.is_zero() { panic!("Attempt to divide by zero"); }
    if dividend < divisor { return (BigInt::zero(), dividend) }
    // x = dividend
    // y = divisor
    // b = 64 (size of a "digit")
    let mut quotient = BigInt::zero();
    let mut lambda = 0;
    let t = divisor.size();
    if divisor.chunks[t] < u64::MAX / 2 {
        while divisor.chunks[t] << lambda < u64::MAX / 2 {
            lambda += 1;
        }
        divisor <<= lambda;
        dividend <<= lambda;
    }
    let n = dividend.size();
    // if y has only 1 "digit", then do long division directly
    if t == 0 {
        let divisor_digit = divisor.chunks[0] as u128;
        let mut remainder = 0;
        let mut current = 0;
        for (i, chunk) in dividend.chunks.iter().enumerate().rev().skip(N - n - 1) {
            current = (remainder << 64) + *chunk as u128;
            quotient.chunks[i] = (current / divisor_digit) as u64;
            remainder = current % divisor_digit;
        }
        return (quotient, BigInt::from(remainder >> lambda));
    }
    // step 2, align and then subtract y from x until x >= aligned
    let mut aligned = divisor.clone();
    for _ in 0..(n - t) {
        aligned <<= 64;
    }
    while dividend >= aligned {
        quotient.chunks[n - t] += 1;
        dividend -= aligned;
    }
    let one = BigInt::from(1);
    let mut x_3digit;
    let mut y_2digit;
    let mut q_u128;
    let mut q_digit;
    // step 3
    for i in ((t + 1)..=n).rev() {
        q_digit = BigInt::zero();
        // step 3.1
        if dividend.chunks[i] == divisor.chunks[t] {
            q_digit.chunks[0] = u64::MAX - 1;
        } else {
            q_u128 = (dividend.chunks[i] as u128) << 64;
            q_u128 += dividend.chunks[i - 1] as u128;
            q_digit = BigInt::from(q_u128 / divisor.chunks[t] as u128);
        }
        // precalc 3digit x and 2digit y
        x_3digit = BigInt::zero();
        x_3digit.chunks[2] = dividend.chunks[i];
        x_3digit.chunks[1] = dividend.chunks[i - 1];
        x_3digit.chunks[0] = dividend.chunks[i - 2];
        y_2digit = BigInt::zero();
        y_2digit.chunks[1] = divisor.chunks[t];
        y_2digit.chunks[0] = divisor.chunks[t - 1];
        // step 3.2
        while q_digit * y_2digit > x_3digit {
            q_digit -= one;
        }
        // move quotient "digit" from temp bigint to its place in quotient
        quotient.chunks[i - t - 1] = q_digit.chunks[0];
        // precalc shifted y
        let mut y_shifted = divisor.clone();
        for _ in 0..(i - t - 1) {
            y_shifted <<= 64;
        }
        // step 3.3 and 3.4
        if dividend >= q_digit  y_shifted {
            dividend -= q_digit  y_shifted;
        } else {
            dividend += y_shifted;
            dividend -= q_digit * y_shifted;
            quotient.chunks[i - t - 1] -= 1;
        }
    }
    // rewind shifts by lambda to get actual remainder
    dividend >>= lambda;
    (quotient, dividend)
}

Умножение

Затем я начал исследовать умножение — вторую по медленности операцию в бенчмарках. По сути, я реализовал тот же алгоритм, что и в книге, так что никакого выигрыша не получил. Хитро перераспределив вычисления в цикле, чтобы избавиться от лишнего BigInt, который использовался для хранения промежуточных результатов, я получил ускорение вдвое.

Так как я уже добавил функцию для вычисления размера (количества занятых блоков) в делении, то воспользовался той же функцией, чтобы выполнять циклы для ненулевых блоков, а также добавил проверки для пропуска итераций цикла, когда один из блоков равен нулю. Хотя это повышает сложность и добавляет ветвлений внутри циклов, так мы всё равно увеличиваем производительность, потому что в большинстве случаев BigInt должен хранить 1024 битов или меньше, то есть быть наполовину пустым. Это дало ускорение ещё в два раза.

Я мог бы реализовать алгоритм Карацубы или даже быстрые преобразования Фурье (FFT), что теоретически могло ещё больше улучшить производительность, но их реализация для моих самодельных BigInt была слишком сложной, а моё умножение и так было уже достаточно быстрым, чтобы не идти по этому пути.

fn bigint_mul(own: BigInt, other: BigInt) -> BigInt {
    let mut result = BigInt::zero();
    let mut intermediate;
    let mut carry;
    let t = own.size();
    let n = other.size();
    if t + n + 1 >= N { panic!("Attempt to multiply with overflow"); }
    for (j, chunk2) in other.chunks.iter().take(n + 1).enumerate() {
        if *chunk2 == 0 { continue; }
        carry = 0;
        for (i, chunk1) in own.chunks.iter().take(t + 1).enumerate() {
            if *chunk1 == 0 && carry == 0 { continue; }
            intermediate = ((chunk1 as u128)  (*chunk2 as u128)) + carry;
            intermediate += result.chunks[i + j] as u128;
            result.chunks[i + j] = intermediate as u64;
            carry = intermediate >> 64;
        }
        result.chunks[t + j + 1] += carry as u64;
    }
    result
}

Сложение и вычитание

Они и так уже были очень быстрыми: меня удивило, что даже моя собственная реализация выполняется примерно за то же время, за которое Rust складывает два нативных числа u128. Я попробовал несколько трюков, но (как и ожидалось) я не так умён, как компилятор и вся его внутренняя магия.

fn bigint_add(own: BigInt, other: BigInt) -> BigInt {
    let mut sum;
    let mut carry = 0;
    let mut sum_overflow;
    let mut carry_overflow;
    let mut result = BigInt::zero();
    let own_iter = own.chunks.iter();
    let other_iter = other.chunks.iter();
    for (i, (chunk1, chunk2)) in own_iter.zip(other_iter).enumerate() {
        (sum, sum_overflow) = chunk1.overflowing_add(*chunk2);
        (sum, carry_overflow) = sum.overflowing_add(carry);
        result.chunks[i] = sum;
        carry = sum_overflow as u64 + carry_overflow as u64;
    }
    if carry != 0 { panic!("Attempt to add with overflow"); }
    result
}

Миллер — Рабин

Я нашёл кучу оптимизаций, которые можно было внести в мою реализацию Миллера — Рабина. Вот, как этот тест выглядел изначально:

 1fn miller_rabin_test(n: BigInt, k: usize) -> PrimeResult {
 2    let zero = BigInt::zero();
 3    let one = BigInt::from(1);
 4    let two = BigInt::from(2);
 5
 6    let mut s = zero;
 7    let mut d = n - one;
 8    while d % two == zero {
 9        d = d / two;
10        s += one;
11    }
12
13    'main_loop: for _ in 0..k {
14        let base = <<rng omitted>>;
15        let mut x = mod_exp(base, d, n);
16        if x == one || x == n - one { continue 'main_loop; }
17
18        while s > zero {
19            x = mod_exp(x, two, n);
20            if x == n - one { continue 'main_loop; }
21            s -= one;
22        }
23        return PrimeResult::Composite;
24    }
25
26    PrimeResult::ProbablePrime
27}

А вот некоторые из внесённых мной серьёзных оптимизаций:

• Вызов mod_exp() в строке 19 больше не требуется, поскольку BigInt достаточно памяти, чтобы выполнять x = (x * x) % n напрямую.

• В строке 15 mod_exp() заменена упрощённой inline-версией, что позволило избавиться от лишней траты ресурсов, связанной с вызовом функций.

• Представление «2» в BigInt в строке 4 было создано только для выполнения проверки на чётность в строке 8. Я написал простую функцию num.is_even(), которой достаточно проверить последний бит, поэтому удалил кучу дополнительных затратных операций деления и лишнее распределение BigInt.

• Аналогично, деление в строке 9 можно заменить операции сдвига d >>= 1. В данном случае замена ещё одной кучи затратных операций деления сдвигами даёт очень большой выигрыш по сравнению с нативным кодом, в котором такое изменение обычно не оправдывает себя.

• В коде присутствует много += one и -= one (где one — это представление «1» в BigInt). Я добавил специальные num.increase() и num.decrease() , которые практически во всех случаях просто выполняют сложение/вычитание последней «цифры» u64, а полное сложение/вычитание BigInt производят только тогда, когда последняя «цифра» равна или 0, или u64::max, то есть в тех редких случаях, когда BigInt действительно требуется для обработки переполнения от прибавления или вычитания единицы.

Эти и другие изменения по отдельности дают всего микросекунды или даже наносекунды выигрыша, но когда они многократно выполняются внутри цикла в тысячах тестов Миллера — Рабина, всё это даёт приличное снижение времени исполнения. По крайней мере, так я думал, пока не провёл бенчмаркинг. Оказалось, что is_even() и d >>=1 с лёгкостью обходят все остальные оптимизации по выгодности, давая выигрыш по 70 тысяч наносекунд каждая! Вот готовый усовершенствованный тест Миллера — Рабина:

fn miller_rabin_test(n: BigInt, k: usize) -> PrimeResult {
    let zero = BigInt::zero();
    let one = BigInt::from(1);
    let n_minus_1 = n.decrease();
    let mut d = n_minus_1;
    let mut s = zero;
    while d.is_even() {
        d >>= 1;
        s = s.increase();
    }
    let mut bytes = [0; (1024 / 16)];
    let mut x;
    let mut base;
    'main_loop: for _ in 0..k {
        rng::insert_random_bytes(&mut bytes).unwrap();
        base = BigInt::from(bytes.as_slice());
        x = one;
        while !d.is_zero() {
            if !d.is_even() { x = (x  base) % n; }
            d = d >> 1;
            base = (base  base) % n;
        }
        if x == one || x == n_minus_1 { continue 'main_loop; }
        while !s.is_zero() {
            x = (x * x) % n;
            if x == n_minus_1 { continue 'main_loop; }
            s = s.decrease();
        }
        return PrimeResult::Composite;
    }
    PrimeResult::ProbablePrime
}

Логика тестирования простоты

Вдохновившись изменениями на этапе 2 (64-битные простые числа), я модифицировал логику проверки на простоту, добавив в начало проверку перебором делителей с использованием заранее вычисленного списка первых 5000 малых простых чисел. По большей части это было бы неэффективно при работе с BigInt, так как перебор делителей требует множества операций деления, то есть чрезвычайно много времени. Трюк заключается в том, что все первые 5000 малых простых чисел умещаются в одну «цифру» (один блок u64). То есть все операции деления внутри перебора делителей будут относиться к добавленному мной особому случаю, в котором можно выполнить всё деление при помощи деления столбиком с использованием u128, отказавшись от затратного алгоритма деления BigInt. Ту же функцию перебора делителей также можно использовать для генерации изначального списка первых 5000 малых простых чисел. В конечном итоге, мы нашли применение для оптимизации перебора делителей этапа 2!

Ещё одно изменение в логике заключается в том, что вместо чтения /dev/urandom для каждой итерации цикла и генерации нового случайного числа мы просто прибавляем 2 к первому случайному числу и получаем нового кандидата. Так как последний бит изменяется так, чтобы быть равным 1, мы знаем, что это нечётное число, то есть прибавляя 2, мы получим следующее нечётное число. Это можно оптимизировать ещё больше, добавив специальную функцию num.increase_by_2(), которая подобно num.increase() будет выполнять полное сложение BigInt только для случаев переполнения, а во всех остальных случаях просто выполняет сложение u64.

Кроме того, это одна из тех задач, которые можно назвать «чрезвычайно параллельными», потому что им не нужна общая память и синхронизация между потоками. Так что можно попросить искать простые числа не один поток CPU, а все шестнадцать, а победит самый быстрый!

Вот те же бенчмарки после всех оптимизаций:

	двоичные	блоки u64	блоки u64 с оптимизацией
a + b и a - b	5537,35 нс	123,57 нс	123,62 нс
a * b	1292283,14 нс	842,32 нс	295,04 нс
a / b и a % b	733446,76 нс	44440,12 нс	831,77 нс
a << b и a >> b	276,85 нс	140,88 нс	126,04 нс
a < b и a > b	2506,02 нс	58,91 нс	58,50 нс
a / 2 (или a << 1)	2638289,48 нс	75121,58 нс	60,89 нс
a % 2 == 0 (или a.is_even())	2447553,14 нс	78400,87 нс	21,65 нс
a - 1 (или a.decrease())	6179,48 нс	103.15ns	67,54 нс

(Все результаты усреднены из 1000 прогонов, измеренных в наносекундах)

Гораздо более быстрые 1024 битов!

Мы наконец добрались до конца этой очень длинной статьи. Давайте соединим всё в функцию для генерации 1024-битных простых чисел:

fn primes_1024bit() -> BigInt {
    const P: usize = 1000;
    let primes = utils::generate_small_primes::<P>();
    let zero = BigInt::zero();
    let mut small_prime = BigInt::zero();
    let mut bytes = [0u8; 1024 / 8];
    rng::insert_random_bytes(&mut bytes).expect("Cannot access /dev/urandom");
    let mut num = BigInt::from(bytes.as_slice())
    num.chunks[(1024 / 64) - 1] |= 0x8000000000000000u64;
    num.chunks[0] |= 1;
    'prime_loop: loop {
        num = num.increase_by_2();
        for i in 0..P {
            small_prime.chunks[0] = primes[i];
            if num % small_prime == zero {
                continue 'prime_loop;
            }
        }
        if miller_rabin_test(num, 10) == PrimeResult::ProbablePrime {
            return num;
        }
    }
}

Вызовем эту функцию в параллельных потоках:

fn run() {
    let (tx, rx) = std::sync::mpsc::channel();
    for _ in 0..16 {
        let thread_tx = tx.clone();
        std::thread::spawn(move || {
            thread_tx.send(primes_1024bit()).unwrap();
        });
    }
    let prime = rx.recv().unwrap();
    prime.print_decimal();
}

И вот результаты:

➜ time cargo run --release
133639768604208228777408136159783586754136713880762782100572086187859339703910900715674773943439684405153138260262492990850200027881950953138966616704637409705491165541761840874200485820151419486204300434469857557841532664407934654743999891926036532834796558113864177048787433702650711105375897281079281724197
cargo run --release  0.58s user 0.01s system 690% cpu 0.086 total

Я наконец-то могу сказать: задача решена!

Для нахождения 1024-битного простого числа в среднем требуется 40 мс. Из-за случайности отдельные вызовы prime_1024bit() могут выполняться от примерно 8 мс до максимум примерно 800 мс, но параллельное исполнение и выбор самого быстрого результата сглаживает эти выбросы. Вот среднее время исполнения при 100 прогонах:

➜ perf stat -r100 ./target/release/primes
--- outputs omitted ---
Performance counter stats for './target/release/primes' (100 runs):
--- other stats omitted ---

0.04109 +- 0.00307 seconds time elapsed  ( +-  7.48% )

В конечном итоге я очень доволен тем, что бросил себе такой вызов. Он вынудил меня выйти из своей зоны комфорта и заставил научиться новому и в программировании, и в том, как изучать новые темы и проводить исследования. В моём todo-списке по-прежнему куча пунктов и идей для улучшений и того, что можно было бы сделать лучше, но мне надо подкопить силы для нового проекта.

Полный код и репозиторий можно найти на github. Можно и не говорить о том, что весь этот код вряд ли криптографически безопасен, но это и не было моей целью.