Math Pro 數學補給站 » 數學軟體 » 用Maxima學密碼學－橢圓曲線

bugmens 發表於 2019-8-6 20:52

用Maxima學密碼學－橢圓曲線

橢圓曲線加密請見連結介紹，本文僅就例2、例3以\(maxima\)示範計算過程。
第7章橢圓曲線密碼系統，[url]https://math.pro/db/attachment.php?aid=5248&k=23e470df1ebb8f1f5ca07d9f3c723d04&t=1565140250[/url]
例子2：在有限體\(F_{23}\)之下，取橢圓曲線\(E\)：\(y^2=x^3+16x+10\)上的兩點\(P=(18,14)\)及\(Q=(5,10)\)。若\(P+Q=R=(x_3,y_3)\)，則\(R=\)?
例子3：承例子2，若\(P+P=2P=R=(x_3,y_3)\)，則\(R=\)?



[color=green]y^2=x^3+ax+b (mod p)
y^2=x^3+16x+10 (mod 23)的橢圓曲線[/color]
[color=red](%i1)[/color]　[color=blue][p,a,b]:[23,16,10];[/color]
[color=red](%o1)[/color]　\(\left[ 23,16,10 \right]\)

[color=green]橢圓曲線點相加[/color]
[color=red](%i2)[/color]
[color=blue]PointAdd(P,Q):=block
([lambda,num,denom,inv_denom,x3,y3],
if P=Q then
  (print("λ=(3*",x_1^2,"+a)/(2*",y_1,")=",
            "(3*",P[1],"^2+",a,")/(2*",P[2],")=",
            lambda: (3*P[1]^2+a)/(2*Q[2]))
  )
else if P[1]=Q[1] and P[2]#Q[2] then
  (print(" P,Q兩點x座標相同,相加結果為無窮遠點"),
   return([inf,inf])
  )
else
  (print("λ=(",y_2,"-",y_1,")/(",x_2,"-",x_1,")=",
            "(",Q[2],"-",P[2],")/(",Q[1],"-",P[1],")=",
            lambda: (P[2]-Q[2])/(P[1]-Q[1]))
  ),
num:num(lambda),/*分子*/
denom:denom(lambda),/*分母*/
print(1/denom,"=",inv_denom:inv_mod(denom,p),"(mod",p,")"),
print("λ=",lambda,"=",num,"*",inv_denom,"=",lambda:mod(num*inv_denom,p),"(mod",p,")"),
print(x_3,"=",λ^2,"-",x_1,"-",x_2,"=",lambda^2,"-",P[1],"-",Q[1],"=",
         x3:lambda^2-P[1]-Q[1],"=",x3:mod(x3,p),"(mod",p,")"),
print(y_3,"=",λ,"(",x_1,"-",x_3,")-",y_1,"=",lambda,"(",P[1],"-",x3,")-",P[2],"=",
         y3:lambda*(P[1]-x3)-P[2],"=",y3:mod(y3,p),"(mod",p,")"),
print("[",x_3,",",y_3,"]=",[x3,y3]),
return([x3,y3])
)$[/color]

[color=green]P=(x1,y1)=(18,14)
Q=(x2,y2)=(5,10)[/color]
[color=red](%i4)[/color]
[color=blue]P:[18,14];
Q:[5,10];[/color]
[color=red](P)[/color]　\(\left[18,14\right]\)
[color=red](Q)[/color]　\(\left[5,14\right]\)

[color=green]計算P+Q=R=(x3,y3)[/color]
[color=red](%i5)[/color]　[color=blue]PointAdd(P,Q);[/color]
\(\displaystyle \lambda=(y_2-y_1)/(x_2-x_1)=(10-14)/(5-18)=\frac{4}{13}\)
\(\displaystyle \frac{1}{13}=16\pmod{23}\)
\(\displaystyle \lambda=\frac{4}{13}=4 \cdot 16=18\pmod{23}\)
\(x_3=\lambda^2-x_1-x_2=324-18-5=301=2\pmod{23}\)
\(y_3=\lambda(x_1-x_3)-y_1=18(18-2)-14=274=21\pmod{23}\)
\(\left[x_3,y_3\right]=\left[2,21\right]\)
[color=red](%o5)[/color]　\(\left[2,21\right]\)

[color=green]計算P+P=2P=R=(x3,y3)[/color]
[color=red](%i6)[/color]　[color=blue]PointAdd(P,P);[/color]
\(\displaystyle \lambda=(3*x_1^2+a)/(2*y_1)=(3*18^2+16)/(2*14)=\frac{247}{7}\)
\(\displaystyle \frac{1}{7}=10\pmod{23}\)
\(\displaystyle \lambda=\frac{247}{7}=247 \cdot 10=9\pmod{23}\)
\(x_3=\lambda^2-x_1-x_2=81-18-18=45=22\pmod{23}\)
\(y_3=\lambda(x_1-x_3)-y_1=9(18-22)-14=-50=19\pmod{23}\)
\(\left[x_3,y_3\right]=\left[22,19\right]\)
[color=red](%o6)[/color]　\(\left[22,19 \right]\)

bugmens 發表於 2019-8-25 22:08

例子1：有限體\(F_{23}\)之下，點\(P=(0,1)\)是橢圓曲線\(E\)：\(y^2=x^3+12x+1\)的生成數。
[table=30%][tr][td]P=(0, 1)
2P = (13, 13)
3P = (5, 5)
4P = (3, 15)
5P = (6, 17)
6P = (19, 2)
7P = (17, 9)
8P = (18, 0)[/td][td]9P = (17, 14)
10P = (19, 21)
11P = (6, 6)
12P = (3, 8)
13P = (5, 18)
14P = (13, 10)
15P = (0, 22)[/td][/tr][/table]



[color=green]y^2=x^3+ax+b (mod p)
y^2=x^3+12x+1 (mod 23)的橢圓曲線[/color]
[color=red](%i1)[/color]　[color=blue][p,a,b]:[23,12,1];[/color]
[color=red](%o1)[/color]　\(\left[ 23,12,1 \right]\)

[color=green]橢圓曲線點相加[/color]
[color=red](%i2)[/color]
[color=blue]PointAdd(P,Q):=block
([lambda,num,denom,inv_denom,x3,y3],
if P=Q then
  (lambda: (3*P[1]^2+a)/(2*Q[2])
  )
else if P[1]=Q[1] and P[2]#Q[2] then
  (print(" P,Q兩點x座標相同,相加結果為無窮遠點"),
   return([inf,inf])
  )
else
  (lambda: (P[2]-Q[2])/(P[1]-Q[1])
  ),
lambda:num(lambda)*inv_mod(denom(lambda),p),
x3:mod(lambda^2-P[1]-Q[1],p),
y3:mod(lambda*(P[1]-x3)-P[2],p),
return([x3,y3])
)$[/color]

[color=green]生成數P=[0,1][/color]
[color=red](%i3)[/color]　[color=blue]P:[0,1];[/color]
[color=red](P)[/color]　\(\left[0,1 \right]\)

[color=green]計算P,2P,3P,...,15P，結果放在L[/color]
[color=red](%i6)[/color]
[color=blue]L:[P]$
for k:2 thru 15 do
  (L:append(L,[PointAdd(last(L),P)])
  )$
L;[/color]
[color=red](%o6)[/color]　\([[0,1],[13,13],[5,5],[3,15],[6,17],[19,2],[17,9],[18,0],[17,14],[19,21],[6,6],[3,8],[5,18],[13,10],[0,22]]\)

[color=green]P,2P,...,15P位置再加上標籤[/color]
[color=red](%i7)[/color]　[color=blue]Label:create_list(append([concat(i,"P")],(L[ i ]+[1,0])),i,1,length(L));[/color]
[color=red](Label)[/color]　\(\matrix{[[1P,1,1],[2P,14,13],[3P,6,5],[4P,4,15],[5P,7,17],[6P,20,2],[7P,18,9],[8P,19,0],\cr
[9P,18,14],[10P,20,21],[11P,7,6],[12P,4,8],[13P,6,18],[14P,14,10],[15P,1,22]]}\)

[color=green]將P,2P,...,15P各點坐標畫出來[/color]
[color=red](%i8)[/color]
[color=blue]draw2d(user_preamble="set size square; set grid; set xtics 1; set ytics 1;",
              grid=true,
              point_type=filled_circle,
              points(L),
              apply(label,Label));[/color]
[color=red](%o8)[/color]　\([gr2d(points,label)]\)
[img]https://math.pro/db/attachment.php?aid=5266&k=8f533f53c25c57ec59b8fe23de9006d6&t=1566741731&noupdate=yes¬humb=yes[/img]

bugmens 發表於 2019-9-1 17:59

前面兩篇文章利用\(maxima\)程式示範橢圓曲線加法(\(P+Q\))和兩倍(\(P+P\))的運算過程，以及藉由\(\pmod{p}\)將橢圓曲線上的每個點全部打散，讓橢圓曲線離散對數問題(\(ECDLP\))變得難以計算。
其實\(maxima\)已有內建的橢圓曲線程式庫，本篇文章將前面兩篇重新用內建橢圓曲線指令計算。


[color=green]要先載入elliptic_curves.mac才能呼叫ec_set_curve,ec_point_p,ec_add[/color]
[color=red](%i1)[/color]　[color=blue]load("elliptic_curves.mac");[/color]
[color=red](%o1)[/color]　[i]C:/maxima-5.43.0/share/maxima/5.43.0/share/contrib/elliptic_curves/elliptic_curves.mac[/i]

[color=green]y^2=x^3+ax+b (mod p)
y^2=x^3+16x+10 (mod 23)的橢圓曲線[/color]
[color=red](%i3)[/color]
[color=blue][p,a,b] : [23,16,10];
ec_set_curve(p,a,b);[/color]
[color=red](%o2)[/color]　\([23,16,10]\)
[color=red](%o3)[/color]　\(true\)

[color=green]P,Q兩點座標[/color]
[color=red](%i5)[/color]
[color=blue]P:[18,14];
Q:[5,10];[/color]
[color=red](P)[/color]　\([18,14]\)
[color=red](Q)[/color]　\([5,10]\)

[color=green]檢查P,Q有沒有在橢圓曲線上[/color]
[color=red](%i7)[/color]
[color=blue]ec_point_p(P);
ec_point_p(Q);[/color]
[color=red](%o6)[/color]　\(true\)
[color=red](%o7)[/color]　\(true\)

[color=green]計算P+Q[/color]
[color=red](%i8)[/color]　[color=blue]ec_add(P,Q);[/color]
[color=red](%o8)[/color]　\([2,21]\)

[color=green]計算P+P[/color]
[color=red](%i9)[/color]　[color=blue]ec_add(P,P);[/color]
[color=red](%o9)[/color]　\([22,19]\)

－－－－－－－－－－－－－－－－－－－－－－
[color=green]要先載入elliptic_curves.mac才能呼叫ec_set_curve,ec_mult[/color]
[color=red](%i1)[/color]　[color=blue]load("elliptic_curves.mac");[/color]
[color=red](%o1)[/color]　[i]C:/maxima-5.43.0/share/maxima/5.43.0/share/contrib/elliptic_curves/elliptic_curves.mac[/i]

[color=green]y^2=x^3+ax+b (mod p)
y^2=x^3+12x+1 (mod 23)的橢圓曲線[/color]
[color=red](%i3)[/color]
[color=blue][p,a,b]:[23,12,1];
ec_set_curve(p,a,b);[/color]
[color=red](%o2)[/color]　\([23,12,1]\)
[color=red](%o3)[/color]　\(true\)

[color=green]生成數P=[0,1][/color]
[color=red](%i4)[/color]　[color=blue]P:[0,1];[/color]
[color=red](P)[/color]　\([0,1]\)

[color=green]用ec_mult(k,P)計算P,2P,3P,...,15P，結果放在L[/color]
[color=red](%i5)[/color]　[color=blue]L:makelist(ec_mult(k,P),k,1,15);[/color]
[color=red](L)[/color]　\([[0,1],[13,13],[5,5],[3,15],[6,17],[19,2],[17,9],[18,0],[17,14],[19,21],[6,6],[3,8],[5,18],[13,10],[0,22]]\)

[color=green]P,2P,...,15P位置再加上標籤[/color]
[color=red](%i6)[/color]　[color=blue]Label:create_list(append([concat(i,"P")],(L[ i ]+[1,0])),i,1,length(L));[/color]
[color=red](Label)[/color]　\(\matrix{[[1P,1,1],[2P,14,13],[3P,6,5],[4P,4,15],[5P,7,17],[6P,20,2],[7P,18,9],[8P,19,0],\cr
[9P,18,14],[10P,20,21],[11P,7,6],[12P,4,8],[13P,6,18],[14P,14,10],[15P,1,22]]}\)

[color=green]將P,2P,...,15P各點坐標畫出來[/color]
[color=red](%i7)[/color]
[color=blue]draw2d(user_preamble="set size square; set grid; set xtics 1; set ytics 1;",
              grid=true,
              point_type=filled_circle,
              points(L),
              apply(label,Label));[/color]
[color=red](%o7)[/color]　\([gr2d(points,label)]\)
[img]https://math.pro/db/attachment.php?aid=5266&k=8f533f53c25c57ec59b8fe23de9006d6&t=1566741731&noupdate=yes%C2%AChumb=yes[/img]

bugmens 發表於 2019-9-7 19:39

上一篇文章以\(y^2=x^3+12x+1 \pmod{23}\)，\(P=[0,1]\)為生成數計算\(1P,2P,\ldots,15P\)，再加上無窮遠點總共16個點。

想要知道有限體上橢圓曲線的點個數，直觀的方法是將\(x=0,1,\ldots,p-1\)代入橢圓曲線方程式\(y^2=x^3+ax+b \pmod{p}\)，以\(jacobi\)符號判斷同餘方程式是否有\(y\)的解，將有\(y\)的解個數加起來就是答案。

關於\(jacobi\)符號請看\(wiki\)介紹
[url]https://zh.wikipedia.org/wiki/%E9%9B%85%E5%8F%AF%E6%AF%94%E7%AC%A6%E5%8F%B7[/url]

但橢圓曲線點個數計算不需要將各個\((x,y)\)點坐標算出來，是利用\(jacobi\)符號判斷是否為二次剩餘。
#\(\displaystyle E(F_{p})=p+1+\sum_{x=0}^{p-1} \left(\frac{x^3+ax+b}{p}\right)\)，其中#\(E(F_p)\)為橢圓曲線點個數，\(\displaystyle \left( \frac{}{}\right)\)為\(jacobi\)符號
但這個方法很沒效率，更有效率的方法請看\(wiki\)介紹。
[url]https://en.wikipedia.org/wiki/Counting_points_on_elliptic_curves[/url]

另外#\(E(F_p)=p+1-t\)，\(t\)定義為\(Frobenius\) \(trace\)，elliptic_curves.mac套件中指令為ec_trace()



[color=green]要先載入elliptic_curves.mac才能呼叫ec_set_curve和ec_trace指令[/color]
[color=red](%i1)[/color]　[color=blue]load("elliptic_curves.mac");[/color]
[color=red](%o1)[/color]　[i]C:/maxima-5.43.0/share/maxima/5.43.0/share/contrib/elliptic_curves/elliptic_curves.mac[/i]

[color=green]要先載入gf.mac才能呼叫msqrt指令(計算二次剩餘的根)[/color]
[color=red](%i2)[/color]　[color=blue]load("gf.mac");[/color]
[color=red](%o2)[/color]　[i]C:/maxima-5.43.0/share/maxima/5.43.0/share/contrib/gf/gf.mac[/i]

