# Estimating the probability of sucess of Pollard's p-1

I'm trying to estimate the probability that Pollard's p-1 factorization in its two-stages variant finds a factor of an RSA modulus product of $$k$$ random $$b$$-bit primes, as a function of the bounds $$B_1$$ and $$B_2$$ used.

That's about $$1-(1-P)^k$$, or about $$k\,P$$ for low $$P$$, where $$P$$ is the probability that a $$b$$-bit prime $$p$$ has $$p-1$$ with its largest prime factor less than $$B_2$$ and its second-largest prime factor less than $$B_1$$.
My (incorrect) idea was that we must have $$P\ge\rho\Bigl(\frac{b-1}{\log_2(B_2)}\Bigr)\rho\Bigl(\frac{\log_2(B_2)}{\log_2(B_1)}\Bigr)$$ where $$\rho$$ is Dickman's function [*], on the assumption that $$p-1$$ behaves like a random even integer, that the probability of the largest prime factor $$q$$ being less than $$B_2$$ is about $$\rho\Bigl(\frac{b-1}{\log_2(B_2)}\Bigr)$$; and then the conditional probability of the second-largest prime factor being less than $$B_1$$ is about $$\rho\Bigl(\frac{\log_2(q)}{\log_2(B_1)}\Bigr)$$, which is at least $$\rho\Bigl(\frac{\log_2(B_2)}{\log_2(B_1)}\Bigr)$$.

Problem is, my estimate is much larger than the one given by gmp-ecm with option -pm1 15e9 [#] which estimates its own probability of success as:

Using B1=15000000000, B2=19849281594647716 (..)
Probability of finding a factor of n digits (assuming one exists):
20      25      30      35      40      45      50      55      60      65
0.84    0.58    0.32    0.15    0.063   0.023   0.0078  0.0024  0.00068 0.00018

For $$n=65$$ decimal digits, that is $$b=216$$, I get $$P$$ at least $$0.0028$$, which is over 15 times more than ECM's estimate $$0.00018$$.

What's the cause of that discrepancy?

[*] Dickman's function $$\rho$$ is defined for $$u\in\mathbb R^+$$ by $$\begin{cases} \rho(u)=1 & \text{if }0\le u\le1\\ \rho(u)=1-\int_1^u\rho(t-1)/t\;dt& \text{if }u>1 \end{cases}$$ As $$B$$ grows, the probability that a random integer about $$B^u$$ has all its prime factors less than $$B$$ converges to $$\rho(u)$$. Sometime it is considered Dickman's function $$F$$ defined by $$F(\alpha)=\rho(1/\alpha)$$. In particular, Richard P. Brent's Some Integer Factorization Algorithms using Elliptic Curves uses $$\rho$$ for that $$F$$.

[#] gmp-ecm with option -pm1 implements Peter L. Montgomery and Alexander Kruppa's Improved Stage 2 to P±1 Factoring Algorithms, in proceedings of ANTS 2008. The first parameter is $$B_1$$. $$B_2$$ is chosen automatically if not specified as second parameter.

• The size of the second-largest factor is not independent from the size of the largest factor, so multiplying probabilities like that doesn't quite work very accurately. Estimating the semi-smoothness probability is a relatively well studied problem, and that's where gmp-ecm gets their approximations from. – Samuel Neves Oct 20 '19 at 13:17

estimate the probability that Pollard's p-1 factorization in its two-stages variant finds a factor of an RSA modulus product of $$k$$ random $$b$$-bit primes, as a function of the bounds $$B_1$$ and $$B_2$$ used

TLDR: GMP-ECM is fine. Modify it for $$b$$ bits rather than n decimal digits, and use $$k\;P$$ where $$P$$ is the low probability shown.

The way GMP-ECM ecm -v -pm1 coumputes $$P$$ is described by Alexander Kruppa in section 5.3.3 of his Doctor's thesis, as $$P=\hat\rho\biggl(\frac{\log(N_\text{eff})}{\log(B_1)},N_\text{eff}\biggr)+\sum_{B_1 where \begin{align}\hat\rho(u,x)&=\rho(u)-\gamma\frac{\rho(u-1)}{\log(x)}\\ \gamma&=\lim_{n\to\infty}\biggl(-\log(n)+\sum_{i=1}^n\frac 1 n\biggr) \approx0.5772156649\\ N_\text{eff}&=N\;e^{-\delta}\\ \delta&=\sum_{q\in\Bbb P}\frac{\log(q)}{(q-1)^2}\approx1.2269688057\\ N&=2^{b-0.5}\text{ is an approximation of the factor} \end{align}

In this, $$\displaystyle\hat\rho\biggl(\frac{\log(N)}{\log(B_1)},N\biggr)$$ estimates the probability that a random integer near $$N$$ has all its primes factors at most $$B_1$$; $$\displaystyle\sum_{B_1 estimates the probability that a random integer near $$N$$ has all its primes factors except the largest at most $$B_1$$ and has that largest factor in $$(B_1,B_2]$$; and $$\delta$$ is a correction to account for the fact that a random prime $$p$$ near $$N$$ has $$p-1$$ divisible by small prime $$q$$ with probability $$1/(q-1)$$, rather than $$1/q$$ for a random integer near $$N$$.

Note: that answer is still missing how to compute the sum for large $$B_2$$. GMP-ECM uses some indirect method when $$q>20000$$, which is always for parameters of practical interest in factoring an RSA modulus. For this reason, I make the question a Community Wiki.

As an illustration: ecm -pm1 35e9 uses $$B_1=35000000000$$, $$B_2=81866328103009612$$, and a 256-bits random prime is found with $$P\approx0.0000119$$. A factor is found in a 3-factors RSA 768-bit modulus product of random 256-bit prime factors for about one in 28000 moduli. Testing a moduli costs about 3.5 hour using one core and 48GiB RAM. Most of the RAM is used <25% of that time, so that one 4-core CPU 64GiB machine can test at a rate over one modulus per hour. If an adversary is content with factoring any one in like a thirty thousand such RSA moduli, Polard's p-1 is a plausible strategy to pull a factor (with expected cost 1170 CPU×days), followed by GNFS to factor a 512-bit composite. This beats GNFS alone, and I believe ECM. This is a setup where it makes sense to generate primes in a way that defeats Polard's p-1 and probably Williams's p+1 too.

Thanks to Paul Zimmermann and Alexander Kruppa for their guidance.

What's the cause of that discrepancy?

Main problem: the question's term $$\displaystyle\frac{\log_2(B_2)}{\log_2(B_1)}$$ is grossly too high. It was supposed to be a lower bound of the probability that the second-highest factor of a large random integer is less than $$B_1$$, knowing that the highest factor is less than $$B_2$$. That was baseless and wrong.