Findprimes

This is a hobby research project of mine. I've been working on it off and on as my copious free time permits. It started off as a dig to the Quantum Computing Industry's claims about Prime Factorization - which is currently still stuck at factoring 21. As in 3 * 7. As a reference point, my Core i9-14900HX laptop can factorize ULLONG_MAX in less than 3 seconds using trial division.

The second - corresponding - part of this project can be found here in my Github repo: primefactors.

Four prime finders over the same range scanner. Each enumerates every prime in a decimal range and prints them in sorted order; they differ in how a candidate is decided and where the work runs.

program	primality method	where it runs
findprimesmp	trial division to sqrt(N)	CPU, pthreads
findprimesmpecm	N-1 Pocklington, N+1 Morrison, and their combined CRT bound, factored with GMP-ECM	CPU, pthreads
findprimesmpecmcgbn	the same three certificates	ECM stage 1 on the GPU via CGBN, curve space split across MPI ranks
findprimesmpecmcgbnmpi	the same three certificates	the same GPU kernels, with MPI splitting the range across ranks and the curve space again inside each slice

Your Linux distro should have GMP-ECM available. Recent Fedora and Ubuntu definitely do.

This was tested on Fedora 41 and 44 with CUDA 12.9 (Ada) and 13.2 (Blackwell). Hardware I used for testing:

• Laptop with an RTX 4080 Ada (sm_89) and a Core(TM) i9-14900HX.
• Intel NUC 13 Extreme with an RTX PRO 4500 Blackwell (sm_120) and a Core(TM) i9-13900K.

No probable-prime tests.

mpz_probab_prime_p() is called nowhere and Miller-Rabin appears nowhere. Every answer is a proof: primality by certificate, compositeness by an exhibited divisor or a failed Fermat or Lucas congruence. That constraint is what the whole design serves — it is why ECM is here at all, since a certificate needs N±1 factored and ECM is what factors it.

The four must agree candidate for candidate. Unrelated machinery over the same input is the strongest correctness signal available here, and verify.py cross-checks all four.

Building

%>> make            # all four
%>> make clean      # cleanup

Requirements, with the versions this was verified against:

	for	verified with
GMP	all four	system libgmp
GMP-ECM (gmp-ecm-devel)	the three ECM builds	7.0.6
CUDA toolkit	the two CGBN builds	12.9
CUDA toolkit	the two CGBN builds	13.2
CGBN headers in /usr/local/include/cgbn	the two CGBN builds	—
Open MPI	the two CGBN builds	5.0.5

The CPU pair builds warning-free under both clang++ and g++ (make CXX=g++), at C++17, C++20 or C++23. The CGBN pair is built by nvcc; GPU_ARCH defaults to sm_89, so change it for a different card.

To run either MPI build under mpirun, put the Open MPI bin directory on PATH:

%>> export PATH=/usr/lib64/openmpi/bin:$PATH
%>> mpirun --oversubscribe -np 4 ./findprimesmpecmcgbn    -s 1000000 -e 1000200 -T 8
%>> mpirun --oversubscribe -np 4 ./findprimesmpecmcgbnmpi -s 1000000 -e 1000200 -T 8

Open MPI 5 launches through PRRTE and finds prterun by PATH lookup, not by its own install prefix. Run without mpirun, either program is a single-rank job.

Usage

Usage: findprimesmpecmcgbnmpi [ -s <range-start> (default 18446744073709551615)]
                    -e <range-end> (required)
       [ -b <number-of-bits> (default 128)]
       [ -T <number-of-threads> (default 4)]
       [ -f <output-file> (default stdout)]
       [ -p (print header at the top)]
       [ -t (print prime discovery time)]
       [ -B <ecm-stage-1-bound> (default 2000, escalates)]
       [ -C <cpu-ecm-curves-per-split> (default 100)]
       [ -G <gpu-ecm-curves-per-launch> (default 4096)]
       [ -N <mpi-ranks> (default: as launched by mpirun)]
       [ -c <ranks-sharing-each-range-slice> (default 1)]
  run under mpirun to spread the range over several hosts:
    mpirun -np 8 --host a,b,c,d ./findprimesmpecmcgbnmpi \
        -e <end> -N 8 -c 2   # 4 range slices, 2 ranks each

That is findprimesmpecmcgbnmpi -h; the other three print the same flags minus the ones they do not have. -B/-C are on all three ECM builds, -G on both CGBN builds, -N/-c only on findprimesmpecmcgbnmpi, and findprimesmp has none of them. Note -h exits 1 like an argument error rather than 0. Primes go to stdout (or -f); progress and the summary line go to stderr.

%>> ./findprimesmpecm -s 1000000 -e 1000200 -T 4 2>/dev/null
1000003
1000033
1000037
...
1000193
1000199

