https://www.cnblogs.com/Kalafinaian/p/7392505.html
橢圓曲線加密中的加法乘法淺析
本文不深入橢圓曲線加密算法的全部知識,只針對橢圓曲線加密中需要用到的加法和乘法計算規則進行淺析。
實際練習中碰到一個比較簡單密碼學的問題,但是涉及到了橢圓曲線加密算法,題目描述如下:
已知橢圓曲線加密Ep(a,b)參數為
p = 15424654874903
a = 16546484
b = 4548674875
G(6478678675,5636379357093)
私鑰為
k = 546768
求公鑰K(x,y)
提示:K=kG
這里需要介紹一下橢圓曲線
一般,橢圓曲線可以用以下二元三階方程的形式來表示:
y² = x³ + ax + b,其中a、b為系數。
它大概的幾何形狀如下圖:

而本文要介紹的加法和乘法,就是基於這樣一個奇怪的幾何圖形來做到的。
橢圓曲線加法(非有限域):
在橢圓曲線上取一點P(Xp,Yp),再取一點Q(Xq,Yq),連接P、Q兩點作一條直線,這條直線將在橢圓曲線上交於第三點G,過G點作垂直於X軸的直線,將過橢圓曲線另一點R(一般是關於X軸對稱的點),R點則被定義為P+Q的結果,既P+Q=R:

當P=Q的情況下,直線將是橢圓曲線在P(Q)點上的切線,而G點是這條切線和曲線的另一個交點,同樣,P+Q=R:

通過上述的圖片和文字描述,已經在幾何圖形上給出了橢圓曲線加法的定義,可是如果要公式化,該如何快速計算呢?
這里只提供快速計算公式,不提供證明,證明可以自己再去解方程組推導一下:
計算P+Q=R
當P!=Q時,兩點縱坐標相減的值與橫坐標相減的值就是直線的斜率:
λ = (Yq - Yp)/(Xq - Xp)
當P=Q,計算過P(Q)點切線的斜率,既橢圓曲線公式兩邊求導相除:
λ = (3Xp² + a)/2Yp
斜率計算之后,對點R的坐標進行計算,公式如下:
Xr = (λ² - Xp - Xq)
Yr = (λ(Xp - Xr) - Yp)
通過上述公式,可以快速計算橢圓曲線上任意兩點的加法和,這里給出加法實現的python代碼:
if P == Q: aaa=(3*pow(P[0],2) + a) bbb=(2*G[1]) k=(aaa/bbb) else: aaa=(Q[1]-P[1]) bbb=(Q[0]-P[0]) k=(aaa/bbb) Rx=(pow(k,2)-P[0] - Q[0]) Ry=(k*(P[0]-Rx) - P[1])
橢圓曲線加法(有限域)
實數范圍上光滑的橢圓曲線在密碼學應用上並不合適,需要進行有限域下的離散化操作才能使用。

