repunit-hunt

Русский · English

repunit-hunt

Гибридный CPU/GPU искатель простых обобщённых репьюнитов

$$R_k(b) = \frac{b^k - 1}{b - 1} = 1 + b + b^2 + \dots + b^{k-1}$$

Конвейер: генерация k → сито по малым делителям (CPU) → trial factoring (GPU, Montgomery 64/96/128 бит) → P-1 (libecm) → PRP (GMP или GWNUM с проверкой Gerbicz–Li).

В репозитории две части: искатель простых репьюнитов (Rust + CUDA + C) и статья (paper/) — статистический анализ найденных последовательностей.

Статья

Каталог paper/ содержит работу «Схема наблюдения в статистике обобщённых репьюнитных простых: точный условный критерий константы Ленстры — Померанса — Вагстаффа».

Текущая версия — 3:

по-русски: paper_ru.pdf, оформление по ГОСТ Р 7.0.7-2021;
по-английски: paper_en.pdf.

Что изменилось по сравнению с версией 1, описано в paper/CHANGES_v3.md.

Коротко: гипотеза ЛПВ предсказывает плотность простых репьюнитов с универсальной константой e^γ ≈ 1,781. Публикуемые с 1993 года оценки лежали примерно на 10 % выше — и это превышение оказалось артефактом схемы наблюдения. Поиск ведётся до незаписанного фронта со случайным числом находок. Положение последней находки достаточно для фронта, и условно на нём N_b − 1 ~ Poisson(λ_b D_b) точно (задача Неймана — Скотта: одно событие на основание). Верная объединённая оценка — (M−B)/S, а не (M−1)/S:

κ̂ = 1,838   точный 95 % ДИ [1,594; 2,109]   p = 0,672 против e^γ

Сама эвристика ЛПВ при конечных n предсказывает интенсивность (κ/ln b)(1 + c_b/t). Константа c_b откалибрована по фактической делимости R_b(p) (в среднем 0,67). С ней при n > 10 оценка равна 1,69–1,71 при любом усечении и ЛПВ не отвергается (p = 0,62–0,73) — § 9 статьи.

Числа и рисунки версии 3 воспроизводят paper/analysis_v2.py, make_results_v3.sh (lpw_second_order.py, cb_empirical.py, cb_theory.py, sim_second_order.py, frontiers.py, stopping_rules.py) при фиксированных зёрнах; каждое десятичное число статьи сверяется с их выводом скриптами check_numbers_v3.py и check_all_numbers_v3.py. Версия 1 (paper_ru_v1.pdf, paper_en_v1.pdf, analysis.py, validate_v1.sh) сохранена как архив; подробности — в paper/README.md.

Структура репозитория

Путь	Что это
`src/`, `native/`, `benches/`, `build.rs`	искатель: сито, trial factoring, P−1, PRP
`config/default.toml`	параметры конвейера с обоснованием каждого значения
`paper/paper_ru.tex`, `paper/paper_ru.pdf`	статья, версия 3, русская
`paper/paper_en.tex`, `paper/paper_en.pdf`	статья, версия 3, английская
`paper/analysis_v2.py`	анализ версии 3 → `results_v2.txt`, `figs_v2/`
`paper/lpw_second_order.py`, `frontiers.py`, `stopping_rules.py`	§ 9, § 4.5, таблица 5 версии 3 → `results_v3.txt`
`paper/paper_ru_v1.tex`, `paper/paper_en_v1.tex`	статья, версия 1 (архив)
`paper/analysis.py`	анализ версии 1 → `results.txt`, `figs/`
`paper/data/`	b-файлы OEIS (лицензия OEIS, см. ниже)
`paper/verify/`	независимый сплошной пересчёт до n < 10⁴ — на нём основано приложение статьи
`paper/verify_32k/`	более поздний прогон, доведённый до k_max = 32000
`paper/METHOD.md`, `paper/FINDINGS.md`	рабочие заметки по методу и находкам (английские версии — `*.en.md`)

Про два каталога проверки. Канонический — verify/: он однороден (k_max = 10000 у всех двадцати оснований) и даёт ровно те числа, что стоят в приложении: 153 подтверждённых члена, 128 из них на основном наборе b ≤ 20. Именно его читают compare_verify.py и check_numbers.py.

verify_32k/ — более поздний прогон, доведённый до k_max = 32000, но неоднородно: у разных оснований граница 10000, 20000 или 32000. Он шире (169 членов) и потому полезен как дополнительное свидетельство, но не подтверждает утверждение вида «проверены все индексы ниже X», и числам статьи не соответствует.

1. Математическая база

1.1 Почему только простые k

Если $d \mid k$, то $R_d(b) \mid R_k(b)$. Значит при составном k репьюнит заведомо составной, и перебирать нужно только простые k.

1.2 Форма делителей

Пусть k простое и q — простой делитель $R_k(b)$. Тогда $b^k \equiv 1 \pmod q$, то есть $\mathrm{ord}_q(b) \mid k$, а k простое, поэтому порядок равен 1 или k.

