橢圓曲線加密中的加法乘法淺析

本文不深入橢圓曲線加密演算法的全部知識，只針對橢圓曲線加密中需要用到的加法和乘法計算規則進行淺析。

實際練習中碰到一個比較簡單密碼學的問題，但是涉及到了橢圓曲線加密演算法，題目描述如下：

已知橢圓曲線加密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[0]-Rx) - G[1]) % mod
    temp=[Rx,Ry]
    #print temp

print temp

參考文獻：

http://blog.51cto.com/11821908/2057726
講解了受限域的曲線下的加法實現計算
https://www.jianshu.com/p/2e6031ac3d50
只講解了無受限域下曲線的加法
https://wenku.baidu.com/view/6f2879cca1c7aa00b52acb5f.html
分數求模原理介紹
https://www.pediy.com/kssd/pediy06/pediy6014.htm
看雪論壇上的詳細介紹，提供了加法運算的驗證集
https://blog.csdn.net/baidu_38271024/article/details/78881031
乘法逆元求解的python實現