2016年9月5日 星期一

橢圓曲線的轉換

參考文章: 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 其實是蒙哥馬利曲線的一個變種

2016年9月4日 星期日

橢圓曲線上奇怪的數學

一個橢圓曲線方程式 E(x, y) 在 x-y 平面滿足二元三次方程式:
            E(x, y) ::  y2 = x3 + ax + b        ---- (0)
假設 P(x1, y1) 是在曲線上的一點, 如果以 P(x1, y1) 當基準點畫一條斜率是 λ 的直線,這條直線用方程式描述:
           -y = λ ( x1 - x ) - y1                   --- (1)
重新整理:
             y =   λ x + c             c  是常數 = y1 - λ x1    --- (2)
如果將式(2)帶入式(0)以就可得出交點,先看等式(0)右邊,如果將通過P的直線帶入時,也就是說:
          ( λ x + c )2 = ( x3 + a x + b)     --- (3)
上式是一個 x 三次方程式,因此 x 會有 3 個根(root),也就是說會有 3 個交點,因為P(x1,y1) 在線上,是已知的根(root), 另外兩交點分別是 Q(x2,y2) 與 R(x3,y3) 滿足曲線方程式:
          y22 = x23 + a x2 + b,  y2 = λ(x2 - x1) + y1   --- (4)
          y32 = x33 + a x3 + b,  y3= λ(x3 - x1) + y1   --- (5)
從方程式(0)來看, 左邊是 y2 , 可以知道當曲線上存在有一點 P(x1, y1), 以同樣的 x1根, 必定有另一點(x1, -y1) 也在曲線上, 因為(-y1)2 = (y1)2 , 如果以座標軸來看, 就是以 y=0 為軸當鏡面的鏡像點, 用 -P = (x1, -y1) 來表示 同樣曲線上的另外兩交點 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 ,當全加一起就會成為零點. 用同樣邏輯可否推廣到其它二元四次  y2 = f(x) 方程式甚至更高次項的方程式與直線交點的根解(root)呢?也就是說將全部根的點全加在一起是否就是零點?待研究...
如果將第(8)式兩邊都加上 R 的鏡像點 -R:
          (P + Q + R) + (-R) = O + (-R)
左邊應用第(7)式讓 R 消滅, 右邊應用第(6)式任何點加上零點還是原來的點就得到:        
           P + Q = -R              --- (9)
方程式(9)描述的是將兩交點 P=(x1, y1)及 Q=(x2, y2), 斜率λ的直線與橢圓曲線其中兩交點, 經過相加運算得出新座標 -R = (x3, -y3) 正好就是最後交點 R(x3, y3) 的鏡像點(x-y 平面上以 y=0 為軸當鏡面), 後續可以推導出座標是:
            x3 = λ2 - x1 - x2           --- (10)
          -y3 = λ(x1 - x3) - y1       --- (11)
也就是說:
          P + Q = (x3, -y3)
已知兩點(x1, y1)及(x2, y2), 相加後的新座標(x3, -y3),還是座落在橢圓曲線上.
λ  是直線的斜率 dy/dx, 對於兩點(x1, y1)及(x2, y2)形成的直線而言:
            λ = (y2 - y1) / (x2 - x1)      其中 x1 != x2
但是當 P = Q 時 x1 = x2, y1 = y2, 導致分母為零無法除,因此不能使用上述 λ 的公式,也就是說當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) = (3x2  + a) / 2y   ---(12)
帶入切點座標(x1, y1)得出斜率 λ':
           λ'(x, y) = (3x12  + a)/ 2y1          --- (13)
