Blum-Blum-Shub

bbs_crack.py

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

from sympy import pollard_rho
from sympy.core.numbers import igcd
from sympy.ntheory import sqrt_mod, nthroot_mod, isprime, factorint
from sympy.ntheory.modular import crt

with open('bbs.txt', 'r') as f:
    xs = [int(x.rstrip()) for x in f.readlines()]

# найдем предположительное n
ys = [xs[i - 1] ** 2 - xs[i] for i in range(1, len(xs))]
n = igcd(*ys)
print(n)


# продолжение вперед
def next_x(x):
    return pow(x, 2, n)


# проверим, что последовательность получена с этим n
for i in range(1, len(xs)):
    assert xs[i] == next_x(xs[i - 1])

# чтобы делать jump-ahead необходимо знать разложение n на p и q

# p = pollard_rho(n)  # took a long time -- for about 17 minutes
# print(p)

p = 1975784453805923
q = n // p
assert isprime(p)
assert isprime(q)

print(p)
print(q)


def jump_ahead(x, m):
    power = pow(2, m, (p - 1) * (q - 1) // 2)
    return pow(x, power, n)


assert jump_ahead(xs[0], 5) == xs[5]


# теперь научимся продолжать последовательность назад
def prev_x(x):
    # извлечем корень в Z/p и Z/q
    xp = [y for y in sqrt_mod(x, p, True) if sqrt_mod(y, p, True)]
    assert len(xp) == 1
    xp = xp[0]
    xq = [y for y in sqrt_mod(x, q, True) if sqrt_mod(y, q, True)]
    assert len(xq) == 1
    xq = xq[0]
    # найдем ответ посредством китайской теоремы об остатках
    return crt([p, q], [xp, xq])[0]


for i in range(1, len(xs)):
    assert xs[i - 1] == prev_x(xs[i])


# научимся делать jump-ahead назад
def jump_ahead_back(x, m):
    power = pow(2, m, (p - 1) * (q - 1) // 2)
    # извлечем корень в Z/p и Z/q
    xp = [y for y in nthroot_mod(x, power, p, True) if sqrt_mod(y, p, True)]
    assert len(xp) == 1
    xp = xp[0]
    xq = [y for y in nthroot_mod(x, power, q, True) if sqrt_mod(y, q, True)]
    assert len(xq) == 1
    xq = xq[0]
    # найдем ответ посредством китайской теоремы об остатках
    return crt([p, q], [xp, xq])[0]


# найдем период последовательности

# число, кратное наименьшему периоду

def is_period(period):
    return jump_ahead(xs[0], period) == xs[0]


t = ((p - 1) // 2 - 1) * ((q - 1) // 2 - 1)

assert is_period(t)

# найдем факторизацию t
fact = factorint(t)

# будем делить на каждый делитель, пока полученное число остается периодом
for d, e in fact.items():
    for _ in range(e):
        if not is_period(t // d):
            break
        t //= d

print(t)