理解"質數有限場域倍乘群的橢圓曲線加密運算"

質數(prime number)有限場域(finite field)倍乘群(multiplicative group),用集合符號表示該質數形成的場域: Zp* = {x | 0 < x < p} 其中 x 是整數, p 是質數
1. 用 mod p (以質數當除數取餘數)的整數經過四則(+-*÷)運算後將映射回質數場域裏面.
2. 場域的質數 p 及原點 0, 不包含在該場域裏面,但在運算中會出現,分別是場域的上限及下限.
3. 場域內的整數 x 與該質數 p 互質,也就除了 1 以外無其它公因數:  GCD(x, p) = 1
4. 假設φ(p)是該場域內與 p 互質的個數,根據 Euler 定理:
              x^φ(p) mod p = 1         其中  φ(p) = p-1  就是集合裏面{ x }與 p 互質的個數
場域內 Zp* 所有整數 x 共有 p-1 個,且都與 p 互質.當 x 倍乘 p-1 次,以 p 當除數再取餘數運算將等於 1
              x^p-1 mod p = 1     ---(1)
兩邊同時乘上 x 就可以得出:
              x^p mod p = x         --- (2)
也就是說場域內的整數 x 倍乘 p 次後,再用 mod p 取餘數就會映射成(還原回)整數本身.這就是為何 ECC 會使用質數的原因? 因為該整數群裡(x)的模倒數(y)可以很容易計算出來.
           x * y  mod p = 1       --- (3)
將第(1)式分解 x
           x^p-1 mod = x * x^p-2 mod p    ---(4)         
比較(3)(4)式得出 x 的模倒數
           y = x^p-2 mod p
也就是說整數的模倒數與質數的關係式是將 x 倍乘 p-2 次後,再以 p 當除數取餘數:
           x^-1  mod p = x modinv p  =   x^p-2 mod p
ECC 用大整數當私鑰,曲線上的點當基底,兩個乘起來形成公鑰後,仍還是曲線上的一點,再指定一個質數 p 去做運算,形成一個整數場域.所有座標內的數值也都被限制在該場域內(mod p).由於曲線上兩點加法產生曲線上另外一點的特性.因此複予了它可以用數學運算法則,當使用通訊協定來加密或解密鑰匙,雖然公開了公鑰及質數 p 還有所使用的橢圓曲線,卻不影響到私鑰的保密性,因為解密出私鑰須耗費相當大的時間或空間(成本).給定一個橢圓曲線的點(公鑰)及質數,破解私鑰,也就是所謂的 ECDLP (橢圓曲線離散對數求解),或許有人認為可以破解 ECDLP,當使用的密碼很小的時候,例如 15, 運算點連加15次就可以找到破解點了(但這只是一種暴力破解法),當密碼是2的 128 次方(也就說 128 個位元)呢?"只要"(暴力)嘗試最多 2^128-1 就一定可以找破解點.但是這個暴力的方法恐怕是好幾?年的事(時間),仔細想一下 2 的 128 次方是一個多大的整數呢?可以用 python 稍微感受一下,它支援大整數,大整數的次方用 ** 來加以運算:
a=2**128
340282366920938463463374607431768211456L
乍看之下好像沒什麼,如果超級電腦每秒運算1百萬個百萬個百萬個百萬指令  = 1e6*1e6*1e6*1e6 = 10²⁴(10的24次方)指令來算,一年共 365*24*60*60 秒,所以一年來運算了指令總共
b=365*24*60*60*1e24
3.1536e+31

將a(等於a¹²⁸)除以 b(等於花一整年的指令數)=花費數(年):

a/b

10790283.070806015

取整數看就好了,總共要花掉 10790283 年的時間啊!假如將密鑰位元數提升一倍,也就是256位元組變成:
c=2**256
c/b
3.671743063080803e+45
一個天文數字 10 的 45 次方,嚇死人!運算時間隨著密碼長度成指數性增加,一旦保密性不夠,再加大密碼的位元數就能因應了!又或者用多部超級電腦同時進行(空間)呢?那恐怕需要計算一下需用多少部電腦來運作,但這樣花的值得嗎?如果不用暴力破解,那就要想辦法解決 ECDLP, 但目前為止還無有效方法可用.反而用亂數的方式隨便猜,機率極低卻有可能猜中,想要破解 ECC 密碼,門都沒有.如果故意放了後門,那就另當別論了.

圓弧點的數學運算

一個半徑為 r 圓的方程式,圓線上的點(x, y)滿足方程式 x² + y² = r² 來描述
或是用 f(x, y) = x² + y² - r² = 0 上所有 P(x, y) 的點畫出來就會形成一個圓, 他的座標就是:
           x = r cos(θ)              0 <  θ < 2Π
           y = r sin(θ)               0 < θ < 2Π
已知一個點 (x, y) 等於是指定一個角度 θ 所在的圓點上,換句話說, θ1 座標是:
          x1 =  r cos(θ1) - > cos(θ1) = x1 / r 
          y1 =  r sin(θ1) - > sin(θ1) = y1 / r
同樣的另外一個點 (x2, y2)角度是θ2,等於知道該點的座標在:
          x2 = r cos(θ2) - > cos(θ2) = x2 / r
          y2 = r sin(θ2) - > sin(θ2) = y2 / r
以此類推當 θ3 = θ1 + θ2 位置在那呢? 當然就在 θ3 所指的圓線上的一點:
           x3 = r cos(θ3))
           y3 = r sin(θ3)