現在將上述的橢圓曲線加法計算公式適當修改,以適應有限域下的計算:
當P!=Q時,兩點縱坐標相減的值與橫坐標相減的值需要與p進行取余操作:
λ = (Yq - Yp)/(Xq - Xp) mod p
當P=Q,計算過P(Q)點切線的斜率,既橢圓曲線公式兩邊求導相除,結果也需要與p進行取余操作:
λ = (3Xp² + a)/2Yp mod p
斜率計算之后,對點R的坐標進行計算,公式如下:
Xr = (λ² - Xp - Xq) mod p Yr = (λ(Xp - Xr) - Yp) mod p
通過比較,有限域下的計算只是對結果進行了取余操作,上述公式看起來已經解決了有限域下的橢圓曲線加法。
但是如果在編寫代碼,計算實際的例子時,有很大可能會得到錯誤的結果,
其根源在於λ = (Yq - Yp)/(Xq - Xp) mod p
或λ = (3Xp² + a)/2Yp mod p
在進行取余計算之前,除數和被除數之前可能並不是一個整除的關系。
如:1/4 mod 23
,如果直接進行處理,將會得到結果0。
但是在分數求模計算中,是如下定義的:
計算a/b(mod n) a/b (mod n)=a*b^-1(mod n) 計算1/b mod n =b^(-1) mod n 就是求y,滿足: yb = 1 mod n y是有限域F(n)上x的乘法逆元素
簡單點說,假設需要求上述的1/4 mod 23
,可以轉化為1*4(-1次方) mod 23
,又可以轉化為1*(4和23的乘法逆元) mod 23
。
而計算乘法逆元,可以通過拓展歐幾里得計算得到,這里對拓展歐幾里得不作展開,只提供一個簡單算法流程描述:
ExtendedEuclid(d,f) 1 (X1,X2,X3):=(1,0,f) 2 (Y1,Y2,Y3):=(0,1,d) 3 if (Y3=0) then return d'=null//無逆元 4 if (Y3=1) then return d'=Y2 //Y2為逆元 5 Q:=X3 div Y3 6 (T1,T2,T3):=(X1-Q*Y1,X2-Q*Y2,X3-Q*Y3) 7 (X1,X2,X3):=(Y1,Y2,Y3) 8 (Y1,Y2,Y3):=(T1,T2,T3) 9 goto 3
得到乘法逆元后,橢圓曲線上的加法運算計算就簡單了,實現Python代碼如下:
#coding:utf-8 #歐幾里得算法求最大公約數 def get_gcd(a, b): k = a // b remainder = a % b while remainder != 0: a = b b = remainder k = a // b remainder = a % b return b #改進歐幾里得算法求線性方程的x與y def get_(a, b): if b == 0: return 1, 0 else: k = a // b remainder = a % b x1, y1 = get_(b, remainder) x, y = y1, x1 - k * y1 return x, y #返回乘法逆元 def yunsle(a,b): #將初始b的絕對值進行保存 if b < 0: m = abs(b) else: m = b flag = get_gcd(a, b) #判斷最大公約數是否為1,若不是則沒有逆元 if flag == 1: x, y = get_(a, b) x0 = x % m #對於Python '%'就是求模運算,因此不需要'+m' #print(x0) #x0就是所求的逆元 return x0 else: print("Do not have!") if P == Q: aaa=(3*pow(P[0],2) + a) bbb=(2*P[1]) if aaa % bbb !=0: val=yunsle(bbb,mod) y=(aaa*val) % mod else: y=(aaa/bbb) % mod else: aaa=(Q[1]-P[1]) bbb=(Q[0]-P[0]) if aaa % bbb !=0: val=yunsle(bbb,mod) y=(aaa*val) % mod else: y=(aaa/bbb) % mod Rx=(pow(k,2)-P[0] - Q[0]) % mod Ry=(k*(P[0]-Rx) - P[1]) % mod
橢圓曲線乘法
簡單介紹完橢圓曲線上定義的加法運算,橢圓曲線上的乘法運算就比較簡單了,因為加法可以退化為加法運算,就像算數上的1*3等價與1+1+1。
假設我們需要求2P,則可以化簡為P+P=2P
同理,當我們需要求3P時,可以化簡為P+2P=3P,其中2P=P+P
最后,我們可以得到規律,當求nP時(n為任意正整數),P+(n-1)P=nP,其中(n-1)P=P+(n-2)P
這樣,通過上述介紹的橢圓曲線加法公式,完全可以進行橢圓曲線的乘法計算
以本文開頭的題目為例,給出Python代碼實現:
#coding:utf-8 #歐幾里得算法求最大公約數 def get_gcd(a, b): k = a // b remainder = a % b while remainder != 0: a = b b = remainder k = a // b remainder = a % b return b #改進歐幾里得算法求線性方程的x與y def get_(a, b): if b == 0: return 1, 0 else: k = a // b remainder = a % b x1, y1 = get_(b, remainder) x, y = y1, x1 - k * y1 return x, y #返回乘法逆元 def yunsle(a,b): #將初始b的絕對值進行保存 if b < 0: m = abs(b) else: m = b flag = get_gcd(a, b) #判斷最大公約數是否為1,若不是則沒有逆元 if flag == 1: x, y = get_(a, b) x0 = x % m #對於Python '%'就是求模運算,因此不需要'+m' #print(x0) #x0就是所求的逆元 return x0 else: print("Do not have!") mod=15424654874903 #mod=23 a=16546484 #a=1 b=4548674875 #b=1 G=[6478678675,5636379357093] #G=[3,10] #次數 k=546768 temp=[6478678675,5636379357093] #temp=[3,10] for i in range(0,k): if i == 0: aaa=(3*pow(G[0],2) + a) bbb=(2*G[1]) if aaa % bbb !=0: val=yunsle(bbb,mod) y=(aaa*val) % mod else: y=(aaa/bbb) % mod else: aaa=(temp[1]-G[1]) bbb=(temp[0]-G[0]) if aaa % bbb !=0: val=yunsle(bbb,mod) y=(aaa*val) % mod else: y=(aaa/bbb) % mod #print y Rx=(pow(y,2)-G[0] - temp[0]) % mod Ry=(y*(G[