#=========================================================# # Rho method to solve n from R=n*Q # # Q,R : Points on the ECC, n : integer # #---------------------------------------------------------# # copy right : Ushiro Yasunori (ISCPC) # # date : 2020/03/16 # #=========================================================# # ecc.py # # set_g(p,r,a,b) : g_p,g_r,g_a,g_b # # inv_mod(a,p) : 1/a (mod p) # # ecc_Add1(A,B) : A+B on ECC # # ecc_Add2(A) : 2*A on ECC # # ecc_Mul(n,A) : n*A on ECC # # y^2=x^3+g_a*x+g_b (mod p) on ECC # # x,y,g_a,g_b:integer, p,r:prime, r:order # #=========================================================# inport numpy as np inport math,sys,time #======= long table size constant ========================# st_bit = 23 #size of t_ptr: st_tabl=2^st_bit #======= global constant =================================# global fin,fout # input,output file global g_p,g_r,g_a,g_b # values on ECC global t_r,s_r # original r,s_r=t_r/g_r global st_size,st_n # table size and number global p_Q,p_R # ECC point: R=n*Q, n:unknown global p_T,p_alp,p_bet # Rho step variable global r_tbl,n_tbl # r-factor global n_sol,n_lp # n and loop list #----- long table for eqaul search -----------------------# st_size = 2 << (st_bit - 1) r_tbl = [0]*30 # r_tbl list n_sol = [0]*30 # R=n*Q, solve n list n_lp = [0]*30 # Rho loop list st_Tab = [] # table for search st_ptr = np.zeros((st_size),dtype=np.int32) #pointer table #=========================================================# # inv_mod(a,p) = 1/a (mod p) # # a:integer p:prime number # #---------------------------------------------------------# # copy right : Ushiro Yasunori (ISCPC) # # date : 2020/03/17 # #=========================================================# def inv_mod(a, p): rn,ro, sn,so = p,a, 0,1 while rn != 0: q = ro//rn un,vn = rn,sn rn = ro - q*rn sn = so - q*sn ro,so = un,vn if (so < 0): so += p return so #=========================================================# # C : ecc_Add2(A) = 2*A on ECC # # C[0] = e^2 - 2*A[0] (mod p) # # C[1] = e*(A[0]-C[0]) - A[1] (mod p) # # d=1/(2*A[1]), e=d*(3*A[0]^2+g_a) (mod p) # #---------------------------------------------------------# # copy right : Ushiro Yasunori (ISCPC) # # date : 2020/01/30 # #=========================================================# def ecc_Add2(A): a2 = 2*A[1] if a2 == 0: return A #special case # compute d,e d = inv_mod(a2, g_p) #d=1/a2 (mod p) h = (3*A[0]*A[0] + g_a) % g_p e = (d*h) % g_p # compute C=[c0,c1] cx = (e*e - 2*A[0]) % g_p cy = (e*(A[0]-cx) - A[1]) % g_p C = [cx, cy] return C #=========================================================# # C : ecc_Add1(A,B) = A+B on ECC # # C[0] = e^2 - A[0] - B[0] (mod p) # # C[1] = e*(A[0]-C[0]) - A[1] (mod p) # # d=1/(B[0]-A[0]), e=d*(B[1]-A[1]) (mod p) # #---------------------------------------------------------# # copy right : Ushiro Yasunori (ISCPC) # # date : 2020/03/19 # #=========================================================# def ecc_Add1(A, B): # special check if A == [0,0]: return B if B == [0,0]: return A if A[0] == B[0]: if A[1] == B[1]: if A[1] == 0: C = [0, 0] else: C = ecc_Add2(A) return C elif (A[1]+B[1]) == g_p: C = [0, 0] return C # compute d,e d = inv_mod(B[0]-A[0], g_p) #d=1/(Bx-Ax) (mod p) e = ( d*(B[1] - A[1]) ) % g_p # compute C=[c0,c1] cx = (e*e - A[0] - B[0]) % g_p cy = (e*(A[0] - cx) - A[1]) % g_p C = [cx, cy] return C #=========================================================# # C : ecc_Mul(n,A) = n*A on ECC # # n: integer, A,C: point in ECC # #---------------------------------------------------------# # copy right : Ushiro Yasunori (ISCPC) # # date : 2020/01/30 # #=========================================================# def ecc_Mul(n, A): W, nw = A, n C = [0, 0] while nw > 0: if (nw&1) == 1: C = ecc_Add1(C, W) W = ecc_Add2(W) nw = nw >> 1 return C #=========================================================# # file read function # # return: integer value # #---------------------------------------------------------# # copy right : Ushiro Yasunori (ISCPC) # # date : 2020/03/19 # #=========================================================# def val_read(name): data = fin.readline() val = data.split() fout.write(name+val[0]+"\n") rv = int(val[0]) return rv #=========================================================# # file open and read data # # g_p,t_r,g_a,g_b : values on ECC # # p_Q,p_R : point on ECC # # return: nv=number of r-factor # #---------------------------------------------------------# # copy right : Ushiro Yasunori (ISCPC) # # date : 2020/03/21 # #=========================================================# def ecc_read(): global g_p,t_r,g_a,g_b # values on ECC global t_Q,t_R,r_tbl # ECC point,r-factor global fin,fout # input,output file t_Q = [0,0] # list t_R = [0,0] # list # file open fin = open("pin.txt", "r") fout = open("pout.txt", "w") fout.write("Rho method to solve n from R=n*Q by python\n") fout.write("----- Input information ---------------\n") # read g_p,t_r,g_a,g_b g_p = val_read(" p=") t_r = val_read(" r=") g_a = val_read(" a=") g_b = val_read(" b=") # read t_Q t_Q[0] = val_read(" Qx=") t_Q[1] = val_read(" Qy=") # read t_Q t_R[0] = val_read(" Rx=") t_R[1] = val_read(" Ry=") # read r-factor data = fin.readline() val = data.split() nv = len(val) rt = 1 for k in range(nv): r_tbl[k] = int(val[k]) rt = rt*r_tbl[k] fin.close() # Check t_Q,t_R Qm = ecc_Mul(t_r, t_Q) Rm = ecc_Mul(t_r, t_R) if Qm != [0,0]: print("Q is not ECC point") sys.exit() if Rm != [0,0]: print("R is not ECC point") sys.exit() # Rho method information fout.write("----- Solving information -------------\n") fout.write(" Search pointer size="+format(st_size)+"\n") return nv #=========================================================# # Initialize in Rho step # # p_T = p_Q + p_R, p_alp = 1, p_bet = 1 # # st_ptr = -1 (all element) : not used flag # #---------------------------------------------------------# # return : Tlim = g_p / val : search upper limit # # val = (l2r - 30)/3, but >=1 # # l2r = log2(g_p) # #---------------------------------------------------------# # copy right : Ushiro Yasunori (ISCPC) # # date : 2020/03/17 # #=========================================================# def rho_init(): global p_T,p_alp,p_bet # Rho step variable global st_ptr,st_n # pointer table,number global fout # output file # set p_T,p_alp,p_bet p_T = ecc_Add1(p_Q,p_R) # T=Q+R on ECC p_alp = 1 p_bet = 1 # set st_n,Tlim st_n = int(math.sqrt(g_r)*2) Tlim = g_p if st_n > st_size: st_n = st_size Tlim = g_p // 4 # set st_ptr for k in range(st_n): # all pointer = -1 st_ptr[k] = -1 return Tlim #=========================================================# # Rho step in Rho method # # ip = p_T[0] (mod 3) # # ip=0: T=2*T, alp=2*alp, bet=2*bet (mod r) # # ip=1: T=T*Q, alp=alp+1 (mod r) # # ip=2: T=T*R, bet=bet+1 (mod r) # #---------------------------------------------------------# # copy right : Ushiro Yasunori (ISCPC) # # date : 2020/03/16 # #=========================================================# def rho_step(): global p_T,p_alp,p_bet # Rho step variable ip = p_T[0] % 3 # ip=p_T[0] (mod 3) if (ip == 0): p_T = ecc_Add2(p_T) # p_T = 2*p_T on ECC p_alp = (2*p_alp) % g_r p_bet = (2*p_bet) % g_r if (ip == 1): p_T = ecc_Add1(p_T,p_Q) # p_T += p_Q on ECC p_alp = (p_alp + 1) % g_r if (ip == 2): p_T = ecc_Add1(p_T,p_R) # p_T += p_R on ECC p_bet = (p_bet + 1) % g_r #=========================================================# # Set list plus Ty # # plas Ty : p_T[1] <= g_p/2 # # return : list = [p_T[0],p_alp,p_bet] # #---------------------------------------------------------# # copy right : Ushiro Yasunori (ISCPC) # # date : 2020/01/30 # #=========================================================# def plus_Ty(): list = [p_T[0], p_alp, p_bet] if p_T[1] > (g_p >> 1): list[1] = g_r - p_alp list[2] = g_r - p_bet return list #=========================================================# # Solve n with T == old-T # # n = (al-Tab[1])/(Tab[2]-bl) (mod g_r) # # Tab[1]=p_alp, Tab[2]=p_bet under p_T[0] > 0 # #---------------------------------------------------------# # copy right : Ushiro Yasunori (ISCPC) # # date : 2020/03/18 # #=========================================================# def sol_equal(loop, al,bl): global fout, n_tbl, n_sol,n_lp Tab = plus_Ty() # get alp,bet : T[0]>0 bval = Tab[2] - bl n = 0 if bval != 0: binv = inv_mod(bval,g_r) # 1/bval (mod g_r) aval = al - Tab[1] n = (aval*binv) % g_r n_sol[n_tbl] = n # set n (solve) n_lp[n_tbl] = loop # set loop #=========================================================# # Solve n with T == 0 # # n = Tab[1]/(g_r-Tab[2]) (mod g_r) # # Tab[1]=p_alp, Tab[2]=p_bet under p_T[0] > 0 # #---------------------------------------------------------# # copy right : Ushiro Yasunori (ISCPC) # # date : 2020/03/18 # #=========================================================# def sol_zero(loop): global fout, n_tbl, n_sol,n_lp Tab = plus_Ty() # get alp,bet : T[0]>0 bval = g_r - p_bet binv = inv_mod(bval,g_r) # 1/bval (mod g_r) n = (p_alp*binv) % g_r n_sol[n_tbl] = n # set n (solve) n_lp[n_tbl] = loop # set loop #=========================================================# # Equal check in Rho method # # if pointer equal: # # if Tx equal: # # sol_equal() # # return 1 #solved by equal # # else: # # change st_Tab # # else: # # change st_ptr # # add st_Tab # # return 0 # #---------------------------------------------------------# # copy right : Ushiro Yasunori (ISCPC) # # date : 2020/01/30 # #=========================================================# def equal_ck(loop): global st_Tab,st_ptr,st_n # check pointer equal pt = p_T[0] % st_n # p_T[0] (mod st_n) stv = st_ptr[pt] Tab = plus_Ty() # get value : T[0]>0 if stv >= 0: # check used table # check Tx equal if p_T[0] == st_Tab[stv][0]: alp = st_Tab[stv][1] # alp in table bet = st_Tab[stv][2] # bet in table sol_equal(loop,alp,bet) # solve equal return 1 # change st_Tab else: st_Tab[stv] = Tab # change st_ptr & add st_Tab else: st_ptr[pt] = len(st_Tab) # point last st_Tab.append(Tab) return 0 #=========================================================# # Rho method in ECC solver # # set Rho table # # initialize Rho step # # loop from solve or limit # # if p_T == zero: sol_zero # #---------------------------------------------------------# # copy right : Ushiro Yasunori (ISCPC) # # date : 2020/03/17 # #=========================================================# def rho_method(): Tlim = rho_init() # Rho step initaial lpmx = int(math.sqrt(g_r)*10) # Rho loop max # Rho step outc = 0 for lp in range(lpmx): outc += 1 rho_step() if p_T == [0,0]: # solve by T=zero sol_zero(lp+1) return 2 if p_T[0] < Tlim: # only feature poin id = equal_ck(lp+1) # check T=old T if id != 0: return 1 # solved by T= old T if outc >= 5000000: print("steps=",(lp+1)//1000000,"M") outc = 0 return 0 # not solve #=========================================================# # Solve n for R=n*Q # # p_Q,p_R : Points on the ECC, n : small integer # #---------------------------------------------------------# # copy right : Ushiro Yasunori (ISCPC) # # date : 2020/03/19 # #=========================================================# def SolRnQ(mlp): global g_r, p_Q,p_R, fout global fout, n_tbl, n_sol,n_lp Rck = p_Q # fout.write(" Q="+format(p_Q)+"\n") # fout.write(" R="+format(p_R)+"\n") for k in range(mlp): if Rck == p_R: n = (k + 1) % g_r n_sol[n_tbl] = n # set n (solve) n_lp[n_tbl] = k + 1 # set loop return n Rck = ecc_Add1(Rck, p_Q) return -1 #=========================================================# # CRT (Chinese Remainder Theorem) with 2 number # #---------------------------------------------------------# # Get c from a=c (mod p), b=c (mod q), gcd(p,q)=1 # # s=1/p (mod q), t=s*(b-a) (mod q) # # c=a+p*t # #---------------------------------------------------------# # copy right : Ushiro Yasunori (ISCPC) # # date : 2020/03/15 # #=========================================================# def CRT_2(a,b, p, q): s = inv_mod(p, q) # s=1/p (mod q) t = ( s*(b - a) ) % q c = a + p*t return c #=========================================================# # CRT (Chinese Remainder Theorem) with n number # #---------------------------------------------------------# # Get c from a[k]=c mod p[k], k=0,1,...,n-1 # #---------------------------------------------------------# # copy right : Ushiro Yasunori (ISCPC) # # date : 2020/03/15 # #=========================================================# def CRT_n(n, A, P): p = P[0] a = A[0] for k in range(n-1): q = P[k+1] b = A[k+1] c = CRT_2(a,b, p,q) p = p*q a = c return c #=========================================================# # main Rho program to solve n from R=n*Q # # Q,R : Points on the ECC, n : integer # #---------------------------------------------------------# # copy right : Ushiro Yasunori (ISCPC) # # date : 2020/03/19 # #=========================================================# nrf = ecc_read() # nrf=number of r-foctor t1 = time.time() tlp = 0 for k in range(nrf): n_tbl = k g_r = r_tbl[k] # r-factor s_r = t_r // g_r p_Q = ecc_Mul(s_r, t_Q) p_R = ecc_Mul(s_r, t_R) if p_Q == [0,0]: n_sol[k] = 0 n_lp[k] = 0 else: if (g_r%2 != 0) and (g_r%3 != 0): id = rho_method() # Rho method else: n = SolRnQ(g_r) # Direct Solve print(k," r=",g_r," loop=",n_lp[k]) ot1 = " r="+format(g_r)+" sol=" ot = ot1+format(n_sol[k])+" loop="+format(n_lp[k])+"\n" fout.write(ot) tlp = tlp + n_lp[k] # compute n from n (mod) fout.write("----- Decoding result n for R=n*Q --------\n") nsol = CRT_n(nrf, n_sol,r_tbl) fout.write("solve n="+format(nsol)+"\n") print("solve n=",nsol) t2 = time.time() fout.write(" solve total loop="+format(tlp)+"\n") uset = "{:.2f}".format(t2-t1) print(" used time=",uset) fout.write(" used time="+uset+" sec\n") fout.close()