再將 θ3 替換成 θ1 + θ2 就可得到:
          x3 = r cos(θ3) = r cos(θ1 + θ2) 
          y3 = r sin(θ3)  = r sin(θ1 + θ2)
根據三角公式將 sin(θ), cos(θ) 換成(x1, y1)與(x2, y2)的運算關係就是:
          x3 = r cos(θ1 + θ2) = r [cos(θ1) cos(θ2) - sin(θ2) cos(θ1)] 
              = r [(x1 / r)(x2 / r) - (y1 / r) (y2 / r)] = (x1 x2 - y1 y2) / r   --- (1)
          y3 = r sin(θ1 + θ2) = r [sin(θ1) cos(θ2) + sin(θ2) cos(θ1)]
                 r [(y1 / r)(x2 / r) + (y2 / r)(x1 / r)] = (x1 y2 + y1 x2) / r   --- (2)
這就是圓座標點加法運算式(Adder),兩圓點相加去算出另一個圓點的座標,用數學記號表示就是:
          P(x1, y1) + Q(x2, y2) = R(x3, y3)     簡寫成 而 P +  Q =  R
因此2倍角座標就很容易推導出來, 因為 (x, y) = (x1, y1) = (x2, y2) 帶入(1)(2)式
         x" = r cos(2θ) = (x² - y²) / r   --- (3)
         y" = r sin(2θ) = 2xy / r             --- (4)
這就是圓座標點加倍運算式(Doubler),簡寫成 2P, 再利用式(1)(2)座標點相加運算式, 將 θ 座標點與兩倍角 2θ 座標點相加就可推導出 3 倍角座標:
         x'" = r cos(3θ) = (x(x² - y²) / r - y 2xy / r) / r = (x³ - 3xy²) / r²
         y'" = r sin(3θ) = (x 2xy / r + y(x² - y²) / r) / r = (3x²y- y³) / r²
以此類推, 用 Double and And 疊帶累加法就可算出 n 倍角的座標,弳度的乘法運算 θ + θ        ... 共 n 個 = n * θ 也是我們熟悉的用法,只要知道一個角度的座標點, 即使不用三角函數也可算出其他角度的座標,再來看弳度 0 的座標,用記號 O 來表示,也就是說 O 的座標位在:
         x0 = r cos(0) =r
         y0 = r sin(0) =0
帶入圓點加法公式去看 P(x, y) + O 是否等於 P(x, y)呢?
           (x*r  - y*0) / r = x
           (x*0 + y*r) / r = y 
確實如此,因此座標(r, 0)的確是中性點 O, 一條線的座標點 f(x, y) = 0, 如果將(x, y)座標同時放大 z 倍時,由原先的座標來看這新的點所形成的圖形會長成什麼樣子呢? 既然新增了一個變數 z, 就把他想成是 z 軸,將(x/z, y/z)映射到(X, Y, Z)來表示,得到新的投射座標,換句話說就是分別用(X/Z,Y/Z)去取代方程式中的(x,y),比方一個圓方程式:

            x² + y² = r²      
將 x = X / Z, y = Y / Z 帶入得到:
            (X / Z)² + (X / Z)² = r²
            X² + Y² = r² Z² = (r Z)² = R²      , R = r Z
可以看到它仍舊是圓的方程式,只不過半徑跟著放大成Z倍.

用 C 語言寫個小程式驗證一下:
// cir.cpp use g++ to compile
// g++  cir.cpp
#include  <iostream>
#include  <math.h>
using namespace std;
class  PointC {
    public:
        double x, y;
        PointC( ) { };
        PointC(double x0, double y0) { x=x0; y=y0; }
};
class  curveC {
        PointC p; // private point Generator
    public:             
        curveC(const PointC &P)         { p = P; } //need to implement PointC() constructor
        curveC(const PointC &P, int n)  { p = P; if ( n > 1 ) scale(n); }
        PointC operator + (const PointC &Another)   { return incPointC(Another);  }
        PointC operator + (const curveC &Another)   { return incPointC(Another.p);}
        PointC operator * (int n)                   { return curveC(p, n).p;      }
        PointC incPointC(const PointC &Another)     { return curveC(p).adder(Another)->p;}
        curveC* adder(const PointC &q) { // self add, call by reference,
            double r = sqrt(p.x * p.x + p.y * p.y);
            double px = p.x;    // copy to prevent update by p.x
            p.x = (px * q.x - p.y * q.y) / r;  // update x
            p.y = (px * q.y + p.y * q.x) / r;   // update y  
            // return curveC(p); // if wrap with class, can use . not -> to obtain member
            return this;
        }
        curveC* doubler() { // self double
            double r = sqrt(p.x * p.x + p.y * p.y);
            double px = p.x;    // copy to prevent update by p.x
            p.x = (p.x * p.x - p.y * p.y) / r ; // update x
            p.y = 2 * px * p.y / r;             // update y
            return this;
        }         
        #define  testBit(n, b)   n & (1L << b)
        curveC* scale(int n)   {   // double-and-add   DAA algorithm
            int k = 32;             // 32 bit integer
            PointC G = p;           // copy to prevent update
            while ( k-- )   if (testBit(n,k))   break;
            while ( k-- ) {         // Double And Add
                doubler();          // self double
                if (testBit(n,k))   adder(G);// self adder with G     
            }
            return this;
        }
        #undef testBit
        curveC* showPoint(void) {
            printf("\tx=%f , y=%f\n", p.x, p.y);
            return this;
        }
};
main() {  
    double degree = 10; 
    double theta  = 2 * 3.14159 * degree / 360;
    double r = 5;      // Circle half diameter = 5
    PointC G = PointC(r * cos(theta), r * sin(theta) );
    curveC c = curveC(G);
    for (int i = 1; i <= 10; i ++) {
        curveC(c * i).showPoint();
        printf("equal%2d\tx=%f , y=%f\n\n",i ,r*cos(theta*i), r*sin(theta*i));     
    }
    c.showPoint();
}