$\mathrm{ord}_q(b) = k$. По малой теореме Ферма $k \mid q-1$; q нечётно, значит $2k \mid q-1$ и $$q = 2mk + 1.$$ Именно эти q перебирает GPU-ядро — вместо всех чисел до предела проверяется лишь каждое $2k$-е.
$\mathrm{ord}_q(b) = 1$ — краевой случай. Тогда $q \mid b-1$, и проверка powmod(b, k, q) == 1 срабатывает всегда, то есть даёт ложный делитель. На самом деле $$R_k(b) = 1 + b + \dots + b^{k-1} \equiv k \pmod q,$$ поэтому $q \mid R_k(b) \iff q \mid k$, а для простого k это значит ровно $k = q$.
$q \mid b$. Тогда $R_k(b) \equiv 1 \pmod q$ — q не делитель никогда.
Делитель обязан быть собственным: $q < R_k(b)$. Иначе простой репьюнит «отсеивается» сам собой: $R_2(10) = 11$, и кандидат $q = 2\cdot1\cdot5+1 = 11$ честно делит его нацело. Проверка нужна только при крошечных k — уже при $(k-1)\log_2 b > 128$ репьюнит больше любого 128-битного q.

Все три ветки реализованы явно и в GPU-ядре (native/cuda/tf_kernel.cu), и в CPU-сите (src/sieve/small.rs), и покрыты тестами.

1.3 Тест на простоту

У $R_k(b)$ нет формы $N \pm 1$ с известной факторизацией, поэтому детерминированное доказательство (N−1/N+1, ECPP) вне рамок этого инструмента. Мы делаем сильный тест Ферма по базе 2 плюс раунды Miller–Rabin — то есть получаем PRP-кандидатов, как PFGW/LLR. Базы берутся только из диапазона $[2, N-2]$: при $a \equiv 0, \pm 1 \pmod N$ раунд не несёт информации, а $a = N$ (например $N = 3$, $a = 3$) объявил бы простое составным — так терялись бы малые репьюниты $R_2(2) = 3$ и $R_2(10) = 11$. Финальное доказательство — отдельным инструментом (Primo/ECPP).

1.4 Счёт по модулю b^k − 1

R_k(b) делит b^k − 1, поэтому

$$x \equiv 3^{E} \pmod{b^k-1} \quad\Longrightarrow\quad x \bmod R_k = 3^{E} \bmod R_k.$$

Это позволяет считать по модулю, который на log₂(b−1) бит длиннее, зато имеет специальную форму 1·b^k + (−1). Для неё GWNUM обходится без редукции Барретта, тогда как произвольный модуль требует gwsetup_general_mod. Разница измерена (b = 10):

k	бит	general mod	b^k − 1	ускорение
5 003	16 617	0.008 мс/итер	0.003	2.81×
10 007	33 240	0.015	0.005	2.73×
20 011	66 472	0.042	0.015	2.79×

На боевом числе R₄₉₀₈₁ (163 041 бит) тест ускорился с 16.6 до 4.3 секунды. Если форма не поддержана (слишком большое основание или экзотическое k), код откатывается на gwsetup_general_mod: корректность важнее скорости.

Умножение на базу — только через GWMUL_MULBYCONST. Отдельный вызов gwsmallmul молча портит результат на больших числах: сверка с GMP на 12 итерациях показала совпадение при 66 476 бит, но расхождение при 163 044 и 332 203 бит (roundoff подскакивал ровно до 0.5), при том что чистые квадраты на тех же размерах считались верно. Штатный механизм — константа задаётся через gwsetmulbyconst и применяется внутри самого возведения в квадрат — даёт совпадение с GMP вплоть до 664 396 бит и экономит одну операцию.

1.5 Контроль ошибок в GWNUM-пути

GWNUM считает Fermat-PRP по базе 3 на IBDWT — это в разы быстрее GMP, но FFT работает с плавающей точкой, поэтому нужен контроль. Схема такая:

ошибка округления (gw_get_maxerr) проверяется каждые 128 итераций; при превышении 0.40 счёт перезапускается с большим FFT (до 4 попыток);
каждый положительный вердикт пересчитывается независимым GMP-путём. Статистически это бесплатно — PRP-кандидатов единицы на миллионы проверок, — зато ложное «PRP» из-за сбоя FFT исключено полностью;
обратный случай — ложное «составное», то есть потерянная находка — ловится только повторным счётом, и для этого есть два механизма: double_check_ratio пересчитывает долю «составных» вторым бэкендом прямо в ходе поиска, а --verify перепроверяет уже готовый журнал (раздел 4.1).

Проверки Гербица–Ли здесь нет, и это осознанно. Она применима, когда показатель даёт длинную цепочку квадратов; сам автор GWNUM пишет это прямо (Prime95, commonb.c): «We can do Gerbicz error checking if b=2 and there are a long string of squarings». У нас E = N−1 — произвольная битовая строка, и цепочка имеет вид x ← x²·3^бит. Инвариант Гербица для неё принимает вид d_{j+1} = d_j^(2^L)·3^(S_j), где S_j — сумма L-битных кусков E, и его проверка стоит ~2L операций на каждые L итераций, то есть больше 100% накладных расходов. Дешёвого варианта для произвольного показателя не существует.

2. Архитектура