上述意思說當 Q = P 意味著是將 Q 點逼近橢圓曲線的切點 P(x1, y1)處.
該條切線所形成的另外交點 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)
        λ' = (3x2 + a)/2y
        x' = λ'2 - x - x = λ'2 - 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) 運算 P2, P3, P4, P5, P6 ...的座標, 把他映設到橢圓曲線空間 P, 例如:
     G  = n*P         n >0, n 是整數, multiply 利用上述 DAA 演算法運算 n*P
     G2 = G*G     = n2P,  square:  n*n P
     G3 = (G2)*G = n3P,  square and multiply n2*n P

     G4 = (G2)2    = n4P,  square and square   (n2)2 P

     G5 = (G2)2*G= n5P,  square then square and multiply (n2)2*n P

     G6 = (G2*G) 2 *G = n6P,  square and multiply then square and multiple (n2*n)2*n P

       ...
附錄:  推導 x 與斜率 λ 的關係式:
將直線方程式 (1) y 帶入曲線方程式 E(x,y), 與假設的 3 個根(x1, x2, x3)所形成的三次方程式, 只要比較 x2 項的係數, 就可得到斜率 λ 與 x1, x2, x3 關係式:
        ( λ x + c )2  - (x3 + a x + b ) = (x - x1)(x - x2)(x - x3) = 0
         x3 = λ2 - x1 - x2
-R 的 y 座標是 - y3, 直接帶回直線方程式(1)解出:
        -y3 = λ(x1 - x3) - y1
詳細參考網路上的文章: 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:    y2 = x3 + 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 就可改善判斷錯誤.

2016年9月3日 星期六

測試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

2016年9月2日 星期五

再談 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)傳給使用者要求使用在指數空間編碼,例如將整數aa^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: a^{\varphi (n)} \equiv 1 \pmod{n}, 其中 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);

}


2016年8月29日 星期一

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

橢圓曲線方程式:
   G(x,y) :: y2 = x3  + a*x + b   其中 4a3 + 27b2 != 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), 經由係數比較就可得出 x3 = s2 - x2 - x1的關係式
       曲線 G(x,y) :: y2 = x3  + a*x + b               
       直線 L(x, y) :: y = s*x + c
條件(一) 當 P != Q 時, 使用一般直線斜率 s 的加法運算求出 R:
      s = ( y2 - y1 ) / (x2 - x1 )
      x3 = s2 - x2 - x1
      y3 = s*(x2-x3) -  y2
條件(二) 當 P == Q 時,  x1=x2, 且 y1=y2, 運算斜率 s 時, 導致除零的錯誤, 要改用雙倍切線斜率運算:
      s = ( 3 *x12 +a )/ ( 2y1)
      x3 = s2 - 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(xa, ya)
3. Bob 收到後, 也隨機挑選一個密鑰 b, 算出 ECC 上的一點當自己的公鑰 PB, 同時把要加密的文件訊息 m 用 Alice 的公鑰編碼成 R, 為了能適用 ECC 的加法運算, m 必須映射到橢圓曲線上, 有一種方式是將 m 當成 x 帶入曲線公式解出 y 得出座標(xm,ym), 再計算編碼過的座標 R = m + b*PA,  連同 PB 一起回傳給 Alice, 讓她用私鑰將 m 解讀出來, 通信過程就算公開了 G(xa,ya) 與 G(xb,yb) 的座標點, 也很難猜到 G(xab,yab) 的座標點 , 因此外人解出 m 的機會很渺茫, 只有兩邊通信方知道共用金鑰 ab 座標點是G(xab, yab), 直接卸下(減掉)該因子就能快速還原文件訊息 m:
     座標 PB:: b*G(x, y) = G(xb, yb)
     座標 R:: m + b*PA(x, y) = m + G(xab, yab) = m + G(xba, yba)
4. Alice 收到 PB, R後, 就用自己的私鑰 a 去計算  R - a*PB 得出 m
     座標 (R  -  a*PB) :: [m + G(xab,yab)]  - G(xab,yab) = m + O = m
5.  ECC 的鏡射點相加是為零點 O:
     - G(xab, yab)  + G(xab, yab)  = O