橢圓曲線的轉換

參考文章: https://en.wikipedia.org/wiki/Montgomery_curve
Montgomery curve(蒙哥馬利曲線) 座標軸(x,y)先用(i,j)表示:
           B*j^2 = i^3 + A*i^2 + i
其中座標(0,0)帶進去可以讓上述等式成立,因此它是該曲線上的一個點.當兩邊都除以 B^3 後:
           (j/B) ^2 = (i/B)^3 +(A/B) * (i/B)^2 + (1/B)^2 * (i/B)
將 j/B 替換成 v 並將 i/B 替換成 u, 就變成縮小版的蒙哥馬利曲線(縮小程度是B * B):
            v^2 = u^3 +(A/B) * v^2 + (1/B)^2 * v
如果再把 u替換成 t - A/3B 就變成縮小再左移版的蒙哥馬利曲線:
            v^2 = (t - A/3B) ^3 + (A/B) * (t-A/3B)^2+ (1/B^2) * (t-A/3B)
            v^2 = [ t^3 - 3*(A/3B) t^2 +3(A/3B)^2*t - (A/3B)^3  ]
                     +           (A/B)* [ t^2 - 2*(A/3B) *t  +(A/3B)^2  ]
                                                     +  [  (1/B)^2 *t  -  A/3(B^3) ]
因j為 A,B 都是常數, 常用的 B = 1 , 就可以更簡化成:
            v^2 = t^3 + 0 + (A^2/3  - 2A^2/3 +1 ) t +(A^3/9 - A^3/27 -A/3)
                   = t^3 + (3-A^2)/3 *t + (2*A^2-9*A)/27 = t^3 + a*t +b
            a = (3-A^2)/3
            b= (2*A^2 - 9*A)/27
因此:
           y^2 = x^3 + a * x + b 是 y^2 =x^3 +A*x^2 + x 的一種變形 
座標變換過程:

(x, y) -> (i/B, j/B) -> (u, v) -> (t-A/3B, v) -> (t, v) -> (x,y)
x , y 軸各放大 B 倍, 等同圖形會縮小 B*B=B^2 倍 --- 變換1
原點右移A/3B, 等同圖形向左平移A/3B --- 變換 2
可見得 y^2 = x^3 + ax + b 其實是蒙哥馬利曲線的一個變種

橢圓曲線上奇怪的數學

一個橢圓曲線方程式 E(x, y) 在 x-y 平面滿足二元三次方程式:
            E(x, y) ::  y² = x³ + ax + b        ---- (0)
假設 P(x₁, y₁) 是在曲線上的一點, 如果以 P(x₁, y₁) 當基準點畫一條斜率是 λ 的直線,這條直線用方程式描述:
           -y = λ ( x₁ - x ) - y₁                   --- (1)
重新整理:
             y =   λ x + c             c  是常數 = y₁ - λ x₁    --- (2)
如果將式(2)帶入式(0)以就可得出交點,先看等式(0)右邊,如果將通過P的直線帶入時,也就是說:
          ( λ x + c )² = ( x³ + a x + b)     --- (3)
上式是一個 x 三次方程式,因此 x 會有 3 個根(root),也就是說會有 3 個交點,因為P(x₁,y₁) 在線上,是已知的根(root), 另外兩交點分別是 Q(x₂,y₂) 與 R(x₃,y₃) 滿足曲線方程式:
          y₂² = x₂³ + a x₂ + b,  y₂ = λ(x₂ - x₁) + y₁   --- (4)
          y₃² = x₃³ + a x₃ + b,  y₃= λ(x₃ - x₁) + y₁   --- (5)
從方程式(0)來看, 左邊是 y² , 可以知道當曲線上存在有一點 P(x₁, y₁), 以同樣的 x₁根, 必定有另一點(x₁, -y₁) 也在曲線上, 因為(-y₁)² = (y₁)² , 如果以座標軸來看, 就是以 y=0 為軸當鏡面的鏡像點, 用 -P = (x₁, -y₁) 來表示 同樣曲線上的另外兩交點 Q, R 也必定有其鏡像點座落在相同曲線上分別稱之為 -Q 與 -R,假設有另外一個虛點 O,可能在無窮遠處,可以讓:
           P + O = P                 --- (6)
           P + (-P) = O             --- (7)     
           P + Q + R = O         --- (8)
上述成立的條件就是當 O 是一個零點,有了零點就可以用加或減的數學運算式來描述交點 P, Q, R 的關係式(6,7,8),方程式(6)描述的是: 交點 P 加上零點還是原來的 P 點,方程式(7)描述的是: 交點 P 與其鏡像點 -P 加在一起就會消滅為零. 方程式(8)描述的是一個二元三次方程式與直線的全部交點 P, Q, R ,當全加一起就會成為零點. 用同樣邏輯可否推廣到其它二元四次  y² = f(x) 方程式甚至更高次項的方程式與直線交點的根解(root)呢？也就是說將全部根的點全加在一起是否就是零點？待研究...
如果將第(8)式兩邊都加上 R 的鏡像點 -R:
          (P + Q + R) + (-R) = O + (-R)