src/
├── main.rs         CLI, пул потоков, NUMA-pinning
├── config.rs       TOML-конфиг + валидация
├── pipeline.rs     конвейер и его каналы
├── tuner.rs        адаптивная глубина TF и решение по P-1
├── affinity.rs     CPU affinity, NUMA, huge pages
├── worklog.rs      JSONL-журнал, resume после перезапуска
├── report.rs       итоговый JSON
├── sieve/
│   ├── kbase.rs    сегментное сито Эратосфена (генератор k)
│   └── small.rs    сито по малым делителям q (перебор «со стороны q»)
└── ffi/            безопасные RAII-обёртки над нативным слоем
native/
├── tests/          автономные проверки GPU-арифметики и GWNUM-пути
├── include/        общий ABI (rh_common.h) + заголовки слоёв
├── cuda/
│   ├── tf_kernel.cu  ядра трёх ширин (Mont64/96/128)
│   ├── tf_host.cu    стримы, CUDA Graphs, persistent-буферы
│   └── rh_mont.cuh   Montgomery-арифметика (шаблоны)
└── prp/
    ├── prp_dispatch.c  выбор бэкенда по размеру числа
    ├── prp_gmp.c       GMP-бэкенд (Miller–Rabin)
    ├── prp_gwnum.c     GWNUM/IBDWT + Gerbicz–Li
    ├── pm1_ecm.c       P-1 с сидом 2^(2k)
    ├── rh_arena.c      арена переиспользуемых mpz_t
    └── rh_alloc.c      thread-local пул для аллокаций GMP

Три решения, определяющие производительность:

Ноль аллокаций в горячем цикле. mpz_init вызывается один раз на поток (арена), временные буферы GMP идут через thread-local bump-аллокатор на huge pages.
Батчи по k на GPU. Один launch считает tf_k_batch показателей: при больших k диапазон m слишком узок, чтобы занять устройство одним k. Раскладка индексов такова, что весь варп работает с одним k — дивергенции в цикле powmod нет.
Проверка каждого найденного делителя на CPU. Делитель с GPU не принимается на веру: он перепроверяется через GMP. Непрошедшая проверка логируется как ошибка — это сигнал о баге в ядре или нестабильности карты.

3. Сборка

Нативный слой рассчитан на Linux (mmap, clock_gettime, sched_setaffinity); на других платформах потребуется правка rh_alloc.c и affinity.rs.

Обязательно: Rust ≥ 1.75, C-компилятор, GMP (dev-пакет). Опционально: CUDA Toolkit ≥ 11.0, GWNUM (Prime95 SDK), libecm (GMP-ECM).

# Полная сборка (CUDA + GMP)
cargo build --release

# Только CPU (без GPU-стадии)
cargo build --release --no-default-features

# С GWNUM и P-1
GWNUM_DIR=/opt/gwnum ECM_DIR=/usr/local cargo build --release

Переменные окружения сборки:

Переменная	Смысл
`CUDA_PATH` / `CUDA_HOME`	корень CUDA Toolkit (иначе ищется `nvcc` в PATH)
`RH_CUDA_ARCHS`	список SM через запятую, по умолчанию `70,75,80,86,89,90` + PTX
`GMP_DIR`	префикс своей сборки GMP (например, с `--enable-fat`)
`RH_GMP_STATIC=1`	статическая линковка GMP
`GWNUM_DIR`	корень GWNUM; без него PRP идёт через GMP
`ECM_DIR` / `RH_ECM_SYSTEM=1`	libecm; без них стадия P-1 отключается
`RH_MAXREG`	`-maxrregcount` для nvcc (подбирается по occupancy)
`RH_PORTABLE=1`	без `-march=native` (для дистрибутивных пакетов)

Фичи cargo: cuda, gwnum, pm1 (все включены по умолчанию), numa. Отсутствие библиотеки не ломает сборку — соответствующая стадия просто выключается, а --devices покажет, какие бэкенды доступны.

3.1 Развёрнутое окружение (WSL2 на этой машине)

Windows-хост проект не собирает: нативный слой POSIX-only. Рабочее окружение подготовлено в WSL2 Ubuntu 26.04:

Что	Где / версия
Rust	1.98.0, системно в `/opt/rust` (PATH задаёт `/etc/profile.d/rust.sh`)
gcc / ar / pkg-config	из `build-essential`
lld	обязателен — его требует `[target.x86_64-unknown-linux-gnu]` в `.cargo/config.toml`
GMP	`libgmp-dev` 6.3.0
libecm	`libecm-dev` 7.0.6 — стадия P-1 доступна
CUDA	`nvidia-cuda-toolkit` 12.4.131, `nvcc` в `/usr/bin`
GPU	GeForce GTX 1650, CC 7.5 (sm_75), 14 SM

.cargo/config.toml задаёт для этой машины RH_CUDA_ARCHS=75 и RH_ECM_SYSTEM=1. Секция [env] не перекрывает уже заданные переменные, поэтому для сборки под другое железо достаточно экспортировать свой список: RH_CUDA_ARCHS=80,90 cargo build --release.

wsl -d Ubuntu
cd /mnt/c/Users/<вы>/Downloads/repunit-hunt

# target выносим в ext4: на /mnt/c (drvfs) сборка в разы медленнее
export CARGO_TARGET_DIR=~/rh-target

