Skip to content

Mersenne numbers

Numerisect 0.13.0

\(M_p = 2^p - 1\). This page covers what Numerisect can do with them, and is honest about where each tool stops.

The short answer

Question Tool Practical reach
Is \(M_p\) prime? Lucas–Lehmer (number_theory.gp) exponents in the thousands
Find a small factor of \(M_p\) Mersenne trial factoring, below exponents in the millions
Hunt for additional factors Staged factor hunt, below materialized \(M_p\), currently \(p\le10^6\)
Factor \(M_p\) completely staged hunt or SNFS, when every cofactor can be resolved strongly input-dependent
Is this a known Mersenne prime? the 56-class prime classifier catalogue lookup

The middle row is the one that changes what is possible, and it is the reason this page exists.

Trial factoring, and why it reaches so far

Every prime factor \(q\) of \(M_p\), for odd prime \(p\), satisfies two classical congruences:

\[ q \equiv 1 \pmod{2p}, \qquad q \equiv \pm 1 \pmod 8 . \]

The first follows because the order of \(2\) modulo \(q\) is exactly \(p\), so \(p \mid q-1\) and \(q\) is odd. The second is the condition for \(2\) to be a quadratic residue modulo \(q\). Together they confine every candidate to \(q = 2kp + 1\) with \(q \equiv \pm 1 \pmod 8\). The second condition retains two of the four possible odd residue classes, so it discards one half of that progression before the modular test.

If the odd exponent \(p\) is composite, a divisor \(q\) can have \(d=\operatorname{ord}_q(2)\) equal to any divisor \(d>1\) of \(p\), not necessarily \(p\) itself. Numerisect therefore asks PARI/GP to factor the exponent, enumerate every such order divisor, and search

\[ q=2kd+1,\qquad d\mid p,\quad d>1. \]

This includes the algebraic factors inherited from \(M_d\mid M_p\). For example, \(1603=7\cdot229\), so \(M_7=127\) divides \(M_{1603}\); the search finds 127 at \(d=7,k=9\). A prime-exponent-only \(q=2kp+1\) search would miss it.

The exponent \(p=2\) remains the trivial exception: \(M_2=3\). The specialized search accepts odd prime or composite exponents starting at 3.

Membership is then decided by a single modular exponentiation:

\[ q \mid M_p \iff 2^p \equiv 1 \pmod q . \]

PARI/GP also completes isprime(q) before the application calls a returned candidate a prime factor; the pseudoprime filter is only a fast rejection step.

\(M_p\) is never constructed. Everything happens modulo the candidate. That is the whole point: \(M_{1000151}\) has 301,076 decimal digits, and its factor \(2000303\) is found at \(k = 1\). Materializing that target merely to test a small candidate would waste memory and arithmetic; modular powering makes the work depend on the candidate size and exponent instead. This is the same reason a Mersenne search trial-factors an exponent before committing to a Lucas–Lehmer test.

There is no formula that predicts the smallest factor—and therefore no mathematically correct fixed value of \(k\) that can be inferred from \(p\). The application's Automatic mode handles this honestly: the native search advances \(k\) until it finds the first factor, reaches the 50,000,000 safety ceiling, or consumes the selected time budget. It reports the largest \(k\) actually tested. Manual mode instead checks every candidate through the bound supplied by the user and reports every factor in that finite range. Factors found before a timeout are retained in either mode.

curl --cookie jar --header 'Content-Type: application/json' \
  --data '{"exponent":1000151,"k_limit":5}' \
  http://127.0.0.1:8765/api/factor-lab/mersenne-factors

Results this actually produced

Every factor below was found by the route above and then verified in a separate PARI/GP session, not by the routine that found it. For each one: \(q\) is prime, \(k\) is a genuine integer in \(q = 2kp+1\), \(q \equiv \pm 1 \pmod 8\), and \(2^p \equiv 1 \pmod q\). Where an exponent has several factors, the product of all of them was also checked to divide \(M_p\).