The default range start is 2^64 - 1: these are big-number benchmarks, not fast prime sieves.

How primality is decided

findprimesmp divides by every odd number up to sqrt(N). Correct, and completely impractical past 2^64 — one prime candidate there costs ~37 s.

The ECM builds resolve a candidate in five stages, cheapest first:

1. Cheap compositeness proofs. A small divisor, then one Fermat congruence to base 2. Between them these settle nearly every composite for the price of one modular exponentiation.
2. N-1 (Pocklington). Factor a divisor F of N-1 with F > sqrt(N) — trial division for the small part, ECM for the rest — then for each prime q | F exhibit a base a with a^(N-1) ≡ 1 (mod N) and gcd(a^((N-1)/q) - 1, N) = 1. Every prime factor of N is then 1 mod F and so exceeds sqrt(N); two of those will not fit inside N.
3. N+1 (Morrison). When N-1 will not come apart far enough. Factor a divisor F of N+1 with (F-1)² > N, pick D with Jacobi symbol (D/N) = -1, and exhibit Lucas parameters giving U_{N+1} ≡ 0 (mod N) and gcd(U_{(N+1)/q}, N) = 1. N-1 and N+1 are unrelated numbers, so this is an independent attempt rather than a retry.
4. Combined bound. When neither side reaches sqrt(N) alone, both congruences still hold at once, and the Chinese remainder theorem turns the pair into a lower bound on the smallest prime factor. Each side then needs only about N^(1/4).
5. Trial division, only if no certificate could be produced.

Only half of N±1 ever has to be factored. Insisting on the complete factorisation is the obvious reading of the theorems and it stalls a couple of hundred bits in.

GPU and MPI

ECM stage 1 is a scalar multiplication [s]P on a Montgomery curve mod M, and every curve is an independent walk — one CGBN instance per curve, thousands per launch. Curves use the inversion-free parameterisation (a24 = d, start (2:1)), which matters for more than speed: building a curve needs no modular inverse, so d is just a counter, and a counter is what lets MPI ranks take provably disjoint slices of curve space without exchanging anything.

Both CGBN builds run that identical kernel. What separates them is what MPI divides.

findprimesmpecmcgbn — curves

Every rank scans the whole range and the ranks divide ECM's curve space. Rank r on attempt a takes slot a*nranks + r; an MPI_Allreduce(MIN) elects the lowest-ranked finder and MPI_Bcast hands its factor to everyone, so all ranks leave with the same factor. Scanning the range on every rank is duplicated work by design — it is the cheap half, and replicating it keeps the ranks in lockstep with no communication.

k hosts try k times the curves per attempt. That is a throughput win on ECM's probabilistic curve search, and it is the only thing here that makes a single hard candidate finish sooner.

findprimesmpecmcgbnmpi — the range, then curves inside it

The same kernel and the same certificates, laid out on two axes:

-N 8 -c 2      slice 0        slice 1        slice 2        slice 3
               rank 0 rank 1  rank 2 rank 3  rank 4 rank 5  rank 6 rank 7
               \____________/ \____________/ \____________/ \____________/
                curve split    curve split    curve split    curve split

-N says how many of the ranks mpirun launched take part; -c says how many of them share each range slice. The active ranks are cut into ceil(N/c) groups of c consecutive ranks. Group g owns range slice g — a contiguous, disjoint run of odd candidates — and the ranks inside a group split curve space between themselves exactly as findprimesmpecmcgbn does.

So the two programs are the ends of one dial. -c 1 (the default) is a pure range split; -c N is the pure curve split findprimesmpecmcgbn already was; anything between is a rectangle — ceil(N/c) candidates in flight, c GPUs on each.

The axes are not interchangeable.

	range split (-c 1)	curve split (-c N)
makes one hard candidate finish sooner	no	yes
turns k nodes into k× candidates/second	yes	no
communication during the scan	none	a collective per ECM attempt
Phase B	full thread pool	one thread, strictly in order

That last row is the real cost of the inner axis. Ranks in a group must reach every collective in the same order, so within a group Phase B walks a sorted list on a single thread, one rank decides every CPU-side ECM split for the group, and the GPU-or-CPU choice has to be unanimous (a rank without a GPU would return a different factor and desynchronise everything after it). A group of one has none of those constraints, because it has nobody to agree with.

Groups, on the other hand, never have to agree about anything. They share no candidate, so a cluster with GPUs on some nodes and not others simply runs its groups at different speeds and still prints the same primes.

