# Pollard's p − 1 algorithm
## Overview
* Thuật toán này được John Pollard vào năm 1974
* Thuật toán này mạnh trong việc phân tích thừa số nguyên tố của một **"smooth number"** là số khi ta phân tích ra thành các số nguyên tố nhỏ.

## Base Concepts
Cho $n$ là một hợp số có thừa số nguyên tố là $p$. Theo định lý **Fermat's little**. Với mọi số nguyên $a$ nguyên tố cùng nhau với $n$. Ta có:
$$a^{K(p-1)} \equiv 1 \pmod p \Rightarrow gcd(a^{K(p-1)}-1, n) = p$$
* Ý tưởng làm cho số mũ là một bội của $(p-1)$ bằng cách biến nó thành một số với rất nhiều thừa số nguyên tố. Nói chung chúng ta lấy tích của tất cả các lũy thừa nguyên tố nhỏ hơn một giới hạn B nào đó.
* Bắt đầu với một số nguyên $x$ và liên tục thay thế nó bằng $x^w \bmod n$ khi $w$ chạy qua những lũy thừa nguyên tố đó.
* Kiểm tra ở mỗi giai đoạn xem $gcd(x-1, n)$ nếu như nó khác 1 thì ta nhận đó là một thừa số nguyên tố của $n$.

## Algorithm
* **Input:** một số nguyên n
* **Output:** danh sách các thừa số nguyên tố của n hoặc thất bại.
* Chọn một giới hạn $B$ (các thừa số nguyên tố sẽ nhỏ hơn B).
* Xác định $M = factorial(B)$
* Chọn một số nguyên dương $a$ nguyên tố cùng nhau với $n$, nếu như $n$ lẻ thì ta luôn chọn $a = 2$.
* Tính toán $gcd(a^{M} -1 , n)$
* Nếu:
* $1<g<n$: thì $g$ là một thừa số nguyên tố.
* $g = 1$ chọn $B$ lớn hơn và thử lại.
* $g = n$ chọn $B$ nhỏ hơn và thử lại.

Dưới đây là một đoạn script khởi tạo thuật toán:

from gmpy2 import fac
from math import gcd
from Crypto.Util.number import *
from sympy import true

n = 1403421

a = 2
B = 100

while True:
    if isPrime(n):
        print(f"[+] p factor : {n}")
        print("[+] Finished")
        break
    M = fac(B)
    tmp = pow(a, M, n) - 1
    g = gcd(tmp, n)
    #print(B)
    if g == 1:
        B += 1
    elif g == n:
        B -= 1
    else:
        print(f"[+] p factor : {g}")
        n = n//g
        

[+] p factor : 2073
[+] p factor : 677
[+] Finished

Như vậy với việc xác định được bound thì ta có thể dễ dàng factor được "smooth number".

## Example
Giờ ta sẽ đến với một challenge Crypto như sau:
Ta có source code:

#!/usr/bin/python

from binascii import hexlify
from gmpy2 import *
import math
import os
import sys

if sys.version_info < (3, 9):
    math.gcd = gcd
    math.lcm = lcm

_DEBUG = False

FLAG  = open('flag.txt').read().strip()
FLAG  = mpz(hexlify(FLAG.encode()), 16)
SEED  = mpz(hexlify(os.urandom(32)).decode(), 16)
STATE = random_state(SEED)

def get_prime(state, bits):
    return next_prime(mpz_urandomb(state, bits) | (1 << (bits - 1)))

def get_smooth_prime(state, bits, smoothness=16):
    p = mpz(2)
    p_factors = [p]
    while p.bit_length() < bits - 2 * smoothness:
        factor = get_prime(state, smoothness)
        p_factors.append(factor)
        p *= factor

    bitcnt = (bits - p.bit_length()) // 2

    while True:
        prime1 = get_prime(state, bitcnt)
        prime2 = get_prime(state, bitcnt)
        tmpp = p * prime1 * prime2
        if tmpp.bit_length() < bits:
            bitcnt += 1
            continue
        if tmpp.bit_length() > bits:
            bitcnt -= 1
            continue
        if is_prime(tmpp + 1):
            p_factors.append(prime1)
            p_factors.append(prime2)
            p = tmpp + 1
            break

    p_factors.sort()

    return (p, p_factors)

