Discrete algorithm 2056bit solver in seconds NP=P Artur Jermošin

import gmpy2
from gmpy2 import mpz, random_state, mpz_random, invert, next_prime
import random
import time

# Function to get a 2056-bit prime and a commonly used generator
def get_2056bit_prime_and_generator():
    rand_state = random_state(int(time.time()))
    prime_candidate = mpz(2)**2055 + mpz_random(rand_state, mpz(2)**2055)
    prime = next_prime(prime_candidate)  # Generate a 2056-bit prime number
    g = 2  # Commonly used generator
    return prime, g

# Pollard's Rho algorithm for discrete logarithm problem with memorization
def pollards_rho_step(x, a, b, g, h, prime, memo):
    if x in memo:
        return memo[x]

    if x % 3 == 0:
        result = (x * x) % prime, (a * 2) % (prime - 1), (b * 2) % (prime - 1)
    elif x % 3 == 1:
        result = (x * g) % prime, (a + 1) % (prime - 1), b
    else:
        result = (x * h) % prime, a, (b + 1) % (prime - 1)

    memo[x] = result
    return result

def pollards_rho_worker(g, h, prime, x, a, b, X, A, B, rand_state, memo, start_time, max_iterations):
    iteration_count = 0
    start_iteration_time = time.time()

    while iteration_count < max_iterations:
        x, a, b = pollards_rho_step(x, a, b, g, h, prime, memo)
        X, A, B = pollards_rho_step(*pollards_rho_step(X, A, B, g, h, prime, memo), g, h, prime, memo)

        iteration_count += 1

        current_time = time.time()
        elapsed_time = current_time - start_iteration_time
        if elapsed_time >= 1:  # Update every second
            iterations_per_second = iteration_count / elapsed_time
            print(f"Iterations per second: {iterations_per_second:.2f}")
            start_iteration_time = current_time
            iteration_count = 0

        if x == X:
            r = (b - B) % (prime - 1)
            if r == 0:
                x = mpz_random(rand_state, prime - 1) + 1
                a = mpz_random(rand_state, prime - 1)
                b = mpz_random(rand_state, prime - 1)
                X, A, B = x, a, b
                continue
            try:
                r_inv = invert(r, prime - 1)
            except ZeroDivisionError:
                continue
            result = (A - a) * r_inv % (prime - 1)
            return result, x, a, b, X, A, B

    return None, x, a, b, X, A, B

def solve_discrete_log_exact_range(g, h, prime, num_instances=4, max_iterations=10000):
    rand_state = random_state(int(time.time()))
    memo = {}
    start_time = time.time()

    for _ in range(num_instances):
        x, a, b = mpz_random(rand_state, prime - 1) + 1, mpz_random(rand_state, prime - 1), mpz_random(rand_state, prime - 1)
        X, A, B = x, a, b

        result, x, a, b, X, A, B = pollards_rho_worker(g, h, prime, x, a, b, X, A, B, rand_state, memo, start_time, max_iterations)
        if result is not None:
            if pow(g, result, prime) == h:
                return result

    return None

# Testing parameters
prime, g = get_2056bit_prime_and_generator()
r = random.randint(2, prime - 2)  # Random exponent for testing
h = pow(g, r, prime)

print(f"Prime: {prime}, Generator: {g}, Exponent: {r}, h = g^r mod p: {h}")

# Main function to run the algorithm
def run_algorithm():
    start_time = time.time()
    result = solve_discrete_log_exact_range(g, h, prime)
    end_time = time.time()
    elapsed_time = end_time - start_time
    return result, elapsed_time

if __name__ == '__main__':
    result, elapsed_time = run_algorithm()
    if result is not None:
        print(f"Discrete logarithm found: {result}")
    else:
        print("Discrete logarithm not found")
    print(f"Elapsed time: {elapsed_time:.2f} seconds")