cargo build --release --no-default-features --features "cuda pm1"
cargo test  --release --no-default-features
$CARGO_TARGET_DIR/release/repunit-hunt --devices

GWNUM собран и подключён (Prime95 SDK 30.19, /opt/gwnum-src/extracted). В дистрибутивах его нет, поэтому порядок такой:

curl -O https://www.mersenne.org/download/software/v30/30.19/p95v3019b20.source.zip
unzip p95v3019b20.source.zip -d /opt/gwnum-src/extracted
cd /opt/gwnum-src/extracted/gwnum && make -f make64     # обязательно!

Последний шаг критичен: каталог linux64/ в архиве содержит только предсобранные ассемблерные FFT-модули, без C-части (gwinit2, allocgiant, gwtogiant). Полная библиотека появляется в gwnum/gwnum.a только после make, и build.rs ищет именно её. Ассемблер GWNUM собран без -fPIC, поэтому при подключённом GWNUM бинарь линкуется как non-PIE — build.rs добавляет -no-pie автоматически.

4. Запуск

# Информация об устройствах и бэкендах
./target/release/repunit-hunt --devices

# Поиск с конфигом
./target/release/repunit-hunt --config config/default.toml

# Быстрый прогон без GPU
./target/release/repunit-hunt --base 10 --kmin 3 --kmax 20000 --no-gpu

Флаг	Смысл
`--config <path>`	TOML-конфиг (флаги ниже имеют приоритет)
`--base <b>`	основание репьюнита
`--kmin` / `--kmax`	диапазон показателя, полуинтервал `[kmin, kmax)`
`--threads <n>`	рабочих потоков, `0` = по числу ядер
`--no-gpu`, `--no-pm1`	отключить соответствующую стадию
`--devices`	показать GPU и доступные бэкенды, выйти
`--double-check <доля>`	пересчитывать долю «составных» вторым бэкендом на лету
`--verify`	перепроверить журнал и выйти (см. 4.1)
`--verify-ratio <доля>`	какую долю «составных» пересчитывать при `--verify` (по умолчанию 0.02)

Логи — через env_logger: RUST_LOG=debug ./target/release/repunit-hunt ...

4.1 Перепроверка журнала (double-check)

# выборочно (2% составных) — быстро
./target/release/repunit-hunt --config config/default.toml --verify

# полностью, включая пересчёт каждого составного
./target/release/repunit-hunt --config config/default.toml --verify --verify-ratio 1.0

Что проверяется по каждой записи worklog.jsonl:

Запись	Проверка
`factored`	делит ли записанный q число R_k(b) и собственный ли он (q < N) — дёшево, проверяются все
`PRP`	находка пересчитывается точной арифметикой GMP — всегда
`composite`	вердикт пересчитывается другим бэкендом: только так ловится потерянная находка

Код возврата 1, если найдено хоть одно расхождение, — режим годится для cron и CI. Выборка детерминирована (хеш от base и k), поэтому один и тот же --verify-ratio всегда проверяет одни и те же числа.

Возобновление работы

Каждый закрытый показатель пишется в worklog.jsonl (append-only JSONL, PRP-находки флашатся немедленно). При следующем запуске такие k пропускаются, так что прерванный поиск продолжается с места остановки. Итог дублируется в results.json (атомарная замена через временный файл).

5. Настройка (config/default.toml)

Ключевые параметры:

bitsieve_q_limit — граница CPU-сита. Дёшево снимает большинство кандидатов; поднимать имеет смысл, пока сито строится за секунды.
tf_q_min / tf_q_hard_max_bits — окно trial factoring. Нижняя граница автоматически поднимается до bitsieve_q_limit: ниже уже отработало сито.
tf_adaptive — глубину TF выбирает тюнер: расширять диапазон q выгодно, пока ожидаемая экономия на PRP (≈ t_prp / ln q) превышает стоимость перебора очередной декады.
tf_k_batch — сколько k уходит на GPU одним батчем (≤ 8192, это MAX_K в tf_host.cu).
mr_rounds — дополнительные базы Miller–Rabin поверх базы 2.
prp_backend — auto переключается на GWNUM при bits ≥ gwnum_threshold_bits.
pin_threads, gmp_pool_mb — NUMA-локальность и размер пула на поток.

6. Проверка результата

results.json содержит PRP, а не доказанные простые. Порядок действий после находки:

Перепроверить независимым инструментом: pfgw64 -tc -q"(10^k-1)/9" или LLR.
Доказать простоту: Primo (ECPP) для чисел до ~50 000 знаков.
Отправить в Prime Pages / профильные проекты.

Ложные срабатывания GPU логируются с уровнем error и суммируются в конце прогона — ненулевой счётчик означает, что результатам trial factoring доверять нельзя (баг ядра, разгон, деградация памяти карты).

6.1 Производительность GPU-стадии

Всё измерено на GTX 1650 (14 SM, Turing) под WSL2, native/tests/bench_gpu.cu, батч 256 показателей.

Накладные расходы запуска. Под WSL GPU паравиртуализован, каждый вызов драйвера идёт через границу VM, поэтому запуск через CUDA Graph (один сабмит вместо пяти вызовов) заметно выгоднее:

вариант	мс/launch
обычные async-вызовы, без графа	51.5
граф, пересоздаваемый каждый раз	37.5
граф, построенный один раз	34.4