[color=green]y^2=x^3+ax+b (mod p)
y^2=x^3+12x+1 (mod 23)的橢圓曲線[/color]
[color=red](%i4)[/color]
[color=blue][p,a,b]:[23,12,1];
ec_set_curve(p,a,b);[/color]
[color=red](%o3)[/color]　\([23,12,1]\)
[color=red](%o4)[/color]　[i]true[/i]

[color=green]以x=0,1,...,p-1計算橢圓曲線點個數[/color]
[color=red](%i5)[/color]
[color=blue]EllipticCurveOrder(p,a,b):=block
([count:0,yy,JacobiSymobl],
for x:0 thru p-1 do
  (print("當x=",x,"時,y"^2,"=",x,""^3,"+",a,"*",x,"+",b,"=",
            yy:x^3+a*x+b,"=",yy:mod(yy,p),"(mod ",p,")"),
   JacobiSymbol:jacobi(yy,p),/*判斷y^2=yy(mod p)是否有解*/
   if JacobiSymbol=0 then/*yy整除p時,有一解*/
     (count:count+1,
      print("(x,y)=(",x,",",0,")")
     )
   else if JacobiSymbol=1 then/*用msqrt()計算兩解的值*/
      (count:count+2,
        Msqrt:msqrt(yy,p),
        print("(x,y)=(",x,",",Msqrt[1],"),(",x,",",Msqrt[2],")")
       )
   else if JacobiSymbol=-1 then
      (print("無解"))
  ),
return(count+1)/*無窮遠點也算進來再+1*/
)$[/color]

[color=green]計算橢圓曲線點個數[/color]
[color=red](%i6)[/color]　[color=blue]EllipticCurveOrder(p,a,b);[/color]
當\(x=0\)時,\(y^2=0^3+12*0+1=1=1\pmod{23}\)
\((x,y)=(0,1),(0,22)\)
當\(x=1\)時,\(y^2=1^3+12*1+1=14=14\pmod{23}\)
無解
當\(x=2\)時,\(y^2=2^3+12*2+1=33=10\pmod{23}\)
無解
當\(x=3\)時,\(y^2=3^3+12*3+1=64=18\pmod{23}\)
\((x,y)=(3,15),(3,8)\)
當\(x=4\)時,\(y^2=4^3+12*4+1=113=21\pmod{23}\)
無解
當\(x=5\)時,\(y^2=5^3+12*5+1=186=2\pmod{23}\)
\((x,y)=(5,5),(5,18)\)
當\(x=6\)時,\(y^2=6^3+12*6+1=289=13\pmod{23}\)
\((x,y)=(6,17),(6,6)\)
當\(x=7\)時,\(y^2=7^3+12*7+1=428=14\pmod{23}\)
無解
當\(x=8\)時,\(y^2=8^3+12*8+1=609=11\pmod{23}\)
無解
當\(x=9\)時,\(y^2=9^3+12*9+1=838=10\pmod{23}\)
無解
當\(x=10\)時,\(y^2=10^3+12*10+1=1121=17\pmod{23}\)
無解
當\(x=11\)時,\(y^2=11^3+12*11+1=1464=15\pmod{23}\)
無解
當\(x=12\)時,\(y^2=12^3+12*12+1=1873=10\pmod{23}\)
無解
當\(x=13\)時,\(y^2=13^3+12*13+1=2354=8\pmod{23}\)
\((x,y)=(13,13),(13,10)\)
當\(x=14\)時,\(y^2=14^3+12*14+1=2913=15\pmod{23}\)
無解
當\(x=15\)時,\(y^2=15^3+12*15+1=3556=14\pmod{23}\)
無解
當\(x=16\)時,\(y^2=16^3+12*16+1=4289=11\pmod{23}\)
無解
當\(x=17\)時,\(y^2=17^3+12*17+1=5118=12\pmod{23}\)
\((x,y)=(17,14),(17,9)\)
當\(x=18\)時,\(y^2=18^3+12*18+1=6049=0\pmod{23}\)
\((x,y)=(18,0)\)
當\(x=19\)時,\(y^2=19^3+12*19+1=7088=4\pmod{23}\)
\((x,y)=(19,2),(19,21)\)
當\(x=20\)時,\(y^2=20^3+12*20+1=8241=7\pmod{23}\)
無解
當\(x=21\)時,\(y^2=21^3+12*21+1=9514=15\pmod{23}\)
無解
當\(x=22\)時,\(y^2=22^3+12*22+1=10913=11\pmod{23}\)
無解
[color=red](%o6)[/color]　\(16\)

[color=green]改用jacobi符號計算橢圓曲線點個數[/color]
[color=red](%i7)[/color]　[color=blue]JacobiList:create_list(jacobi(x^3+a*x+b,p),x,0,p-1);[/color]
[color=red](JacobiList)[/color]　\(\left[1,-1,-1,1,-1,1,1,-1,-1,-1,-1,-1,-1,1,-1,-1,-1,1,0,1,-1,-1,-1 \right]\)

[color=green]jacobi結果全部相加[/color]
[color=red](%i8)[/color]　[color=blue]apply("+",JacobiList);[/color]
[color=red](%o8)[/color]　\(-8\)

[color=green]橢圓曲線點個數[/color]
[color=red](%i9)[/color]　[color=blue]p+1+%;[/color]
[color=red](%o9)[/color]　\(16\)

[color=green]利用ec_trace()指令直接得到trace值[/color]
[color=red](%i10)[/color]　[color=blue]t:ec_trace();[/color]
[color=red](t)[/color]　\(8\)

[color=green]橢圓曲線點個數[/color]
[color=red](%i11)[/color]　[color=blue]p+1-t;[/color]
[color=red](%o11)[/color]　\(16\)

bugmens 發表於 2019-9-21 10:08

橢圓曲線加密的安全性倚賴橢圓曲線離散對數問題(\(Elliptic\) \(Curve\) \(Discrete\) \(Logarithm\) \(Problem\)，\(ECDLP\))的困難度上，橢圓曲線離散對數問題敘述如下：

令\(E\)為基於有限體\(F_p\)的橢圓曲線，\(P,Q\)為橢圓曲線上的點坐標(\(P,Q \in E(F_p)\))且\(Q=kP,k \in N\)，當\(p\)很大時，計算\(k\)值非常困難。

解決橢圓曲線離散對數問題的方法會在另一篇文章討論，這裡僅示範暴力搜尋法找出\(k\)值。
elliptic_curves.mac套件中解橢圓曲線離散對數問題的指令為k:ec_log(P,Q);



[color=green]要先載入elliptic_curves.mac才能呼叫ec_set_curve,ec_add,ec_point_order,ec_log[/color]
[color=red](%i1)[/color]　[color=blue]load("elliptic_curves.mac");[/color]
[color=red](%o1)[/color]　[i]C:/maxima-5.43.0/share/maxima/5.43.0/share/contrib/elliptic_curves/elliptic_curves.mac[/i]

[color=green]y^2=x^3+ax+b (mod p)
y^2=x^3+12x+1 (mod 23)的橢圓曲線[/color]
[color=red](%i3)[/color]
[color=blue][p,a,b]:[23,12,1];
ec_set_curve(p,a,b);[/color]
[color=red](%o2)[/color]　\([23,12,1]\)
[color=red](%o3)[/color]　[i]ture[/i]

[color=green]計算橢圓曲線點個數[/color]
[color=red](%i4)[/color]　[color=blue]order:p+1-ec_trace();[/color]
[color=red](order)[/color]　16

[color=green]生成數P=[0,1][/color]
[color=red](%i5)[/color]　[color=blue]P:[0,1];[/color]
[color=red](P)[/color]　\([0,1]\)

[color=green]計算P的階數n，每個點坐標有不同的階數,但都是橢圓曲線點個數的因數
例如：
P=[0,1],階數16,16P=無窮遠點
P=[13,13],階數8,8P=無窮遠點
P=[5,5],階數16,16P=無窮遠點
P=[3,15],階數4,4P=無窮遠點[/color]
[color=red](%i6)[/color]　[color=blue]n:ec_point_order(P,order);[/color]
[color=red](n)[/color]　16

[color=red](%i7)[/color]　[color=blue]Q:[17,14];[/color]
[color=red](Q)[/color]　\([17,14]\)

[color=green]計算2P,3P,...,kP，看何時等於Q[/color]
[color=red](%i9)[/color]
[color=blue]kP: P$
for k:2 thru n do
  (print(k," P=",kP:ec_add(P,kP)),
   if kP=Q then
      (return(k))
  );[/color]
\(2P=[13,13]\)
\(3P=[5,5]\)
\(4P=[3,15]\)
\(5P=[6,17]\)
\(6P=[19,2]\)
\(7P=[17,9]\)
\(8P=[18,0]\)
\(9P=[17,14]\)
[color=red](%o9)[/color]　9

[color=green]直接呼叫ec_log計算k[/color]
[color=red](%i10)[/color]　[color=blue]k:ec_log(P,Q);[/color]
[color=red](k)[/color]　9

bugmens 發表於 2019-10-4 09:25

關於maxima的Galois體指令可以參考以下pdf檔
C:\maxima-5.43.0\share\maxima\5.43.0\share\contrib\gf\gf_manual.pdf



[color=green]\(F[2]/(x^4+x+1)\)的Galois體[/color]
[color=red](%i1)[/color]　[color=blue]gf_set_data(2,4);[/color]
[color=red](%o1)[/color]　Structure [GF-DATA]

[color=green]\(F[2]/(x^4+x+1)\)的Galois體相關訊息[/color]
[color=red](%i2)[/color]　[color=blue]gf_info();[/color]
characteristic = 2
reduction polynomial = x^4+x+1
primitive element = x
nr of elements = 16
nr of units = 15
nr of primitive elements = 8
[color=red](%o2)[/color]　false

[color=green]\(F[2]/(x^4+x+1)\)不可約多項式為\(x^4+x+1\)，元素個數16個[/color]
[color=red](%i3)[/color]　[color=blue][p,ReductionPoly,g,n,m]:gf_infolist();[/color]
[color=red](%o3)[/color]　\([2,x^4+x+1,x,16,15]\)

[color=green]\(F[2]/(x^4+x+1)\)各元素運算過程[/color]
[color=red](%i5)[/color]
[color=blue]gpower:1$
for i:1 thru n-1 do
  (printList:["g"^i,"=(",gpower,")*(",g,")=",gpower:expand(gpower*g)],
    for j:1 thru hipow(gpower,x) do
      (if coeff(gpower,x^j)>=2 then/*係數超過2，取2同餘*/
         (printList:append(printList,["=",gpower:polymod(gpower,2)]),
          return()
         )
      ),
    /*次方超過4次方，取x^4+x+1同餘*/
    if hipow(gpower,x)>=hipow(ReductionPoly,x) then
      (printList:append(printList,["=",gpower:gf_exp(g,i)])),
    apply(print,printList)
  )$[/color]
\(g=(1)*(x)=x\)
\(g^2=(x)*(x)=x^2\)
\(g^3=(x^2)*(x)=x^3\)
\(g^4=(x^3)*(x)=x^4=x+1\)
\(g^5=(x+1)*(x)=x^2+x\)
\(g^6=(x^2+x)*(x)=x^3+x^2\)
\(g^7=(x^3+x^2)*(x)=x^4+x^3=x^3+x+1\)
\(g^8=(x^3+x+1)*(x)=x^4+x^2+x=x^2+1\)
\(g^9=(x^2+1)*(x)=x^3+x\)
\(g^{10}=(x^3+x)*(x)=x^4+x^2=x^2+x+1\)
\(g^{11}=(x^2+x+1)*(x)=x^3+x^2+x\)
\(g^{12}=(x^3+x^2+x)*(x)=x^4+x^3+x^2=x^3+x^2+x+1\)
\(g^{13}=(x^3+x^2+x+1)*(x)=x^4+x^3+x^2+x=x^3+x^2+1\)
\(g^{14}=(x^3+x^2+1)*(x)=x^4+x^3+x=x^3+1\)
\(g^{15}=(x^3+1)*(x)=x^4+x=1\)

[color=green]轉二進位副程式[/color]
[color=red](%i6)[/color]
[color=blue]bits(n) := block
([l :""],
while n # 0 do
  (l : concat(mod(n, 2),l),
   n : floor(n / 2)),
   return(l)
  )$[/color]

[color=green]16個元素的多項式、二進位、十進位表示法[/color]
[color=red](%i7)[/color]
[color=blue]addrow(matrix(["次方","多項式","二進位","十進位"]),
             matrix(["g"^"0","0","0","0"]),
       genmatrix(lambda([i,j],base10:gf_p2n(gf_exp(g,i)),
                                     if j=1 then "g"^i
                             else if j=2 then gf_exp(g,i)
                             else if j=3 then bits(base10)
                             else if j=4 then base10),n-1,4));[/color]
[color=red](%o7)[/color]
\(\left[\matrix{次方&多項式&二進位&十進位\cr
g^0&0&0&0\cr
g&x&10&2\cr
g^2&x^2&100&4\cr
g^3&x^3&1000&8\cr
g^4&x+1&11&3\cr
g^5&x^2+x&110&6\cr
g^6&x^3+x^2&1100&12\cr
g^7&x^3+x+1&1011&11\cr
g^8&x^2+1&101&5\cr
g^9&x^3+x&1010&10\cr
g^{10}&x^2+x+1&111&7\cr
g^{11}&x^3+x^2+x&1110&14\cr
g^{12}&x^3+x^2+x+1&1111&15\cr
g^{13}&x^3+x^2+1&1101&13\cr
g^{14}&x^3+1&1001&9\cr
g^{15}&1&1&1}\right]\)

[color=green]兩元素a,b[/color]
[color=red](%i9)[/color]
[color=blue]a:x^3+x;
b:x^3+x^2+1;[/color]
[color=red](a)[/color]　\(x^3+x\)
[color=red](b)[/color]　\(x^3+x^2+1\)

[color=green]相加[/color]
[color=red](%i10)[/color]　[color=blue]print("g"^gf_index(a),"+g"^gf_index(b),"=(",a,")+(",b,")=",a+b,"=",add:polymod(a+b,2),"=g"^gf_index(add))$[/color]
\(g^9+g^{13}=(x^3+x)+(x^3+x^2+1)=2x^3+x^2+x+1=x^2+x+1=g^{10}\)

[color=green]相加指令為gf_add()[/color]
[color=red](%i11)[/color]　[color=blue]gf_add(a,b);[/color]
[color=red](%o11)[/color]　\(x^2+x+1\)

[color=green]相減[/color]
[color=red](%i12)[/color]　[color=blue]print("g"^gf_index(a),"-g"^gf_index(b),"=(",a,")-(",b,")=",a-b,"=",sub:polymod(a-b,2),"=g"^gf_index(sub))$[/color]
\(g^9-g^{13}=(x^3+x)-(x^3+x^2+1)=-x^2+x-1=x^2+x+1=g^{10}\)

[color=green]相減指令為gf_sub()[/color]
[color=red](%i13)[/color]　[color=blue]gf_sub(a,b);[/color]
[color=red](%o13)[/color]　\(x^2+x+1\)

[color=green]相乘[/color]
[color=red](%i14)[/color]
[color=blue]print("g"^gf_index(a),"*g"^gf_index(b),"=(",a,")*(",b,")=",expand(a*b),"=",
        mult:polymod(remainder(a*b,ReductionPoly),2),"=g"^gf_index(mult))$[/color]
\(g^9*g^{13}=(x^3+x)*(x^3+x^2+1)=x^6+x^5+x^4+2x^3+x=x^3+x+1=g^7\)

[color=green]相乘指令為gf_mult()[/color]
[color=red](%i15)[/color]　[color=blue]gf_mult(a,b);[/color]
[color=red](%o15)[/color]　\(x^3+x+1\)