exponent \(p\) decimal digits of \(M_p\) factors found
2,000,003 602,061 160000241, 8924785387159
30,000,001 9,030,901 1380000047
50,000,017 15,051,505 1131900384847, 3615901229407, 355682820932119
70,000,027 21,072,108 9520003673
400,000,009 120,412,001 10876800244729
600,000,001 180,617,998 27600000047, 2418000004031
999,999,001 301,029,695 357999642359

The last row is the one to look at. \(M_{999999001}\) has 301,029,695 decimal digits. Representing such a value is already substantial, while applying a general-purpose factorization to it is not realistic. Testing twenty million candidates, up to roughly \(4 \times 10^{16}\), took 10.7 seconds in the recorded local run because the target is never built. That timing describes one machine and build, not a portable benchmark.

Two known Mersenne prime exponents, 6,972,593 and 20,996,011, were searched as controls and correctly yielded nothing.

How the search is made fast

PARI/GP searched this progression correctly but at about 1.1 million candidates a second, using generic arbitrary-precision arithmetic on one core. A compiled helper, numerisect-mfactor, now does the scanning for any finite range and sustains about 1.3 billion. Three things account for it:

  1. The progression is sieved first. \(q = 2kd+1\) is divisible by a small prime \(r\) exactly when \(k \equiv -(2d)^{-1} \pmod r\), an arithmetic progression in \(k\), so a sieve removes those \(k\) with no modular exponentiation at all. About 96% of the range goes this way at the default bound.
  2. Candidates below \(2^{64}\) use 64-bit arithmetic with a 128-bit intermediate, instead of arbitrary precision. This is where nearly all the work lands.
  3. Every core is used. The surviving tests are independent.

Above \(2^{64}\) the helper falls back to GMP. That path is slower but correct, and the report counts how many candidates needed it, so a caller can see which ran.

PARI/GP keeps both jobs it should keep. It factors the exponent and enumerates the order divisors before the scan, and it confirms afterwards that every reported \(q\) is prime and genuinely divides \(2^d - 1\). The helper is a scanner, not an authority: \(2^d \equiv 1 \pmod q\) makes \(q\) a divisor and says nothing about it being prime. Nothing reaches a report on the scanner's word alone.

Automatic mode, which stops at the first factor and preserves findings when its budget expires, still runs in PARI/GP. If no C compiler is available the whole search falls back there too, and the answers are the same.

The optional GPU accelerator

The modular exponentiations are perfectly independent, which is what a GPU is for. numerisect_mfactor_kernel.cu runs them on the device using Montgomery multiplication, which replaces a 128-bit division per multiply with two 64-bit multiplies and a shift. Candidates the device cannot cover are counted as deferred rather than skipped, because silently dropping them would turn an untested range into an apparent absence of factors; Numerisect sends such ranges to the C helper instead.

Building it needs the CUDA toolkit. When nvcc is present, Numerisect compiles numerisect_mfactor_cuda.cu into data/tools/bin/numerisect-mfactor-cuda on first use, and rebuilds it whenever the host program or the kernel file it includes is newer than the binary. Without nvcc or without a device, the C helper does all the work and the answers are the same.

Verified on hardware. On an RTX 5090 (compute capability 12.0) the device returns the published factorizations for every case tested, matches PARI/GP's prime divisors over the same ranges, and agrees exactly with the compiled CPU helper. The suite skips these tests cleanly where no CUDA toolkit is installed.

What it is worth, measured: about twelve times the C helper, once the sieve moved onto the device.

wall clock, \(k \le 4 \times 10^9\)
C helper, 24 cores 23.3 s
CUDA build, RTX 5090 2.0 s

Getting there took two corrections. The first GPU build sieved on the host, gathered the survivors and copied them across, and measured slower than the C helper: the work per candidate is only about thirty 64-bit Montgomery squarings, which does not pay for moving that candidate over PCIe. Moving the sieve onto the device made it slower still, because one thread per prime is catastrophically unbalanced. In a 67-million-entry segment the thread holding \(r = 3\) performs 22 million serialized writes while a thread near the sieve bound performs about seventy.