e = 0x10001

while True:
    p, p_factors = get_smooth_prime(STATE, 1024, 16)
    if len(p_factors) != len(set(p_factors)):
        continue
    # Smoothness should be different or some might encounter issues.
    q, q_factors = get_smooth_prime(STATE, 1024, 17)
    if len(q_factors) != len(set(q_factors)):
        continue
    factors = p_factors + q_factors
    if e not in factors:
        break

if _DEBUG:
    import sys
    sys.stderr.write(f'p = {p.digits(16)}\n\n')
    sys.stderr.write(f'p_factors = [\n')
    for factor in p_factors:
        sys.stderr.write(f'    {factor.digits(16)},\n')
    sys.stderr.write(f']\n\n')

    sys.stderr.write(f'q = {q.digits(16)}\n\n')
    sys.stderr.write(f'q_factors = [\n')
    for factor in q_factors:
        sys.stderr.write(f'    {factor.digits(16)},\n')
    sys.stderr.write(f']\n\n')

n = p * q

m = math.lcm(p - 1, q - 1)
d = pow(e, -1, m)

c = pow(FLAG, e, n)

print(f'n = {n.digits(16)}')
print(f'c = {c.digits(16)}')

**output.txt**

n = 4f7aa864f662a42a92220e372f5ff25a142aef26106a0dbdf573a66594966ac5dd03848745bb6a80402cad7ac6f2bf93f9ed840edd9c157dfd5d265ce2403e155a29666df8f9b98167ad2452e5a63fd0b7b14ffe966db60c6e2c65b0f602f5c22eb030c0335187759909abd4df622118c23463bcc42650e0a7761257452bf40069ca50dbe0c922d8823a9dcc4231b3952d31d1e977cb520528c6a450405f2a2ee6134db8c61ceb4478a647b0469712cc4f3d1369ef3dfd3d876a2c77bac5a149ccf3723a6e8c3ba1deb0675f25def8da9de2b3ac8b3e38d5ac5c9736b9af087b3fc53450136428e07d58fbc00f6609a4cc14eb0a13a7e76056a241256e03e95d
c = e41a61908eb48b85dc78975c288e62a271b1f237fdc958162727d2930b9af850e908137655c5955a078ff1aa63f5509fbaf79d179d24d209a061c36e0709437b8d2641f41d354bdea062084ea3be8637ed1c4bd8cf63d16c942976dd9d6188fc5e419afae17493d7cdb93d84052637d15e7fa1f852f4f5d786c86bfd024df0dfcf8431e7230cfbbce76a1835b178020ef839af42c377706918a50aac56f79285d743f4a177425eb00eaeb2bebe99343911ab653fe64bb61e140153b113f8554fe29561756fafc7460683d59dd3ee50eb48b718443b9f49e663b6dd02b0a15297468ec30a4f487e328103cdbc59d1d66fc4f03ef75ae45d6ce2035fdfaeb86b7

Đọc source code ở trên ta nhận thấy các **"smooth prime"** được tạo từ các số nguyên tố 16 bit và 17 bit nên ta dễ dàng tìm được Bound của nó là $B_p = 2^{16} + 1$ và $B_q = 2^{17} + 1$ Ta sẽ sử dụng thuật toán **Pollard's p-1** để factor N.
Tại vì N chỉ có 2 thừa số nguyên tố nên ta chỉ cần sử dụng một Bound thì đã factor ra được $p$ và $q$.


Khi đã factor được thì mọi chuyện trở nên dễ dàng rồi. :monkey:

---END---