左邊應用第(7)式讓 R 消滅, 右邊應用第(6)式任何點加上零點還是原來的點就得到:        
           P + Q = -R              --- (9)
方程式(9)描述的是將兩交點 P=(x₁, y₁)及 Q=(x₂, y₂), 斜率λ的直線與橢圓曲線其中兩交點, 經過相加運算得出新座標 -R = (x₃, -y₃) 正好就是最後交點 R(x₃, y₃) 的鏡像點(x-y 平面上以 y=0 為軸當鏡面), 後續可以推導出座標是:
            x₃ = λ² - x₁ - x₂           --- (10)
          -y₃ = λ(x₁ - x₃) - y₁       --- (11)
也就是說:
          P + Q = (x₃, -y₃)
已知兩點(x₁, y₁)及(x₂, y₂), 相加後的新座標(x₃, -y₃),還是座落在橢圓曲線上.
λ  是直線的斜率 dy/dx, 對於兩點(x₁, y₁)及(x₂, y₂)形成的直線而言:
            λ = (y₂ - y₁) / (x₂ - x₁)      其中 x₁ != x₂
但是當 P = Q 時 x₁ = x₂, y₁ = y₂, 導致分母為零無法除,因此不能使用上述 λ 的公式,也就是說當P, Q 相同時要改成 Q 逼近 P, 讓兩交點(Q => P)形成一點, 通過該點的直線就是切線, 此切線會與橢圓曲線相交在 R' 點,它的鏡像 -R'  = P + Q = P + P = 2P, 用橢圓曲線方程式 E(x,y) 把函式 y = f(x), 作一次偏微分就可得出橢圓曲線的切線方程式:
           λ' (x, y) = dy/dx = (∂E/∂x) / (∂E/∂y) = (3x²  + a) / 2y   ---(12)
帶入切點座標(x₁, y₁)得出斜率 λ':
           λ'(x, y) = (3x₁²  + a)/ 2y₁          --- (13)
上述意思說當 Q = P 意味著是將 Q 點逼近橢圓曲線的切點 P(x₁, y₁)處.
該條切線所形成的另外交點 R', 其鏡像點 -R 是:
           -R' = P + Q@P = 2P
也就是說與切線相交的第3交點 R' 正好就是兩倍切點的鏡像點,  用 - 2P 表示:
            R' = - 2P