[color=green]反元素[/color]
[color=red](%i16)[/color]　[color=blue]print("1/g"^gf_index(b),"=",1/b,"=",inv_b:gf_inv(b),"=g"^gf_index(inv_b))$[/color]
\(\displaystyle 1/g^{13}=\frac{1}{x^3+x^2+1}=x^2=g^2\)

[color=green]反元素指令為gf_inv()[/color]
[color=red](%i17)[/color]　[color=blue]gf_inv(b);[/color]
[color=red](%o17)[/color]　\(x^2\)

[color=green]相除[/color]
[color=red](%i18)[/color]
[color=blue]print("g"^gf_index(a),"/g"^gf_index(b),"=",a/b,"=(",a,")*(",inv_b:gf_inv(b),")=",
        expand(a*inv_b),"=",div:polymod(remainder(a*inv_b,ReductionPoly),2),"=g"^gf_index(div))$[/color]
\(\displaystyle g^9/g^{13}=\frac{x^3+x}{x^3+x^2+1}=(x^3+x)*(x^2)=x^5+x^3=x^3+x^2+x=g^{11}\)

[color=green]相除指令為gf_div()[/color]
[color=red](%i19)[/color]　[color=blue]gf_div(a,b);[/color]
[color=red](%o19)[/color]　\(x^3+x^2+x\)

bugmens 發表於 2019-10-11 12:44

二元橢圓曲線加密請見連結介紹，本文僅就例4、例5以\(maxima\)示範計算過程。
[url]https://math.pro/db/attachment.php?aid=5248&k=23e470df1ebb8f1f5ca07d9f3c723d04&t=1565140250[/url]

※橢圓曲線在\(F_{2^n}\)下的運算規則
例子4：在有限體\(F_{2^4}\)之下，取橢圓曲線\(y^2+xy=x^3+g^8x^2+g^2\)上的兩點\(P=(g^3,g^9)\)及\(Q=(g,g^5)\)，其中\(g=(0010)\)為\(F_{2^4}\)的生成數，且不可約多項式為\(f(x)=x^4+x+1\)。若\(P+Q=R=(x_3,y_3)\)，則\(R=\)?
例子5：承例子4，若\(P+P=2P=R=(x_3,y_3)\)，則\(R=\)?



[color=green]F[2]/(x^4+x+1) 的Galois體[/color]
[color=red](%i1)[/color]　[color=blue]gf_set_data(2,4);[/color]
[color=red](%o1)[/color]　Structure [GF-DATA]

[color=green]F[2]/(x^4+x+1) 不可約多項式為x^4+x+1，元素個數16個[/color]
[color=red](%i2)[/color]　[color=blue][p,ReductionPoly,g,n,m]:gf_infolist();[/color]
[color=red](%o2)[/color]　\(\left[2,x^4+x+1,x,16,15 \right]\)

[color=green]產生g^1,g^2,...,g^15[/color]
[color=red](%i3)[/color]　[color=blue]g:create_list(gf_exp(g,i),i,1,n-1);[/color]
[color=red](g)[/color]　\(\left[x,x^2,x^3,x+1,x^2+x,x^3+x^2,x^3+x+1,x^2+1,x^3+x,x^2+x+1,x^3+x^2+x,x^3+x^2+x+1,x^3+x^2+1,x^3+1,1\right]\)

[color=green]y^2+xy=x^3+ax^2+b (mod 2)(mod x^4+x+1)的橢圓曲線[/color]
[color=red](%i5)[/color]
[color=blue]a:g[8];
b:g[2];[/color]
[color=red](a)[/color]　\(x^2+1\)
[color=red](b)[/color]　\(x^2\)

[color=green]分式的分子和分母係數超過2以(mod 2)化簡[/color]
[color=red](%i6)[/color]
[color=blue]PolyMod2(poly,PrintList):=block
([numpoly:num(poly),denompoly:denom(poly),numpolymod2,denompolymod2],
numpolymod2:polymod(numpoly,2),/*分子多項式對2同餘*/
denompolymod2:polymod(denompoly,2),/*分母多項式對2同餘*/
if denompoly=1 then/*分母為1,是多項式*/
    (if numpoly#numpolymod2 then
        (PrintList:append(PrintList,["=",poly:numpolymod2]))
    )
else
    (if numpoly#numpolymod2 or denompoly#denompolymod2 then
       (PrintList:append(PrintList,["=",poly:numpolymod2/denompolymod2]))
    ),
return([poly,PrintList])
)$[/color]

[color=green]分式的分子和分母超過4次方以(mod x^4+x+1)化簡[/color]
[color=red](%i7)[/color]
[color=blue]PolyModGF(poly,PrintList):=block
([numpoly:num(poly),denompoly:denom(poly),numpolymodGF,denompolymodGF],
numpolymodGF:polymod(remainder(numpoly,ReductionPoly),2),/*分子多項式對x^4+x+1同餘*/
denompolymodGF:polymod(remainder(denompoly,ReductionPoly),2),/*分母多項式對x^4+x+1同餘*/
if denompoly=1 then/*分母為1,是多項式*/
   (if numpoly#numpolymodGF then
        (PrintList:append(PrintList,["=",poly:numpolymodGF]))
   )
else/*分母不為1是分式,針對分子和分母分開對x^4+x+1同餘*/
  (if numpoly#numpolymodGF or denompoly#denompolymodGF then
     (PrintList:append(PrintList,["=",poly:numpolymodGF/denompolymodGF]))
  ),
return([poly,PrintList])
)$[/color]

[color=green]針對分式的分母求反元素[/color]
[color=red](%i8)[/color]
[color=blue]Inverse(poly,PrintList):=block
([numpoly:num(poly),denompoly:denom(poly)],
if denompoly#1 then/*分母不為1才需要計算反元素*/
  (PrintList:append(PrintList,["=(",numpoly,")(",gf_inv(denompoly),")"]),
   PrintList:append(PrintList,["=",poly:expand(numpoly*gf_inv(denompoly))]),
   [poly,PrintList]:PolyMod2(poly,PrintList),
   [poly,PrintList]:PolyModGF(poly,PrintList)
  ),
return([poly,PrintList])
)$[/color]

[color=green]橢圓曲線點相加
注意：
多項式係數超過2以(mod 2)化簡
超過4次方以(mod x^4+x+1)化簡[/color]
[color=red](%i9)[/color]
[color=blue]PointAdd(P,Q):=block
([R:[0,0],lambda,inv_P1,PrintList],
if P=[inf,inf] then (return(Q))
else if Q=[inf,inf] then (return(P))
else if P=Q then
  (if P[1]=0 then
     (print(" P,Q兩點x座標為0,相加結果為無窮遠點"),
       return([inf,inf])
     ),
   /*計算λ=(P1^2+P2)/P1*/
   PrintList:["λ=",(x_1^2+y_1)/x_1,"=",lambda:(P[1]^2+P[2])/P[1]],
   if P[1]^2#expand(P[1]^2) then/*若P[1]超過1項,將平方展開過程列出來*/
     (PrintList:append(PrintList,["=",lambda:expand(P[1]^2+P[2])/P[1]])),
   [lambda,PrintList]:PolyMod2(lambda,PrintList),/*若係數大於2,同餘2*/
   [lambda,PrintList]:PolyModGF(lambda,PrintList),/*若超過4次方,同餘x^4+x+1*/
   [lambda,PrintList]:Inverse(lambda,PrintList),/*若λ為分式,將分母取反元素後和分子相乘*/
   apply(print,PrintList),/*將λ計算過程印出來*/
   /*計算R1=λ^2+λ+a*/
   print(x_3,"=λ"^2,"+λ+a=(",lambda,")"^2,"+(",lambda,")+(",a,")"),
   PrintList:[x_3,"=(",expand(lambda^2),")+",lambda+a,
                         "=",R[1]:expand(lambda^2)+lambda+a],
   [R[1],PrintList]:PolyMod2(R[1],PrintList),/*若係數大於2,同餘2*/
   [R[1],PrintList]:PolyModGF(R[1],PrintList),/*若超過4次方,同餘x^4+x+1*/
   apply(print,PrintList),/*將R1計算過程印出來*/
   /*計算R2=P1^2+λR1+R1*/
   print(y_3,"=",x_1^2,"+λ",x_3,"+",x_3,"=(",P[1],")"^2,"+(",lambda,")(",R[1],")+(",R[1],")"),
   PrintList:[y_3,"=(",expand(P[1]^2),")+(",expand(lambda*R[1]),")+",R[1],
                          "=",R[2]:expand(P[1]^2+lambda*R[1])+R[1]],
   [R[2],PrintList]:PolyMod2(R[2],PrintList),/*若係數大於2,同餘2*/
   [R[2],PrintList]:PolyModGF(R[2],PrintList),/*若超過4次方,同餘x^4+x+1*/
   apply(print,PrintList)/*將R2計算過程印出來*/
  )
else
  (if gf_add(P[1],Q[1])=0 then
      (print(" P,Q兩點x座標相同,相加結果為無窮遠點"),
       return([inf,inf])
      ),
   /*計算λ=(P2+Q2)/(P1+Q1)*/
   PrintList:["λ=",(y_1+y_2)/(x_1+x_2),"=",lambda:(P[2]+Q[2])/(P[1]+Q[1])],
   [lambda,PrintList]:PolyMod2(lambda,PrintList),/*若係數大於2,同餘2*/
   [lambda,PrintList]:PolyModGF(lambda,PrintList),/*若超過4次方,同餘x^4+x+1*/
   [lambda,PrintList]:Inverse(lambda,PrintList),/*若λ為分式,將分母取反元素後和分子相乘*/
   apply(print,PrintList),/*將λ計算過程印出來*/
   /*計算R1=λ^2+λ+P1+Q1+a*/
   print(x_3,"=λ"^2,"+λ+",x_1,"+",x_2,"+a=(",lambda,")"^2,"+(",lambda,")+(",P[1],")+(",Q[1],")+(",a,")"),
   PrintList:[x_3,"=(",expand(lambda^2),")+",lambda+P[1]+Q[1]+a,
                         "=",R[1]:expand(lambda^2)+lambda+P[1]+Q[1]+a],
   [R[1],PrintList]:PolyMod2(R[1],PrintList),/*若係數大於2,同餘2*/
   [R[1],PrintList]:PolyModGF(R[1],PrintList),/*若超過4次方,同餘x^4+x+1*/
   apply(print,PrintList),/*將R1計算過程印出來*/
   /*計算R2=λ(P1+R1)+R1+P2*/
   print(y_3,"=λ(",x_1,"+",x_3,")+",x_3,"+",y_1,"=(",lambda,")(",P[1]+R[1],")+(",R[1],")+(",P[2],")"),
   PrintList:[y_3,"=(",expand(lambda*(P[1]+R[1])),")+",R[1]+P[2],
                          "=",R[2]:expand(lambda*(P[1]+R[1]))+R[1]+P[2]],
   [R[2],PrintList]:PolyMod2(R[2],PrintList),/*若係數大於2,同餘2*/
   [R[2],PrintList]:PolyModGF(R[2],PrintList),/*若超過4次方,同餘x^4+x+1*/
   apply(print,PrintList)/*將R2計算過程印出來*/
   ),
print("(",x_3,",",y_3,")=(",R[1],",",R[2],")=(",if R[1]=0 then 0 else "g"^gf_index(R[1]),",",
                                                                      if R[2]=0 then 0 else "g"^gf_index(R[2]),")"),
return(R)
)$[/color]

[color=green]P=(x1,y1)=(g^3,g^9)
Q=(x2,y2)=(g^1,g^5)[/color]
[color=red](%i11)[/color]
[color=blue]P:[g[3],g[9]];
Q:[g[1],g[5]];[/color]
[color=red](P)[/color]　\(\left[x^3,x^3+x \right]\)
[color=red](Q)[/color]　\(\left[x,x^2+x \right]\)

[color=green]計算P+Q=R=(x3,y3)[/color]
[color=red](%i12)[/color]　[color=blue]R:PointAdd(P,Q);[/color]
\(\displaystyle \lambda=\frac{y_2+y_1}{x_2+x_1}=\frac{x^3+x^2+2x}{x^3+x}=\frac{x^3+x^2}{x^3+x}=(x^3+x^2)(x^3+x^2)=x^6+2x^5+x^4=x^6+x^4=x^3+x^2+x+1\)
\(x_3=\lambda^2+\lambda+x_1+x_2+a=(x^3+x^2+x+1)^2+(x^3+x^2+x+1)+(x^3)+(x)+(x^2+1)\)
\(x_3=(x^6+2x^5+3x^4+4x^3+3x^2+2x+1)+2x^3+2x^2+2x+2=x^6+2x^5+3x^4+6x^3+5x^2+4x+3=x^6+x^4+x^2+1=x^3+x\)
\(y_3=\lambda(x_1+x_3)+x_3+y_1=(x^3+x^2+x+1)(2x^3+x)+(x^3+x)+(x^3+x)\)
\(y_3=(2x^6+2x^5+3x^4+3x^3+x^2+x)+2x^3+2x=2x^6+2x^5+3x^4+5x^3+x^2+3x=x^4+x^3+x^2+x=x^3+x^2+1\)
\((x_3,y_3)=(x^3+x,x^3+x^2+1)=(g^9,g^{13})\)
[color=red](R)[/color]　\(\left[x^3+x,x^3+x^2+1 \right]\)

[color=green]計算P+P=2P=R=(x3,y3)[/color]
[color=red](%i13)[/color]　[color=blue]R:PointAdd(P,P);[/color]
\(\displaystyle λ=\frac{y_1+x_1^2}{x_1}=\frac{x^6+x^3+x}{x^3}=\frac{x^2+x}{x^3}=(x^2+x)(x^3+x^2+x+1)=x^5+2x^4+2x^3+2x^2+x=x^5+x=x^2\)
\(x_3=λ^2+λ+a=(x^2)^2+(x^2)+(x^2+1)\)
\(x_3=(x^4)+2x^2+1=x^4+2x^2+1=x^4+1=x\)
\(y_3=x_1^2+λx_3+x_3=(x^3)^2+(x^2)(x)+(x)\)
\(y_3=(x^6)+(x^3)+x=x^6+x^3+x=x^2+x\)
\((x_3,y_3)=(x,x^2+x)=(g,g^5)\)
[color=red](R)[/color]　\(\left[x,x^2+x \right]\)

bugmens 發表於 2020-1-29 07:17

在\(F[2]/(x^4+x+1)\)Galois體之下，點\(P=(g^3,g^9)=(x^3,x^3+x)\)是橢圓曲線\(E\)：\(y^2+xy=x^3+g^8x^2+g^2\)的生成數。

\(1P=(g^3,g^9)=(x^3,x^3+x)\)
\(2P=(g^1,g^5)=(x,x^2+x)\)
\(3P=(g^9,g^{13})=(x^3+x,x^3+x^2+1)\)
\(4P=(g^8,g^2)=(x^2+1,x^2)\)
\(5P=(g^7,g^8)=(x^3+x+1,x^2+1)\)
\(6P=(1,0)=(1,0)\)
\(7P=(g^{13},g^{12})=(x^3+x^2+1,x^3+x^2+x+1)\)
\(8P=(0,g^1)=(0,x)\)
\(9P=(g^{13},g^1)=(x^3+x^2+1,x)\)
\(10P=(1,1)=(1,1)\)
\(11P=(g^7,g^{11})=(x^3+x+1,x^3+x^2+x)\)
\(12P=(g^8,1)=(x^2+1,1)\)
\(13P=(g^9,g^{10})=(x^3+x,x^2+x+1)\)
\(14P=(g^1,g^2)=(x,x^2)\)
\(15P=(g^3,g^1)=(x^3,x)\)



