<jupyter_start><jupyter_text>Analyse pseudo-primes<jupyter_code># Fermat strorng pseudoprimes base a
fs2_pp_mask = [fermat_strong_test(i, 2) for i in composites.tolist()]
fs2_pp = composites[fs2_pp_mask]
# This is A001262 sequence (https://oeis.org/A001262)
fs2_pp
# estimate probability that odd composite number is pseudoprime:
len(fs2_pp), len(fs2_pp) / len(composites)
# Lucas Selfridge strong pseudoprimes
# oesis sequence: https://oeis.org/A217255 
ls_pp_mask = [lucas_selfridge_test(i) for i in composites.tolist()]
ls_pp = composites[ls_pp_mask]
ls_pp
len(ls_pp), len(ls_pp) / len(comps)
# first intersection is of order 10^10000 !
set(ls_pp) & set(fs2_pp)
# Fermat pseudo primes with bases 3 and 5 
fs3_pp_mask = [fermat_strong_test(i, 3) for i in composites.tolist()]
fs3_pp = composites[fs3_pp_mask]
fs3_pp
len(fs3_pp)
fs5_pp_mask = [fermat_strong_test(i, 5) for i in composites.tolist()]
fs5_pp = composites[fs5_pp_mask]
fs5_pp
len(fs5_pp)
print(len(set(fs2_pp) & set(fs3_pp)))
set(fs2_pp) & set(fs3_pp)

# No intersection up to 10^9
set(fs2_pp) & set(fs3_pp) & set(fs5_pp)<jupyter_output><empty_output><jupyter_text>Fermat pseudo primes vs Fermat strong pseudo primes<jupyter_code>f2w_pp_mask = [fermat_test(i, 2) for i in composites.tolist()]
f2w_pp = composites[f2w_pp_mask]

len(f2w_pp)
f3w_pp_mask = [fermat_test(i, 2) for i in composites.tolist()]
f3w_pp = composites[f3w_pp_mask]

fw_pp = set(f2w_pp) & set(f3w_pp)
print(len(fw_pp))
fw_pp
2 ** 340 % 341 
d = 340 // 4
2 ** (d - 1) % d, d
inv2 = powmod(2, n - 2, n)

def div2(x, n):
    if n % 2 == 1:
        x += n
    return x // 2 % n

def lucas_double(u_k, v_k, k):
    # computes U_k, V_k -> U_2k, V_2k
    return u_k * v_k % n, (v_k * v_k + -2 * powmod(q, k, n)) % n

def lucas_sum(u_k, v_k, u_m, v_m):
    u_km = div2(u_k * v_m + u_m * v_k, n)
    v_km = (v_k * v_m + D * u_k * u_m)  % n
    return u_km, v_km
lucas_sum(u[1], v[1], u[2], v[2])
list(zip(range(len(u)), u, v))<jupyter_output><empty_output><jupyter_text>Analyze factorization methods<jupyter_code>from algorithms.number_theory import Pollard_rho_factor
from random import randint
from math import gcd, sqrt

from algorithms.number_theory import powmod

"""Count number of steps during factorization"""

def Pollard_rho_Floyd(n, x0=2, c=1):
    """Pollard's Rho method
    Attempts to find a divisor using pollard rho method
    with Floyd cycle finding algorithm.

    Returns either any prime factor of n or n itself.

    If n = p*q then algorithm exits approximately in O(p^(1/4))
    """

    def f(x):
        return (x * x + c) % n

    x, y, g = x0, x0, 1
    steps = 0
    while g == 1:
        steps += 1
        x = f(x)
        y = f(f(y))
        g = gcd(abs(x - y), n)
    return g, steps


def Pollard_rho_factor(n, check_prime=False):
    if check_prime and Miller_Rabin_test(n):
        return n

    nstesps, maxsteps = 0, 0
    while True:
        x0 = randint(2, n - 1)
        c = randint(1, n - 1)
        g, s = Pollard_rho_Floyd(n, x0, c)
        nstesps, maxsteps = nstesps + s, max(maxsteps, s)
        if g != n:
            return g, sqrt(sqrt(g)), nstesps, maxsteps

def Pollard_pm1(n, primes, max_B=1000000):
    """Pollard's p - 1 method
    Attempts to find some B-powersmooth factor of n

    Fcator p is B-powersmooth if p - 1 = p1^d1 * ... * pn^dn
    and max(p1, ..., pn) < B

    Parameters
    ----------
    n : int
        integer to factorize
    primes : list
        sorted list of primes < max_B
    max_B : int
        maximal powersmoothness of extracted factor
    """
    B = 10
    g = 1
    sm_steps, b_steps = 0, 0
    while B < max_B and g < n:
        a = randint(2, n - 2)
        g = gcd(a, n)
        b_steps += 1
        if g != 1:
            return g, sm_steps, b_steps
        for p in primes:
            if p >= B:
                break
            sm_steps += 1
            pd = 1  # p^d
            while pd * p <= B:
                pd *= p
            a = powmod(a, pd, n)
            g = gcd(a - 1, n)
            if g != 1 and g != n:
                return g, sm_steps, b_steps
        B *= 2
    return 1, sm_steps, b_steps
Pollard_rho_factor(27)
comps = composites[1000000:2000000]
factors = [Pollard_pm1(n, primes) for n in comps]
ok = list(i % d[0] == 0 and i != d[0] for i, d in zip(comps, factors))
np.all(ok) # all factorized correctly
factors2 = [Pollard_rho_factor(n) for n in comps]
ok = list(i % d[0] == 0 and i != d[0] for i, d in zip(comps, factors2))
ok = np.all(ok)
ok
# number of steps for Pollard p-1
sns.distplot(np.array(factors)[:, 1])
np.array(factors)[:, 1].max()
# number of steps for Pollard rho
sns.distplot(np.array(factors2)[:, 2])
np.array(factors2)[:, 2].max()
# number of steps for Pollard rho
sns.distplot(np.array(factors2)[:, 3])
np.array(factors2)[:, 3].max()
Pollard_rho_factor(primes[100500]*primes[333111])
Pollard_rho_factor(primes[100501]*primes[233111])
Pollard_rho_factor(primes[130501]*primes[233141])<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12949206986548886, 'fraction_numerical': 0.043164023288496285, 'mean_word_length': 2.9919871794871793, 'pattern_counts': {'":': 0, '<': 18, '<?xml version=': 0, '>': 15, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '121170', 'license': 'MIT', 'n_tokens_mistral': 2178, 'n_tokens_neox': 2018, 'n_words': 550, 'path': 'benchmark/primes.ipynb', 'repo_name': 'vadimadr/python-algorithms', 'size': '4981'}