#=========================================================# # Rho method to solve n from R=n*Q # # Q,R : Points on the ECC, n : integer # #---------------------------------------------------------# # copy right : Ushiro Yasunori (ISCPC) # # date : 2020/01/30 # #=========================================================# # 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 ========================# ct_bit = 6 #size of t_c,t_Q,t_R: ct_tabl=2^ct_bit st_bit = 20 #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 ct_size,ct_mask # table size for Rno step global st_size,st_mask # table size for solved check global p_Q,p_R # ECC point: Q=n*R, n:unknown global t_c,t_Q,t_R # table in Rho step global p_T,p_alp,p_bet # Rho step variable t_c = [] # value table for Rho step t_Q = [] # Q pointfor Rho step t_R = [] # R pointfor Rho step #----- long table for eqaul search -----------------------# ct_size = 2 << (ct_bit - 1) ct_mask = ct_size - 1 st_size = 2 << (st_bit - 1) st_mask = st_size - 1 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/01/30 # #=========================================================# 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/01/30 # #=========================================================# 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]: 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 open and read data # # g_p,g_r,g_a,g_b : values on ECC # # Q,R : Points on the ECC, R = n*Q on ECC # #---------------------------------------------------------# # copy right : Ushiro Yasunori (ISCPC) # # date : 2020/01/30 # #=========================================================# def ecc_read(): global g_p,g_r,g_a,g_b # values on ECC global p_Q,p_R # ECC point: Q=n*R, n:unknown global fin,fout # input,output file # file open fin = open("in.txt", "r") fout = open("out.txt", "w") fout.write("Rho method to solve n from R=n*Q by python\n") fout.write("----- Input information ---------------\n") # read g_p,g_r data = fin.readline() val = data.split() outs = " p=" + val[0] + " r=" + val[1] + "\n" fout.write(outs) g_p = int(val[0]) g_r = int(val[1]) # read g_a,g_b data = fin.readline() val = data.split() outs = " a=" + val[0] + " b=" + val[1] + "\n" fout.write(outs) g_a = int(val[0]) g_b = int(val[1]) # read p_Q data = fin.readline() val = data.split() outs = " Q=(" + val[0] + "," + val[1] +")\n" fout.write(outs) p_Q = [int(val[0]), int(val[1])] # read p_R data = fin.readline() val = data.split() outs = " R=(" + val[0] + "," + val[1] +")\n" fout.write(outs) p_R = [int(val[0]), int(val[1])] fin.close() # Check p_Q,p_R Qm = ecc_Mul(g_r, p_Q) Rm = ecc_Mul(g_r, p_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("----- Rho method information ----------\n") slp = format( int(math.sqrt(math.pi*g_r/2)) ) fout.write(" Standard Rho steps:sqrt(πr/2)="+slp+"\n") fout.write(" Rho table size="+format(ct_size)+"\n") fout.write(" Search pointer size="+format(st_size)+"\n") #=========================================================# # Set table in Rho step # # set t_c,t_Q,t_R # # t_c[k] = 7^k (mod g_r), k=0,1,...,ct_size-1 # # t_Q[k]=t_c[k]*p_Q, t_R[k]=t_c[k]*p_R # #---------------------------------------------------------# # copy right : Ushiro Yasunori (ISCPC) # # date : 2020/01/30 # #=========================================================# def tabl_set(): global t_c,t_Q,t_R # table in Rho step cv = 1 for k in range(ct_size): t_c.append(cv) Q = ecc_Mul(cv,p_Q) t_Q.append(Q) R = ecc_Mul(cv,p_R) t_R.append(R) cv = cv*3 % g_r #=========================================================# # 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/01/30 # #=========================================================# def rho_init(): global p_T,p_alp,p_bet # Rho step variable global st_ptr # pointer table 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_ptr for k in range(st_size): # all pointer = -1 st_ptr[k] = -1 # set Tlim l2r = int(math.log2(g_p)) sft = (l2r - 30) // 3 if sft < 0: val = 1 else: val = 2 << sft # val=2^sft Tlim = g_p // val fout.write(" Average search interval="+format(val)+"\n") fout.write("----- Solved by Rho information -------\n") return Tlim #=========================================================# # Rho step in Rho method # # ip = p_T[0] (mod ct_size) # # if p_T[0] < g_p/2: # # p_T += t_Q[ip] on ECC, p_alp += t_c[ip] (mod g_r) # # else: # # p_T += t_R[ip] on ECC, p_bet += t_c[ip] (mod g_r) # #---------------------------------------------------------# # copy right : Ushiro Yasunori (ISCPC) # # date : 2020/01/30 # #=========================================================# def rho_step(): global p_T,p_alp,p_bet # Rho step variable p2 = g_p >> 2 # p2=g_p/2 ip = p_T[0] & ct_mask # ip=p_T[0] (mod ct_size) cv = t_c[ip] if p_T[0] < p2: p_T = ecc_Add1(p_T,t_Q[ip]) # p_T += t_Q[ip] on ECC p_alp = (p_alp + cv) % g_r else: p_T = ecc_Add1(p_T,t_R[ip]) # p_T += t_R[ip] on ECC p_bet = (p_bet + cv) % 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/01/30 # #=========================================================# def sol_equal(loop, al,bl): global fout 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 nstr = format(n) lstr = format(loop) + "\n" outs = " solve: n=" + nstr + " Rho steps=" + lstr else: outs = " can not solve for alpha = beta\n" fout.write(outs) #=========================================================# # 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/01/30 # #=========================================================# def sol_zero(loop): Tab = plus_Ty() # get alp,bet : T[0]>0 bval = g_r - Tab[2] binv = inv_mod(bval,g_r) # 1/bval (mod g_r) n = (Tab[1]*binv) % g_r nstr = format(n) lstr = format(loop) + "\n" outs = " solve: n=" + nstr + " Rho steps=" + lstr fout.write(outs) #=========================================================# # 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 # check pointer equal pt = p_T[0] & st_mask # p_T[0] (mod st_size) 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/01/30 # #=========================================================# def rho_method(): tabl_set() # set Rho table 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 >= 1000000: print("steps=",(lp+1)//1000000,"M") outc = 0 return 0 # not solve #=========================================================# # 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/01/30 # #=========================================================# t1 = time.time() ecc_read() # file data read id = rho_method() # Rho method usize = format(len(st_Tab))+"\n" fout.write(" used search table size="+usize) if id == 0: fout.write(" Unsolved, insufficient number of iterations") t2 = time.time() uset = "{:.2f}".format(t2-t1) fout.write(" used time="+uset+" sec\n") fout.close()