[color=green]\(F[2]/(x^4+x+1)\)的Galois體[/color]
[color=red](%i1)[/color]　[color=blue]gf_set_data(2,4);[/color]
[color=red](%o1)[/color]　Structure [GF-DATA]

[color=green]\(F[2]/(x^4+x+1)\)不可約多項式為\(x^4+x+1\)，元素個數16個[/color]
[color=red](%i2)[/color]　[color=blue][p,ReductionPoly,g,n,m]:gf_infolist();[/color]
[color=red](%o2)[/color]　\([2,x^4+x+1,x,16,15]\)

[color=green]產生\(g^1,g^2,\ldots,g^{15}\)[/color]
[color=red](%i3)[/color]　[color=blue]g:create_list(gf_exp(g,i),i,1,n-1);[/color]
[color=red](g)[/color]　\([x,x^2,x^3,x+1,x^2+x,x^3+x^2,x^3+x+1,x^2+1,x^3+x,x^2+x+1,x^3+x^2+x,x^3+x^2+x+1,x^3+x^2+1,x^3+1,1]\)

[color=green]\(y^2+xy=x^3+ax^2+b \pmod{2}\pmod{x^4+x+1}\)的橢圓曲線[/color]
[color=red](%i5)[/color]
[color=blue]a:g[8];
b:g[2];[/color]
[color=red](a)[/color]　\(x^2+1\)
[color=red](b)[/color]　\(x^2\)

[color=green]橢圓曲線點相加[/color]
[color=red](%i6)[/color]
[color=blue]PointAdd(P,Q):=block
([R:[0,0],lambda],
if P=[inf,inf] then (return(Q))
else if Q=[inf,inf] then (return(P))
else if P=Q then
  (if P[1]=0 then
     (print(" P,Q兩點x座標為0,相加結果為無窮遠點"),
       return([inf,inf])
     ),
  lambda:gf_add(P[1],gf_div(P[2],P[1])),
  R[1]:gf_add(gf_add(gf_add(gf_mult(lambda,lambda)),lambda),a),
  R[2]:gf_add(gf_add(gf_mult(P[1],P[1]),gf_mult(lambda,R[1])),R[1])
  )
else
  (if gf_add(P[1],Q[1])=0 then
      (print(" P,Q兩點x座標相同,相加結果為無窮遠點"),
       return([inf,inf])
      ),
   lambda:gf_div(gf_add(P[2],Q[2]),gf_add(P[1],Q[1])),
   R[1]:gf_add(gf_add(gf_add(gf_add(gf_mult(lambda,lambda),lambda),P[1]),Q[1]),a),
   R[2]:gf_add(gf_add(gf_mult(lambda,gf_add(P[1],R[1])),R[1]),P[2])
  ),
return(R)
)$[/color]

[color=green]生成數\(P=(g^3,g^9)\)[/color]
[color=red](%i7)[/color]　[color=blue]P:[g[3],g[9]];[/color]
[color=red](P)[/color]　\([x^3,x^3+x]\)

[color=green]計算\(P,2P,3P,\ldots,15P\)，結果放在\(L\)[/color]
[color=red](%i10)[/color]
[color=blue]L:[P]$
for k:2 thru 15 do
  (L:append(L,[PointAdd(last(L),P)])
  )$
L;[/color]
[color=red](%o10)[/color]　\(\matrix{[[x^3,x^3+x],[x,x^2+x],[x^3+x,x^3+x^2+1],[x^2+1,x^2],[x^3+x+1,x^2+1],[1,0],[x^3+x^2+1,x^3+x^2+x+1], \cr
[0,x],[x^3+x^2+1,x],[1,1],[x^3+x+1,x^3+x^2+x],[x^2+1,1],[x^3+x,x^2+x+1],[x,x^2],[x^3,x]]}\)

[color=green]\(P,2P,\ldots,15P\)轉換成g次方坐標[/color]
[color=red](%i11)[/color]　[color=blue]Position:create_list(append([gf_index(L[ i ][1]),gf_index(L[ i ][2])]),i,1,length(L));[/color]
[color=red](Position)[/color]　\([[3,9],[1,5],[9,13],[8,2],[7,8],[0,false],[13,12],[false,1],[13,1],[0,0],[7,11],[8,0],[9,10],[1,2],[3,1]]\)

[color=green]gf_index(0)回傳false,以-1替換[/color]
[color=red](%i12)[/color]　[color=blue]Position:subst(-1,false,Position);[/color]
[color=red](Position)[/color]　\([[3,9],[1,5],[9,13],[8,2],[7,8],[0,-1],[13,12],[-1,1],[13,1],[0,0],[7,11],[8,0],[9,10],[1,2],[3,1]]\)

[color=green]\(P,2P,\ldots,15P\)位置再加上標籤[/color]
[color=red](%i13)[/color]　[color=blue]Label:create_list(append([concat(i,"P")],(Position[ i ]+[1,0])),i,1,length(Position));[/color]
[color=red](Label)[/color]　\(\matrix{[[1P,4,9],[2P,2,5],[3P,10,13],[4P,9,2],[5P,8,8],[6P,1,-1],[7P,14,12],[8P,0,1],\cr
[9P,14,1],[10P,1,0],[11P,8,11],[12P,9,0],[13P,10,10],[14P,2,2],[15P,4,1]]}\)

[color=green]將\(P,2P,\dots,15P\)各點坐標畫出來[/color]
[color=red](%i14)[/color]
[color=blue]draw2d(user_preamble="set size square; set grid;
                               set xtics ('0' -1,'1' 0,'g^1' 1,'g^2' 2,'g^3' 3,'g^4' 4,'g^5' 5,'g^6' 6,'g^7' 7,'g^8' 8,'g^9' 9,'g^{10}' 10,'g^{11}' 11,'g^{12}' 12,'g^{13}' 13,'g^{14}' 14);
                               set ytics ('0' -1,'1' 0,'g^1' 1,'g^2' 2,'g^3' 3,'g^4' 4,'g^5' 5,'g^6' 6,'g^7' 7,'g^8' 8,'g^9' 9,'g^{10}' 10,'g^{11}' 11,'g^{12}' 12,'g^{13}' 13,'g^{14}' 14);",
              grid=true,
              point_type=filled_circle,
              points(Position),
              apply(label,Label));[/color]
[color=red](%o14)[/color]　[gr2d(points,label)]
[img]https://math.pro/db/attachment.php?aid=5331&k=7f1ed5571206a8422e410e26857e1fa7&t=1580253365&noupdate=yes¬humb=yes[/img]

bugmens 發表於 2020-2-9 22:16

[color=green]要先載入elliptic_curves.mac才能呼叫ec2_set_curve,ec2_add[/color]
[color=red](%i1)[/color]　[color=blue]load("elliptic_curves.mac");[/color]
[color=red](%o1)[/color]　C:/maxima-5.43.0/share/maxima/5.43.0/share/contrib/elliptic_curves/elliptic_curves.mac

[color=green]\(y^2+xy=x^3+ax^2+b\pmod{2}\pmod{x^4+x+1}\)的橢圓曲線
\(a=g^8=x^2+1=101_2=5_{10}\)，\(b=g^2=x^2=100_2=4_{10}\)[/color]
[color=red](%i3)[/color]
[color=blue][ReductionPoly,a,b]:[x^4+x+1,5,4];
ec2_set_curve(ReductionPoly,a,b);[/color]
[color=red](%o2)[/color]　\([x^4+x+1,5,4]\)
[color=red](%o3)[/color]　[i]true[/i]

[color=green]\(4P=[g^8,g^2]=[x^2+1,x^2]=[101_2,100_2]=[5_{10},4_{10}]\)
正確答案應為\(8P=4P+4P=[0,g^1]=[0,x]=[0_2,10_2]=[0_{10},2_{10}]\)
卻得到錯誤答案\([0,5]\)[/color]
[color=red](%i4)[/color]　[color=blue]ec2_add([5,4],[5,4]);[/color]
[color=red](%o4)[/color]　\([0,5]\)

－－－－－－－－－－－－－－－－－－－－－
\(8P=4P+4P=[0,x]=[0_2,10_2]=[0_{10},2_{10}]\)計算過程為
\(\displaystyle \lambda= \frac{y_1+x_1^2}{x_1}= \frac{(x^2+1)^2+x^2}{x^2+1}=\frac{x^4+3x^2+1}{x^2+1}=\frac{x^4+x^2+1}{x^2+1}=\frac{x^2+x}{x^2+1}=( x^2+x )( x^3+x+1 ) = x^5+x^4+x^3+2x^2+x = x^5+x^4+x^3+x = x^3+x^2+x+1\)
\(x_3 =\lambda^2 +\lambda+a=( x^3+x^2+x+1 )^2 +( x^3+x^2+x+1 )+( x^2+1 )\)
\(x_3 =( x^6+2x^5+3x^4+4x^3+3x^2+2x+1 )+ x^3+2x^2+x+2 = x^6+2x^5+3x^4+5x^3+5x^2+3x+3 = x^6+x^4+x^3+x^2+x+1 = 0 \)
\(y_3 = x_1^2 +\lambda x_3 + x_3 =( x^2+1 )^2 +( x^3+x^2+x+1 )( 0 )+( 0 ) \)
\(y_3 =( x^4+2x^2+1 )+( 0 )+ 0 = x^4+2x^2+1 = x^4+1 = x \)
\(( x_3 , y_3 )=( 0 , x )=( 0 , g )\)

－－－－－－－－－－－－－－－－－－－－－
開啟ec2_add副程式原始碼
檔案路徑C:\maxima-5.43.0\share\maxima\5.43.0\share\contrib\elliptic_curves\ec2.lisp

(defun ec2-add (pt qt)
  (cond
    ((null pt) qt)
    ((null qt) pt)
    [color=red]((equal pt qt) (ec2-double pt))[/color][color=blue]因為pt和qt相等，由ec2-double副程式計算[/color]
    ((equal pt (ec2-neg qt)) nil)
    (t (let* ((*f2-red* *ec2-red*)
              (xp (car pt))
              (yp (cadr pt))
              (xq (car qt))
              (yq (cadr qt))
              (xpq (logxor xp xq))
              (m (f2-times (logxor yp yq) (f2-inv xpq)))
              (xx (logxor (f2-times m m) m xpq *ec2-a*))
              (yy (logxor (f2-times m (logxor xp xx)) xx yp)) )
        (list xx yy) ))))
;;
(defun ec2-double (pt)
  (cond
    ((null pt) nil)
    ((= 0 (car pt)) nil)
    (t (let* ((*f2-red* *ec2-red*)
              (x (car pt))
              (y (cadr pt))
              (x2 (f2-times x x))
              (xi (f2-inv x))
              [color=red](xx (logxor x2 (f2-times *ec2-b* (f2-times xi xi))))[/color][color=blue]\(\displaystyle xx=x^2+\frac{b}{x^2}=x^2+b\cdot xi^2=0_2=0_{10}\)[/color]
              [color=red](tt (logxor x (f2-times y xi)))[/color][color=blue]\(\displaystyle tt=x+\frac{y}{x}=x+y\cdot xi=1111_2=15_{10}\)[/color]
              (yy (logxor x2 [color=red](f2-times tt xx)[/color] xx)) )[color=blue]但計算\(tt\cdot xx=15_{10} \cdot 0_{10}=7_{10}\)不為0[/color]
        (list xx yy) ))))
－－－－－－－－－－－－－－－－－－－－－
開啟f2-times副程式原始碼
檔案路徑在C:\maxima-5.43.0\share\maxima\5.43.0\src\numth.lisp

(defun f2-times (a b)
  (declare (optimize (speed 3) (safety 0)))
  [color=red](if (zerop b) (return-from f2-times 0))[/color][color=blue]多加這行，當\(b=0\)時直接回傳0[/color]
  (let* ((ilen (integer-length b))
         (a1 (ash a (1- ilen)))
         (ab a1) )
    (do ((i (- ilen 2) (1- i)) (k 0))
        ((< i 0) (f2-red ab))
        (declare (fixnum i k))
      (decf k)
      (when (logbitp i b)
        (setq a1 (ash a1 k)
              ab (logxor ab a1)
              k 0 )))))

－－－－－－－－－－－－－－－－－－－－－
重新執行範例程式
[color=green]要先載入elliptic_curves.mac才能呼叫ec2_set_curve,ec2_add[/color]
[color=red](%i1)[/color]　[color=blue]load("elliptic_curves.mac");[/color]
[color=red](%o1)[/color]　C:/maxima-5.43.0/share/maxima/5.43.0/share/contrib/elliptic_curves/elliptic_curves.mac

[color=green]讓更正的程式碼能重新載入[/color]
[color=red](%i2)[/color]　[color=blue]load("numth.lisp");[/color]
[color=red](%o2)[/color]　C:/maxima-5.43.0/share/maxima/5.43.0/src/numth.lisp

[color=green]\(y^2+xy=x^3+ax^2+b\pmod{2}\pmod{x^4+x+1}\)的橢圓曲線
\(a=g^8=x^2+1=101_2=5_{10}\)，\(b=g^2=x^2=100_2=4_{10}\)[/color]
[color=red](%i4)[/color]
[color=blue][ReductionPoly,a,b]:[x^4+x+1,5,4];
ec2_set_curve(ReductionPoly,a,b);[/color]
[color=red](%o3)[/color]　\([x^4+x+1,5,4]\)
[color=red](%o4)[/color]　[i]true[/i]

