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

質數(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);

}

理解橢圓曲線密碼學(Elliptic Curve Cryptography)

橢圓曲線方程式:
   G(x,y) :: y² = x³  + a*x + b   其中 4a³ + 27b² != 0
P Q R 是在橢圓曲線的3個點, 並且都同在一個直線上:
         P + Q + R = O (橢圓曲線零點)
上式並不是算術運算的加法, 只是一個術語.
         因此 P + Q = - R
負數代表的是鏡像(mirror), 因為橢圓曲線對於 Y 軸而言是一個對稱曲線, 因此可以想像成 R 與 -R 是分別在 Y 正負軸上面等距離的點
已知 P=G(x1,y1), Q=G(x2,y2) 求解 R=G(x3,y3), 帶入上述橢圓曲線G(x,y)及直線方程式L(x,y), 經由係數比較就可得出 x₃ = s² - x₂ - x₁的關係式
       曲線 G(x,y) :: y² = x³  + a*x + b               
       直線 L(x, y) :: y = s*x + c
條件(一) 當 P != Q 時, 使用一般直線斜率 s 的加法運算求出 R:
      s = ( y2 - y1 ) / (x2 - x1 )
      x3 = s² - x2 - x1
      y3 = s*(x2-x3) -  y2
條件(二) 當 P == Q 時,  x1=x2, 且 y1=y2, 運算斜率 s 時, 導致除零的錯誤, 要改用雙倍切線斜率運算:
      s = ( 3 *x1² +a )/ ( 2y1)
      x3 = s² - 2x1
      y3 = s*(x1  -x3) - y1
條件(三) n(n>1)倍G運算, 使用雙倍運算及累加運算(DAA),演算法用2進位法:
      (1) 找出 n 當中最高位元是 1 時, 並開始初始化 R
               R = G ; // 初始化 R 等於 G(x,y)
               k = 要處理的位元數                  
      (2) while( k-- > 0 ) { // 累進運算
               R= optDouble(R);// 先將點作兩倍運算, 用上述條件(二)的雙倍運算
               if(  n.testBit(k) ) R=optINC(R, G); //當n測試到位元 1 時執行上述條件(一)的加法運算

           } //  重複步驟2直到 k = 0 為止
ECC 的對稱加密通信演算法:
當有兩端要通訊時, 兩邊各選擇一個整數(通常是大整數)當成私鑰, 將它乘上 G就形成 public key例如 pubA = a * G, pubB =b * G, 兩邊交換後, A 將收到的 pubB 乘上自己的私鑰,B收到pubA也乘上自己的私鑰, 結果 a*pubB = a*b*G = b *a * G = b * pubA 就是共用金鑰( a*b=b*a)所對應 ECC 上的一點, 把該點的資訊當作 AES 加解密之密碼, 傳輸的資料經 AES 加密過, 對外人來說是一串亂碼, 就算中間 pubA 或 pubB 被擷取, 也很難算出 a, b , a*b 來. 因為 a, b 對局外人來講是未知(大整)數, 只有通信方兩端知道自己的私鑰(a, b)及共同金鑰(ab)所對應的點ab *G(x,y), 實際上兩邊通信方也很難推算出共用金鑰 ab, 因此 ECC 可以讓通信保密

實際上觀念很簡單, 使用小整數來說:
1. Alice 與 Bob 共同約定使用相同ECC編碼,橢圓曲線(EC)上的點假設是 G(x,y)當作初始點.
2. Alice 使用 3 當作密鑰, Alice 將 3*G(x,y) 映設到新的橢圓曲線上一點假設它是 G(xp,yp), 實際上G(xp,yp)就是 Alice 的公開鑰匙, 將它傳給 Bob.
3.  Bob 使用 5 當作密鑰, 並將 5*G(x,y) 映設到新的橢圓曲線上一點假設它是 G(xq,yq), 實際上G(xq,yq)就是 Bob 的公開鑰匙, 將它傳給 Alice.
4. Alice 和 Bob 分別用自己的私鑰乘上所收到點的資訊並重新映設到新的橢圓曲線上一點,例如 Alice收到G(xq,yq), 而 Bob收到G(xp,yp), 實際上乘上自己的鑰匙後得到:
            Alice: 3*G(xq,yq) = 3*5*G(x,y) = 15*G(x,y)
            Bob:   5*G(xp,yp) = 5*3*G(x,y) = 15*G(x,y)
因此 3*5 稱為可以稱之共同金鑰, 而15*G(x,y) 是 ECC 上的一個點, 雖然兩邊很難算出 15 這個數字, 但可以利用該點共同座標當作訊息替後續要傳送的資料加以編碼(例如用 SHA-256 轉成亂數當成 AES-256 的密碼將訊息編碼). 重頭到尾, 外人無法得知 3, 5, 15 這個數字, 外人只能察覺到兩人在交換 ECC 上點的資訊(公開鑰匙), 如果要反推 3,5,15, 便需要從頭 G(x,y)一一去嘗試, 一旦試到 15 才能找到答案.但密碼學裡很重要的觀念是'數'的量很龐大,當我們使用超級大整數時例如 256 bit來將私鑰編碼時就有 2的256次方個組合, 這根本是一個無法想像的數字.而 ECC 點的運算,必需從頭到尾一一加以運算, 就算知道起點和終點,仍然很難分解出該私鑰點所對應的之相對位置. 更遑論是大整數的運算了. 當駭客無法破解密碼時,唯一可以用的方法是暴力窮舉法一一測試, 但由於是大整數的關係, 要一直不斷嘗試, 真的非常困難. 這就是 ECC 相當安全的重要原因.