Reporting. Each group's root packs its primes as length-prefixed records — a 32-bit word count, then that many little-endian limbs — and one MPI_Gatherv on the whole job hands them to rank 0, which merges them into the sorted set its own slice filled and prints the union. Nothing on the wire is fixed-width, so 128-bit primes and 16,000-bit primes travel through the same code. Non-root ranks of a group send nothing (a group replicates its slice, so the root already speaks for them), and rank 0 sends nothing (its primes are already in the set).

Because the slices are disjoint and a std::set re-sorts the union, the printed list is byte-identical at every -N, -c and -T — verified against findprimesmpecm across -np/-N 1–8 and -c 1–4 at 20, 64, 96, 128, 256 and 384 bits, including the shapes that do not divide evenly (-np 3 -c 2, -np 6 -c 4) and the degenerate ones: -c above the rank count, -N above the launched world, and a range with fewer candidates than groups. The 384-bit sweep is the one that exercises everything at once — it is where N-1 starts failing and the N+1 Morrison certificate has to carry a candidate, and where the GPU does real work (14 splits in 36 launches at -c 1).

Both flags only select from what mpirun launched. A job's world size is fixed at launch, so -N splits MPI_COMM_WORLD down and says so when asked for more than it has; it cannot conjure ranks up. (Growing the job would mean MPI_Comm_spawn, which places children on the local host unless separately handed a host list — so it would not reach other nodes, which is the point.) Reaching several nodes is mpirun's job:

%>> mpirun -np 8 --host a,b,c,d ./findprimesmpecmcgbnmpi -e <end> -N 8 -c 2

On a single GPU, more ranks is slower on either axis — ranks queue on the one device, and with -c > 1 they also replicate the scan. Both designs target one rank per GPU across hosts; that configuration is not measured here, only that the distribution is correct.

Performance

10^18 + 0..1000 — 23 primes, identical output from both:

	findprimesmp	findprimesmpecm
wall, -T 8	70.5 s	0.04 s
wall, -T 32	33.5 s	0.04 s

Scaling, same 1000-wide window at -T 8:

	10^12	10^15	10^18
findprimesmp	0.07 s	1.83 s	70.47 s
findprimesmpecm	0.00 s	0.01 s	0.03 s

Trial division tracks sqrt(N) as advertised. The certificates are flat in N because their cost follows how hard N±1 is to factor.

Bit Width

There is no width cap in any of the ECM builds. What bounds the answer is the certificate needing a factored divisor of N±1 to reach sqrt(N) — so the ceiling depends on how N±1 factors, not on how big N is:

N-1	prime	width	wall
smooth by construction	6223*2^16384+1	16,397 bits	4.1 s
arbitrary	nextprime(2^768)	769 bits	280 s
arbitrary	nextprime(2^1024)	1,025 bits	not certified in 900 s

For Proth-form primes k*2^n+1, N-1 = k*2^n is factored by inspection and width barely registers. For an arbitrary prime, roughly half of N-1 has to come apart, and that wall arrives just under 1,024 bits. A range scan costs what its unluckiest candidate costs.

Going further needs ECPP, which GMP-ECM does not provide.

Testing

%>> pip3 install sympy
%>> python3 verify.py          # needs sympy

For each binary: 8 ranges × 10 thread counts against a computed reference, asserting both that the output matches and that every -T agrees — then a cross-check between the binaries. The ECM builds additionally get ranges at 2^64, 2^96, 2^256 and 2^384.

verify.py drives both CGBN builds as a single MPI rank, which exercises neither program's distribution. That is covered separately, and every shape must produce byte-identical output:

#!/bin/bash

export PATH=/usr/lib64/openmpi/bin:$PATH

for np in 1 2 4; do
    mpirun --oversubscribe -np $np ./findprimesmpecmcgbn -s ... -e ... -T 8
done

for np in 1 2 4 8; do
    for c in 1 2 4; do
        mpirun --oversubscribe -np $np ./findprimesmpecmcgbnmpi \
            -s ... -e ... -T 8 -N $np -c $c
    done
done

For findprimesmpecmcgbnmpi that is the check that earns its keep: a group whose ranks fall out of step pair an Allreduce with the wrong cofactor, and a Gatherv with a miscomputed displacement drops a slice — both show up at once as a wrong or missing prime. Worth including in the sweep: -c larger than the rank count (clamped to one slice), -N larger than the world (diagnosed, runs with what it has), and a range with fewer candidates than groups (empty slices must retire quietly, not error).

Two intentional quirks the reference model encodes: only odd candidates are tested, so 2 is never reported, and 1 is reported as prime. Both are preserved deliberately — all four programs agree, and the cross-check between them is worth more than the tidier answer.

License

MIT License - see LICENSE, which is the verbatim MIT text. Copyright (C) 2019-2026 Stefan Teleman.

AI

I used Claude Opus 5.0 for some of the very difficult parts of the primality proofs and to test / verify / find bugs.