[color=green]\(4P=[g^8,g^2]=[x^2+1,x^2]=[101_2,100_2]=[5_{10},4_{10}]\)
得到正確答案\(8P=4P+4P=[0,g^1]=[0,x]=[0_2,10_2]=[0_{10},2_{10}]\)[/color]
[color=red](%i5)[/color]　[color=blue]ec2_add([5,4],[5,4]);[/color]
[color=red](%o5)[/color]　\([0,2]\)
－－－－－－－－－－－－－－－－－－－－－
結論：將[url=https://math.pro/db/attachment.php?aid=5355&k=3917bf84843323d579d5fc645e4dbd2b&t=1581257128]numth.lisp[/url]取代原本的檔案C:\maxima-5.43.0\share\maxima\5.43.0\src\numth.lisp

bugmens 發表於 2020-2-10 20:55

以\(maxima\)內建程式庫計算橢圓曲線加法(\(P+Q\))、兩倍(\(P+P=2P\))、畫出\(1P,2P,\ldots,15P\)散佈圖。
要注意的是\(maxima\)程式庫以十進位表示有限體\(F_{2^n}\)元素，以下列出有限體\(F_{2^4}\)各元素。
\(\left[\matrix{次方&多項式&二進位&十進位\cr
g^0&0&0&0\cr
g&x&10&2\cr
g^2&x^2&100&4\cr
g^3&x^3&1000&8\cr
g^4&x+1&11&3\cr
g^5&x^2+x&110&6\cr
g^6&x^3+x^2&1100&12\cr
g^7&x^3+x+1&1011&11\cr
g^8&x^2+1&101&5\cr
g^9&x^3+x&1010&10\cr
g^{10}&x^2+x+1&111&7\cr
g^{11}&x^3+x^2+x&1110&14\cr
g^{12}&x^3+x^2+x+1&1111&15\cr
g^{13}&x^3+x^2+1&1101&13\cr
g^{14}&x^3+1&1001&9\cr
g^{15}&1&1&1}\right]\)

※橢圓曲線在\(F_{2^n}\)下的運算規則
例子4：在有限體\(F_{2^4}\)之下，取橢圓曲線\(y^2+xy=x^3+g^8x^2+g^2\)上的兩點\(P=(g^3,g^9)\)及\(Q=(g,g^5)\)，其中\(g=(0010)\)為\(F_{2^4}\)的生成數，且不可約多項式為\(f(x)=x^4+x+1\)。若\(P+Q=R=(x_3,y_3)\)，則\(R=\)?
例子5：承例子4，若\(P+P=2P=R=(x_3,y_3)\)，則\(R=\)?


[color=green]要先載入elliptic_curves.mac才能呼叫ec2_set_curve,ec2_point_p,ec2_add[/color]
[color=red](%i1)[/color]　[color=blue]load("elliptic_curves.mac");[/color]
[color=red](%o1)[/color]　[i]C:/maxima-5.43.0/share/maxima/5.43.0/share/contrib/elliptic_curves/elliptic_curves.mac[/i]

[color=green]將[url=https://math.pro/db/attachment.php?aid=5355&k=3917bf84843323d579d5fc645e4dbd2b&t=1581257128]numth.lisp[/url]取代原本的檔案C:\maxima-5.43.0\share\maxima\5.43.0\src\numth.lisp
讓更正的程式碼能重新載入[/color]
[color=red](%i2)[/color]　[color=blue]load("numth.lisp");[/color]
[color=red](%o2)[/color]　[i]C:/maxima-5.43.0/share/maxima/5.43.0/src/numth.lisp[/i]

[color=green]\(y^2+xy=x^3+ax^2+b\pmod{2}\pmod{x^4+x+1}\)的橢圓曲線
\(a=g^8=x^2+1=101_2=5_{10}\)，\(b=g^2=x^2=100_2=4_{10}\)[/color]
[color=red](%i4)[/color]
[color=blue][ReductionPoly,a,b]:[x^4+x+1,5,4];
ec2_set_curve(ReductionPoly,a,b);[/color]
[color=red](%o3)[/color]　\([x^4+x+1,5,4]\)
[color=red](%o4)[/color]　[i]true[/i]

[color=green]P,Q兩點座標
\(P=[g^3,g^9]=[x^3,x^3+x]=[1000_2,1010_2]=[8_{10},10_{10}]\)
\(Q=[g^1,g^5]=[x,x^2+x]=[10_2,110_2]=[2_{10},6_{10}]\)[/color]
[color=red](%i6)[/color]
[color=blue]P:[8,10];
Q:[2,6];[/color]
[color=red](P)[/color]　\([8,10]\)
[color=red](Q)[/color]　\([2,6]\)

[color=green]檢查\(P,Q\)有沒有在橢圓曲線上[/color]
[color=red](%i8)[/color]
[color=blue]ec2_point_p(P);
ec2_point_p(Q);[/color]
[color=red](%o7)[/color]　[i]true[/i]
[color=red](%o8)[/color]　[i]true[/i]

[color=green]計算\(P+Q\)
\(P+Q=[x^3+x,x^3+x^2+1]=[1010_2,1101_2]=[10_{10},13_{10}]\)[/color]
[color=red](%i9)[/color]　[color=blue]ec2_add(P,Q);[/color]
[color=red](%o9)[/color]　\([10,13]\)

[color=green]計算\(P+P=2P\)
\(P+P=2P=[x,x^2+x]=[10_2,110_2]=[2_{10},6_{10}]\)[/color]
[color=red](%i10)[/color]　[color=blue]ec2_add(P,P);[/color]
[color=red](%o10)[/color]　\([2,6]\)

－－－－－－－－－－－－－－－－－－－－－－
在\(F[2]/(x^4+x+1)\)Galois體之下，點\(P=(g^3,g^9)=(x^3,x^3+x)\)是橢圓曲線\(E\)：\(y^2+xy=x^3+g^8x^2+g^2\)的生成數。

\(1P=(g^3,g^9)=(x^3,x^3+x)=(1000_2,1010_2)=(8_{10},10_{10})\)
\(2P=(g^1,g^5)=(x,x^2+x)=(10_2,110_2)=(2_{10},6_{10})\)
\(3P=(g^9,g^{13})=(x^3+x,x^3+x^2+1)=(1010_2,1101_2)=(10_{10},13_{10})\)
\(4P=(g^8,g^2)=(x^2+1,x^2)=(101_2,100_2)=(5_{10},4_{10})\)
\(5P=(g^7,g^8)=(x^3+x+1,x^2+1)=(1011_2,101_2)=(11_{10},5_{10})\)
\(6P=(1,0)=(1,0)=(1_2,0_2)=(1_{10},0_{10})\)
\(7P=(g^{13},g^{12})=(x^3+x^2+1,x^3+x^2+x+1)=(1101_2,1111_2)=(13_{10},15_{10})\)
\(8P=(0,g^1)=(0,x)=(0_2,10_2)=(0_{10},2_{10})\)
\(9P=(g^{13},g^1)=(x^3+x^2+1,x)=(1101_2,10_2)=(13_{10},2_{10})\)
\(10P=(1,1)=(1,1)=(1_2,1_2)=(1_{10},1_{10})\)
\(11P=(g^7,g^{11})=(x^3+x+1,x^3+x^2+x)=(1011_2,1110_2)=(11_{10},14_{10})\)
\(12P=(g^8,1)=(x^2+1,1)=(101_2,1_2)=(5_{10},1_{10})\)
\(13P=(g^9,g^{10})=(x^3+x,x^2+x+1)=(1010_2,111_2)=(10_{10},7_{10})\)
\(14P=(g^1,g^2)=(x,x^2)=(10_2,100_2)=(2_{10},4_{10})\)
\(15P=(g^3,g^1)=(x^3,x)=(1000_2,10_2)=(8_{10},2_{10})\)



[color=green]要先載入elliptic_curves.mac才能呼叫ec2_set_curve,ec2_mult[/color]
[color=red](%i1)[/color]　[color=blue]load("elliptic_curves.mac");[/color]
[color=red](%o1)[/color]　[i]"C:/maxima-5.43.0/share/maxima/5.43.0/share/contrib/elliptic_curves/elliptic_curves.mac"[/i]

[color=green]將[url=https://math.pro/db/attachment.php?aid=5355&k=3917bf84843323d579d5fc645e4dbd2b&t=1581257128]numth.lisp[/url]取代原本的檔案C:\maxima-5.43.0\share\maxima\5.43.0\src\numth.lisp
讓更正的程式碼能重新載入[/color]
[color=red](%i2)[/color]　[color=blue]load("numth.lisp");[/color]
[color=red](%o2)[/color]　[i]C:/maxima-5.43.0/share/maxima/5.43.0/src/numth.lisp[/i]

[color=green]\(F[2](x^4+x+1)\)的Galois體[/color]
[color=red](%i3)[/color]　[color=blue]gf_set_data(2,4);[/color]
[color=red](%o3)[/color]　[i]Structure [GF-DATA][/i]

[color=green]\(y^2+xy=x^3+ax^2+b\pmod{2}\pmod{x^4+x+1}\)的橢圓曲線
\(a=g^8=x^2+1=101_2=5_{10}\)，\(b=g^2=x^2=100_2=4_{10}\)[/color]
[color=red](%i5)[/color]
[color=blue][ReductionPoly,a,b]:[x^4+x+1,5,4];
ec2_set_curve(ReductionPoly,a,b);[/color]
[color=red](%o4)[/color]　\([x^4+x+1,5,4]\)
[color=red](%o5)[/color]　[i]true[/i]

[color=green]\(P=[g^3,g^9]=[x^3,x^3+x]=[1000_2,1010_2]=[8_{10},10_{10}]\)[/color]
[color=red](%i6)[/color]　[color=blue]P:[8,10];[/color]
[color=red](P)[/color]　\([8,10]\)

[color=green]用ec2_mult(k,P)計算\( P,2P,3P,\ldots,15P\)，結果放在L[/color]
[color=red](%i7)[/color]　[color=blue]L:makelist(ec2_mult(k,P),k,1,15);[/color]
[color=red](L)[/color]　\([[8,10],[2,6],[10,13],[5,4],[11,5],[1,0],[13,15],[0,2],[13,2],[1,1],[11,14],[5,1],[10,7],[2,4],[8,2]]\)

)[color=green]\( P,2P,\ldots,15P \)轉換成g次方坐標[/color]
[color=red](%i8)[/color]　[color=blue]Position:create_list(append([gf_index(gf_n2p(L[ i ][1])),gf_index(gf_n2p(L[ i ][2]))]),i,1,length(L));[/color]
[color=red](Position)[/color]　\([[3,9],[1,5],[9,13],[8,2],[7,8],[0,false],[13,12],[false,1],[13,1],[0,0],[7,11],[8,0],[9,10],[1,2],[3,1]]\)

[color=green]gf_index(0)回傳false,以-1替換[/color]
[color=red](%i9)[/color]　[color=blue]Position:subst(-1,false,Position);[/color]
[color=red](Position)[/color]　\(\)

[color=green]\(P,2P,\ldots,15P\)位置再加上標籤[/color]
[color=red](%i10)[/color]　[color=blue]Label:create_list(append([concat(i,"P")],(Position[ i ]+[1,0])),i,1,length(Position));[/color]
[color=red](Label)[/color]　\(\matrix{[[1P,4,9],[2P,2,5],[3P,10,13],[4P,9,2],[5P,8,8],[6P,1,-1],[7P,14,12],[8P,0,1],\cr
[9P,14,1],[10P,1,0],[11P,8,11],[12P,9,0],[13P,10,10],[14P,2,2],[15P,4,1]]}\)

[color=green][/color]
[color=red](%i12)[/color]
[color=blue]draw2d(user_preamble="set size square; set grid;
                               set xtics ('0' -1,'1' 0,'g^1' 1,'g^2' 2,'g^3' 3,'g^4' 4,'g^5' 5,'g^6' 6,'g^7' 7,'g^8' 8,'g^9' 9,'g^{10}' 10,'g^{11}' 11,'g^{12}' 12,'g^{13}' 13,'g^{14}' 14);
                               set ytics ('0' -1,'1' 0,'g^1' 1,'g^2' 2,'g^3' 3,'g^4' 4,'g^5' 5,'g^6' 6,'g^7' 7,'g^8' 8,'g^9' 9,'g^{10}' 10,'g^{11}' 11,'g^{12}' 12,'g^{13}' 13,'g^{14}' 14);",
              grid=true,
              point_type=filled_circle,
              points(Position),
              apply(label,Label));[/color]
[color=red](%o12)[/color]　[i][gr2d(points,label)][/i]
[img]https://math.pro/db/attachment.php?aid=5331&k=7f1ed5571206a8422e410e26857e1fa7&t=1580253365&noupdate=yes%C2%AChumb=yes[/img]

bugmens 發表於 2020-2-22 09:31

上一篇文章在\(F[2]/(x^4+x+1)\)Galois體之下，點\(P=(g^3,g^9)=(x^3,x^3+x)\)是橢圓曲線\(E\)：\(y^2+xy=x^3+g^8x^2+g^2\)的生成數計算\(1P,2P,\ldots,15P\)，再加上無窮遠點總共16個點。

想要知道有限體上橢圓曲線的點個數，直觀的方法是將\(x=0,1,g^1,g^2,\dots,g^{14}\)代入橢圓曲線方程式\(y^2+xy=x^3+ax^2+b\)，解出\(y\)的解，將有\(y\)的解個數加起來就是答案。

將方程式\(y^2+xy=x^3+ax^2+b\)同除\(x^2\)得到\(\displaystyle \left(\frac{y}{x}\right)^2+\left(\frac{y}{x}\right)=x+a+\frac{b}{x^2}\)
令\(\displaystyle z=\frac{y}{x}\)，\(\displaystyle \beta=x+a+\frac{b}{x^2}\)得到\(z^2+z=\beta\)，\(z^2+z-\beta=0\)，\(z^2+z+\beta=0\)
利用[url=https://books.google.com.tw/books?id=0_vegzgyqGMC&lpg=PP1&hl=zh-TW&pg=PA26#v=onepage&q&f=false]這裡的方法[/url](\(maxima\)指令為f2_solve_quadratic(\(\beta\),n))解出\(z\)，再計算\(y=z\cdot x\)解出\(y的解\)。



[color=green]要先載入elliptic_curves.mac才能呼叫ec2_set_curve,f2_solve_quadratic[/color]
[color=red](%i1)[/color]　[color=blue]load("elliptic_curves.mac");[/color]
[color=red](%o1)[/color]　[i]C:/maxima-5.43.0/share/maxima/5.43.0/share/contrib/elliptic_curves/elliptic_curves.mac[/i]

[color=green]將[url=https://math.pro/db/attachment.php?aid=5355&k=3917bf84843323d579d5fc645e4dbd2b&t=1581257128]numth.lisp[/url]取代原本的檔案C:\maxima-5.43.0\share\maxima\5.43.0\src\numth.lisp
讓更正的程式碼能重新載入[/color]
[color=red](%i2)[/color]　[color=blue]load("numth.lisp");[/color]
[color=red](%o2)[/color]　[i]C:/maxima-5.43.0/share/maxima/5.43.0/src/numth.lisp[/i]

[color=green]\(F[2](x^4+x+1)\)的Galois體[/color]
[color=red](%i3)[/color]　[color=blue]gf_set_data(2,4);[/color]
[color=red](%o3)[/color]　[i]Structure [GF-DATA][/i]

[color=green]\(F[2]/(x^4+x+1)\)不可約多項式為\(x^4+x+1\)，元素個數16個[/color]
[color=red](%i4)[/color]　[color=blue][p,ReductionPoly,g,n,m]:gf_infolist();[/color]
[color=red](%o4)[/color]　\([2,x^4+x+1,x,16,15]\)

[color=green]產生\(g^1,g^2,\ldots,g^{15}\)[/color]
[color=red](%i5)[/color]　[color=blue]g:create_list(gf_exp(g,i),i,1,n-1);[/color]
[color=red](g)[/color]　\([x,x^2,x^3,x+1,x^2+x,x^3+x^2,x^3+x+1,x^2+1,x^3+x,x^2+x+1,x^3+x^2+x,x^3+x^2+x+1,x^3+x^2+1,x^3+1,1]\)

[color=green]\(y^2+xy=x^3+ax^2+b\pmod{2}\pmod{x^4+x+1}\)的橢圓曲線[/color]
[color=red](%i8)[/color]
[color=blue]a:g[8];
b:g[2];
ec2_set_curve(ReductionPoly,gf_p2n(a),gf_p2n(b));[/color]
[color=red](a)[/color]　\(x^2+1\)
[color=red](b)[/color]　\(x^2\)
[color=red](%o8)[/color]　[i]true[/i]

[color=green]分式的分子和分母係數超過2以\(\pmod{2}\)化簡[/color]
[color=red](%i9)[/color]
[color=blue]PolyMod2(poly,PrintList):=block
([numpoly:num(poly),denompoly:denom(poly),numpolymod2,denompolymod2],
numpolymod2:polymod(numpoly,2),/*分子多項式對2同餘*/
denompolymod2:polymod(denompoly,2),/*分母多項式對2同餘*/
if denompoly=1 then/*分母為1,是多項式*/
    (if numpoly#numpolymod2 then
        (PrintList:append(PrintList,["=",poly:numpolymod2]))
    )
else
    (if numpoly#numpolymod2 or denompoly#denompolymod2 then
       (PrintList:append(PrintList,["=",poly:numpolymod2/denompolymod2]))
    ),
return([poly,PrintList])
)$[/color]

[color=green]分式的分子和分母超過4次方以\(\pmod{x^4+x+1}\)化簡[/color]
[color=red](%i10)[/color]
[color=blue]PolyModGF(poly,PrintList):=block
([numpoly:num(poly),denompoly:denom(poly),numpolymodGF,denompolymodGF],
numpolymodGF:polymod(remainder(numpoly,ReductionPoly),2),/*分子多項式對x^4+x+1同餘*/
denompolymodGF:polymod(remainder(denompoly,ReductionPoly),2),/*分母多項式對x^4+x+1同餘*/
if denompoly=1 then/*分母為1,是多項式*/
   (if numpoly#numpolymodGF then
        (PrintList:append(PrintList,["=",poly:numpolymodGF]))
   )
else/*分母不為1是分式,針對分子和分母分開對x^4+x+1同餘*/
  (if numpoly#numpolymodGF or denompoly#denompolymodGF then
     (PrintList:append(PrintList,["=",poly:numpolymodGF/denompolymodGF]))
  ),
return([poly,PrintList])
)$[/color]

[color=green]針對分式的分母求反元素[/color]
[color=red](%i11)[/color]
[color=blue]Inverse(poly,PrintList):=block
([numpoly:num(poly),denompoly:denom(poly)],
if denompoly#1 then/*分母不為1才需要計算反元素*/
  (PrintList:append(PrintList,["=(",numpoly,")(",gf_inv(denompoly),")"]),
   PrintList:append(PrintList,["=",poly:expand(numpoly*gf_inv(denompoly))]),
   [poly,PrintList]: PolyMod2(poly,PrintList),
   [poly,PrintList]: PolyModGF(poly,PrintList)
  ),
return([poly,PrintList])
)$[/color]

[color=green]用大寫\(X,Y\)代表橢圓曲線，\(Y^2+XY=X^3+aX^2+b\)
用大寫\(\displaystyle Z=\frac{Y}{X}\)的變數變換
用小寫\(x\)代表Galois體運算符號[/color]
[color=red](%i12)[/color]
[color=blue]CountingPoints(ReductionPoly,a,b):=block
([count:1,Points:[[inf,inf]]/*將無窮遠點先放入*/,Degree,beta,PrintList,X,Y:[0,0],Z:[0,0]],
print("當X=0時,Y"^2,"+XY=X"^3,"+aX"^2,"+b=>Y"^2,"=b"),
Degree:hipow(ReductionPoly,x),
print("解方程式Y"^2,"=b,得Y=b^2^",(Degree-1),"=(",b,")^2^",(Degree-1),"=(",b,")"^2^(Degree-1),"=",Y[1]:gf_exp(b,2^(Degree-1))),
print("(X,Y)=(0,",Y[1],")=(0,","g"^gf_index(Y[1]),")"),
Points:append(Points,[[0,gf_exp(b,2^(Degree-1))]]),/*當X=0時Y=b^2^(Degree-1)*/
count:count+1,/*當X=0時,Y只有一解(重根)*/
for i:1 thru m do/*m為gf_infolist()除了0之外的Galois體的元素個數*/
   (print("－－－－－－－－－－－"),
    print("當X=",X:g[ i ],"時"),
    PrintList:["β=X+a+","b"/"X"^2,"=(",X,")+(",a,")+",b/X^2,"=",beta:X+a+b/X^2],
    if b/X^2#ratsimp(b/X^2) then
      (PrintList:append(PrintList,["=",beta:X+a+expand(b/X^2)])),
    if denom(beta:ratsimp(beta))#1 then/*把帶分式換成假分式*/
      (PrintList:append(PrintList,["=",beta])),
    [beta,PrintList]: PolyMod2(beta,PrintList),/*若係數大於2,同餘2*/
    [beta,PrintList]: PolyModGF(beta,PrintList),/*若超過4次方,同餘x^4+x+1*/
    [beta,PrintList]:Inverse(beta,PrintList),/*若λ為分式,將分母取反元素後和分子相乘*/
    apply(print,PrintList),
    print("解方程式Z"^2,"+Z+β=Z"^2,"+Z+(",beta,")=0"),
    Z:f2_solve_quadratic(gf_p2n(beta),Degree),
    if Z=false then
      (print("Z無解")
      )
    else
      (print("Z=",Z:[gf_n2p(Z[1]),gf_n2p(Z[2])]),
       PrintList:["Y1=Z1*X=(",Z[1],")(",X,")=",Y[1]:expand(Z[1]*X)],
       [Y[1],PrintList]: PolyMod2(Y[1],PrintList),/*若係數大於2,同餘2*/
       [Y[1],PrintList]: PolyModGF(Y[1],PrintList),/*若超過4次方,同餘x^4+x+1*/
       [Y[1],PrintList]:Inverse(Y[1],PrintList),/*若λ為分式,將分母取反元素後和分子相乘*/
       apply(print,PrintList),
       PrintList:["Y2=Z2*X=(",Z[2],")(",X,")=",Y[2]:expand(Z[2]*X)],
       [Y[2],PrintList]: PolyMod2(Y[2],PrintList),/*若係數大於2,同餘2*/
       [Y[2],PrintList]: PolyModGF(Y[2],PrintList),/*若超過4次方,同餘x^4+x+1*/
       [Y[2],PrintList]:Inverse(Y[2],PrintList),/*若λ為分式,將分母取反元素後和分子相乘*/
       apply(print,PrintList),
       print("(X,Y)=(",X,",",Y[1],")=(","g"^gf_index(X),",",if Y[1]=0 then 0 else "g"^gf_index(Y[1]),")和",
                           "(",X,",",Y[2],")=(","g"^gf_index(X),",",if Y[2]=0 then 0 else "g"^gf_index(Y[2]),")"),
       Points:append(Points,[[X,Y[1]],[X,Y[2]]]),
       count:count+2
      )
   ),
return([count,Points])
)$[/color]

[color=green]以\(x=0,1,g^1,\ldots,g^{14}\)計算橢圓曲線點個數[/color]
[color=red](%i13)[/color]　[color=blue]CountingPoints(ReductionPoly,a,b);[/color]
當\(X=0\)時,\(Y^2+XY=X^3+aX^2+b\)=>\(Y^2=b\)
解方程式\(Y^2=b\),得\(Y={b^2}^{3}={(x^2)^2}^{3}=(x^2)^{8}=x\)
\((X,Y)=(0,x)=(0,g)\)
－－－－－－－－－－－
當\(X= x \)時
\(\displaystyle \beta=X+a+\frac{b}{X^2}=(x)+(x^2+1)+1=x^2+x+2=x^2+x\)
解方程式\(Z^2+Z+\beta=Z^2+Z+(x^2+x)=0\)
\(Z=[x+1,x]\)
\(Y1=Z1*X=(x+1)(x)=x^2+x\)
\(Y2=Z2*X=(x)(x)=x^2\)
\((X,Y)=(x,x^2+x)=(g,g^5)\)和\((x,x^2)=(g,g^2)\)
－－－－－－－－－－－
當\(X=x^2\)時
\(\displaystyle \beta=X+a+\frac{b}{X^2}=(x^2)+(x^2+1)+\frac{1}{x^2}=2x^2+\frac{1}{x^2}+1=\frac{2x^4+x^2+1}{x^2}=\frac{x^2+1}{x^2}=(x^2+1)(x^3+x^2+1)\)
\(=x^5+x^4+x^3+2x^2+1=x^5+x^4+x^3+1=x^3+x^2\)
解方程式\(Z^2+Z+\beta=Z^2+Z+(x^3+x^2)=0\)
\(Z\)無解
－－－－－－－－－－－
當\(X= x^3\)時
\(\displaystyle \beta=X+a+\frac{b}{X^2}=(x^3)+(x^2+1)+\frac{1}{x^4}=x^3+x^2+\frac{1}{x^4}+1=\frac{x^7+x^6+x^4+1}{x^4}=\frac{x^2+1}{x+1}=(x^2+1)(x^3+x^2+x)\)
\(=x^5+x^4+2x^3+x^2+x=x^5+x^4+x^2+x=x+1\)
解方程式\(Z^2+Z+\beta=Z^2+Z+( x+1 )=0\)
\(Z=[x^3+x^2+1,x^3+x^2]\)
\(Y1=Z1*X=(x^3+x^2+1)(x^3)=x^6+x^5+x^3=x\)
\(Y2=Z2*X=(x^3+x^2)(x^3)=x^6+x^5=x^3+x\)
\((X,Y)=(x^3,x)=(g^3,g)\)和\((x^3,x^3+x)=(g^3,g^9)\)
－－－－－－－－－－－
當\(X=x+1\)時
\(\displaystyle \beta=X+a+\frac{b}{X^2}=(x+1)+(x^2+1)+\frac{x^2}{(x+1)^2}=\frac{x^2}{(x+1)^2}+x^2+x+2=\frac{x^2}{x^2+2x+1}+x^2+x+2=\frac{x^4+3x^3+6x^2+5x+2}{x^2+2x+1}\)
\(\displaystyle =\frac{x^4+x^3+x}{x^2+1}=\frac{x^3+1}{x^2+1}=(x^3+1)(x^3+x+1)=x^6+x^4+2x^3+x+1=x^6+x^4+x+1=x^3+x^2\)
解方程式\(Z^2+Z+\beta=Z^2+Z+(x^3+x^2)=0\)
\(Z\)無解
－－－－－－－－－－－
當\(X=x^2+x\)時
\(\displaystyle \beta=X+a+\frac{b}{X^2}=(x^2+x)+(x^2+1)+\frac{x^2}{(x^2+x)^2}=\frac{x^2}{(x^2+x)^2}+2x^2+x+1=\frac{x^2}{x^4+2x^3+x^2}+2x^2+x+1\)
\(\displaystyle =\frac{2x^4+5x^3+5x^2+3x+2}{x^2+2x+1}=\frac{x^3+x^2+x}{x^2+1}=(x^3+x^2+x)(x^3+x+1)=x^6+x^5+2x^4+2x^3+2x^2+x=x^6+x^5+x=x^3\)
解方程式\(Z^2+Z+\beta=Z^2+Z+(x^3)=0\)
\(Z\)無解
－－－－－－－－－－－
當\(X=x^3+x^2\)時
\(\displaystyle \beta=X+a+\frac{b}{X^2}=(x^3+x^2)+(x^2+1)+\frac{x^2}{(x^3+x^2)^2}=\frac{x^2}{(x^3+x^2)^2}+x^3+2x^2+1=\frac{x^2}{x^6+2x^5+x^4}+x^3+2x^2+1\)
\(\displaystyle =\frac{x^7+4x^6+5x^5+3x^4+2x^3+x^2+1}{x^4+2x^3+x^2}=\frac{x^7+x^5+x^4+x^2+1}{x^4+x^2}=\frac{x^3+x+1}{x^2+x+1}=(x^3+x+1)(x^2+x)\)
\(=x^5+x^4+x^3+2x^2+x=x^5+x^4+x^3+x=x^3+x^2+x+1\)
解方程式\(Z^2+Z+\beta=Z^2+Z+(x^3+x^2+x+1)=0\)
\(Z\)無解
－－－－－－－－－－－
當\(X=x^3+x+1\)時
\(\displaystyle \beta=X+a+\frac{b}{X^2}=(x^3+x+1)+(x^2+1)+\frac{x^2}{(x^3+x+1)^2}=\frac{x^2}{(x^3+x+1)^2}+x^3+x^2+x+2\)
\(\displaystyle =\frac{x^2}{x^6+2x^4+2x^3+x^2+2x+1}+x^3+x^2+x+2=\frac{x^9+x^8+3x^7+6x^6+5x^5+9x^4+8x^3+6x^2+5x+2}{x^6+2x^4+2x^3+x^2+2x+1}\)
\(\displaystyle =\frac{x^9+x^8+x^7+x^5+x^4+x}{x^6+x^2+1}=\frac{x+1}{x^3+1}=(x+1)(x)=x^2+x\)
解方程式\(Z^2+Z+\beta=Z^2+Z+(x^2+x)=0\)
\(Z=[x+1,x]\)
\(Y1=Z1*X=(x+1)(x^3+x+1)=x^4+x^3+x^2+2x+1=x^4+x^3+x^2+1=x^3+x^2+x\)
\(Y2=Z2*X=(x)(x^3+x+1)=x^4+x^2+x=x^2+1\)
\((X,Y)=(x^3+x+1,x^3+x^2+x)=(g^7,g^{11})\)和\((x^3+x+1,x^2+1)=(g^7,g^8)\)
－－－－－－－－－－－
當\(X=x^2+1\)時
\(\displaystyle \beta=X+a+\frac{b}{X^2}=(x^2+1)+(x^2+1)+\frac{x^2}{(x^2+1)^2}=\frac{x^2}{(x^2+1)^2}+2x^2+2=\frac{x^2}{x^4+2x^2+1}+2x^2+2=\frac{2x^6+6x^4+7x^2+2}{x^4+2x^2+1}=\frac{x^2}{x^4+1}=x\)
解方程式\(Z^2+Z+\beta=Z^2+Z+(x)=0\)
\(Z=[x^3+x,x^3+x+1]\)
\(Y1=Z1*X=(x^3+x)(x^2+1)=x^5+2x^3+x=x^5+x=x^2\)
\(Y2=Z2*X=(x^3+x+1)(x^2+1)=x^5+2x^3+x^2+x+1=x^5+x^2+x+1=1\)
\((X,Y)=(x^2+1,x^2)=(g^8,g^2)\)和\((x^2+1,1)=(g^8,1)\)
－－－－－－－－－－－
當\(X=x^3+x\)時
\(\displaystyle \beta=X+a+\frac{b}{X^2}=(x^3+x)+(x^2+1)+\frac{x^2}{(x^3+x)^2}=\frac{x^2}{(x^3+x)^2}+x^3+x^2+x+1=\frac{x^2}{x^6+2x^4+x^2}+x^3+x^2+x+1\)
\(\displaystyle =\frac{x^7+x^6+3x^5+3x^4+3x^3+3x^2+x+2}{x^4+2x^2+1}=\frac{x^7+x^6+x^5+x^4+x^3+x^2+x}{x^4+1}=\frac{x^3+x^2}{x}=(x^3+x^2)(x^3+1)=x^6+x^5+x^3+x^2=x^2+x\)
解方程式\(Z^2+Z+\beta=Z^2+Z+(x^2+x)=0\)
\(Z=[x+1,x]\)
\(Y1=Z1*X=(x+1)(x^3+x)=x^4+x^3+x^2+x=x^3+x^2+1\)
\(Y2=Z2*X=(x)(x^3+x)=x^4+x^2=x^2+x+1\)
\((X,Y)=(x^3+x,x^3+x^2+1)=(g^9,g^{13})\)和\((x^3+x,x^2+x+1)=(g^9,g^{10})\)
－－－－－－－－－－－
當\(X=x^2+x+1\)時
\(\displaystyle \beta=X+a+\frac{b}{X^2}=(x^2+x+1)+(x^2+1)+\frac{x^2}{(x^2+x+1)^2}=\frac{x^2}{(x^2+x+1)^2}+2x^2+x+2=\frac{x^2}{x^4+2x^3+3x^2+2x+1}+2x^2+x+2\)
\(\displaystyle =\frac{2x^6+5x^5+10x^4+11x^3+11x^2+5x+2}{x^4+2x^3+3x^2+2x+1}=\frac{x^5+x^3+x^2+x}{x^4+x^2+1}=\frac{x^3}{x^2+x}=(x^3)(x^2+x+1)=x^5+x^4+x^3=x^3+x^2+1\)
解方程式\(Z^2+Z+\beta=Z^2+Z+(x^3+x^2+1)=0\)
\(Z\)無解
－－－－－－－－－－－
當\(X=x^3+x^2+x\)時
\(\displaystyle \beta=X+a+\frac{b}{X^2}=(x^3+x^2+x)+(x^2+1)+\frac{x^2}{(x^3+x^2+x)^2}=\frac{x^2}{(x^3+x^2+x)^2}+x^3+2x^2+x+1\)
\(\displaystyle =\frac{x^2}{x^6+2x^5+3x^4+2x^3+x^2}+x^3+2x^2+x+1=\frac{x^7+4x^6+8x^5+11x^4+10x^3+7x^2+3x+2}{x^4+2x^3+3x^2+2x+1}\)
\(\displaystyle =\frac{x^7+x^4+x^2+x}{x^4+x^2+1}=\frac{x^3+x^2+x}{x^2+x}=(x^3+x^2+x)(x^2+x+1)= x^5+2x^4+3x^3+2x^2+x=x^5+x^3+x=x^3+x^2\)
解方程式\(Z^2+Z+\beta=Z^2+Z+(x^3+x^2)=0\)
\(Z\)無解
－－－－－－－－－－－
當\(X=x^3+x^2+x+1\)時
\(\displaystyle \beta=X+a+\frac{b}{X^2}=(x^3+x^2+x+1)+(x^2+1)+\frac{x^2}{(x^3+x^2+x+1)^2}=\frac{x^2}{(x^3+x^2+x+1)^2}+x^3+2x^2+x+2\)
\(\displaystyle =\frac{x^2}{x^6+2x^5+3x^4+4x^3+3x^2+2x+1}+x^3+2x^2+x+2=\frac{x^9+4x^8+8x^7+14x^6+18x^5+18x^4+16x^3+11x^2+5x+2}{x^6+2x^5+3x^4+4x^3+3x^2+2x+1}\)
\(\displaystyle =\frac{x^9+x^2+x}{x^6+x^4+x^2+1}=\frac{x^3+x^2}{x^3+x}=(x^3+x^2)(x^3+x^2)=x^6+2x^5+x^4=x^6+x^4=x^3+x^2+x+1\)
解方程式\(Z^2+Z+\beta=Z^2+Z+(x^3+x^2+x+1)=0\)
\(Z\)無解
－－－－－－－－－－－
當\(X=x^3+x^2+1\)時
\(\displaystyle \beta=X+a+\frac{b}{X^2}=(x^3+x^2+1)+(x^2+1)+\frac{x^2}{(x^3+x^2+1)^2}=\frac{x^2}{(x^3+x^2+1)^2}+x^3+2x^2+2\)
\(\displaystyle =\frac{x^2}{x^6+2x^5+x^4+2x^3+2x^2+1}+x^3+2x^2+2=\frac{x^9+4x^8+5x^7+6x^6+10x^5+6x^4+5x^3+7x^2+2}{x^6+2x^5+x^4+2x^3+2x^2+1}\)
\(\displaystyle =\frac{x^9+x^7+x^3+x^2}{x^6+x^4+1}=\frac{x^3+x^2+1}{x^3+x^2+x}=(x^3+x^2+1)(x+1)=x^4+2x^3+x^2+x+1=x^4+x^2+x+1=x^2\)
解方程式\(Z^2+Z+\beta=Z^2+Z+(x^2)=0\)
\(Z=[x^3+1,x^3]\)
\(Y1=Z1*X=(x^3+1)(x^3+x^2+1)=x^6+x^5+2x^3+x^2+1=x^6+x^5+x^2+1=x^3+x^2+x+1\)
\(Y2=Z2*X=(x^3)(x^3+x^2+1)=x^6+x^5+x^3=x\)
\((X,Y)=(x^3+x^2+1,x^3+x^2+x+1)=(g^{13},g^{12})\)和\((x^3+x^2+1,x)=(g^{13},g)\)
－－－－－－－－－－－
當\(X=x^3+1\)時
\(\displaystyle \beta=X+a+\frac{b}{X^2}=(x^3+1)+(x^2+1)+\frac{x^2}{(x^3+1)^2}=\frac{x^2}{(x^3+1)^2}+x^3+x^2+2=\frac{x^2}{x^6+2x^3+1}+x^3+x^2+2\)
\(\displaystyle =\frac{x^9+x^8+4x^6+2x^5+5x^3+2x^2+2}{x^6+2x^3+1}=\frac{x^9+x^8+x^3}{x^6+1}=\frac{x^2+x+1}{x^3+x^2+1}=(x^2+x+1)(x^2)=x^4+x^3+x^2=x^3+x^2+x+1\)
解方程式\(Z^2+Z+\beta=Z^2+Z+(x^3+x^2+x+1)=0\)
\(Z\)無解
－－－－－－－－－－－
當\(X=1\)時
\(\displaystyle \beta=X+a+\frac{b}{X^2}=(1)+(x^2+1)+x^2=2x^2+2=0\)
解方程式\(Z^2+Z+\beta=Z^2+Z+(0)=0\)
\(Z=[0,1]\)
\(Y1=Z1*X=(0)(1)=0\)
\(Y2=Z2*X=(1)(1)=1\)
\((X,Y)=(1,0)=(1,0)\)和\((1,1)=(1,1)\)
[color=red](%o13)[/color]　\(\matrix{[16,[[\infty,\infty],[0,x],[x,x^2],[x,x^2+x],[x^3,x^3+x],[x^3,x],[x^3+x+1,x^2+1],[x^3+x+1,x^3+x^2+x],[x^2+1,x^2],\cr
[x^2+1,1],[x^3+x,x^2+x+1],[x^3+x,x^3+x^2+1],[x^3+x^2+1,x],[x^3+x^2+1,x^3+x^2+x+1],[1,0],[1,1]]]}\)

bugmens 發表於 2020-3-3 22:17

將上一篇文章濃縮成一個副程式



[color=green]要先載入elliptic_curves.mac才能呼叫ec2_set_curve,f2_solve_quadratic[/color]
[color=red](%i1)[/color]　[color=blue]load("elliptic_curves.mac");[/color]
[color=red](%o1)[/color]　[i]C:/maxima-5.43.0/share/maxima/5.43.0/share/contrib/elliptic_curves/elliptic_curves.mac[/i]

[color=green]將[url=https://math.pro/db/attachment.php?aid=5355&k=3917bf84843323d579d5fc645e4dbd2b&t=1581257128]numth.lisp[/url]取代原本的檔案C:\maxima-5.43.0\share\maxima\5.43.0\src\numth.lisp
讓更正的程式碼能重新載入[/color]
[color=red](%i2)[/color]　[color=blue]load("numth.lisp");[/color]
[color=red](%o2)[/color]　[i]C:/maxima-5.43.0/share/maxima/5.43.0/src/numth.lisp[/i]

[color=green]計算橢圓曲線點個數和坐標副程式[/color]
[color=red](%i3)[/color]
[color=blue]CountingPoints(ReductionPoly,a,b):=block
([count:1,Points:[[inf,inf]]/*將無窮遠點先放入*/,p,g,n,m,Degree,beta,X,Z:[0,0]],
gf_set_data(2,ReductionPoly),/*F[2]/(x^4+x+1)的Galois體*/
[p,ReductionPoly,g,n,m]:gf_infolist(),
ec2_set_curve(ReductionPoly,gf_p2n(a),gf_p2n(b)),/*設定橢圓曲線參數*/
Degree:hipow(ReductionPoly,x),/*x^4+x+1的最高次方Degree=4*/
Points:append(Points,[[0,gf_exp(b,2^(Degree-1))]]),/*當X=0時Y=b^2^(Degree-1)*/
count:count+1,/*當X=0時,Y只有一解(重根)*/
for i:1 thru m do/*m為gf_infolist()除了0之外的Galois體的元素個數*/
   (X:gf_exp(g,i),/*X坐標*/
    beta:gf_add(gf_add(X,a),gf_div(b,gf_mult(X,X))),/*計算β=X+a+b/X^2*/
    Z:f2_solve_quadratic(gf_p2n(beta),Degree),/*解方程式Z^2+Z+β=0*/
    if Z#false then/*若Z有解*/
      (Z:[gf_n2p(Z[1]),gf_n2p(Z[2])],
       Points:append(Points,[[X,gf_mult(X,Z[1])],[X,gf_mult(X,Z[2])]]),/*計算Y=X*Z得到Y坐標*/
       count:count+2/*得到2個點坐標*/
      )
   ),
return([count,Points])
)$[/color]

[color=green]計算\(Y^2+XY=X^3+aX^2+b\pmod{2}\pmod{x^4+x+1}\)橢圓曲線的點個數和坐標
\(a=x^2+1,b=x^2\)[/color]
[color=red](%i4)[/color]　[color=blue][count,Points]:CountingPoints(x^4+x+1,x^2+1,x^2);[/color]
[color=red](%o4)[/color]　\(\matrix{[16,[[\infty,\infty],[0,x],[x,x^2],[x,x^2+x],[x^3,x^3+x],[x^3,x],[x^3+x+1,x^2+1],[x^3+x+1,x^3+x^2+x],[x^2+1,x^2],\cr
[x^2+1,1],[x^3+x,x^2+x+1],[x^3+x,x^3+x^2+1],[x^3+x^2+1,x],[x^3+x^2+1,x^3+x^2+x+1],[1,0],[1,1]]]}\)

[color=green]計算以不同的不可分解多項式p(x)的\(Y^2+XY=X^3+aX^2+b\pmod{2}\pmod{p(x)}\)橢圓曲線點個數
\(a=0,b=1\)[/color]
[color=red](%i7)[/color]
[color=blue]a:0$
b:1$
for k:1 thru 10 do
  (gf_set_data(2,k),
   [p,ReductionPoly,g,n,m]:gf_infolist(),
   print("Y"^2,"+XY=X"^3,"+",a,"X"^2,"+",b,"(mod 2) (mod",ReductionPoly,")",
           ",橢圓曲線點個數#E(F",2,""^k,")=",CountingPoints(ReductionPoly,a,b)[1])
  )$[/color]
\(Y^2+XY=X^3+0X^2+1 \pmod{2} \pmod{x}\),橢圓曲線點個數#\(E(F_{2})=4\)
\(Y^2+XY=X^3+0X^2+1 \pmod{2} \pmod{x^2+x+1}\),橢圓曲線點個數#\(E(F_{2^2})=8\)
\(Y^2+XY=X^3+0X^2+1 \pmod{2} \pmod{x^3+x+1}\),橢圓曲線點個數#\(E(F_{2^3})=4\)
\(Y^2+XY=X^3+0X^2+1 \pmod{2} \pmod{x^4+x+1}\),橢圓曲線點個數#\(E(F_{2^4})=16\)
\(Y^2+XY=X^3+0X^2+1 \pmod{2} \pmod{x^5+x^2+1}\),橢圓曲線點個數#\(E(F_{2^5})=44\)
\(Y^2+XY=X^3+0X^2+1 \pmod{2} \pmod{x^6+x+1}\),橢圓曲線點個數#\(E(F_{2^6})=56\)
\(Y^2+XY=X^3+0X^2+1 \pmod{2} \pmod{x^7+x+1}\),橢圓曲線點個數#\(E(F_{2^7})=116\)
\(Y^2+XY=X^3+0X^2+1 \pmod{2} \pmod{x^8+x^4+x^3+x+1}\),橢圓曲線點個數#\(E(F_{2^8})=288\)
\(Y^2+XY=X^3+0X^2+1 \pmod{2} \pmod{x^9+x+1}\),橢圓曲線點個數#\(E(F_{2^9})=508\)
\(Y^2+XY=X^3+0X^2+1 \pmod{2} \pmod{x^{10}+x^3+1}\),橢圓曲線點個數#\(E(F_{2^{10}})=968\)

[color=green]計算以不同的不可分解多項式p(x)的\(Y^2+XY=X^3+aX^2+b\pmod{2}\pmod{p(x)}\)橢圓曲線點個數
\(a=1,b=1\)[/color]
[color=red](%i10)[/color]
[color=blue]a:1$
b:1$
for k:1 thru 10 do
  (gf_set_data(2,k),
   [p,ReductionPoly,g,n,m]:gf_infolist(),
   print("Y"^2,"+XY=X"^3,"+",a,"X"^2,"+",b,"(mod 2) (mod",ReductionPoly,")",
           ",橢圓曲線點個數#E(F",2,""^k,")=",CountingPoints(ReductionPoly,a,b)[1])
  )$[/color]
\(Y^2+XY=X^3 + 1 X^2 + 1 \pmod{2} \pmod{x}\),橢圓曲線點個數#\(E(F_{2} )=2\)
\(Y^2+XY=X^3 + 1 X^2 + 1 \pmod{2} \pmod{x^2+x+1}\),橢圓曲線點個數#\(E(F_{2^2})=8\)
\(Y^2+XY=X^3 + 1 X^2 + 1 \pmod{2} \pmod{x^3+x+1}\),橢圓曲線點個數#\(E(F_{2^3})=14\)
\(Y^2+XY=X^3 + 1 X^2 + 1 \pmod{2} \pmod{x^4+x+1}\),橢圓曲線點個數#\(E(F_{2^4})=16\)
\(Y^2+XY=X^3 + 1 X^2 + 1 \pmod{2} \pmod{x^5+x^2+1}\),橢圓曲線點個數#\(E(F_{2^5})=22\)
\(Y^2+XY=X^3 + 1 X^2 + 1 \pmod{2} \pmod{x^6+x+1}\),橢圓曲線點個數#\(E(F_{2^6})=56\)
\(Y^2+XY=X^3 + 1 X^2 + 1 \pmod{2} \pmod{x^7+x+1}\),橢圓曲線點個數#\(E(F_{2^7})=142\)
\(Y^2+XY=X^3 + 1 X^2 + 1 \pmod{2} \pmod{x^8+x^4+x^3+x+1}\),橢圓曲線點個數#\(E(F_{2^8})=288\)
\(Y^2+XY=X^3 + 1 X^2 + 1 \pmod{2} \pmod{x^9+x+1}\),橢圓曲線點個數#\(E(F_{2^9})=518\)
\(Y^2+XY=X^3 + 1 X^2 + 1 \pmod{2} \pmod{x^{10}+x^3+1}\),橢圓曲線點個數#\(E(F_{2^{10}})=968\)

bugmens 發表於 2020-3-22 08:15

美國國家標準暨技術研究院(National Institute of Standards and Technology，NIST)在[url=https://nvlpubs.nist.gov/nistpubs/FIPS/NIST.FIPS.186-4.pdf]FIPS.186-4[/url]文件中推薦數個二元橢圓曲線參數
二元橢圓曲線還分為pseudo-random curve(\(E：y^2+xy=x^3+x^2+b\))和koblitz curve(\(E：y^2+xy=x^3+ax^2+1，a=0\)或1)兩類。
其中koblitz curve可由二階遞迴數列算出在不同\(GF(F_{q^n})\)，\(q=2\)下橢圓曲線點個數，計算方式如下(第107頁)。
1.計算在\(GF(F_q)\)下的橢圓曲線點個數#\(E(F_q)\)。
2.求出frobenius trace \(t=q+1-\)#\(E(F_q)\)。
3.二階遞迴數列\(V_n\)為
　\(V_0=2,V_1=t\)
　\(V_n=V_1V_{n-1}-qV_{n-2}\)
4.在\(GF(F_{q^n})\)的橢圓曲線點個數為#\(E(F_{q^n})=q^n+1-V_n\)

以下就koblitz curve不同\(a\)值進行公式推導[table=98%][tr][td]當\(a=0\)時
\(E：y^2+xy=x^3+0x^2+1\)[/td][td]當\(a=1\)
\(E：y^2+xy=x^3+1x^2+1\)[/td][/tr]
[tr][td]1.在\(GF(F_2)\)下的橢圓曲線點個數
\(E(F_2)=\{\;(\infty,\infty),(0,1),(1,0),(1,1) \}\;\)
#\(E(F_2)=4\)
2.frobenius trace \(t\)
\(t=q+1-\)#\(E(F_2)=2+1-4=-1\)
3.二階遞迴數列\(V_n\)
\(V_0=2,V_1=t=-1\)
\(V_n=V_1V_{n-1}-qV_{n-2}=-V_{n-1}-2V_{n-2}\)
4.在\(GF(F_{2^n})\)的橢圓曲線點個數
#\(E(F_{2^n})=2^n+1-V_n\)
[/td]
[td]1.在\(GF(F_2)\)下的橢圓曲線點個數
\(E(F_2)=\{\;(\infty,\infty),(0,1) \}\;\)
#\(E(F_2)=2\)
2.frobenius trace \(t\)
\(t=q+1-\)#\(E(F_2)=2+1-2=1\)
3.二階遞迴數列\(V_n\)
\(V_0=2,V_1=t=1\)
\(V_n=V_1V_{n-1}-qV_{n-2}=V_{n-1}-2V_{n-2}\)
4.在\(GF(F_{2^n})\)的橢圓曲線點個數
#\(E(F_{2^n})=2^n+1-V_n\)
[/td][/tr][/table]

再解出\(V_n\)一般式
[table=98%][tr][td]求二階遞迴數列\(V_n+V_{n-1}+2V_{n-2}=0\)一般式
特徵方程式為\(x^2+x+2=0\)，得二根\(\displaystyle x={-1\pm \sqrt{7}i \over 2}\)
設一般式\(\displaystyle V_n=c_1 \left({-1-\sqrt{7}i}\over 2\right)^n+c_2 \left({-1+\sqrt{7}i}\over 2\right)^n\)
將\(V_0=2,V_1=-1\)代入得到\(\cases{\displaystyle c_1+c_2=2 \cr c_1 \left({-1-\sqrt{7}i}\over 2\right)+c_2 \left({-1+\sqrt{7}i}\over 2\right)=-1}\)
解出\(c_1=c_2=1\)
一般式\(\displaystyle V_n=\left({-1-\sqrt{7}i}\over 2\right)^n+ \left({-1+\sqrt{7}i}\over 2\right)^n\)[/td][td]求二階遞迴數列\(V_n-V_{n-1}+2V_{n-2}=0\)一般式
特徵方程式為\(x^2-x+2=0\)，得二根\(\displaystyle x={1\pm \sqrt{7}i \over 2}\)
設一般式\(\displaystyle V_n=c_1 \left({1-\sqrt{7}i}\over 2\right)^n+c_2 \left({1+\sqrt{7}i}\over 2\right)^n\)
將\(V_0=2,V_1=1\)代入得到
\(\cases{\displaystyle c_1+c_2=2 \cr c_1 \left({1-\sqrt{7}i}\over 2\right)+c_2 \left({1+\sqrt{7}i}\over 2\right)=1}\)
解出\(c_1=c_2=1\)
一般式\(\displaystyle V_n=\left({1-\sqrt{7}i}\over 2\right)^n+ \left({1+\sqrt{7}i}\over 2\right)^n\)[/td][/tr][/table]



[color=green]當\(a=0\)時\(E：y^2+xy=x^3+0x^2+1\)的二階遞迴數列[/color]
[color=red](%i3)[/color]
[color=blue]V[0]:2;
V[1]:-1;
V[n]:=V[1]*V[n-1]-2*V[n-2];[/color]
[color=red](V[0])[/color]　2
[color=red](V[1])[/color]　\(-1\)
[color=red](V[n])[/color]　\(V_n:=V_1V_{n-1}-2V_{n-2}\)

[color=green]利用二階遞迴數列計算\(GF(F_{2^2})\)到\(GF(F_{2^{10}})\)橢圓曲線點個數，符合上一篇執行結果[/color]
[color=red](%i4)[/color]
[color=blue]for n:2 thru 10 do
  (print("V[",n,"]=V[1]*V[",n-1,"]-2V[",n-2,"]=(",V[1],")*(",V[n-1],")-2*(",V[n-2],")=",V[n],
           "，橢圓曲線點個數#E(F 2"^n,")=2"^n,"+1-V[",n,"]=",2^n+1-V[n])
  );[/color]
\(V[2]=V[1]*V[1]-2V[0]=(-1)*(-1)-2*(2)=-3\)，橢圓曲線點個數#\(E(F_{2^2})=2^2+1-V[2]=8\)
\(V[3]=V[1]*V[2]-2V[1]=(-1)*(-3)-2*(-1)=5\)，橢圓曲線點個數#\(E(F_{2^3})=2^3+1-V[3]=4\)
\(V[4]=V[1]*V[3]-2V[2]=(-1)*(5)-2*(-3)=1\)，橢圓曲線點個數#\(E(F_{2^4})=2^4+1-V[4]=16\)
\(V[5]=V[1]*V[4]-2V[3]=(-1)*(1)-2*(5)=-11\)，橢圓曲線點個數#\(E(F_{2^5})=2^5+1-V[5]=44\)
\(V[6]=V[1]*V[5]-2V[4]=(-1)*(-11)-2*(1)=9\)，橢圓曲線點個數#\(E(F_{2^6})=2^6+1-V[6]=56\)
\(V[7]=V[1]*V[6]-2V[5]=(-1)*(9)-2*(-11)=13\)，橢圓曲線點個數#\(E(F_{2^7})=2^7+1-V[7]=116\)
\(V[8]=V[1]*V[7]-2V[6]=(-1)*(13)-2*(9)=-31\)，橢圓曲線點個數#\(E(F_{2^8})=2^8+1-V[8]=288\)
\(V[9]=V[1]*V[8]-2V[7]=(-1)*(-31)-2*(13)=5\)，橢圓曲線點個數#\(E(F_{2^9})=2^9+1-V[9]=508\)
\(V[10]=V[1]*V[9]-2V[8]=(-1)*(5)-2*(-31)=57\)，橢圓曲線點個數#\(E(F_{2^{10}})=2^{10}+1-V[10]=968\)
[color=red](%o4)[/color]　done

[color=green]當\(a=0\)時，NIST建議以下4種koblitz橢圓曲線[/color]
[color=red](%i5)[/color]　[color=blue]Degree:[233,283,409,571];[/color]
[color=red](Degree)[/color]　\(\left[233,283,409,571 \right]\)

[color=green]計算橢圓曲線點個數[/color]
[color=red](%i6)[/color]
[color=blue]for d in Degree do
  (print("K-",d,"橢圓曲線"),
   print("V[",d,"]=(",(-1-sqrt(7)*%i)/2,")"^d,"+(",(-1+sqrt(7)*%i)/2,")"^d,"=",V[d]),
   print("橢圓曲線點個數#E(F 2"^d,")=2"^d,"+1-V[",d,"]=",2^d+1-V[d])
  );[/color]
\(K-233\)橢圓曲線
\( \displaystyle V[\;233]\;=\left({-\sqrt{7}\%i-1 \over 2}\right)^{233}+\left({\sqrt{7}\%i-1\over 2}\right)^{233}=-137381546011108235394987299651366779\)
橢圓曲線點個數#\(E(F_{2^{233}})=2^{233}+1-V[\;233]\;=13803492693581127574869511724554051042283763955449008505312348098965372\)
\(K-283\)橢圓曲線
\( \displaystyle V[\;283]\;=\left({-\sqrt{7}\%i-1\over 2}\right)^{283}+\left({\sqrt{7}\%i-1\over 2}\right)^{283}=7777244870872830999287791970962823977569917\)
橢圓曲線點個數#\(E(F_{2^{283}})=2^{283}+1-V[\;283]\;=\)
\(15541351137805832567355695254588151253139246935172245297183499990119263318817690415492\)
\(K-409\)橢圓曲線
\( \displaystyle V[\;409]\;=\left({-\sqrt{7}\%i-1\over 2}\right)^{409}+\left({\sqrt{7}\%i-1\over 2}\right)^{409}=10457288737315625927447685387048320737638796957687575791173829\)
橢圓曲線點個數#\(E(F_{2^{409}})=2^{409}+1-V[\;409]\;=\)
\(1322111937580497197903830616065542079656809365928562438569297580091522845156996764202693033831109832056385466362470925434684\)
\(K-571\)橢圓曲線
\( \displaystyle V[\;571]\;=\left({-\sqrt{7}\%i-1\over 2}\right)^{571}+\left({\sqrt{7}\%i-1\over 2}\right)^{571}=-148380926981691413899619140297051490364542574180493936232912339534208516828973111459843\)
橢圓曲線點個數#\(E(F_{2^{571}})=2^{571}+1-V[\;571]\;=\)
\(7729075046034516689390703781863974688597854659412869997314470502903038284579120849072535914090826847338826851203301405845094699896266469247718729686468370014222934741106692\)
[color=red](%o6)[/color]　done

[color=green]刪除\(V_n\)數列定義[/color]
[color=red](%i7)[/color]　[color=blue]kill(V);[/color]
[color=red](%o7)[/color]　done

[color=green]當\(a=1\)時\(E：y^2+xy=x^3+1x^2+1\)的二階遞迴數列[/color]
[color=red](%i10)[/color]
[color=blue]V[0]:2;
V[1]:1;
V[n]:=V[1]*V[n-1]-2*V[n-2];[/color]
[color=red](V[0])[/color]　2
[color=red](V[1])[/color]　1
[color=red](V[n])[/color]　\(V_n:=V_1V_{n-1}-2V_{n-2}\)

[color=green]利用二階遞迴數列計算\(GF(F_{2^2})\)到\(GF(F_{2^{10}})\)橢圓曲線點個數，符合上一篇執行結果[/color]
[color=red](%i11)[/color]
[color=blue]for n:2 thru 10 do
  (print("V[",n,"]=V[1]*V[",n-1,"]-2V[",n-2,"]=(",V[1],")*(",V[n-1],")-2*(",V[n-2],")=",V[n],
           "，橢圓曲線點個數#E(F 2"^n,")=2"^n,"+1-V[",n,"]=",2^n+1-V[n])
  );[/color]
\(V[2]=V[1]*V[1]-2V[0]=(1)*(1)-2*(2)=-3\)，橢圓曲線點個數#\(E(F_{2^2})=2^2+1-V[2]=8\)
\(V[3]=V[1]*V[2]-2V[1]=(1)*(-3)-2*(1)=-5\)，橢圓曲線點個數#\(E(F_{2^3})=2^3+1-V[3]=14\)
\(V[4]=V[1]*V[3]-2V[2]=(1)*(-5)-2*(-3)=1\)，橢圓曲線點個數#\(E(F_{2^4})=2^4+1-V[4]=16\)
\(V[5]=V[1]*V[4]-2V[3]=(1)*(1)-2*(-5)=11\)，橢圓曲線點個數#\(E(F_{2^5})=2^5+1-V[5]=22\)
\(V[6]=V[1]*V[5]-2V[4]=(1)*(11)-2*(1)=9\)，橢圓曲線點個數#\(E(F_{2^6})=2^6+1-V[6]=56\)
\(V[7]=V[1]*V[6]-2V[5]=(1)*(9)-2*(11)=-13\)，橢圓曲線點個數#\(E(F_{2^7})=2^7+1-V[7]=142\)
\(V[8]=V[1]*V[7]-2V[6]=(1)*(-13)-2*(9)=-31\)，橢圓曲線點個數#\(E(F_{2^8})=2^8+1-V[8]=288\)
\(V[9]=V[1]*V[8]-2V[7]=(1)*(-31)-2*(-13)=-5\)，橢圓曲線點個數#\(E(F_{2^9})=2^9+1-V[9]=518\)
\(V[10]=V[1]*V[9]-2V[8]=(1)*(-5)-2*(-31)=57\)，橢圓曲線點個數#\(E(F_{2^{10}})=2^{10}+1-V[10]=968\)
[color=red](%o11)[/color]　done

[color=green]當\(a=1\)時，NIST建議以下1種koblitz橢圓曲線[/color]
[color=red](%i12)[/color]　[color=blue]d:163;[/color]
[color=red](d)[/color]　163

[color=green]計算橢圓曲線點個數[/color]
[color=red](%i15)[/color]
[color=blue]print("K-",d,"橢圓曲線")$
print("V[",d,"]=(",(1-sqrt(7)*%i)/2,")"^d,"+(",(1+sqrt(7)*%i)/2,")"^d,"=",V[d])$
print("橢圓曲線點個數#E(F 2"^d,")=2"^d,"+1-V[",d,"]=",2^d+1-V[d])$[/color]
\(K-163\)橢圓曲線
\(\displaystyle V[\;163]\;=\left( {1-\sqrt{7}\%i \over 2} \right)^{163}+\left( {\sqrt{7}\%i+1 \over 2} \right)^{163}=-4845466632539410776804317\)
橢圓曲線點個數#\(E(F_{2^{163}})=2^{163}+1-V[\;163]\;=11692013098647223345629483507196896696658237148126\)

bugmens 發表於 2020-3-28 14:56

橢圓曲線加密的安全性倚賴橢圓曲線離散對數問題(Elliptic Curve Discrete Logarithm Problem，ECDLP)的困難度上，橢圓曲線離散對數問題敘述如下：

令\(E\)為基於二元體\(F_{2^n}\)的橢圓曲線，\(P,Q\)為橢圓曲線上的點坐標(\(P,Q\in E(F_{2^n})\))且\(Q=kP,k\in N\)，當\(n\)很大時，計算\(k\)值非常困難。

解決橢圓曲線離散對數問題的方法會在另一篇文章討論，這裡僅示範暴力搜尋法找出\(k\)值。
可惜在C:\maxima-5.43.0\share\maxima\5.43.0\share\contrib\elliptic_curves\ec2.lisp程式碼中沒有解橢圓曲線離散對數問題的指令



[color=green]要先載入elliptic_curves.mac才能呼叫ec2_set_curve,ec2_add[/color]
[color=red](%i1)[/color]　[color=blue]load("elliptic_curves.mac");[/color]
[color=red](%o1)[/color]　[i]C:/maxima-5.43.0/share/maxima/5.43.0/share/contrib/elliptic_curves/elliptic_curves.mac[/i]

[color=green]\(y^2+xy=x^3+ax^2+b\pmod{2}\pmod{x^4+x+1}\)的橢圓曲線
\(a=g^8=x^2+1=101_2=5_{10}\)，\(b=g^2=x^2=100_2=4_{10}\)[/color]
[color=red](%i3)[/color]
[color=blue][ReductionPoly,a,b]:[x^4+x+1,5,4];
ec2_set_curve(ReductionPoly,a,b);[/color]
[color=red](%o2)[/color]　\([x^4+x+1,5,4]\)
[color=red](%o3)[/color]　[u]true[/u]

[color=green]\(P,Q\)兩點坐標
\(P=[g^3,g^9]=[x^3,x^3+x]=[1000_2,1010_2]=[8_{10},10_{10}]\)
\(Q=[g^{13},g^1]=[x^3+x^2+1,x]=[1101_2,10_2]=[13_{10},2_{10}]\)[/color]
[color=red](%i5)[/color]
[color=blue]P:[8,10];
Q:[13,2];[/color]
[color=red](P)[/color]　\([8,10]\)
[color=red](Q)[/color]　\([13,2]\)

[color=green]計算\(2P,3P,\ldots,kP\)，看何時等於\(Q\)[/color]
[color=red](%i7)[/color]
[color=blue]kP: P$
for k:2 thru 15 do
  (print(k," P=",kP:ec2_add(P,kP)),
   if kP=Q then
      (return(k))
  );[/color]
\(2P=[2,6]\)
\(3P=[10,13]\)
\(4P=[5,4]\)
\(5P=[11,5]\)
\(6P=[1,0]\)
\(7P=[13,15]\)
\(8P=[0,2]\)
\(9P=[13,2]\)
[color=red](%o7)[/color]　9

查看完整版本: 用Maxima學密碼學－橢圓曲線