ECC 使用的 ElGamal 非對稱加解密通信演算法:
1. 兩人約定用 ECC 的點 G(x, y) 當作通訊協定
2. Alice 挑選一個私鑰  a, 算出 ECC 上的一點座標 PA 傳給 Bob, PA 就是 Alice 公鑰:
     座標 PA:: a*G(x, y) = G(x_a, y_a)
3. Bob 收到後, 也隨機挑選一個密鑰 b, 算出 ECC 上的一點當自己的公鑰 PB, 同時把要加密的文件訊息 m 用 Alice 的公鑰編碼成 R, 為了能適用 ECC 的加法運算, m 必須映射到橢圓曲線上, 有一種方式是將 m 當成 x 帶入曲線公式解出 y 得出座標(x_m,y_m), 再計算編碼過的座標 R = m + b*PA,  連同 PB 一起回傳給 Alice, 讓她用私鑰將 m 解讀出來, 通信過程就算公開了 G(x_a,y_a) 與 G(x_b,y_b) 的座標點, 也很難猜到 G(x_ab,y_ab) 的座標點 , 因此外人解出 m 的機會很渺茫, 只有兩邊通信方知道共用金鑰 ab 座標點是G(x_ab, y_ab), 直接卸下(減掉)該因子就能快速還原文件訊息 m:
     座標 PB:: b*G(x, y) = G(x_b, y_b)
     座標 R:: m + b*PA(x, y) = m + G(x_ab, y_ab) = m + G(x_ba, y_ba)
4. Alice 收到 PB, R後, 就用自己的私鑰 a 去計算  R - a*PB 得出 m
     座標 (R  -  a*PB) :: [m + G(x_ab,y_ab)]  - G(x_ab,y_ab) = m + O = m
5.  ECC 的鏡射點相加是為零點 O:
     - G(x_ab, y_ab)  + G(x_ab, y_ab)  = O
6. 鏡射點, 實際上只要反轉 y 座標為負值:
     - G(x_ab, y_ab) = G(x_ab, - y_ab)
   m = R - G(x_ab, y_ab) = R + G(x_ab, - y_ab) 利用ECC 的加法運算就可算出座標值 (x_m, y_m), 其中的 x_m 就是解出來的訊息

Android 讀取加速度, 亮度, 接近等三種感應器

// 主程式: MainActivity.java

package com.example.mint.accelerometer;

import android.hardware.Sensor;
import android.hardware.SensorEvent;
import android.hardware.SensorEventListener;
import android.hardware.SensorManager;
import android.os.Bundle;
import android.support.v7.app.AppCompatActivity;
import android.widget.TextView;

public class MainActivity extends AppCompatActivity  implements SensorEventListener {
    SensorManager sensorMan;
    TextView textview;
    String accMeter, lightMeter, proximityMeter;
    public void onAccuracyChanged(Sensor arg0, int arg1) {
    }
    public void onSensorChanged(SensorEvent event) {
        switch( event.sensor.getType() ) {
            case Sensor.TYPE_ACCELEROMETER:
                accMeter = "\n加速感應器: x=" + event.values[0] + "    y=" + event.values[1] + "    z=" + event.values[2];
                break;
            case Sensor.TYPE_PROXIMITY:
                proximityMeter = "\n接近感應器:" + event.values[0];
                break;
            case Sensor.TYPE_LIGHT:
                lightMeter = "\n亮度感應器:" + event.values[0];
                break;
        }
        textview.setText(accMeter+proximityMeter+lightMeter);
    }
    @Override    public void onCreate(Bundle savedInstanceState) {
        super.onCreate(savedInstanceState);
        setContentView(R.layout.activity_main);
        textview = (TextView)findViewById(R.id.meterInfo);
        sensorMan = (SensorManager) getSystemService(SENSOR_SERVICE);
    }
    protected void onResume() {
        super.onResume();
        sensorMan.registerListener(this, sensorMan.getDefaultSensor(Sensor.TYPE_ACCELEROMETER), SensorManager.SENSOR_DELAY_NORMAL);
        sensorMan.registerListener(this, sensorMan.getDefaultSensor(Sensor.TYPE_PROXIMITY), SensorManager.SENSOR_DELAY_NORMAL);
        sensorMan.registerListener(this, sensorMan.getDefaultSensor(Sensor.TYPE_LIGHT), SensorManager.SENSOR_DELAY_NORMAL);
    }

    @Override    protected void onPause() {
        super.onPause();
        sensorMan.unregisterListener(this);
    }

}

// activity_main.xml

<?xml version="1.0" encoding="utf-8"?>
<RelativeLayout xmlns:android="http://schemas.android.com/apk/res/android"
    xmlns:tools="http://schemas.android.com/tools"
    android:layout_width="match_parent"
    android:layout_height="match_parent"
    android:paddingBottom="@dimen/activity_vertical_margin"
    android:paddingLeft="@dimen/activity_horizontal_margin"
    android:paddingRight="@dimen/activity_horizontal_margin"
    android:paddingTop="@dimen/activity_vertical_margin"
    tools:context="com.example.mint.accelerometer.MainActivity">
    <TextView
        android:layout_width="wrap_content"
        android:layout_height="wrap_content"
        android:text="Hello World!"
        android:id="@+id/meterInfo" />
</RelativeLayout>