6. 鏡射點, 實際上只要反轉 y 座標為負值:
     - G(xab, yab) = G(xab, - yab)
   m = R - G(xab, yab) = R + G(xab, - yab) 利用ECC 的加法運算就可算出座標值 (xm, ym), 其中的 xm 就是解出來的訊息


2016年8月23日 星期二

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>

2016年8月8日 星期一

使用 gulp.js 讓流程自動化

gulp 中文語意是:吞嚥, 我覺得翻譯成 '嗑' 比較貼切. 它就是將東西吞下去. 當我們寫了許多Javascript 要將它整合成一個檔案並將它醜化,需要通過工具程式及一些流程來完成. 當然用Makefile 就可做到, 但寫個 Makefile 很麻煩,要考慮檔案的相關性,用 gulp.js 很直覺就能直接使用,且他是用 Javascript 語法, 對於寫 Javascript 來說相當便利.但 gulp.js 可不只這樣, 他還可以讓各個流程同時進行.先簡單以上述例子來練習. 先用 sudo npm -g 來安裝 gulp 及 requirejs , linux mint 會將它放在 /usr/lib/node_modules 目錄裏面:
                       sudo npm install gulp requirejs -g
再切換到專案目錄,安裝所需要的工具模組(例如串接檔案成單一個, 去醜化等):
                       cd /home/mint/myproject
                       npm install gulp-concat gulp-uglify
接著寫一個 gulpfile.js 讓 gulp 根據文件內容去執行.下面例子是將 1.js, 2.js, 3.js 結合成單一檔案 my.js, gulp.task 用來自訂一個流程, 第一個參數是自定的名稱, 第2個參數是匿名函數,實現
流程化, 在流程化裏面可以直接套用 gulp 的函數諸如 .src .pipe .dest 等等. 通過管線(.pipe)就可以將輸出/入檔案相互串接,達到想要完成的目的(流程).

         // gulpfile.js
         var gulp = require('gulp');
         var concat = require('gulp-concat');
         var uglify = require('gulp-uglify');
         gulp.task('combine', function() {
                  // .src( string_array ) -> .pipe(function( )  ) -> .pipe( function( ) ) -> .pipe( function( ) )
         gulp.src(['1.js', '2.js', '3.js']).pipe(concat('my.js')).pipe(uglify()).pipe(gulp.dest('build/'));
         });

最後要執行該流程時,將流程名稱提供給 gulp, 鍵入以下命令:
         gulp combine
在 build/ 目錄內就會生成一個醜化過後的檔案  my.js
對於其他檔案例如 index.htm 想要複製到 build 目錄裏面時, 也可以自定一個流程將檔案複製過去:
         gulp.task('copy', function() {
                  gulp.src('index.htm').pipe(gulp.dest('build/'));
         });
將所有的流程集合到一個字串陣列,並命名一個 gulp 內定名稱 'default', 就可以讓上述兩個流程同時進行:
          gulp.task('default', ['combine', 'copy']);
完整的內容如下:
         // gulpfile.js
         var gulp = require('gulp');
         var concat = require('gulp-concat');
         var uglify = require('gulp-uglify');
         gulp.task('default', ['combine', 'copy']);
         gulp.task( 'combine', function() {
                gulp.src(['1.js', '2.js', '3.js']).pipe(concat('my.js')).pipe(uglify()).pipe(gulp.dest('build/'));
         });
         gulp.task('copy', function() {
                  gulp.src('index.htm').pipe(gulp.dest('build/'));
         });
要執行上述並行流程並不需要再提供流程名稱, gulp 會先找到 qulpfile.js, 再到檔案裏面找尋 default 這個流程名稱,解析它並自動執行,因此只要鍵入以下命令就可:
          gulp

另外在 linux 下面透過 cat 其實就可以將檔案合併, 透夠管線, 一條指令就可以達到上述同樣功能了:
          cat 1.js 2.js 3.js | ugilfyjs > my.js