要注意的是 2P 雖位在橢圓曲線上, 但不在該條切線上, 橢圓曲線與切線的交點除了P, Q@P外, 另一交點是 -2P, 而運算自己加上自己時會轉移(transform)到另外一個橢圓曲線點 2P, 假設三交點分別是切點 P = Q= (x, y), R'(x', y'), 帶入公式(12)(10)(11)
        λ' = (3x² + a)/2y
        x' = λ'² - x - x = λ'² - 2x         ---(14)
       -y' = λ'(x - x') - y                     --- (15)
一旦有了基礎座標加法及2倍座標的運算公式, 利用 Double And Add (DAA)演算法, 可以導出倍乘(multiply)座標 3P, 4P, 5P, 6P..., 而這些座標點全都位在橢圓曲線上:
　  3P = 2P + P   double and and
      4P = 2(2P)     double and double
      5P =(2(2P) + P    double and double and add
      6P =(2(2P + P)) + P   double (double and add) and add
      ...
甚至稍微修改一下 DAA 演算法, 更進一步變成 Square And Multiply (SAM)演算法, 自行定義出指數(power) 運算 P², P³, P⁴, P⁵, P⁶ ...的座標, 把他映設到橢圓曲線空間 P, 例如:
     G  = n*P         n >0, n 是整數, multiply 利用上述 DAA 演算法運算 n*P
     G² = G*G     = n²P,  square:  n*n P
     G³ = (G²)*G = n³P,  square and multiply n²*n P

     G⁴ = (G²)²    = n⁴P,  square and square   (n²)² P

     G⁵ = (G²)²*G= n⁵P,  square then square and multiply (n²)²*n P

     G⁶ = (G²*G) ² *G = n⁶P,  square and multiply then square and multiple (n²*n)²*n P

       ^...
附錄:  推導 x 與斜率 λ 的關係式:
將直線方程式 (1) y 帶入曲線方程式 E(x,y), 與假設的 3 個根(x₁, x₂, x₃)所形成的三次方程式, 只要比較 x² 項的係數, 就可得到斜率 λ 與 x₁, x₂, x₃ 關係式:
        ( λ x + c )²  - (x³ + a x + b ) = (x - x₁)(x - x₂)(x - x₃) = 0
         x₃ = λ² - x₁ - x₂
-R 的 y 座標是 - y₃, 直接帶回直線方程式(1)解出:
        -y₃ = λ(x₁ - x₃) - y₁
詳細參考網路上的文章: http://peter.havercan.net/geometry/elliptic-curve-visualization.html
用 c++ 寫個程式驗證一下:
// ecc.cpp
#include <iostream>
using namespace std;
class  pointEC {
    public:
        double x, y;
        pointEC( ) { };
        pointEC(double x0, double y0) { x=x0; y=y0; }
};
class  curveEC { // curve:    y² = x³ + ax + b
        int isPoint(){ return ((p.y * p.y - p.x * p.x * p.x - a * p.x -b) < 1e-9) ? 1:0; }
        double a = 3;
        double b = -12;    
        pointEC p; // private point Generator
        pointEC addpointEC(const pointEC &p, const pointEC &q) { // P + Q = -R
            double m,fx;
            m = (p.y - q.y)/(p.x - q.x);// slope         
            fx = m * m - p.x - q.x;
            return pointEC(fx, m * (p.x - fx) - p.y);// return new object
        }
        void adder(const pointEC &q) { // self adder
            double m;
            double px;
            px = p.x;
            m = (p.y - q.y)/(p.x - q.x);// slope
            p.x = m * m - px - q.x;     // x of intersect point
            p.y = m * (px - p.x) - p.y; // y of intersect mirror point
        }
        void doubler() { // self double
            double m;
            double px;
            px = p.x;
            m=(3*p.x * p.x + a )/(2*p.y); // tangent@(x1,y1)
            p.x = m * m - 2 * p.x;          // x of intersect point
            p.y = m * (px - p.x) - p.y;     // y of intersect mirror point
        }      
        #define  testBit(n, b)   n & (1L << b)
        void scale(int n)   {   // double-and-add   DAA algorithm
            int k = 32;             // 32 bit integer
            pointEC G = p;          // copy to prevent update
            while ( k-- )   if (testBit(n,k))   break;
            while ( k-- ) {         // Double And Add
                doubler();          // self double
                if (testBit(n,k)) adder(G);// self adder with G    
            }
        }   
        #undef testBit
    public:     
        curveEC (const pointEC &G)        { p = G; } //need to implement pointEC() constructor
        curveEC (const curveEC &C)        { p = C.p;}
        curveEC (const pointEC &G, int n) { p = G  ; if (n > 1)  scale(n); }
        curveEC (const curveEC &C, int n) { p = C.p; if (n > 1)  scale(n); }
        curveEC operator + (const pointEC &Another)   { return curveEC(addpointEC(p,Another));   }
        curveEC operator + (const curveEC &Another)   { return curveEC(addpointEC(p,Another.p)); }
        curveEC operator * (int n)                    { return curveEC(p, n);                    }     
        void showPoint(string msg) {
            if( isPoint() ) printf("Yes\t");
            else            printf("No\t");
            cout << msg;
            printf("\tx=%8.4f , y=%8.4f\n", p.x, p.y);
        }
};
int main() {
    curveEC  c1 = curveEC(pointEC(4, 8));
    curveEC  c2 = curveEC(c1, 2);
    curveEC  c5 = curveEC(c1, 5);
    for ( int i = 1; i < 30 ; i++ ) {
        printf("c%-3d ",i);
        curveEC( c1 * i ).showPoint(":");
    }
    printf("Test for curve\n");
    c1.showPoint("\tc1");
    (c2 * 3).showPoint("\tc2*3");
    (c1 + c2 *4 + c5 *3).showPoint("\tc1+c2*4+c5*3");
}

輸出結果:

c1   Yes    :    x=  4.0000 , y=  8.0000
c2   Yes    :    x=  2.1602 , y= -2.1355
c3   Yes    :    x= 24.1877 , y=-119.2119
c4   Yes    :    x= 11.5206 , y= 39.3905
c5   Yes    :    x=  1.9013 , y=  0.7598
c6   Yes    :    x=  6.0004 , y=-14.9012
c7   Yes    :    x=121.0615 , y=1332.1461
c8   Yes    :    x=  2.8896 , y=  4.5603
c9   Yes    :    x=  2.7062 , y= -3.9922
c10  Yes    :    x= 79.2098 , y=-705.1261
c11  Yes    :    x=  6.6953 , y= 17.5559
c12  Yes    :    x=  1.8750 , y= -0.4660
c13  Yes    :    x=  9.9976 , y=-31.8947
c14  Yes    :    x= 30.2486 , y=166.6002
c15  Yes    :    x=  2.2599 , y=  2.5142
c16  Yes    :    x=  3.6788 , y= -6.9873
c17  No    :    x=2169.3136 , y=-101037.7420
c18  Yes    :    x=  4.3684 , y=  9.1905
c19  Yes    :    x=  2.0761 , y= -1.7824
c20  Yes    :    x= 19.7784 , y=-88.2287
c21  Yes    :    x= 13.4167 , y= 49.4305
c22  Yes    :    x=  1.9404 , y=  1.0616
c23  Yes    :    x=  5.4084 , y=-12.7446
c24  Yes    :    x=207.5399 , y=2989.9734
c25  Yes    :    x=  3.0989 , y=  5.2016
c26  Yes    :    x=  2.5459 , y= -3.4841
c27  Yes    :    x= 55.8274 , y=-417.3156
c28  Yes    :    x=  7.5175 , y= 20.8659
c29  Yes    :    x=  1.8612 , y= -0.1770
Test for curve
Yes        c1    x=  4.0000 , y=  8.0000
Yes        c2*3    x=  6.0004 , y=-14.9012
Yes        c1+c2*4+c5*3    x=207.5399 , y=2989.9734
上述紅色是因浮點運算誤差所造成, 這點可能是一個奇特點, 因為數值變化很大. 修改 isPoint() 的 1e-9 改成 1e-6 就可改善判斷錯誤.

測試RSA 演算法的兩個質數(p及q),當加解密p或q時是否仍正確

用Euler theorem 證明 RSA 演算法時, 有個缺陷:當想要加解密的整數值若剛好等於所挑選的兩個質數p或q, 因為 n=p*q, 所以GCD(p,n)=p, 且GCD(q,n)=q, 並非是 1,造成Euler theorem並不適用. 用Java寫個程式驗證看看,在 1~100 整數裏面, 任取兩個質數來測試,到底加解密是否仍然可行,底下是原始碼,若加解密過程出現不相等情況就會印出錯誤訊息,經測試結果顯示,任何選取兩個小於 100 的質數乘積所形成的整數空間Z={x|x < (p*q)},當將它用RSA演算法來將該所選取的質數拿去加解密,仍然是適用的.或許可以直接就認定任取兩個質數(p,q)下面式子是成立的:
    p^(ed) mod n = p -- (1)
    q^(ed) mod n = q -- (2)
    其中 n = p*q,   φ = (p-1) * (q-1), (e*d) mod φ =1, 1 < e < φ
另外測試程式顯示結果可以發現 e 是有機會等於 d 的, 如果用這組key (e,d) 加解密到底安不安全的?我也不知道

 // RSA algorithm: // ( p, q) are two prime number // step1: n = p*q // step2: φ = (p-1) * (q-1) // step3: choose e : 1 < e < φ, GCD(e,φ) =1 to guarantee d exist // step4: calculate d: e*d mod n = 1 // step5: encode and decode // step5a: encode p and decode through: ( p^e )^d mod(n) to see if equals to p // step5b: encode q and decode through: ( q^e)^d mod(n) to see if equals to q import java.io.*; import java.util.List; import java.util.ArrayList; import java.math.BigInteger; import java.util.Base64; class TestRSA { private static boolean isPrime(int n) { if ( n % 2 == 0 ) return false; for(int i = 3; i * i <= n; i += 2) { if(n%i==0) return false; } return true; } private static boolean coprime(int a, int m) { int r; while ( m > 0 ) { r = a % m; a = m; m = r; } if( a==1 ) return true; else return false; } private static int modinv(int d, int n) { int a,q,r,tf; int tt = 0; int ts = 1; a = n; while ( d > 0 ) { q = a / d; tf = tt - q *ts ; tt = ts; ts = tf; r = a % d; a = d; d = r; } if (a != 1) return 0; return ( tt + n) % n; } public static void main(String[ ] args ) { List < Integer > prime = new ArrayList < Integer > ( ); for(int i = 1; i < 100; i ++ ) { if( isPrime(i) ) prime.add(i); }
for(int pp : prime ) { boolean test = true; BigInteger p = new BigInteger(""+pp); for(int qq:prime ) { if( pp==qq ) continue ; BigInteger q = new BigInteger(""+qq); int nn = pp*qq; // step 1: n = p*q BigInteger n = new BigInteger(""+nn); int φ = (pp-1)*(qq-1); // step 2: φ = (p-1) * (q-1) for(int ee=2; ee <φ; ee ++){// step 3: choose e such that 1 < e < φ and gcd(e,φ) =1
if( coprime(ee, φ) ) { BigInteger e = new BigInteger(""+ee); int dd = modinv(ee, φ); // step 4: d = modinv(e,φ); BigInteger d = new BigInteger(""+dd); BigInteger py = p.modPow(e,n); // step 5a: py = p^e mod n , px =py^d mod n BigInteger px = py.modPow(d,n); // decode py BigInteger qy = q.modPow(e,n); // step 5b: pq= q^e mod n , qx =pq^d mod n BigInteger qx = qy.modPow(d,n); // decode qy String text="(p,q) mod n =("+pp+","+qq+") mod "+nn+":enc-dec:("+px.toString()+","+qx.toString()+")\t(e,d) mod φ =("+ee+","+dd+") mod "+φ; if( ee==dd ) text = text +"\tsame (e=d)"; if( ! ( q.toString().equals(qx.toString()) && p.toString().equals(px.toString()) ) ) { text = text +"\tError!"; test = false; System.out.println(text); } } // key pair exist

} // for all the key pair (e,d) need to test } // for all the prime number need to test if( test ) System.out.println(""+pp+"->全正確"); else System.out.println(""+pp+"->錯誤!!"); }// end of test loop



} // end of main }// end of TestRSA

再談 RSA 演算法

• n = pq, where p and q are distinct primes.
• phi, φ(n) = (p-1)(q-1)
• e < n such that gcd(e, phi)=1
• d = e-1 mod phi.
• c = me mod n

m = cd mod n

參考影片:https://www.youtube.com/watch?v=QSlWzKNbKrU 
1. 取兩個質數p及q越大越好,假設(p,q)=(3,11)
2. 令 n = p *q = 3 * 11 = 33
3. 根據定理, 在 n 的整數空間中與n互質的個數等於(p-1)*(q-1),該值就是phi(n),等於就是刪掉3的倍數及11的倍數,因此 phi(n)=2*10=20個:
         1         2         3    
12345678901234567890123456789012
VV VV VV V  VV VV VV  V VV VV VV   
如果質數大於1時,通用公式為:(pi^ei -pi^ei')*...的總乘積,i是質數的個數
phi(3*11)=(3^(1)-3^(0))*(11^(1)-11^(0))=(3-1)*(11-1)=20
4. 選擇e < 33,且gcd(e,20)=1,也就是說與phi(n)互質的整數,例如3為例:

       因為 gcd(3,20) =1 故必定存在 d*3  = 1 mod 20
5. 求出 d = e^-1 mod phi(33)
    因為 phi(33) = 20 , e=3

而  (3*7) mod 20 =1 

    所以 d = 7
6. (33,3)就是公鑰,7就是私鑰,伺服器將(33,3)傳給使用者要求使用在指數空間編碼,例如將整數a用a^e mod n 來編碼,以a=4為例,公鑰e=4,n=33:
   4^3 mod 33 = 64 mod 33 = 31
   將4編碼成31,傳回伺服器
7.伺服器同樣使用指數解碼,但改用d來解碼,指數d=7,n=33
   31^7 mod 33 = (-2)^7 mod 33 = -128 mod 33 = -4*33+4 mod 33
   = 4 mod 33 --> 解碼等於 4,正好就是使用者傳的資料!
8.Euler's theorem: 其中 gcd(a,n)=1
  如果p是質數,與p互質(coprime)的個數1..p-1全與p互質=phi(p)= p-1,
  因此 a^(phi(p))=a^(p-1) = 1 mod p,兩邊同乘以 a
        a^p = a mod p
  上式就是Fermat's little theorem定理:當p是質數時對任何整數a而言
        a^p=a mod p

9.來看編碼/解碼程序中 e,d,n,phi(n)關係:
  指數e,d是由phi(n)整數空間推導出來,因為gcd(e,phi(n))=1,因此d必定存在,且互為倒數,換句話說:
  e*d = 1 mod phi(n) => ed=1+phi(n)*k, k 是整數
而 a^(ed) = a^(1+phi(n)*k)= a * a^(phi(n)*k)
上述尤拉定理 a^phi(n) = 1 mod n,將上式映射到 mod n 空間就會得到:
a^ed mod n = a * a^(phi(n)*k) mod n = a *1^k mod n = a mod n
詳細參考:https://zh.wikipedia.org/wiki/RSA%E5%8A%A0%E5%AF%86%E6%BC%94%E7%AE%97%E6%B3%95

• 將一個整數應映射到指數空間,加密時使用指數e:
     y= f(a) = a^e mod n
  解密時使用同樣的函數映射到指數空間,但改用指數d:
     c=f(y) = y^d mod n = a^(e*d) mod n = a mod n
  
• 有點要注意的是尤拉定理有個條件 gcd(a,n)=1 才會成立,而n=p*q是兩個質數乘積,只有當a=p或q兩點時會讓gcd(a,n)!=1,因此幾乎可以說a^phi(n)=1 mod n 是成立的除了 a=p, 及 a=q 之外.如果再將 a=p,a=q 單獨處理證明
• p^(ed) mod n = p 且 q^(ed) mod n = q 就可

10. 公鑰交換, Alice 與 Bob 約定好基底 r
     Alice 選擇私鑰 a, 計算公鑰  Pa = r^a mod n, 並將 Pa 傳給 Bob
     Bob 選擇私鑰 b, 計算公鑰  Pb = r^b mod n, 並將 Pb 傳給 Alice
     兩人收到後再用私鑰用同樣函式計算出 key
     Alice 收到 r^b mod n 計算出  (r^b)^a mod n  = r^(ba) mod n
     Bob 收到 r^a mod n 計算出  (r^a)^b mod n = r^(ab) mod n
     之後使用 a*b 來當 key 加密之後要傳輸的資料

11 這樣資料密鑰交換有沒問題呢? 如果駭客監聽到 Pa 或 Pb,而且他也知道r及n,想要破解a或b
             Pa = r^a mod n
             Pb = r^b mod n
兩條方程式2個變數(a,b),因為對稱的關係,一條方程式各有一個變數,看似好像有解,但是很難解
單純看一個方程式就好,因為它是mod n的運算關係式,已知 p, r, n想要求出 a,也就是說:
            p = r^a mod n  --- (1)
因為 p 是餘數,且包含r,a,n等全都是整數,因此可以假設另一個整數k可以讓下列關係式成立:
             r^a = p + k*n  --- (2)
上式經 mod n 運算後, 最後一項(k*n)會被消去成0,就等同 r^a mod n = p,因為係數全都是整數的關係,如果想要解方程式可以利用嘗試錯誤法勉強可以將上式(1)湊出答案.也就是說讓a從 1, 2,3 ... 一直不停的帶進方程式裏面,再求出兩邊的答案,當兩邊等式成立時就找到答案了.但萬一a,p,n及r都是很大的一個整數的話,比如說 1024 位元組的一個數字,大約是2^1024 次方數量級的一個大整數的話,可是要花費很長一段時間的.


12 ElGamal 加解密演匴法:
     1.Alice選一個質數=p及產生器=g,並挑選私鑰=a {2,3,...,p-2},用指數函數計算出h,公鑰就是(p,g,h)可以公開傳給Bob:
                 h = g^a mod p                --- (1)
     2.Bob利用Alice傳來的公鑰(p,g,h)將明文m加密,但他先隨機挑選一個亂數:
                  k      條件式: 滿足1 < k < p-1 且  gcd(k,p-1)=1
        也就是說k與p-1必需互質無公因數.同樣使用指數函數計算,並將明文m編碼成(r,s)兩個數傳回給Alice:
                 r  =        g^k  mod p       --- (2)
                 s = m *  h^k  mod p      --- (3)
     3. Alice收到Bob所傳來的的(r,s)後,用私鑰a並用公式計算就可解碼:
                 m ≡ s *  (r^−a)  mod p   --- (4)
                 (-a) = p - 1 - a                     (5)

     4. 理由是將(2),(3)式將帶入(4)式就可一窺之迷:
                    (4)                         (2)             (3)        
                 s * r^-a  mod p =  m  * h^k  *  (g^k)^-a   mod p,   再把(1)式 h = g^a 帶入
                 = m *  (g^a )^k   * (g^k)^-a   mod p  可以看到後面兩項互相消去了,得到明文 m
                 = m mod p

     5. 在第4式中有一個難題是 r^-a mod p = (r^a)^-1 mod p 是什麼? 它就是指數r^a的模p倒數.當r,a,p已知,因為 p 是質數,要求出(r^a)^-1 mod p並不難,但我可以用更快方式算出來:
             Euler's theorem說: 如果n是質數, Zn 空間與n互質的個數phi(n) = n-1個
                                    a^phi(n) = 1 mod n , 且GCD(a, n)= 1, 將 phi(n) 用 n-1 替換就得到:
                                    a^(n-1) = 1  mod n, 其中 n 是質數
             相同來看  r^p-1 ≡ 1 mod p , 其中 p 是質數, gcd(r,p)=1
              (r^a)*(r^-a) mod p = 1 mod p = r^(p-1) mod p
               r^(a + -a ) mod p = 1 mod p  = r^(p-1) mod p
              (a)+(-a) = p-1  -->     就可以很簡單推導出
                                (-a) = p - 1 - a   --- (5)

     6. 範例:選一個質數p=809,p-1=808=8*101,q=101,
                  (a) 計算 order=q=101的generator g=16  ???
                  (b) Alice 選擇密鑰 a =68, 計算 h = g^a mod p = 16^68 mod 809 = 46
                  (c) (p,g,h) = (809, 16, 46) 傳給 Bob, 密鑰a=68要保密
                  (d) Bob將明文m=100,選擇亂數k=89加密成:
                        r = g^k mod p=16^89 mod 809  = 342
                        s = m * h^k mod p = 100 * 46^89 mod 809 = 745
                  (f) Alice 收到r=342, s=745,
                       要計算 s*r^(-a) mod p 先計算出 r^(-a) mod p, a=68,  p=809,
                       而 (-a) = p -1 -a = 809 -1 -68 = 740
                       而  s=745, r = 342
                       最後 s * r^(-a) mod p = 745 * 342^740  mod p = 100 mod p
13 術語:
    multiplicative group Zp*
      Zp*  是刪除無 modulus inverse 的整數空間群組,也就是說裏面的整數元素與p全都互質
     例如 Z3的整數空間集合元素共有4 個 = {0,1,2,3}但  Z3* = {1,2}
     當 p 是質數時 Zp = {0} | Zp*  ,  Zp* ={1,2,3,...,p-1}
   order: 集合順序
     例如 |Z3| = 4 , 但 |Z3*|=2
     當元素可以分解成質數的乘積時 Z = a^n * b^p * c^q ...
     |Zp*| = ( a^n - a^n-1) * (b^p - b^p-1) * (c^q - c^q-1)
     例如 6 = 2 *3,  Zp={ 0 ,1 ,2, 3, 4, 5} =6
     |Zp*| = (2-1) * (3-1) = 1*2=2 ={1,5}
     當元素是質數時 |Zp*| = p -1 = phi(p)
   
     映射到指數空間並取餘數:令 Z11*={1,2,3,4,5,6,7,8,9,10}, f(a) = a^Zp* mod 11
     令 a = 2
          a^1 mod 11 = 2
          a^2 mod 11 = 4
          a^3 mod 11 = 8
          a^4 mod 11 = 5
          a^5 mod 11 = 10
          a^6 mod 11 = 9
          a^7 mod 11 = 7
          a^8 mod 11 = 3
          a^9 mod 11 = 6
          a^10 mod 11 = 1
          a^11 mod 11 = 2
    ord(2) = 10 個順序:  當 2^k mod 11 =1 時,最小的整數 k=10 就是 order
  primitive element:generator
  cyclic group: 當group裏面包含有最大 order 的元素時稱之為 cyclic group,該元素就稱為generator,例如 2^i 就可以產生 10 個元素, 最大order元素也包含在內,因此 2是 generator.
當 a 屬於 cycle group時 a^|Zp*| = 1 且 ord(a) 可以被 |Zp*|整除
     費曼定理 a^p mod p = a mod p -> a^(p-1) = 1 mod p
         p-1 = |Zp*|  --> a^|Zp*| =1         
        ord(1)  =1
        ord(2) = 10      --> generator
        ord(3) = 5
        ord(4) = 5
        ord(5) =5
        ord(6) = 10     --> generator
        ord(7) = 10     --> generator
        ord(8) = 10     --> generator
        ord(9) = 5
        ord(10) = 2
      可能的 order ={1,2,5,10} 必定可以被 k = 10 整除
       Z47* 當 a=5 時 |Z47*|=46, 因此 5 在Z47是 generator
     cyclic group 可以產生漂亮的 DPL, 而 generator 可以產生所有的元素.
                       






// 有趣的演算法: Double and Add  及 Square And Multiply
int testBit(int a, int b) {
        if(  a& (1 << b)  ) return 1;
else return 0;;
}
main()
{
int i, b;
int n;
int a=0xa2;
int k;
int e,p,c;

// double-and-add   DAA algorithm
b=0;
for( k=8; k>0 ;k-- ) {
if( testBit(a,k) ) {
b=1; // generator
break;
}
}
while( k-- ) {
b += b; // double
if( testBit(a,k)  ) b++; // add
}
// square-and-multiply   SAM algorithm e^p
e=5;
p=3;
for( k=8; k>0 ;k-- ) {
if( testBit(p,k) ) {
c = e; // generator
break;
}
}
while( k-- ) {
c *= c; // square
if( testBit(p,k)  ) c *= e; // multiply
}

printf("a=%d  b=%d\n",a,b);
printf("%d^%d=%d\n",e,p,c);

}