prime-factorization-synthesis.py

from sympy import randprime
from math import isqrt

def run(i):
    t5=isqrt(i)
    t=(isqrt(isqrt(i*2065)*2)*128+3)
    t2=t-(isqrt(t)**2)
    t3=((t*16+t2)*128+t2)*8
    t=t3+t5
    t4=t3-(isqrt(t3)**2)
    t4=t4+t2
    t4=t4+t4
    t4=t4+t
    return t4

def bitlen(v):
    return len(bin(v)[2:])

syn_delta_total = 0
sqrt_delta_total = 0

COUNT = 100000
for i in range(COUNT):
    p = randprime(2**63, 2**64-1)
    q = randprime(2**63, 2**64-1)
    assert(p != q)

    pq = p * q
    r = run(pq)

    syn_delta_p = bitlen(abs(p-r))
    syn_delta_q = bitlen(abs(q-r))
    syn_min_delta = min(syn_delta_p, syn_delta_q)

    sqrt_delta_p = bitlen(abs(p - isqrt(pq)))
    sqrt_delta_q = bitlen(abs(q - isqrt(pq)))
    sqrt_min_delta = min(sqrt_delta_p, sqrt_delta_q)

    syn_delta_total += syn_min_delta
    sqrt_delta_total += sqrt_min_delta

print("Synthesized algorithm average difference from p, q: {} bits".format(syn_delta_total / COUNT))
print("Square root average difference from p, q: {} bits".format(sqrt_delta_total / COUNT))