import random
import math

def inverse(x, p):
    """
    Calculate the modular inverse of x ( mod p )    
    """
    inv1 = 1
    inv2 = 0
    n=1
    while p != 1 and p!=0:        
        quotient = x // p        
        inv1, inv2 = inv2, inv1 - inv2 * quotient
        x, p = p, x % p        
        n = n+1
    
    return inv2

def dblpt(pt, p):
    """
    Calculate pt+pt = 2*pt
    """
    if pt is None:
        return None
    (x,y)= pt
    if y==0:
        return None 
    
    slope= 3*pow(x,2,p)*pow(2*y,p-2,p)  
    xsum= pow(slope,2,p)-2*x     
    ysum= slope*(x-xsum)-y 
    
    return (xsum%p, ysum%p)

def addpt(p1,p2, p):
    """
    Calculate p1+p2
    """
    if p1 is None or p2 is None:
        return None
    (x1,y1)= p1
    (x2,y2)= p2
    if x1==x2:
        return dblpt(p1, p)        
        
    # calculate (y1-y2)/(x1-x2)  modulus p    
    slope=(y1-y2)*pow(x1-x2,p-2,p)
    xsum= pow(slope,2,p)-(x1+x2)  
    ysum= slope*(x1-xsum)-y1 
    
    return (xsum%p, ysum%p)

def ptmul(pt,a, p):
    """
    Calculate pt*a
    """
    scale= pt    
    acc=None
    while a:        
        if a&1:
            if acc is None:
                acc= scale                 
            else:     
                acc= addpt(acc,scale, p)             
        scale= dblpt(scale, p)
        a >>= 1
        
    return acc

def ptdiv(pt,a,p,n):  
    """
    Calculate pt/a
    """   
    divpt=inverse(a, n)%n 
    
    return ptmul(pt, divpt, p)


def getuncompressedpub(compressed_key):
    """
    returns uncompressed public key
    """
    y_parity = int(compressed_key[:2]) - 2    
    x = int(compressed_key[2:], 16)
    a = (pow(x, 3, p) + 7) % p
    y = pow(a, (p+1)//4, p)    
    if y % 2 != y_parity:
        y = -y % p   
        
    return (x,y)

def compresspub(uncompressed_key):
    """
    returns compressed public key
    """
    (x,y)=uncompressed_key
    y_parity = y&1 
    head='02'
    if y_parity ==1:
        head='03'    
    compressed_key = head+'{:064x}'.format(x) 
    
    return compressed_key


#secp256k1 constants
Gx=0x79BE667EF9DCBBAC55A06295CE870B07029BFCDB2DCE28D959F2815B16F81798
Gy=0x483ADA7726A3C4655DA4FBFC0E1108A8FD17B448A68554199C47D08FFB10D4B8
n=0xfffffffffffffffffffffffffffffffebaaedce6af48a03bbfd25e8cd0364141
p = 2**256 - 2**32 - 977
g= (Gx,Gy)


rangePower=129
rb=2**rangePower
re=2**(rangePower+1)-1
print ("Bit range 2^",math.log2(re-rb))
print ("Begin > 0x%x"%rb,">DEC ",rb)
print ("End   > 0x%x"%re,">DEC ",re)

#MANUAL

compressed_key='03633CBE3EC02B9401C5EFFA144C5B4D22F87940259634858FC7E59B1C09937852'
point=getuncompressedpub(compressed_key)
divisor = pow(2,20)
print("")

rbdiv = rb//divisor
rediv = re//divisor
print ("Div Begin > 0x%x"%rbdiv,">DEC ",rbdiv)
print ("Div End   > 0x%x"%rediv,">DEC ",rediv)
print ("Div Bit range 2^",math.log2(rediv-rbdiv))
newpub=ptdiv(point,divisor,p,n)
(partGx,partGy)=ptdiv(g,divisor,p,n)
idx=random.randrange(0,divisor-1)
#Randomly generate position in range 0..divisor-1
print("pos > 0x%x"%idx)
if idx<divisor:
    if idx==0:
        (searchpubx,searchpuby)=newpub
        (fracX, fracY)=(partGx,partGy)
    else:   
        (fracX, fracY)=ptmul((partGx,partGy),idx,p)        
        (searchpubx,searchpuby) = addpt(newpub,(fracX,p-fracY), p)
        (fracX, fracY) = addpt((fracX, fracY),(partGx,partGy), p)   
    print("searchpub > ",compresspub((searchpubx,searchpuby)))
else:
    print("idx>=divisor")