Исходный код захватывал граф через stream capture, и его сигнатура включала m_start — а конвейер увеличивает m на каждом запуске, так что граф пересоздавался всегда. Теперь граф строится вручную, а между запусками обновляются только параметры узла ядра (cudaGraphExecKernelNodeSetParams). Проверяется это не по времени (разброс частот под WSL до 20%), а прямым счётчиком: rh_gpu_graph_builds() показывает 1 пересборку на 208 запусков.

Фильтр малых простых оказался главным тормозом. У Turing нет аппаратного целочисленного деления, и 54 операции q % p на кандидата стоили дороже, чем powmod, который они экономят. Деление заменено умножением на обратный по модулю 2⁶⁴:

$$p \mid q \iff (q \cdot p^{-1} \bmod 2^{64}) \le \lfloor (2^{64}-1)/p \rfloor$$

Константы генерирует native/tests/gen_small_primes.py (там же самопроверка).

число простых в фильтре	54	32	16	12	8
деление, мс/launch (k~10⁵)	36.1	21.3	—	—	13.1
умножение, мс/launch (k~10⁵)	12.5	10.9	9.9	9.4	9.0
умножение, мс/launch (k~10⁶)	13.8	12.1	11.0	10.5	10.2

Кривая плоская в диапазоне 8..16, по умолчанию берётся 12. Суммарно GPU-стадия ускорена с 36.1 до 10.5 мс на запуск — в 3.4 раза.

Буфер попаданий. Исходные 16 записей на запуск оказались малы: в начале диапазона m (где q лишь немного больше границы CPU-сита) плотность делителей высокая, и батч из сотен показателей даёт десятки попаданий. На длинном прогоне это дало 8 переполнений — а каждый потерянный делитель означает, что кандидат вместо мгновенного отсева уходит на полноценный PRP-тест в сотню тысяч бит. Ёмкость поднята до 256 записей (6 КиБ), а поле lost в rh_tf_result_t теперь сообщает точное число не поместившихся попаданий, а не просто факт переполнения. Эффект виден на эталонном прогоне (b=10, k<5000): GPU снимает 43 кандидата вместо 31, и на PRP уходит 275 вместо 287 — то есть найденные, но потерянные ранее делители вернулись.

7. Что проверено

Проверка	Результат
`cargo test` (сито, краевые случаи, границы сегментов)	8/8
b=10, k < 3000, CPU	`[2, 19, 23, 317, 1031]` — совпадает со списком известных простых репьюнитов
b=10, k < 3000, с GPU	тот же список; GPU снял 11 кандидатов, ложных срабатываний 0
b=2, k < 700	`[2,3,5,7,13,17,19,31,61,89,107,127,521,607]` — ровно показатели Мерсенна
Каждая запись `worklog.jsonl`	делитель действительно делит R_k; ни одно простое не отсеяно
Ветвление ядра против честного деления	50 388 комбинаций (b,k,q), 0 расхождений
P-1 (принудительный режим)	снял 66 кандидатов, все 66 делителей проверены — настоящие
Montgomery-арифметика на GPU	134 вектора против эталона: Mont64 43/43, Mont96 50/50, Mont128 41/41
Ядра `rh_tf_k64/k96/k128`	каждое находит настоящий делитель репьюнита в узком окне по m
Два потока на одну карту (`gpu_devices=[0,0]`)	результат не меняется, ложных срабатываний нет
Компиляция GWNUM-бэкенда	с настоящими заголовками Prime95 SDK — чисто
GWNUM на живых числах	вердикты совпали с GMP на 6 размерах (1 050…69 761 бит) и на всех 353 кандидатах диапазона k=3000…6000
R_49081(10), 163 041 бит	GWNUM нашёл PRP за 11.9 с, GMP-верификация подтвердила
Ускорение GWNUM против GMP	2.9× (6.6k бит) → 8.1× (70k бит); на диапазоне 3000–6000 — 4.4×
Откат при отсутствии GWNUM	`prp_backend="gwnum"` без библиотеки считает на GMP с предупреждением, кандидат не теряется
`cargo clippy`	без замечаний
Double-check на лету	179 «составных» пересчитано вторым бэкендом, 0 расхождений
`--verify` на чистом журнале	119 делителей + 5 PRP + 179 составных, 0 расхождений, код 0
`--verify` на испорченном журнале	пойманы все три подделки: фальшивый делитель, ложный PRP и потерянная находка; код 1
Длинный прогон, k = 5000…60000	итог оптимизаций: 1:18:48 → 14:28 (в 5.4 раза), находка R_49081(10) на месте, ошибок 0
Перепроверка журнала	3376 делителей, находка подтверждена, 50 пересчитанных составных — расхождений нет
Двенадцать оснований, b = 2…13, k < 600	списки PRP совпали с независимым эталоном для каждого b, включая вырожденные случаи b=4, 8, 9

Воспроизвести нативные проверки: bash native/tests/run_tests.sh (векторы перегенерируются python native/tests/gen_vectors.py native/tests/mont_vectors.h).

Что было найдено проверками и исправлено

