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

Question

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.

Answer 1

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<q\le B_2\\\quad q\in\Bbb P}\hat\rho\biggl(\frac{\log(N_\text{eff}/q)}{\log(B_1)},N_\text{eff}/q\biggr)/q$$ 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<q\le B_2\\\quad q\in\Bbb P}\hat\rho\biggl(\frac{\log(N/q)}{\log(B_1)},N/q\biggr)/q$ 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.