Splitting the work by (prime, chunk) rather than by prime fixed it. A small prime's work now spreads across every chunk instead of landing on one thread, and writes stay inside a narrow window. Nothing but the prime table, once, and the handful of hits ever crosses the bus.

Past \(2^{64}\): two-limb Montgomery

Single-limb Montgomery needs \(q < 2^{63}\), which at an exponent near \(10^9\) is reached around \(k = 4.6\) billion. The device used to stop there and hand the rest to the C helper, which tests wide candidates with GMP at roughly a hundredth the rate.

The same algorithm now runs on two limbs, covering \(q < 2^{127}\). At that exponent the bound moves to \(k\) near \(8 \times 10^{28}\), which is past anything the k ceiling allows, so in practice the device covers the whole range. Both kernels run over each segment and their results are merged.

The arithmetic is CIOS Montgomery multiplication for a two-limb odd modulus. One detail is worth recording: \(R \bmod n\) is built by doubling 1 exactly 128 times, not by reducing \(2^{128}\) with repeated subtraction, which runs about \(2^{128}/n\) times and does not terminate in practice for a small \(n\).

Measured across the boundary, where the C helper must use GMP for a third of the candidates:

\(k \le\) C helper CUDA build
\(10^{10}\) 55.2 s 6.4 s 8.7x
\(2 \times 10^{10}\) 121.3 s 13.9 s 8.7x

Both agree exactly on the factors found, and the device defers nothing.

Both scanners are checked against each other in the test suite, and against the published factorizations.

A hit is a divisor, not a prime factor

\(2^d \equiv 1 \pmod q\) makes \(q\) a divisor of \(2^d - 1\) and says nothing about \(q\) being prime. Composite divisors genuinely occur: \(2047 = 23 \cdot 89\) divides \(2^{11} - 1\), and the kernel's divisibility test is correctly true for it. The pipeline sieves those out beforehand and PARI/GP confirms primality afterwards.

Two facts handed to GMP-ECM rather than rediscovered

The staged hunt gives the unresolved cofactor to P−1, P+1 and ECM. Two properties of \(M_p\) are known before any of them starts, so they are passed as parameters:

  • The cofactor divides \(2^p-1\). GMP-ECM's -base2 reduces modulo that form instead of performing a general division.
  • Every prime factor \(q\) satisfies \(q = 2kp+1\), so \(2p \mid q-1\). For P−1 the group order is \(q-1\), so -go 2*p preloads a known divisor of it and \(B_1\) no longer has to reach \(p\) itself.

The second is worth a measurement. The factor \(2000303 = 2\cdot 1\cdot 1000151+1\) of \(M_{1000151}\) has \(q-1 = 2 \cdot 1000151\), whose large prime is \(1000151\). Plain P−1 at \(B_1=1000\) cannot find it; with -go 2*1000151 it falls out in stage 1 at the same bound.

P+1 gets no group order, because \(q+1\) need not be divisible by \(2p\), and ECM gets none because a curve order is not known in advance. Numerisect refuses -go for ECM rather than passing a number nobody can justify.

Finding nothing is inconclusive

For prime \(p\), an exhausted \(k\) range means no factor of the form \(2kp+1\) exists below the bound searched. For composite \(p\), it means every selected \(q=2kd+1\) order progression was exhausted through that bound. Neither outcome is evidence that \(M_p\) is prime.

Finding one or several trial factors is still not a complete factorization. Dividing them out can leave a cofactor almost as large as \(M_p\): for \(M_{87083}\), the displayed eight-digit factor leaves a cofactor of roughly 26,207 decimal digits. Completely factoring a cofactor of that size is not currently practical. Numerisect therefore says factors found, never all factors, unless an engine has actually resolved and verified every cofactor.

Staged factor hunt and factor inventory