Mont96 давала неверный результат на всех входах: 96-битный REDC делал сдвиг на 128 бит вместо 96 (терялись младшие 32 бита) и «телепортировал» перенос через лимб. В Mont128 та же цепочка случайно корректна, потому что там младшие лимбы зануляются по построению. Обе ширины теперь ведут цепочку переносов аппаратно (add.cc → addc.cc → addc.cc → addc).
Делитель мог совпасть с самим числом: R_2(10) = 11 отсеивался «делителем» 11. Теперь требуется собственный делитель q < R_k(b) — в сите, в ядре и в rh_prp_verify_factor.
Базы Miller–Rabin вне [2, N−2]: для N = 3 раунд с базой 3 объявлял простое составным, терялись R_2(2) = 3 и R_2(10) = 11.
Показатель мог исчезнуть из результатов при ЛЮБОЙ ошибке PRP. pipeline.rs на Err лишь писал строку в лог и не создавал записи в журнале: кандидат не «составной», не «PRP», а просто отсутствует — заметить пропажу можно только сверкой с внешним списком. Найдено расширенным пересчётом: b = 2, k = 1279 (M_1279, известное простое Мерсенна) пропал из прогона --kmax 20000, оставив ERROR ... PRP k=1279: internal error. Тот же k, посчитанный отдельно, проходит без ошибок всеми бэкендами — сбой проявлялся только под многопоточной нагрузкой. Теперь при ошибке кандидат досчитывается на GMP (не зависит ни от FFT, ни от состояния GWNUM), а если не вышло и это — пишется запись "status":"failed", чтобы показатель остался виден.
Оценка размера числа завышалась вдвое для b = 2. В prp_dispatch.c стояло est_bits = (k-1)·(floor(log2 b) + 1); для b = 2 это (k-1)·2 вместо (k-1)·1. Следствие: измеренный порог переключения на GWNUM (2500 бит) фактически применялся к b = 2 как 1250 бит — то есть ровно в той зоне, где GWNUM по собственным замерам проекта медленнее GMP (при 1017 бит — в 10 раз), и где его настройка на слишком малом числе и давала сбой. Заменено на честный log2. Это и есть первопричина пропажи M_1279.
Потерянная находка в GWNUM-пути: roundoff не проверялся после цикла. Ошибка округления контролировалась раз в 128 итераций, а финальное значение gw_get_maxerr() только записывалось в статистику и не сравнивалось с порогом. Скачок на одной из последних 127 итераций проходил незамеченным, и вердикт «составное» возвращался по заведомо испорченному остатку — то есть находка терялась молча. Найдено сплошным пересчётом последовательностей OEIS (paper/verify_sequences.sh): R_3181(23) объявлялся составным, хотя это PRP — член A204940. Диагноз: maxerr = 0.5000 (максимум из возможных), первое превышение порога — на итерации 14383 из 14384, то есть после последней периодической проверки на 14336-й. С margin = 1 FFT растёт с 768 до 1024, maxerr падает до 0.0015, вердикт верен. Частота: 1 случай на 96 пар (b, k) в скане b ∈ {3..26}, k ∈ {1009..10007}. Добавлена финальная проверка порога с возвратом RH_ERR_FFT_ERROR, что запускает штатный перезапуск с увеличенным FFT.
Кандидат исчезал, если GWNUM не укладывался в roundoff за 4 попытки. prp_dispatch.c возвращал ошибку, а pipeline.rs на ошибке лишь пишет лог и не создаёт записи в журнале — показатель пропадал из результатов без следа. Теперь такой случай досчитывается на GMP: медленно, но верно.
Фиктивная проверка Гербица в GWNUM-пути сравнивала неверное равенство и при этом никогда не запускалась (условие срабатывало после L² ≈ 4·10⁶ итераций — на два порядка больше практических размеров). Заменена схемой из раздела 1.4.
flatten() на итераторе строк журнала мог зациклиться на повторяющейся ошибке чтения — заменён на map_while.
rh_gmp_pool_reset() убран из Rust-обёртки: арена держит буферы из пула, и сброс вершины испортил бы живые числа.

Как проверялся сам механизм перепроверки. Мало показать, что на чистом журнале расхождений нет, — нужно убедиться, что механизм вообще способен их увидеть. Поэтому в журнал намеренно вносились три подделки: делитель q=13 для R_7(10) (не делит), R_19(10) помеченный составным (на деле простое) и R_5(10) помеченный как PRP (на деле 41·271). Все три пойманы, каждая со своим диагнозом, код возврата 1.

Итог оптимизаций: три прогона на одном диапазоне

k = 5000…60000, числа от 16 000 до 200 000 бит, одна и та же машина:

	прогон 1	прогон 2	прогон 3
время	1:18:48	1:03:30	14:28
снято ситом	1915	2298	2298
снято GPU	165	278	1078
PRP-тестов	3308	2812	2012
память	715 МиБ	744 МиБ	1241 МиБ
R_49081	✓	✓	✓ (8.3 с)

Между прогоном 1 и 2: ускоренное GPU-ядро, увеличенный буфер попаданий, калибровка тюнера. Между 2 и 3: счёт по модулю b^k−1, порог GWNUM 2500, исправленная модель стоимости PRP, все 16 логических ядер.

