p = 0xfffffffffffffffffffffffffffffffffffffffffffffffffffffffefffffc2f

n = 0xfffffffffffffffffffffffffffffffebaaedce6af48a03bbfd25e8cd0364141

E = EllipticCurve(GF(p), [0, 7])

G = E.point( (0x79be667ef9dcbbac55a06295ce870b07029bfcdb2dce28d959f2815b16f81798,0x483ada7726a3c4655da4fbfc0e1108a8fd17b448a68554199c47d08ffb10d4b8))   # Base point

## Input
y = 0x483ada7726a3c4655da4fbfc0e1108a8fd17b448a68554199c47d08ffb10d4b8

## Field parameters
# Field modulus
p = 0xfffffffffffffffffffffffffffffffffffffffffffffffffffffffefffffc2f
# Cube root of 1
beta = 60197513588986302554485582024885075108884032450952339817679072026166228089408
 
## Actual code
xcubed = (y*y - 7) % p
print ("xcubed = 0x%x" % xcubed)

x = pow(xcubed, (p + 2) / 9, p)
print ("x1 = 0x%x" % x)
print ("x2 = 0x%x" % (x * beta % p))
print ("x3 = 0x%x" % (x * beta * beta % p))
x1=x * beta % p
x2=(x * beta * beta % p)
P=E.point((x,y))
P2=E.point((x1,y))
P3=E.point((x2,y))

assert G*78074008874160198520644763525212887401909906723592317393988542598630163514318 == P
assert G*37718080363155996902926221483475020450927657555482586988616620542887997980018 == P2