一個橢圓曲線方程式 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 就可改善判斷錯誤.
沒有留言:
張貼留言