Topic: Theoretical minimum # of logic operations to perform privkey->pubkey(secp256k1) (Read 1279 times)

#code by hakabit supraspecially for gmaxwell
#algorithm source: https://dspace.cvut.cz/bitstream/handle/10467/77227/F3-DP-2018-Zhumaniezov-Alisher-priloha-Testing.py

import random

def bin_ext_gcd(a, b):

    global cc
    cc = cc + 1

    if a == 0:
        #print(a, b, 0, 1)
        return [0, 1]

    if b == 0:
        #print(a, b, 1, 0)
        return [1, 0]

    if a == b:
        #print(a, b, 1, 0)
        return [1, 0]

    if (a % 2 == 0) & (b % 2 == 0):
        x, y = bin_ext_gcd(a // 2, b // 2)
        #print(a, b, x, y)
        return [x, y]

    if a % 2 == 0:
        x, y = bin_ext_gcd(a // 2, b)
        if x % 2 == 0:
            x, y = x // 2, y
        else:
            x, y = (x + b) // 2, y - a // 2
        #print(a, b, x, y)
        return [x, y]

    if b % 2 == 0:
        x, y = bin_ext_gcd(a, b // 2)
        if y % 2 == 0:
            x, y = x, y // 2
        else:
            x, y = x - b // 2, (y + a) // 2
        #print(a, b, x, y)
        return [x, y]

    if a > b:
        x, y = bin_ext_gcd((a - b) // 2, b)
        if x % 2 == 0:
            x, y = x // 2, y - x // 2
        else:
            x, y = (x + b) // 2, y - (a - b) // 2 - (x + b) // 2
        #print(a, b, x, y)
        return [x, y]

    if a < b:
        x, y = bin_ext_gcd(a, (b - a) // 2)
        if y % 2 == 0:
            x, y = x - y // 2, y // 2
        else:
            x, y = x - (b - a) // 2 - (y + a) // 2, (y + a) // 2
        #print(a, b, x, y)
        return [x, y]

cc = 0

p = 2**256 - 2**32 - 977


#cc=0
#a=12
#res=bin_ext_gcd(p,a)
#print ("initial test, should be in=12, out=48246703848881748093154577086953294938862493610683568349773993336628681113193")
#print (res,cc)



print ("test# calls# in [out, mod symm out]")
for i in range(100):

    #print(i)
    ii = random.getrandbits(256)
    ii = ii % p

    if (ii < 2):
        ii = ii + 2

    cc = 0
    res = bin_ext_gcd(p, ii)
    print(i,cc,ii,res)

#end

#code by hakabit supraspecially for gmaxwell
#algorithm source: http://cacr.uwaterloo.ca/hac/about/chap14.pdf, page 608, section "14.61 Algorithm Binary extended gcd algorithm"

import random
#print(random.getrandbits(256))

def ext_binary_gcd(x, y):

    g = 1

    iterx = 0
    iter1 = 0
    iter2 = 0
    iter3 = 0

    while (x % 2 == 0) and (y % 2 == 0) :    # yeah, funny cycle when one number is prime :-)

        iter1 = iter1 + 1

        x = x // 2
        y = y // 2
        g = 2 * g

    u = x
    v = y
    A = 1
    B = 0
    C = 0
    D = 1

    while (u != 0) :

        iterx = iterx + 1

        #step4
        while (u % 2 == 0) :

            iter2 = iter2 + 1

            u = u // 2

            if (A % 2 == 0) and (B % 2 == 0) :

                A = A // 2
                B = B // 2

            else :

                A = (A + y) // 2
                B = (B - x) // 2

        while (v % 2 == 0) :

            iter3 = iter3 + 1

            v = v // 2

            if (C % 2 == 0) and (D % 2 == 0) :

                C = C // 2
                D = D // 2

            else :

                C = (C + y) // 2
                D = (D - x) // 2

        if (u >= v) :

            u = u - v
            A = A - C
            B = B - D

        else :

            v = v - u
            C = C - A
            D = D - B

    #a = C
    #b = D

    pp = 2 ** 256 - 2 ** 32 - 977
    b = D % pp   #BLM-rig

    #return (a, b, g*v)
    return (b, iterx, iter1, iter2, iter3, iter2 + iter3)

p = 2**256 - 2**32 - 977

#a=12
#res=ext_binary_gcd(p,a)
#print ("initial test, should be in=12, out=48246703848881748093154577086953294938862493610683568349773993336628681113193")
#print ("in (out, iterbig, iter1, iter2, iter3): ",a,res)



print("test# in (out, iterx, iter1, iter2, iter3, itersum): ")
for i in range(1000):

    #print(i)
    ii = random.getrandbits(256)
    ii = ii % p

    if (ii < 2):
        ii = ii + 2

    res = ext_binary_gcd(p, ii)
    print(i, ii, res)

#end

def inv(a, m):
    if a < 0 or m <= a: a = a % m
    c, d = a, m
    uc, vc, ud, vd = 1, 0, 0, 1
    while c != 0:
        q, c, d = divmod(d, c) + (c,)
        uc, vc, ud, vd = ud - q*uc, vd - q*vc, uc, vc
    if ud > 0: return ud
    else: return ud + m

def to_jacobian((Xp, Yp)):
    """
    Convert point to Jacobian coordinates

    :param (Xp,Yp,Zp): First Point you want to add
    :return: Point in Jacobian coordinates
    """
    return (Xp, Yp, 1)

def from_jacobian((Xp, Yp, Zp), P):
    """
    Convert point back from Jacobian coordinates

    :param (Xp,Yp,Zp): First Point you want to add
    :param P: Prime number in the module of the equation Y^2 = X^3 + A*X + B (mod p)
    :return: Point in default coordinates
    """
    z = inv(Zp, P)
    return ((Xp * z**2) % P, (Yp * z**3) % P)

def jacobian_add((Xp, Yp, Zp), (Xq, Yq, Zq), A, P):
    """
    Add two points in elliptic curves

    :param (Xp,Yp,Zp): First Point you want to add
    :param (Xq,Yq,Zq): Second Point you want to add
    :param P: Prime number in the module of the equation Y^2 = X^3 + A*X + B (mod p)
    :param A: Coefficient of the first-order term of the equation Y^2 = X^3 + A*X + B (mod p)
    :return: Point that represents the sum of First and Second Point
    """
    if not Yp:
        return (Xq, Yq, Zq)
    if not Yq:
        return (Xp, Yp, Zp)
    U1 = (Xp * Zq ** 2) % P
    U2 = (Xq * Zp ** 2) % P
    S1 = (Yp * Zq ** 3) % P
    S2 = (Yq * Zp ** 3) % P
    if U1 == U2:
        if S1 != S2:
            return (0, 0, 1)
        return jacobian_double((Xp, Yp, Zp), A, P)
    H = U2 - U1
    R = S2 - S1
    H2 = (H * H) % P
    H3 = (H * H2) % P
    U1H2 = (U1 * H2) % P
    nx = (R ** 2 - H3 - 2 * U1H2) % P
    ny = (R * (U1H2 - nx) - S1 * H3) % P
    nz = (H * Zp * Zq) % P
    return (nx, ny, nz)

def fast_add(a,b, A, P):
    return from_jacobian(jacobian_add(to_jacobian(a), to_jacobian(b), A, P), P)

P = 2**256 - 2**32 - 977
A=0

a=39006303077722201472019215118849942292522572644014512222183524570138876820125,67456145645462716033455438706152384917972612258004579542185540136844047420439

b=4032983015753143990395647783770666587927265353624430905763286836981504199392,44353125519324157186344456159742269880631179110473143840214086765587351124293

g2=fast_add(a,b,A,P)

print g2

Gx = 0x79BE667EF9DCBBAC55A06295CE870B07029BFCDB2DCE28D959F2815B16F81798

Gy= 0x483ADA7726A3C4655DA4FBFC0E1108A8FD17B448A68554199C47D08FFB10D4B8

P = 0xFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFEFFFFFC2F
  
#N=2**256 - 2**32 - 977
 

def powm(n, e, mod)
    fail ArgumentError, 'negative exponent' if e < 0
    prod = 1
    base = n % mod
    until e.zero?
      prod = (prod * base) % mod if e.odd? #e.odd can give error on newest ruby versions
      e >>= 1
      base = (base * base) % mod
    end
    prod
  end
   
    def add(x1,y1,x2,y2)   
        px, py = x1,y1
        qx, qy = x2,y2
       
   if px == qx && py==qy then       
        lam = (3 * px * px) * powm(2 * py, P - 2, P)     #can be ignored in estimations
   else
        lam = (qy - py) * powm(qx - px, P - 2, P)
   end
   
        rx = lam*lam - px - qx
        ry = lam * (px - rx) - py

      [rx % P, ry % P]
 end
   
        
privkey=0x18E14A7B6A307F426A94F8114701E7C8E774E7F9A47E2C2035DB29A206321725

          tx=Gx     
          ty=Gy
    
       rx=0
       ry=0    
    
for i in 0..255

             if privkey[i]==1 then
          
               if rx == 0 && ry==0  then                             
                  rx,ry=tx,ty
               else
                   rx,ry = add(rx, ry,tx,ty)
               end         
                           
             end
             
            tx,ty = add(tx, ty,tx,ty)        
      
   end
    
    
puts rx,ry