pow-pow-solve.sage

from Crypto.Util.number import long_to_bytes, bytes_to_long
from hashlib import shake_128

# Utils & Arguments
def encode(x):
    return long_to_bytes(x, 256)

def H(g, h):
    return bytes_to_long(shake_128(encode(g) + encode(h)).digest(int(16)))

n = 20074101780713298951367849314432888633773623313581383958340657712957528608477224442447399304097982275265964617977606201420081032385652568115725040380313222774171370125703969133604447919703501504195888334206768326954381888791131225892711285554500110819805341162853758749175453772245517325336595415720377917329666450107985559621304660076416581922028713790707525012913070125689846995284918584915707916379799155552809425539923382805068274756229445925422423454529793137902298882217687068140134176878260114155151600296131482555007946797335161587991634886136340126626884686247248183040026945030563390945544619566286476584591

# Find prime in [2, 10000000)
MAXV = 10000000
is_prime = [True for _ in range(MAXV)]
is_prime[0] = is_prime[1] = False

for i in range(2, MAXV):
    if is_prime[i]:
        for j in range(i * i, MAXV, i):
            is_prime[j] = False      
primes = [ i for i in range(MAXV) if is_prime[i]]

# Let M be 2 * 3 * 5 * ... * 9999991 (primes in [2, 10000000))
M = 1
for i in primes:
    M *= i

# Find b that
#     b  =  2^(2^k)
#     h  =  1
#     g  =  b^M
#     m  =  H(g, h)
# s.t.
#     m | M
#
# then we can get pi by
#     pi =  b^(-(M // m) * r)
# where r = 2^(2^64) % m
h = 1
b = Mod(2, n)
g = b^M

idx = 0
while True:
    idx += 1
    if idx % 10000 == 0:
        print(idx)
    m = H(int(g), int(h))
    if M % m == 0:
        break
    b *= b
    g *= g

# Get pi
assert(b^M == g)
m = H(int(g), int(h))
assert(M % m == 0)
r = pow(2, 2^64, m)
pi = (b^(-M // m))^(r)
assert((pi^m) * (g^r) == 1)

# Print g, h, pi
print(f'g = {g}')
print(f'h = {h}')
print(f'pi = {pi}')