Christoph Heinrichs and Fabian Januszewski developed an open-source CUDA implementation of the Self-Initializing Quadratic Sieve that factors semiprimes (up to 512 bits) on GPU hardware. Which started as a proof-of-concept project we implemented an integer factorization pipeline which runs entirely on GPU, including sieving, large prime processing, linear algebra stage (Block-Wiedemann), and square-root step, achieving SOTA performance for integers in the 100 digit range, e.g. RSA-100. We validated our implementation factoring RSA-110 & RSA-120 on consumer grade hardware and RSA-130 and RSA-140 on clusters of A100 and H100 GPUs respectively.

Implementation details:

• six-stage pipeline, ~63,000 lines of code, ~140 CUDA kernels
• GPU Architectures: Tested on NVIDIA GPUs from Turing (2018) to Blackwell (2025)
• Current version: 1.0.4, first public release 03.06.2026 (v1.0.0)

We maintain two github repositories, one for cuda-mpqs including block-wiedemann and one for block-wiedemann only.

• License: LGPL-3.0-only + NVIDIA CUDA Toolkit exception
• Source code cuda-mpqs (CH+FJ): https://github.com/drjanosch42/cuda-mpqs
• Block-Wiedemann Linear-algebra solver (FJ): https://github.com/drjanosch42/block-wiedemann

Performance

For RSA-100 we achieve SOTA performance. On a consumer RTX 5070 Ti GPU cuda-mpqs factors RSA-100 in < 2 minutes, on a data-center class H100 SXM GPU in < 1 minute.

The sieve dominates wall-clock time on every platform (≈77–91%). Runtimes refer to v1.0.4.

Scaling

Runtime for RSA-100 class moduli scales on ‘large’ GPUs anti-proportionally with CUDA core count. This also applies to the sieve and linear algebra stages individually and is visualized in the graphs on the right hand side which is based on v1.0.0 timings. Since then we saw a speed improvement in v1.0.4 across the board.

Factorizations

So far we factored RSA-100, RSA-110, RSA-120, RSA-130, RSA-140. Factoring RSA-140 took less than 7 hours wall time on a 4-node cluster with 16 H100-SMX GPUs, including ~40% oversieving. We emphasize that we only implemented a single large prime variant at this stage and maintained an LP contribution <60%.

Under the hood

The sieve is a three-kernel pipeline

polynomial generation → meta-sieve → sieve-and-scan

running over a hypercube of factor-base primes with Gray-code `b`-traversal. It supports a double-buffered, zero-synchronization batch mode that overlaps sieving with post-processing, plus CUDA-graph capture for low-launch-overhead replay. A single large-prime variant pairs partial relations through a GPU slab hash table with a fully asynchronous hot loop (no `cudaStreamSynchronize`).

Optional experimental matrix preprocessing follows a CADO-NFS-style expanded-matrix merge/filter design with a fully GPU-resident packed backend.

Block Wiedemann runs in three sub-stages

Krylov sequence → Berlekamp–Massey/`lingen` → reconstruction

over an SpMM engine that auto-selects among 10 distinct GF(2) sparse-matrix-vector formats and matching specialized CUDA kernels per matrix row block.

• All-GPU pipeline. Every stage — sieve, post-processing, large-prime matching, matrix build, linear algebra, and square root — runs on the GPU. The same code scales from a 25-watt Jetson Orin Nano to a multi-node GPU cluster.
• 512-bit arithmetic on the GPU. A custom `uint512` type implements full add/sub/mul/div/mod with Montgomery (CIOS) modular multiplication. Every operation is annotated `__host__ __device__`, so the identical code path validates on the CPU.
• No merge tree. The GPU matrix backend packs each entry into a single `uint32_t` as `(column_index << 8) | exponent` and carries `sqrt_Q` products straight through merges via Montgomery multiplication — eliminating the classical merge tree.
• 10 ways to multiply a sparse matrix. The Block Wiedemann engine measures each matrix and auto-selects from 10 GF(2) SpMM formats, from dense bitslice to Golomb-coded streams.
• Distributed sieve, one code path. Solo, cluster-coordinator, and cluster-worker all run the same sieve loop; the distributed mode is injected through a `nullptr`-by-default data-tap, so solo runs pay zero cluster overhead. Workers stream relations to the coordinator over a length-prefixed, CRC32-checked TCP wire protocol.
• Self-tuning. A four-stage autotuner with a convex sieve parameter optimizer and a persistent history database lets repeated runs at a given size start from previously discovered optima.