Любопытно, что прогноз «подешевевший PRP уменьшит выгоду от trial factoring» не подтвердился: GPU стал снимать вчетверо больше. Обе стороны экономики изменились одновременно, и правильная модель позволила тюнеру пользоваться отсевом там, где прежняя, занижавшая стоимость PRP вдвое, этого не делала.

Что показал длинный прогон

Час с четвертью непрерывной работы на диапазоне k = 5000…60000 (числа от 16 000 до 200 000 бит), все стадии включены:

стадия	снято кандидатов
CPU-сито (q ≤ 2²⁴)	1915 из 5388 (36%)
GPU trial factoring	165
P-1	0
PRP	3308 проверено, 1 найден

Память стабильна: 683 → 693 МиБ за 79 минут, пик 715 МиБ. Утечек нет.
Тюнер калибруется правильно: измеренная им скорость GPU (0.30 G/s) совпала с независимым замером бенчмарка (340 млн кандидатов/с), а оценка стоимости PRP (20 с на 163k бит) — с фактическими 16.7 с.
P-1 не запустился ни разу, и это верное решение. При B1 = 10⁵ этап стоит ~3.7·10⁵ модульных умножений, а весь PRP-тест числа в 163 041 бит — 1.6·10⁵. То есть P-1 здесь дороже теста, который должен экономить, и адаптивный режим честно его отключает. Стадия окупается лишь когда B1 много меньше размера числа — на миллионах бит, а не на сотнях тысяч.

Настройка стадий по измерениям

Три параметра были выставлены не «на глаз», а по замерам этого прогона.

Калибровка GPU при старте. В тюнере стояла вшитая оценка 8 млрд кандидатов/с — «как на топовой карте». Реальность GTX 1650 — 0.30 млрд, то есть промах в 27 раз, и пока не накопится статистика, глубина TF выбиралась наугад. Теперь GPU-поток делает один пробный запуск (~10 мс) и сообщает тюнеру настоящую скорость; стартовые константы заодно занижены — ошибка в сторону мелкого TF безопаснее, чем в сторону бесполезно глубокого.

Предел CPU-сита масштабируется по диапазону. Сито — самая дешёвая стадия, но его стоимость не зависит от k, а отдача — зависит: чем дороже PRP-тест, тем выгоднее снять кандидата заранее. Замер (k = 5000…15000):

предел q	2²⁴	2²⁶	2²⁷	2²⁸
время	0.2 с	1.0 с	2.25 с	5.4 с
отсеяно	30.1%	32.3%	33.1%	33.8%
память	179 МиБ	256 МиБ	275 МиБ	519 МиБ

Потолок поднят до 2²⁷ (дальше отдача падает, а память удваивается), но фактический предел считается как k_max · 2000 и лишь затем ограничивается конфигом: на коротком прогоне сито до 2²⁷ строилось бы дольше, чем идёт весь поиск.

Модель стоимости PRP была неверной — и это стоило дороже всего. В тюнере стояло t ∝ bits^1.15, то есть почти линейная зависимость. Подгонка по 2756 реальным замерам из журнала (числа 20 000…199 310 бит) даёт показатель 2.24:

	модель в коде	по замерам
формула	`bits^1.15`	`bits^2.24`
средняя ошибка	6.11 с	1.03 с
прогноз для 163 041 бит	9.7 с	17.8 с (факт 16.6)
прогноз для 200 000 бит	12.3 с	28.2 с

Показатель 2.24 физически обоснован: R_k(b) не имеет специальной формы, поэтому GWNUM работает через gwsetup_general_mod — IBDWT-умножение плюс редукция Барретта. Итерация стоит O(n log n), итераций n, итого ~n²·log n.

Последствие ошибки было прямым: тюнер вдвое занижал выгоду от trial factoring и почти не пользовался GPU. Замер на одном поддиапазоне (k = 9000…18000):

модель	снято ситом	снято GPU	ушло на PRP	доля GPU
`bits^1.15`	454	3	490	0.6%
`bits^2.24`	413	172	362	32.2%

Порог GWNUM был завышен впятеро. Стояло 10 000 бит, а измеренная точка пересечения с GMP — около 2000:

бит	1017	1994	3010	5017	9966
ускорение GWNUM	0.10×	1.06×	1.95×	3.21×	4.31×

Весь диапазон 2000…10 000 бит считался на GMP, теряя там от двух до четырёх раз. Порог опущен до 2500.

P-1 выключен по умолчанию. Прежняя эвристика брала вероятность успеха «≈4%» с потолка. Честный подсчёт в модульных умножениях: stage 1 стоит ≈ 2.9·B1, весь PRP-тест — ≈ bits. При B1 = 10⁵ и числе в 163 041 бит это 2.9·10⁵ против 1.6·10⁵, то есть этап дороже того, что экономит. Вероятность теперь оценивается функцией Дикмана ρ(ln m / ln B1) по таблице, а не выдумывается. P-1 начинает окупаться на миллионах бит — там его и включать.

Число потоков и виртуализация

Замер на i7-10700 (8 физических ядер, 16 логических, AVX2, L3 16 МиБ), диапазон k = 9000…18000:

потоков	4	8	12	16	20	24
время	41.1 с	22.7 с	17.0 с	15.8 с	15.8 с	15.8 с
мкс/бит	6.5	7.2	7.8	9.6	—	—