The second form on Factor integers → Mersenne numbers joins the native operations that are useful after trial factoring:

  1. PARI/GP exhaustively searches the selected finite \(k\) range.
  2. PARI/GP constructs \(M_p\), verifies that each reported divisor divides it, removes every occurrence, and records exact multiplicities.
  3. GMP-ECM runs Pollard \(p-1\), Williams \(p+1\), and an elliptic-curve campaign on the exact remaining cofactor, with the selected CPU count exposed to its native OpenMP runtime.
  4. After every discovery, PARI/GP reconciles the entire inventory again and attempts a bounded rigorous primality proof for the cofactor.

The web response shows only a bounded cofactor preview. The automatically saved report contains the exact complete cofactor, even when it has tens of thousands of digits. The status vocabulary is deliberate:

  • proven_prime means PARI/GP completed a rigorous proof;
  • composite means the cofactor is definitely composite;
  • probable_prime means screening succeeded but the proof budget expired;
  • unknown means even bounded screening did not finish;
  • complete prime factorization: yes appears only when all reported divisors and the final cofactor are rigorously prime.

P−1, P+1, and ECM are factor-discovery methods, not exhaustive searches. Finishing all three without a new divisor leaves the result incomplete. For a cofactor of practical NFS size, it can subsequently be submitted to the ordinary SIQS/NFS/SNFS workflow. A 26,207-digit cofactor such as the one left by the known factor of \(M_{87083}\) is far beyond a realistic complete NFS factorization.

curl --cookie jar --header 'Content-Type: application/json' \
  --data '{"exponent":87083,"trial_k_limit":100000,"trial_seconds":60,"stage_seconds":60}' \
  http://127.0.0.1:8765/api/factor-lab/mersenne-hunt

Complete factorization

A Mersenne number exposes algebraic structure that a special number field sieve can use: \(M_p = 2^p - 1\) gives the symbolic relation \(x^p-1\) at \(x=2\). The special-form analyser reports that relation and a difficulty indicator; the selected NFS engine still owns production polynomial selection and parameters:

2^1061-1  ->  homogeneous, n = 2^1061 - 1
              SNFS suitable: yes
              SNFS polynomial: x^1061 - 1
              SNFS difficulty: 319

Hand the number to the ordinary factoring pipeline to actually run it. A displayed special-form relation is advice, not a claim that the factorization is practical or that an engine-ready polynomial has already been selected.

The analyser used to hang here

Detecting the algebraic split called factor() on \(\Phi_n(b)\), and \(\Phi_n(b)\) is the whole input whenever \(n\) is prime. Asking for a special-form report on \(2^{1061}-1\) therefore asked PARI to factor a 320-digit number inside a metadata routine, which exhausted its budget and returned nothing at all. The split is now attempted only below a size cap, and above it is reported as not attempted, which is inconclusive rather than a claim that no algebraic factor exists.

Primality

nt_lucas_lehmer implements the Lucas–Lehmer test: \(M_p\) is prime exactly when \(s_{p-2} \equiv 0 \pmod{M_p}\) for \(s_0 = 4\), \(s_{i+1} = s_i^2 - 2\). This one is a proof, not a probable-prime test. It constructs \(M_p\), so its reach is governed by memory and time rather than by the arithmetic, and it is far slower than trial factoring: on a large exponent, search for a factor first.

pl_lucas_lehmer_riesel covers the related \(k \cdot 2^n - 1\) family.

Lucas–Lehmer's nonzero final residue proves that \(M_p\) is composite but does not identify a divisor. In the application, a composite result therefore links directly to Factor integers → Mersenne numbers, where the separate progression search tests \(q=2kp+1\) candidates without constructing \(M_p\).

Routes

POST /api/factor-lab/mersenne-factors   trial-factor over every odd order divisor d of p
POST /api/factor-lab/mersenne-hunt      trial, P-1, P+1, ECM, and exact reconciliation
POST /api/factor-lab/special-form       recognise the form, report the SNFS polynomial
POST /api/number-theory/lucas-lehmer    prove primality of M_p