Por que a fatoração de produtos de primos próximos é muito mais lenta do que produtos de primos diferentes

Esta é uma questão puramente acadêmica, sem nenhuma consideração prática. Isto não é dever de casa, eu abandonei o ensino médio há muito tempo. Estou apenas curioso, e não consigo dormir bem sem saber o porquê.

Eu estava brincando com Python. Decidi fatorar números inteiros grandes e medir o tempo de execução das chamadas para cada entrada.

Usei vários números e descobri que alguns demoram muito mais para serem fatorados do que outros.

Decidi então investigar mais a fundo, escrevi rapidamente uma função de peneira de primos para gerar primos para teste. Descobri que um produto de um par de primos moderadamente grandes (dois primos de quatro dígitos) leva muito mais tempo para ser fatorado do que um produto de um primo muito grande (seis dígitos+) e um primo pequeno (<=três dígitos).

A princípio, pensei que minha primeira função simples para teste era ineficiente, e esse é realmente o caso. Então, escrevi uma segunda função que extraía números primos diretamente de uma lista pré-gerada de números primos. A segunda função era de fato mais eficiente, mas estranhamente ela exibe o mesmo padrão.

Aqui estão alguns números que usei:

13717421 == 3607 * 3803
13189903 == 3593 * 3671
56267023 == 7187 * 7829
65415743 == 8087 * 8089

12345679 == 37 * 333667
38760793 == 37 * 1047589
158202851 == 151 * 1047701
762312571 == 727 * 1048573

Código:

import numpy as np
from itertools import cycle

def factorize(n):
    factors = []
    while not n % 2:
        factors.append(2)
        n //= 2

    i = 3
    while i**2 <= n:
        while not n % i:
            factors.append(i)
            n //= i
        i += 2
    
    return factors if n == 1 else factors + [n]

TRIPLE = ((4, 2), (9, 6), (25, 10))
WHEEL = ( 4, 2, 4, 2, 4, 6, 2, 6 )
def prime_sieve(n):
    primes = np.ones(n + 1, dtype=bool)
    primes[:2] = False
    for square, double in TRIPLE:
        primes[square::double] = False
    
    wheel = cycle(WHEEL)
    k = 7
    while (square := k**2) <= n:
        if primes[k]:
            primes[square::2*k] = False
        
        k += next(wheel)
    
    return np.flatnonzero(primes)
    
PRIMES = list(map(int, prime_sieve(1048576)))
TEST_LIMIT = PRIMES[-1] ** 2

def factorize_sieve(n):
    if n > TEST_LIMIT:
        raise ValueError('Number too large')

    factors = []
    for p in PRIMES:
        if p**2 > n:
            break
        while not n % p:
            factors.append(p)
            n //= p
    
    return factors if n == 1 else factors + [n]

Resultado do teste:

In [2]: %timeit factorize(13717421)
279 μs ± 4.29 μs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

In [3]: %timeit factorize(12345679)
39.6 μs ± 749 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

In [4]: %timeit factorize_sieve(13717421)
64.1 μs ± 688 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

In [5]: %timeit factorize_sieve(12345679)
12.6 μs ± 146 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)

In [6]: %timeit factorize_sieve(13189903)
64.6 μs ± 964 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

In [7]: %timeit factorize_sieve(56267023)
117 μs ± 3.88 μs per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

In [8]: %timeit factorize_sieve(65415743)
130 μs ± 1.38 μs per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

In [9]: %timeit factorize_sieve(38760793)
21.1 μs ± 232 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

In [10]: %timeit factorize_sieve(158202851)
21.4 μs ± 385 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

In [11]: %timeit factorize_sieve(762312571)
22.1 μs ± 409 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

Como você pode ver claramente, a fatoração de dois primos médios leva, em média, muito mais tempo do que dois extremos. Por que é esse o caso?