Оптимум — ровно число ЛОГИЧЕСКИХ ядер, дальше плато. Каждый отдельный тест при этом медленнее (9.6 против 7.2 мкс/бит), но суммарная пропускная способность выше: Hyper-Threading хорошо скрывает задержки памяти на FFT-нагрузке. Значение по умолчанию threads = 0 как раз и означает «по числу логических ядер» — менять его вручную не стоит.

Про WSL2. Процессорные инструкции там выполняются нативно, поэтому GWNUM работает на полной скорости, а PRP — это почти всё время работы. Виртуализация бьёт только по вызовам CUDA-драйвера, и это измерено: без CUDA Graphs запуск стоил 51.5 мс против 34.4 мс с ними, то есть накладные составляли треть. Но после ускорения PRP втрое GPU-стадия даёт около процента общего времени (15.68 с против 15.80 с), так что полный уход из виртуализации принёс бы единицы процентов — ценой портирования rh_alloc.c (mmap/madvise), affinity.rs (sched_setaffinity) и clock_gettime под MSVC.

Куда осмысленнее смотреть на процессор: у i7-10700 нет AVX-512, а GWNUM на нём заметно быстрее.

Достижимая глубина trial factoring

Множитель m передаётся в ядро как uint64_t, поэтому перебором достижимо только

$$q = 2mk + 1 < 2^{65},k.$$

Для k ≈ 10³ это ~2⁷⁵, для k ≈ 10⁶ — ~2⁸⁵. Практический вывод: 128-битная ветка (q ≥ 2⁹⁵) включается лишь при k > 2³⁰, то есть в обычном поиске работают Mont64 и Mont96. Если тюнер запросит глубину выше достижимой, диапазон по m молча ограничивается — см. clamp в gpu_worker.

Лицензия

Репозиторий лицензирован раздельно.

Исходный код искателя (src/, native/, benches/, build.rs, config/) — на выбор MIT или Apache-2.0, как и объявлено в Cargo.toml.

Статья и анализ (paper/: текст, рисунки, скрипты, производные данные) — CC BY 4.0.

Исключение: paper/data/*.txt — b-файлы, загруженные из OEIS. Они включены для воспроизводимости без обращения к сети и остаются на условиях OEIS (CC BY-NC-SA 4.0), а не на условиях лицензий этого репозитория. Загрузить их заново из первоисточника: bash paper/fetch_data.sh.

Показатели простых Мерсенна получены из открытого отчёта PrimeNet проекта GIMPS.

Цитирование

Докучаев Т. С. Схема наблюдения в статистике обобщённых репьюнитных простых:
точный условный критерий константы Ленстры — Померанса — Вагстаффа.
Версия 3. 2026. Препринт архиворг.ру AX-135915.
https://github.com/Wilps93/repunit-hunt/releases/tag/paper-v3
ORCID: 0009-0006-0510-5225

Navigation Menu

Repository files navigation

repunit-hunt

Статья

Структура репозитория

1. Математическая база

1.1 Почему только простые k

1.2 Форма делителей

1.3 Тест на простоту

1.4 Счёт по модулю b^k − 1

1.5 Контроль ошибок в GWNUM-пути

2. Архитектура

3. Сборка

3.1 Развёрнутое окружение (WSL2 на этой машине)

4. Запуск

4.1 Перепроверка журнала (double-check)

Возобновление работы

5. Настройка (config/default.toml)

6. Проверка результата

6.1 Производительность GPU-стадии

7. Что проверено

Итог оптимизаций: три прогона на одном диапазоне

Что показал длинный прогон

Настройка стадий по измерениям

Число потоков и виртуализация

Достижимая глубина trial factoring

Лицензия

Цитирование

About

Releases

Packages

Contributors

Languages

Name		Name	Last commit message	Last commit date
.cargo		.cargo
benches		benches
config		config
native		native
paper		paper
src		src
.gitattributes		.gitattributes
.gitignore		.gitignore
Cargo.lock		Cargo.lock
Cargo.toml		Cargo.toml
LICENSE-APACHE		LICENSE-APACHE
LICENSE-MIT		LICENSE-MIT
README.en.md		README.en.md
README.md		README.md
build.rs		build.rs
rust-toolchain.toml		rust-toolchain.toml

Latest commit

History

Folders and files

Repository files navigation

repunit-hunt

Статья

Структура репозитория

1. Математическая база

1.1 Почему только простые k

1.2 Форма делителей

1.3 Тест на простоту

1.4 Счёт по модулю b^k − 1

1.5 Контроль ошибок в GWNUM-пути

2. Архитектура

3. Сборка

3.1 Развёрнутое окружение (WSL2 на этой машине)

4. Запуск

4.1 Перепроверка журнала (double-check)

Возобновление работы

5. Настройка (config/default.toml)

6. Проверка результата

6.1 Производительность GPU-стадии

7. Что проверено

Итог оптимизаций: три прогона на одном диапазоне

Что показал длинный прогон

Настройка стадий по измерениям

Число потоков и виртуализация

Достижимая глубина trial factoring

Лицензия

Цитирование

About

Topics

Resources

Stars

Watchers

Forks

Releases

Packages

Contributors

Languages