問題敘述
令\(N=p^8q^5=80171134270603235454533\)是未知分解的整數，\(r>s\)和\(gcd(r,s)=1\)，給定\(N\)為輸入，能在多項式時間\(\log N\)和\(r=\Omega(\log^3 max(p,q))\)情況下回復質因數
步驟1：將\(r\)和\(s\)表示成\(\cases{r=u\cdot \alpha+a\cr s=u\cdot \beta+b}\)
設矩陣\(M=\left[\matrix{\lfloor\;r^{1/3}\rfloor\;&-s\cr0&r}\right]\)，經LLL化簡後短的非零向量\(v=(\lfloor\;r^{1/3}\rfloor\;\cdot\alpha,\gamma)\)，其中\(\gamma=-s\cdot\alpha+r\cdot\beta\)，\(\beta\in\mathbb{Z}\) 從\(v\)向量得到\(\alpha,\beta\)值，符合\(\alpha>0\)且\(0\le\beta\le\alpha\) Case 1：\(\beta=0\)或(\(\beta\ne0\)和\(\displaystyle\lfloor\;\frac{r}{\alpha}\rfloor\;\le\frac{s}{\beta}\)) 　\(\displaystyle u=\lfloor\;\frac{r}{\alpha}\rfloor\;\) 　\(a=r-u\cdot\alpha\) 　\(b=s-u\cdot\beta\) Case 2：\(\beta\ne0\)和\(\displaystyle\lfloor\;\frac{r}{\alpha}\rfloor\;>\frac{s}{\beta}\) 　\(\displaystyle u=\lceil\;\frac{r}{\alpha}\rceil\;\) 　\(a=r-u\cdot\alpha\) 　\(b=s-u\cdot\beta\)	設矩陣\(M=\left[\matrix{\lfloor\;8^{1/3}\rfloor\;&-5\cr0&8}\right]=\left[\matrix{2&-5\cr0&8}\right]\)，經LLL化簡後\(B=\left[\matrix{2&3\cr4&-2}\right]\) 原本短的非零向量\(v=(2,3)\)無法求得合適的\(u,\alpha,\beta,a,b\) 改用\(v=(4,-2)=(\lfloor\;r^{1/3}\rfloor\;\cdot \alpha,\gamma)\)，其中\(\gamma=-s\cdot\alpha+r\cdot\beta\)，\(\beta\in\mathbb{Z}\) \(4=\lfloor\;8^{1/3}\rfloor\;\cdot \alpha\)，得\(\alpha=2\) \(-2=-5\cdot 2+8\cdot\beta\)，得\(\beta=1\) 符合\(\alpha>0\)且\(0\le\beta\le\alpha\) \(\displaystyle \lfloor\;\frac{r}{\alpha}\rfloor\;=\lfloor\;\frac{8}{2}\rfloor\;=4\)，\(\displaystyle \frac{s}{\beta}=\frac{5}{1}=5\) 符合Case 1:\(\beta=0\)或\((\beta\ne0\)和\(\displaystyle \lfloor\;\frac{r}{\alpha}\rfloor\;\le\frac{s}{\beta})\) 　\(\displaystyle u=\lfloor\;\frac{r}{\alpha}\rfloor\;=\lfloor\;\frac{8}{2}\rfloor\;=4\) 　\(a=r-u\cdot\alpha=8-4\cdot2=0\) 　\(b=s-u\cdot\beta=5-4\cdot1=1\) \(\cases{8=4\cdot2+0\cr5=4\cdot1+1}\)
步驟2：根據\(a,b\ge0\)(Case1)使用BDH方法
\(a,b\ge0\)(Case1) \(\displaystyle N=p^rq^s=p^{u\cdot \alpha+a}\cdot q^{u\cdot\beta+b}=(p^{\alpha}q^{\beta})^u\cdot p^aq^b=P^u\cdot Q\)，其中\(P=p^{\alpha}q^{\beta}\)和\(Q=p^aq^b\)	\(N=p^8q^5=P^4\cdot Q\)，其中\(P=p^2q^1=p^2q\)和\(Q=p^0q^1=q\)
步驟2-1.設同餘方程式，計算參數\(X\)
設同餘方程式\(f(x)=(V+x)^u\pmod{P^u}\) 其中\(V\)是使得\(P=V+x_0\)的整數和\(V\)的高位元和\(P\)相同，能找到上限\(X= P\cdot Q^{-1/u}\)較小的解\(x_0\)，其中\(|\;x_0|\;< P\cdot Q^{-1/u}\)	\(V=\lfloor\;N^{1/4}\rfloor\;=532113\) \(P=p^2q=61^2\cdot 53=197213\)，\(Q=q=53\) \(X=P\cdot Q^{-1/4}=73091\) 但這樣\(V\)和\(X\)所需的矩陣\(M\)要非常大才能找到較小的解\(x_0\) 設同餘方程式\(f(x)=(V+x)^4\pmod{P^4}\) 取\(V=197000\)，\(V\)的高位元和\(P\)相同，上限\(X=1000\) 只要維度9的矩陣\(M\)經LLL化簡就能找到較小的解\(x_0=213\)，其中\(|\;x_0|\;<1000\)
步驟2-2.產生矩陣\(M\)
\(g_{i,k}(xX)=N^{m-k}(xX)^i f(xX)^k\) \(g_{i,k}(xX)\)，\(i=0,1,\ldots,u-1\)和\(k=0,1,\ldots,m-1\) \(g_{j,m}(xX)\)，\(j=0,1,\ldots,d-um-1\) 注意到對所有\(i,k\)，\(g_{i,k}(x_0)\equiv 0\pmod{P^{um}}\) 注意到對所有\(j,m\)，\(g_{j,m}(x_0)\equiv 0\pmod{P^{um}}\)	原本矩陣\(M\)需要\(d=2\cdot u\cdot(u+1)=2\cdot 4\cdot5=40\)，改取維度\(d=9\)減少LLL執行時間 設\(P^c>Q\)，取\(c=1\) \(\displaystyle m=\lfloor\;\frac{d}{u+c}-\frac{1}{2}\rfloor\;=\lfloor\;\frac{9}{4+1}-\frac{1}{2}\rfloor\;=1\)太小，改取\(m=2\) \(i=0,k=0,g_{0,0}(xX)=N^2\cdot1\cdot1=N^2=\bbox[border:1px solid black]{N^2}\) \(i=1,k=0,g_{1,0}(xX)=N^2(xX)\cdot1=N^2Xx=\bbox[border:1px solid black]{N^2X}x\) \(i=2,k=0,g_{2,0}(xX)=N^2(xX)^2\cdot1=N^2X^2x^2=\bbox[border:1px solid black]{N^2X^2}x^2\) \(i=3,k=0,g_{3,0}(xX)=N^2(xX)^3\cdot1=N^2X^3x^3=\bbox[border:1px solid black]{N^2X^3}x^3\) \(i=0,k=1,g_{0,1}(xX)=N\cdot1\cdot f(xX)=N(V+Xx)^4\) \(=NV^4+4NV^3Xx+6NV^2X^2x^2+4NVX^3x^3+\bbox[border:1px solid black]{NX^4}x^4\) \(i=1,k=1,g_{1,1}(xX)=N(xX)f(xX)=NXx(V+Xx)^4\) \(=NV^4Xx+4NV^3X^2x^2+6NV^2X^3x^3+4NVX^4x^4+\bbox[border:1px solid black]{NX^5}x^5\) \(i=2,k=1,g_{2,1}(xX)=N(xX)^2f(xX)=NX^2x^2(V+Xx)^4\) \(=NV^4X^2x^2+4NV^3X^3x^3+6NV^2X^4x^4+4NVX^5x^5+\bbox[border:1px solid black]{NX^6}x^6\) \(i=3,k=1,g_{3,1}(xX)=N(xX)^3f(xX)=NX^3x^3(V+Xx)^4\) \(=NV^4X^3x^3+4NV^3X^4x^4+6NV^2X^5x^5+4NVX^6x^6+{NX^7}x^7\) \(j=0,g_{0,2}(xX)=1\cdot f(xX)^2=(V+Xx)^8\) \(=V^8+8V^7Xx+28V^6X^2x^2+56V^5X^3x^3+70V^4X^4x^4+56V^3X^5x^5+28V^2X^6x^6+8VX^7x^7+\bbox[border:1px solid black]{X^8}x^8\)
\(M=\matrix{&\matrix{1&　x　&　x^2&　x^3&　x^4&　x^5&　x^6&　x^7&　x^8}\cr \matrix{g_{0,0}(xX)\cr g_{1,0}(xX)\cr g_{2,0}(xX)\cr g_{3,0}(xX)\cr g_{0,1}(xX)\cr g_{1,1}(xX)\cr g_{2,1}(xX)\cr g_{3,1}(xX)\cr g_{0,2}(xX)}&\left[\matrix{N^2&&&&&&&&\cr 0&N^2X&&&&&&&\cr 0&0&N^2X^2&&&&&&\cr 0&0&0&N^2X^3&&&&&\cr *&*&*&*&NX^4&&&&\cr 0&*&*&*&*&NX^5&&&\cr 0&0&*&*&*&*&NX^6&&\cr 0&0&0&*&*&*&*&NX^7&\cr *&*&*&*&*&*&*&*&X^8}\right]}\) *代表非零數字
步驟2-3.經LLL化簡後的短向量產生不需要同餘\(P^{um}\)的方程式，得到公鑰\(N\)的因數\(P\)和\(Q\)
矩陣\(M\)經LLL化簡為\(B\) \(B=LLL(M)\) lattice經LLL化簡後第一列\(b_1\)為整個lattice中較短向量 所形成的方程式不需要再同餘\(P^{um}\)	本範例第一列\(b_1\)無法得到正確答案 lattice經LLL化簡後第二列\(b_2\)為整個lattice中較短向量 \(b_2=[520814821077069868175677184373000000000000,\) 　\(-2430603060877548431773082696764000000000000,\) 　\(-68074195511300390211470172818000000000000,\) 　\(-817458291481540650461827996000000000000,\) 　\(-5507602630239396764545467000000000000,\) 　\(-22691467064000000000000000000000000,\) 　\(-57592556000000000000000000000000,\) 　\(-83528000000000000000000000000,\) 　\(-53000000000000000000000000]\) 產生不需要同餘\(P^{um}=P^8\)的方程式 \(h(x)=520814821077069868175677184373000000000000\) 　\(\displaystyle-2430603060877548431773082696764000000000000\left(\frac{x}{1000}\right)\) 　\(\displaystyle-68074195511300390211470172818000000000000\left(\frac{x}{1000}\right)^2\) 　\(\displaystyle-817458291481540650461827996000000000000\left(\frac{x}{1000}\right)^3\) 　\(\displaystyle-5507602630239396764545467000000000000\left(\frac{x}{1000}\right)^4\) 　\(\displaystyle-22691467064000000000000000000000000\left(\frac{x}{1000}\right)^5\) 　\(\displaystyle-57592556000000000000000000000000\left(\frac{x}{1000}\right)^6\) 　\(\displaystyle-83528000000000000000000000000\left(\frac{x}{1000}\right)^7\) 　\(\displaystyle-53000000000000000000000000\left(\frac{x}{1000}\right)^8\) \(h(x)=-53(-213+x)(197000+x)^4(394213+x)(77701967369+394000x+x^2)\) 得到正整數解\(x=213\) 得到因數\(P=V+x=197000+213=197213\) \(N\)可分解成\(80171134270603235454533=197213^4\cdot53\) \(P=197213,Q=53\)
步驟3.從\(P,Q\)求出\(N\)的質因數\(p,q\)
已知\(P=p^{\alpha}q^{\beta}\)，\(Q=p^aq^b\)，\(\gamma=a\beta-b\alpha\)，則\(\cases{p=Q^{\displaystyle\frac{\beta}{\gamma}} \cdot P^{\displaystyle\frac{-b}{\gamma}}\cr q=Q^{\displaystyle\frac{-\alpha}{\gamma}} \cdot P^{\displaystyle\frac{a}{\gamma}}}\)	\(\alpha=2,\beta=1,a=0,b=1,\gamma=a\beta-b\alpha=-2\) \(\cases{p=Q^{\displaystyle\frac{\beta}{\gamma}} \cdot P^{\displaystyle\frac{-b}{\gamma}}=53^{\displaystyle\frac{1}{-2}}\cdot 197213^{\displaystyle\frac{-1}{-2}}=61\cr q=Q^{\displaystyle\frac{-\alpha}{\gamma}} \cdot P^{\displaystyle\frac{a}{\gamma}}=53^{\displaystyle\frac{-2}{-2}}\cdot 197213^{\displaystyle\frac{0}{-2}}=53}\)

請下載LLL.zip，解壓縮後將LLL.mac放到C:\maxima-5.49.0\share\maxima\5.49.0\share目錄下
要先載入LLL.mac才能使用LLL指令
(%i1)　load("LLL.mac");
(%o1)　C:/maxima-5.49.0/share/maxima/5.49.0/share/LLL.mac

要因數分解的公鑰\(N\)
(%i2)　N:80171134270603235454533;
(%o2)　\(80171134270603235454533\)

\(N=p^rq^s\)，因數\(p\)的次方\(r\)，因數\(q\)的次方\(s\)
(%i4)
r:8;
s:5;
(%o3)　\(8\)
(%o4)　\(5\)

將\(r,s\)表示成\(\cases{r=u\alpha+a\cr s=u\beta+b}\)
(%i10)
u:4;
alpha:2;
beta:1;
a:0;
b:1;
gamma:a*beta-b*alpha;
(%o5)　\(4\)
(%o6)　\(2\)
(%o7)　\(1\)
(%o8)　\(0\)
(%o9)　\(1\)
(%o10)　\(-2\)

假設\(Q< P^c\)
(%i11)　c:1;
(%o11)　\(1\)

矩陣維度\(d\)，按照BDH方法維度\(d=40\)，改成\(d=9\)但解的上限\(X\)也變小
(%i13)
d:2*u*(u+1);
d:9;
(%o12)　\(40\)
(%o13)　\(9\)

參數\(m\)，原本\(m=1\)無法得到正確答案，改取\(m=2\)
(%i15)
m:floor(d/(u+c)-1/2);
m:2;
(%o14)　\(1\)
(%o15)　\(2\)

希望能找到\(|\;x|\;<X\)，\(f(x)\equiv0\pmod{P^u}\)
(%i16)　X:1000;
(%o16)　\(1000\)

\(V\)的高位元和\(P\)相同，當作\(P\)的近似值
(%i17)　V:197000;
(%o17)　\(197000\)

同餘方程式\(f(x)=(V+x)^u \equiv0\pmod{P^u}\)
(%i18)　fx: ('V+x)^u;
(%o18)　\((x+V)^4\)

將\(x\)以\(Xx\)代替，\(f(Xx)=(V+Xx)^u \equiv0\pmod{P^u}\)
(%i19)　fXx:subst(x=x*'X,fx);
(%o19)　\((Xx+V)^4\)

以x升冪排序顯示
(%i20)　powerdisp:true;
(%o20)　true

\(g(xX)\)多項式
(%i21)　gxX:[];
(%o21)　\([]\)

產生\(g_{i,k}(xX)=N^{m-k}(Xx)^if^k(xX)\)多項式，\(i=0,\ldots,u-1\)，\(k=0,\ldots,m-1\)
(%i22)
for k:0 thru m-1 do
  (for i:0 thru u-1 do
     (print("i=",i,",k=",k,",g",i,",",k,"(xX)=N"^(m-k),"*","(xX)"^i,"*","f(xX)"^k,"=",gik:'N^(m-k)*(x*'X)^i*fXx^k,"=",expand(gik)),
      gxX:append(gxX,[gik])
     )
  );
\(i=0,k=0,g0,0(xX)=N^2*1*1=N^2=N^2\)
\(i=1,k=0,g1,0(xX)=N^2*(xX)*1=N^2Xx=N^2Xx\)
\(i=2,k=0,g2,0(xX)=N^2*(xX)^2*1=N^2X^2x^2=N^2X^2x^2\)
\(i=3,k=0,g3,0(xX)=N^2*(xX)^3*1=N^2X^3x^3=N^2X^3x^3\)
\(i=0,k=1,g0,1(xX)=N*1*f(xX)=N(V+Xx)^4=NV^4+4NV^3Xx+6NV^2X^2x^2+4NVX^3x^3+NX^4x^4\)
\(i=1,k=1,g1,1(xX)=N*(xX)*f(xX)=NXx(V+Xx)^4=NV^4Xx+4NV^3X^2x^2+6NV^2X^3x^3+4NVX^4x^4+NX^5x^5\)
\(i=2,k=1,g2,1(xX)=N*(xX)^2*f(xX)=NX^2x^2(V+Xx)^4=NV^4X^2x^2+4NV^3X^3x^3+6NV^2X^4x^4+4NVX^5x^5+NX^6x^6\)
\(i=3,k=1,g3,1(xX)=N*(xX)^3*f(xX)=NX^3x^3(V+Xx)^4=NV^4X^3x^3+4NV^3X^4x^4+6NV^2X^5x^5+4NVX^6x^6+NX^7x^7\)
(%o22)　done

產生\(g_{j,m}(xX)=x^jf^m(xX)\)多項式，\(j=0,\ldots,d-mu-1\)
(%i23)
for j:0 thru d-m*u-1 do
  (print("j=",j,",g",j,",",m,"(xX)=","(xX)"^j,"*f(xX)"^m,"=",gim: (x*'X)^j*fXx^m,"=",expand(gim)),
   gxX:append(gxX,[gim])
  );
\(j=0,g0,2(xX)=1*f(xX)^2=(V+Xx)^8=V^8+8V^7Xx+28V^6X^2x^2+56V^5X^3x^3+70V^4X^4x^4+56V^3X^5x^5+28V^2X^6x^6+8VX^7x^7+X^8x^8\)
(%o23)　done

全部的\(g(xX)\)多項式
(%i24)　gxX;
(%o24)　\([N^2,N^2Xx,N^2X^2x^2,N^2X^3x^3,N(V+Xx)^4,NXx(V+Xx)^4,NX^2x^2(V+Xx)^4,NX^3x^3(V+Xx)^4,(V+Xx)^8]\)

\(x^1,\ldots,x^{d-1}\)
(%i25)　xpower:create_list(x^i,i,1,d-1);
(%o25)　\([x,x^2,x^3,x^4,x^5,x^6,x^7,x^8]\)

取\(g(xX)\)多項式係數(常數項在最後一行)
(%i26)　M:augcoefmatrix(gxX,xpower);
(%o26)　\(\left[\matrix{0&0&0&0&0&0&0&0&N^2\cr
N^2X&0&0&0&0&0&0&0&0\cr
0&N^2X^2&0&0&0&0&0&0&0\cr
0&0&N^2X^3&0&0&0&0&0&0\cr
4NV^3X&6NV^2X^2&4NVX^3&NX^4&0&0&0&0&NV^4\cr
NV^4X&4NV^3X^2&6NV^2X^3&4NVX^4&NX^5&0&0&0&0\cr
0&NV^4X^2&4NV^3X^3&6NV^2X^4&4NVX^5&NX^6&0&0&0\cr
0&0&NV^4X^3&4NV^3X^4&6NV^2X^5&4NVX^6&NX^7&0&0\cr
8V^7X&28V^6X^2&56V^5X^3&70V^4X^4&56V^3X^5&28V^2X^6&8VX^7&X^8&V^8}\right]\)

將常數項移到第一行
(%i27)　M:addcol(col(M,d),submatrix(M,d));
(%o27)　\(\left[\matrix{N^2&0&0&0&0&0&0&0&0\cr
0&N^2X&0&0&0&0&0&0&0\cr
0&0&N^2X^2&0&0&0&0&0&0\cr
0&0&0&N^2X^3&0&0&0&0&0\cr
NV^4&4NV^3X&6NV^2X^2&4NVX^3&NX^4&0&0&0&0\cr
0&NV^4X&4NV^3X^2&6NV^2X^3&4NVX^4&NX^5&0&0&0\cr
0&0&NV^4X^2&4NV^3X^3&6NV^2X^4&4NVX^5&NX^6&0&0\cr
0&0&0&NV^4X^3&4NV^3X^4&6NV^2X^5&4NVX^6&NX^7&0\cr
V^8&8V^7X&28V^6X^2&56V^5X^3&70V^4X^4&56V^3X^5&28V^2X^6&8VX^7&X^8}\right]\)

將\(N=80171134270603235454533,V=197213,X=1000\)代入矩陣\(M\)
(%i28)　M:ev(M,[N=N,V=V,X=X]);
(%o28)　\(\left[\matrix{6427410770235092574143943139305425635110248089&0&0&0&0&0&0&0&0\cr
0&6427410770235092574143943139305425635110248089000&0&0&0&0&0&0&0\cr
0&0&6427410770235092574143943139305425635110248089000000&0&0&0&0&0&0\cr
0&0&0&6427410770235092574143943139305425635110248089000000000&0&0&0&0&0\cr
120748830390373400001175677184373000000000000&2451752901327378680226917303236000000000000&18668169299447045788529827182000000000000&63174853805235349538172004000000000000&80171134270603235454533000000000000&0&0&0&0\cr
0&120748830390373400001175677184373000000000000000&2451752901327378680226917303236000000000000000&18668169299447045788529827182000000000000000&63174853805235349538172004000000000000000&80171134270603235454533000000000000000&0&0&0\cr
0&0&120748830390373400001175677184373000000000000000000&2451752901327378680226917303236000000000000000000&18668169299447045788529827182000000000000000000&63174853805235349538172004000000000000000000&80171134270603235454533000000000000000000&0&0\cr
0&0&0&120748830390373400001175677184373000000000000000000000&2451752901327378680226917303236000000000000000000000&18668169299447045788529827182000000000000000000000&63174853805235349538172004000000000000000000000&80171134270603235454533000000000000000000000&0\cr
2268453123948987361000000000000000000000000&92119923815187304000000000000000000000000&1636648392655612000000000000000000000000&16615719722392000000000000000000000000&105429693670000000000000000000000000&428140888000000000000000000000000&1086652000000000000000000000000&1576000000000000000000000000&1000000000000000000000000}\right]\)

LLL化簡
(%i29)　B: LLL(M);
(%o29)　\(\left[\matrix{2268453123948987361000000000000000000000000&92119923815187304000000000000000000000000&1636648392655612000000000000000000000000&16615719722392000000000000000000000000&105429693670000000000000000000000000&428140888000000000000000000000000&1086652000000000000000000000000&1576000000000000000000000000&1000000000000000000000000\cr
520814821077069868175677184373000000000000&-2430603060877548431773082696764000000000000&-68074195511300390211470172818000000000000&-817458291481540650461827996000000000000&-5507602630239396764545467000000000000&-22691467064000000000000000000000000&-57592556000000000000000000000000&-83528000000000000000000000000&-53000000000000000000000000\cr
119574028217671068321357761887635110248089&-1120941543841003168053234143016000000000000&2618519389228227254415838318708000000000000&39977022196844180948953767576000000000000&287653869286346057041819502000000000000&1202647754392000000000000000000000000&3052405468000000000000000000000000&4426984000000000000000000000000&2809000000000000000000000000\cr
-241520619786717295077270256278597856642032&1006233591109786490302081721964050089220000&791289620910282234486801424154248089000000&-1025261299321869239983536773032000000000000&295501763176439064773555173486000000000000&985527115520611944767267972000000000000000&1714483878272561581454944000000000000000000&6207338904904808000000000000000000000000&3938666817833000000000000000000000000\cr
146111802828159872978714311180467622752981&-337364242705630748590601841116939841131000&-899218919995713091129794664660248089000000&-2698981263760596036060949910100000000000000&-3289903627281272224831697418325000000000000&-1009722386745134958108177997000000000000000&-1714560548593106805454944000000000000000000&-6207450101888920000000000000000000000000&-3938737374295000000000000000000000000\cr
61367641518645839679950446956111346617966&-57506707444030967438209001123160337309000&-542612286676673217701391991592248089000000&-3026902524618602725326328336604000000000000&2590332438445698909444733467517000000000000&-961953111791700696190903414000000000000000&-1714408988253107429454944000000000000000000&-6207230289938232000000000000000000000000&-3938597899707000000000000000000000000\cr
15067948572848914227471586284476558432910&-141458761630094005226905014839668317681000&335539551509200464614265223688160566000000&-94901188582953634568900965644267000000000&780356555960664279505166852075000000000000&-2333400442619330164545029877000000000000000&1004715999866625459367336000000000000000000&3758047457599007913741000000000000000000000&1172394756143066893000000000000000000000000\cr
-203445131193957092545439571131643721426395&1241053695755686421865190490359020206491000&-1100111943534622115960142187693663193000000&-1156898039491760619610370540935110000000000&450136837007525240038159764818000000000000&-1635079766422598279193272752000000000000000&-199868442363979553097839000000000000000000&846535899130441456508000000000000000000000&-3865583546108321988000000000000000000000000\cr
-233825991062131604151461386720904468484216&929448995310335846542698455837482677551000&1028704530536093855735204722429708747000000&-1376103751198842125959550070103822000000000&829839950910549569092835862550000000000000&2286782227716711908322343822000000000000000&-3321284626741256170546210000000000000000000&4165942424452207105995000000000000000000000&-760406771341748937000000000000000000000000}\right]\)

找出\(N\)的因數\(P,Q\)
(%i30)
for i:1 thru d do
  (print("第",i,"列向量B[",i,"]=",B[ i ]),
   print("產生不需要同餘P"^"um","=p"^(u*m),"的方程式h(x)"),
   printList:["h(x)=",B[ i ][1]],
   for j:2 thru d do
    (if B[ i ][j]>=0 then printList:append(printList,["+"]),/*若係數為正則補印+號*/
     printList:append(printList,[B[ i ][j],"(",x/X,")"^(j-1)])
    ),
   apply(print,printList),/*再用apply(print,)將全部內容印在同一行*/
   print("h(x)=",hx:sum(B[ i ][j+1]*(x/X)^j,j,0,d-1),"=",factor(hx)),
   print("正整數解為",posIntRoot:sublist(solve(hx,x),lambda([x],integerp(rhs(x)) and rhs(x)>0))),
   if length(posIntRoot)>0 then/*若有正整數解*/
     (for root in posIntRoot do
        (x:rhs(root),
         print("當正整數解x=",x,"時，可能的因數P=V+x=",V,"+",x,"=",P: V+x),
         if mod(N,P^u)=0 then
            (print("N可被P"^u,"整除，N可分解成",N,"=",P,""^u,"*",Q:N/(P^u)),
             i:d/*將i設為d，直接結束for迴圈*/
            )
         else
           (print("N無法被P"^u,"整除"))
        )
     ),
   print("－－－－－－－－－")
  );
第1列向量\(B[1]=[2268453123948987361000000000000000000000000,92119923815187304000000000000000000000000,\)
　\(1636648392655612000000000000000000000000,16615719722392000000000000000000000000,\)
　\(105429693670000000000000000000000000,428140888000000000000000000000000,\)
　\(1086652000000000000000000000000,1576000000000000000000000000,1000000000000000000000000]\)
產生不需要同餘\(P^um=p^8\)的方程式\(h(x)\)
\(h(x)=2268453123948987361000000000000000000000000\)
\(\displaystyle+92119923815187304000000000000000000000000\left(\frac{x}{1000}\right)\)
\(\displaystyle+1636648392655612000000000000000000000000\left(\frac{x}{1000}\right)^2\)
\(\displaystyle+16615719722392000000000000000000000000\left(\frac{x}{1000}\right)^3\)
\(\displaystyle+105429693670000000000000000000000000\left(\frac{x}{1000}\right)^4\)
\(\displaystyle+428140888000000000000000000000000\left(\frac{x}{1000}\right)^5\)
\(\displaystyle+1086652000000000000000000000000\left(\frac{x}{1000}\right)^6\)
\(\displaystyle+1576000000000000000000000000\left(\frac{x}{1000}\right)^7\)
\(\displaystyle+1000000000000000000000000\left(\frac{x}{1000}\right)^8\)
\(h(x)=2268453123948987361000000000000000000000000\)
　\(+92119923815187304000000000000000000000x\)
　\(+1636648392655612000000000000000000x^2\)
　\(+16615719722392000000000000000x^3\)
　\(+105429693670000000000000x^4\)
　\(+428140888000000000x^5\)
　\(+1086652000000x^6+1576000x^7+x^8\)
\(=(197000+x)^8\)
正整數解為\([]\)
－－－－－－－－－
第2列向量\(B[2]=[520814821077069868175677184373000000000000,-2430603060877548431773082696764000000000000,\)
　\(-68074195511300390211470172818000000000000,-817458291481540650461827996000000000000,\)
　\(-5507602630239396764545467000000000000,-22691467064000000000000000000000000,\)
　\(-57592556000000000000000000000000,-83528000000000000000000000000,-53000000000000000000000000]\)
產生不需要同餘\(P^um=p^8\)的方程式\(h(x)\)
\(h(x)=520814821077069868175677184373000000000000\)
\(\displaystyle-2430603060877548431773082696764000000000000\left(\frac{x}{1000}\right)\)
\(\displaystyle-68074195511300390211470172818000000000000\left(\frac{x}{1000}\right)^2\)
\(\displaystyle-817458291481540650461827996000000000000\left(\frac{x}{1000}\right)^3\)
\(\displaystyle-5507602630239396764545467000000000000\left(\frac{x}{1000}\right)^4\)
\(\displaystyle-22691467064000000000000000000000000\left(\frac{x}{1000}\right)^5\)
\(\displaystyle-57592556000000000000000000000000\left(\frac{x}{1000}\right)^6\)
\(\displaystyle-83528000000000000000000000000\left(\frac{x}{1000}\right)^7\)
\(\displaystyle-53000000000000000000000000\left(\frac{x}{1000}\right)^8\)
\(h(x)=520814821077069868175677184373000000000000\)
　\(-2430603060877548431773082696764000000000x\)
　\(-68074195511300390211470172818000000x^2\)
　\(-817458291481540650461827996000x^3\)
　\(-5507602630239396764545467x^4\)
　\(-22691467064000000000x^5\)
　\(-57592556000000x^6\)
　\(-83528000x^7-53x^8\)
\(=-53(-213+x)(197000+x)^4(394213+x)(77701967369+394000x+x^2)\)
正整數解為\([x=213]\)
當正整數解\(x=213\)時，可能的因數\(P=V+x=197000+213=197213\)
\(N\)可被\(P^4\)整除，\(N\)可分解成\(80171134270603235454533=197213^4\cdot53\)
－－－－－－－－－
(%o30)　done

\(N=P^uQ\)的\(P,Q\)值
(%i32)
P;
Q;
(%o31)　\(197213\)
(%o32)　\(53\)

\(N=p^rq^s\)的質因數\(p,q\)
(%i34)
p: Q^(beta/gamma)*P^(-b/gamma);
q: Q^(-alpha/gamma)*P^(a/gamma);
(%o34)　\(61\)
(%o35)　\(53\)

問題敘述
令\(N=p^8q^3=28540809637096203437\)是未知分解的整數，\(r>s\)和\(gcd(r,s)=1\)，給定\(N\)為輸入，能在多項式時間\(\log N\)和\(r=\Omega(\log^3 max(p,q))\)情況下回復質因數
步驟1：將\(r\)和\(s\)表示成\(\cases{r=u\cdot \alpha+a\cr s=u\cdot \beta+b}\)
設矩陣\(M=\left[\matrix{\lfloor\;r^{1/3}\rfloor\;&-s\cr0&r}\right]\)，經LLL化簡後短的非零向量\(v=(\lfloor\;r^{1/3}\rfloor\;\cdot\alpha,\gamma)\)，其中\(\gamma=-s\cdot\alpha+r\cdot\beta\)，\(\beta\in\mathbb{Z}\) 從\(v\)向量得到\(\alpha,\beta\)值，符合\(\alpha>0\)且\(0\le\beta\le\alpha\) Case 1：\(\beta=0\)或(\(\beta\ne0\)和\(\displaystyle\lfloor\;\frac{r}{\alpha}\rfloor\;\le\frac{s}{\beta}\)) 　\(\displaystyle u=\lfloor\;\frac{r}{\alpha}\rfloor\;\) 　\(a=r-u\cdot\alpha\) 　\(b=s-u\cdot\beta\) Case 2：\(\beta\ne0\)和\(\displaystyle\lfloor\;\frac{r}{\alpha}\rfloor\;>\frac{s}{\beta}\) 　\(\displaystyle u=\lceil\;\frac{r}{\alpha}\rceil\;\) 　\(a=r-u\cdot\alpha\) 　\(b=s-u\cdot\beta\)	設矩陣\(M=\left[\matrix{\lfloor\;r^{1/3}\rfloor\;&-s\cr0&r}\right]=\left[\matrix{2&-3\cr0&8}\right]\)，經LLL化簡後\(B=\left[\matrix{2&-3\cr4&2}\right]\) 原本短的非零向量\(v=(2,-3)\)無法求得合適的\(u,\alpha,\beta,a,b\) 改用\(v=(4,2)=(\lfloor\;r^{1/3}\rfloor\;\cdot \alpha,\gamma)\)，其中\(\gamma=-s\cdot\alpha+r\cdot\beta\)，\(\beta\in\mathbb{Z}\) \(4=\lfloor\;8^{1/3}\rfloor\;\cdot \alpha\)，得\(\alpha=2\) \(2=-3\cdot 2+8\cdot\beta\)，得\(\beta=1\) 符合\(\alpha>0\)且\(0\le\beta\le\alpha\) \(\displaystyle \lfloor\;\frac{r}{\alpha}\rfloor\;=\lfloor\;\frac{8}{2}\rfloor\;=4\)，\(\displaystyle \frac{s}{\beta}=\frac{3}{1}=3\) 符合Case 2:\(\beta\ne0\)和\(\displaystyle \lfloor\;\frac{r}{\alpha}\rfloor\;>\frac{s}{\beta})\) 　\(\displaystyle u=\lceil\;\frac{r}{\alpha}\rceil\;=\lceil\;\frac{8}{2}\rceil\;=4\) 　\(a=r-u\cdot\alpha=8-4\cdot2=0\) 　\(b=s-u\cdot\beta=3-4\cdot1=-1\) \(\cases{8=4\cdot2+0\cr3=4\cdot1-1}\)
步驟2：根據\(a,b\le0\)(Case2)使用Coppersmith方法
\(a,b\le0\)(Case2) \(\displaystyle N=p^rq^s=p^{u\cdot\alpha+a}q^{u\cdot\beta+b}=\frac{(p^{\alpha}q^{\beta})^u}{p^{-a}q^{-b}}=\frac{P^u}{Q}\)，其中\(P=p^{\alpha}q^{\beta}\)和\(Q=p^{-a}q^{-b}\)	\(\displaystyle N=p^8q^3=\frac{P^4}{Q}\)，其中\(P=p^2q^1=p^2q\)和\(Q=p^0q^1=q^1\)
步驟2-1.設同餘方程式，計算參數\(X\)
設同餘方程式\(f(x)=(X\cdot t+x)^u \equiv0\pmod{N}\) 其中\(X=\lfloor\;N^{1/u}\rfloor\;\)和\(|\;x_0|\;<X\) 針對\(t\)進行窮舉，其中\(t\)的範圍為\(\displaystyle0\le t\le \frac{P}{X}\le\frac{2P}{N^{1/u}}=2Q^{1/u}\)	\(X=\lfloor\;N^{1/4}\rfloor\;=73091\)，\(0\le t\le 2Q^{1/4}=5\) \(P=X\cdot t+x_0\)，\(x_0=P-X\cdot t\) \(x_0=197213-73091\cdot0=197213\) \(x_0=197213-73091\cdot1=124122\) \(x_0=197213-73091\cdot2=51031\) \(x_0=197213-73091\cdot3=-22060\)，當\(t=3\)時才有較小的解\(x_0\) \(x_0=197213-73091\cdot4=-95151\) \(x_0=197213-73091\cdot5=-168242\) 但這樣\(X\)所需的矩陣\(M\)要非常大才能找到較小的解\(x_0\) 設同餘方程式\(f(x)=(X\cdot t+x)^4\pmod{N}\) 取\(V=X\cdot t=197000\)，Coppersmith方法計算上限\(X=316\) 只要維度8的矩陣\(M\)經LLL化簡就能找到較小的解\(x_0=213\)，其中\(|\;x_0|\;<316\)
步驟2-2.產生矩陣\(M\)
定義三角\((hk)\times(hk)\)矩陣\(M=(m_{i,j})\)，\(m_{i,j}=e_{i,j}X^{j-1}\) \(e_{i,j}\)是\(q_{u,v}(x)=N^{(h-1-v)}x^u(p(x))^v\)的\(\displaystyle x^{j-1}\)項係數 其中\(\displaystyle v=\lfloor\;\frac{i-1}{k}\rfloor\;\)和\(u=(i-1)-kv\) 注意到對所有\(u,v\ge 0\)，\(q_{u,v}(x_0)\equiv 0\pmod{N^{h-1}}\)	\(u=0,v=0,q_{0,0}(xX)=N=\bbox[border:1px solid black]{N}\) \(u=1,v=0,q_{1,0}(xX)=NXx=\bbox[border:1px solid black]{NX}x\) \(u=2,v=0,q_{2,0}(xX)=N(Xx)^2=\bbox[border:1px solid black]{NX^2}x^2\) \(u=3,v=0,q_{3,0}(xX)=N(Xx)^3=\bbox[border:1px solid black]{NX^3}x^3\) \(u=0,v=1,q_{0,1}(xX)=(V+Xx)^4\) \(=V^4+4V^3Xx+6V^2X^2x^2+4VX^3x^3+\bbox[border:1px solid black]{X^4}x^4\) \(u=1,v=1,q_{1,1}(xX)=(Xx)(V+Xx)^4\) \(=V^4Xx+4V^3X^2x^2+6V^2X^3x^3+4VX^4x^4+\bbox[border:1px solid black]{X^5}x^5\) \(u=2,v=1,q_{2,1}(xX)=(Xx)^2(V+Xx)^4\) \(=V^4X^2x^2+4V^3X^3x^3+6V^2X^4x^4+4VX^5x^5+\bbox[border:1px solid black]{X^6}x^6\) \(u=3,v=1,q_{3,1}(xX)=(Xx)^3(V+Xx)^4\) \(=V^4X^3x^3+4V^3X^4x^4+6V^2X^5x^5+4VX^6x^6+\bbox[border:1px solid black]{X^7}x^7\)
\(M=\matrix{&\matrix{1&　x&　x^2&x^3　&x^4&x^5&　x^6&x^7}\cr \matrix{q_{0,0}(xX)\cr q_{1,0}(xX)\cr q_{2,0}(xX)\cr q_{3,0}(xX)\cr q_{0,1}(xX)\cr q_{1,1}(xX)\cr q_{2,1}(xX)\cr q_{3,1}(xX)}&\left[\matrix{N&&&&&&&\cr 0&NX&&&&&&\cr 0&0&NX^2&&&&&\cr 0&0&0&NX^3&&&&\cr *&*&*&*&X^4&&&\cr 0&*&*&*&*&X^5&&\cr 0&0&*&*&*&*&X^6&\cr 0&0&0&*&*&*&*&X^7}\right]}\) *代表非零數字
步驟2-3.經LLL化簡後的短向量產生不需要同餘\(N^{h-1}\)的方程式，得到公鑰\(N\)的因數\(P\)和\(Q\)
矩陣\(M\)經LLL化簡為\(B\) \(B=LLL(M)\) lattice經LLL化簡後第一列\(b_1\)為整個lattice中較短向量 所形成的方程式不需要再同餘\(N^{h-1}\)	lattice經LLL化簡後第一列\(b_1\)為整個lattice中較短向量 \(b_1=[-214382400364144710,\) 　\(262262951620456068,\) 　\(970699742880291792,\) 　\(-1149094400899629312,\) 　\(-786698634624195328,\) 　\(278086338099347456,\) 　\(557584281975848960,\) 　\(314636844829229056]\) 產生不需要同餘\(N^{h-1}=N\)的方程式 \(r(x)=-214382400364144710\) 　\(\displaystyle+262262951620456068\left(\frac{x}{316}\right)\) 　\(\displaystyle+970699742880291792\left(\frac{x}{316}\right)^2\) 　\(\displaystyle-1149094400899629312\left(\frac{x}{316}\right)^3\) 　\(\displaystyle-786698634624195328\left(\frac{x}{316}\right)^4\) 　\(\displaystyle+278086338099347456\left(\frac{x}{316}\right)^5\) 　\(\displaystyle+557584281975848960\left(\frac{x}{316}\right)^6\) 　\(\displaystyle+314636844829229056\left(\frac{x}{316}\right)^7\) \(r(x)=(x-213)^3(x^4+1199x^3+718310x^2+226574471x+22184534430)\) 得到整數解\(x=213\) 得到因數\(P=V+x=197000+213=197213\) \(N\)可分解成\(\displaystyle28540809637096203437=\frac{197213^4}{53}\) \(P=197213,Q=53\)
步驟3.從\(P,Q\)求出\(N\)的質因數\(p,q\)
已知\(P=p^{\alpha}q^{\beta}\)，\(Q=p^{-a}q^{-b}\)，\(\gamma=-a\beta+b\alpha\)，則\(\cases{p=Q^{\displaystyle\frac{\beta}{\gamma}} \cdot P^{\displaystyle\frac{b}{\gamma}}\cr q=Q^{\displaystyle\frac{-\alpha}{\gamma}} \cdot P^{\displaystyle\frac{-a}{\gamma}}}\)	\(\alpha=2,\beta=1,a=0,b=-1,\gamma=-a\beta+b\alpha=-2\) \(\cases{p=Q^{\displaystyle\frac{\beta}{\gamma}} \cdot P^{\displaystyle\frac{b}{\gamma}}=53^{\displaystyle\frac{1}{-2}} \cdot 197213^{\displaystyle\frac{-1}{-2}}=61\cr q=Q^{\displaystyle\frac{-\alpha}{\gamma}} \cdot P^{\displaystyle\frac{a}{\gamma}}=53^{\displaystyle\frac{-2}{-2}} \cdot 197213^{\displaystyle\frac{0}{-2}}=53}\)

請下載LLL.zip，解壓縮後將LLL.mac放到C:\maxima-5.49.0\share\maxima\5.49.0\share目錄下
要先載入LLL.mac才能使用Coppersmith_Howgrave指令
(%i1)　load("LLL.mac");
(%o1)　C:/maxima-5.49.0/share/maxima/5.49.0/share/LLL.mac

要因數分解的公鑰\(N\)
(%i2)　N:28540809637096203437;
(%o2)　\(28540809637096203437\)

\(N=p^rq^s\)，因數\(p\)的次方\(r\)，因數\(q\)的次方\(s\)
(%i4)
r:8;
s:3;
(%o3)　\(8\)
(%o4)　\(3\)

將\(r,s\)表示成\(\cases{r=u\alpha+a\cr s=u\beta+b}\)
(%i10)
u:4;
alpha:2;
beta:1;
a:0;
b:-1;
gamma:-a*beta+b*alpha;
(%o5)　\(4\)
(%o6)　\(2\)
(%o7)　\(1\)
(%o8)　\(0\)
(%o9)　\(-1\)
(%o10)　\(-2\)

\(V\)的高位元和\(P\)相同，當作\(P\)的近似值
(%i11)　V:197000;
(%o11)　\(197000\)

同餘方程式\(f(x)=(V+x)^u\equiv0\pmod{N}\)
(%i12)　fx: ('V+x)^u;
(%o12)　\((x+V)^4\)

Coppersmith_Howgrave方法副程式,參數\(h=2\)
(%i13)　x:Coppersmith_Howgrave(fx,N,2);
參數\(h=2\)
\(p(x)\)最高次方\(k=4\)
\(\displaystyle X=ceiling(\frac{1}{\sqrt{2}}(hk)^{-1/(hk-1)}N^{(h-1)/(hk-1)})=ceiling(\frac{1}{\sqrt{2}}8^{-1/7}619081497^{1/7})=316\)
\(q_{uv}=N^{h-1-v}x^up(x)^v=[28540809637096203437,28540809637096203437x,28540809637096203437x^2,28540809637096203437x^3,\)
\(x^4+4Vx^3+6V^2x^2+4V^3x+V^4,x(x^4+4Vx^3+6V^2x^2+4V^3x+V^4),x^2(x^4+4Vx^3+6V^2x^2+4V^3x+V^4),x^3(x^4+4Vx^3+6V^2x^2+4V^3x+V^4)]\)
用\(316x\)取代\(x\),得到\(q_{uv}=[28540809637096203437,9018895845322400286092x,2849971087121878490405072x^2,900590863530513602968002752x^3,\)
\(9971220736x^4+24864942848000x^3+23251869024000000x^2+9663751472000000000x+1506138481000000000000,\)
\(316x(9971220736x^4+24864942848000x^3+23251869024000000x^2+9663751472000000000x+1506138481000000000000),\)
\(99856x^2(9971220736x^4+24864942848000x^3+23251869024000000x^2+9663751472000000000x+1506138481000000000000),\)
\(31554496x^3(9971220736x^4+24864942848000x^3+23251869024000000x^2+9663751472000000000x+1506138481000000000000)]\)
產生矩陣\(M=\left[\matrix{28540809637096203437&0&0&0&0&0&0&0\cr
0&9018895845322400286092&0&0&0&0&0&0\cr
0&0&2849971087121878490405072&0&0&0&0&0\cr
0&0&0&900590863530513602968002752&0&0&0&0\cr
1506138481000000000000&9663751472000000000&23251869024000000&24864942848000&9971220736&0&0&0\cr
0&475939759996000000000000&3053745465152000000000&7347590611584000000&7857321939968000&3150905752576&0&0\cr
0&0&150396964158736000000000000&964983566988032000000000&2321838633260544000000&2482913733029888000&995686217814016&0\cr
0&0&0&47525440674160576000000000000&304934807168218112000000000&733701008110331904000000&784600739637444608000&314636844829229056}\right]\)
LLL化簡\(B=\left[\matrix{-214382400364144710&262262951620456068&970699742880291792&-1149094400899629312&-786698634624195328&278086338099347456&557584281975848960&314636844829229056\cr
-95647068587766726&182979578697288624&-11740466593266768&93960183407089280&326727727726846464&-1172332296114931712&-373382331680256000&1258547379316916224\cr
563806567324327383&-1713882702679529004&741868428215794752&1343812760321008064&-686844257598941440&533048178880539648&-742781918489255936&-314636844829229056\cr
339672331849247613&-1051715069006529168&836490663580953744&-511700036361113088&1179476536175185408&-648320914932780032&557584281975848960&-943910534487687168\cr
380209352163547656&-229641609525161168&809752168407502960&-1587976988878818624&93307054123140096&254958689875439616&-1300366200465104896&-629273689658458112\cr
12884067955383654&1107969310258392544&-1103177827036703360&630163920939967104&-1502628534289998592&-1230204982072495104&742781918489255936&-629273689658458112\cr
-2380570610455190403&901304238417789672&2592981941171391968&1408390870111434112&460521577533715200&777917668536230912&995686217814016&-629273689658458112\cr
16254669409485033140&9971967498989572660&6766766811505476528&3583338340143789952&2952987490748172800&3784414259565920256&1674744218363174912&1573184224146145280}\right]\)
產生不需要同餘\(N\)的方程式
\(r(x)=-214382400364144710\)
　\(+262262951620456068\left(\frac{x}{316}\right)\)
　\(+970699742880291792\left(\frac{x}{316}\right)^2\)
　\(-1149094400899629312\left(\frac{x}{316}\right)^3\)
　\(-786698634624195328\left(\frac{x}{316}\right)^4\)
　\(+278086338099347456\left(\frac{x}{316}\right)^5\)
　\(+557584281975848960\left(\frac{x}{316}\right)^6\)
　\(+314636844829229056\left(\frac{x}{316}\right)^7\)
\(r(x)=x^7+560x^6+88256x^5-78896923x^4-36416186172x^3+9720995662557x^2+829946049431823x-214382400364144710\)
\(=(x-213)^3(x^4+1199x^3+718310x^2+226574471x+22184534430)\)
整數解為\([x=213]\)
(%o13)　\([x=213]\)

\(\displaystyle N=\frac{P^u}{Q}\)的\(P\)值
(%i14)　P:V+rhs(x[1]);
(%o14)　\(197213\)

\(\displaystyle N=\frac{P^u}{Q}\)的\(Q\)值
(%i15)　Q: (P^u)/N;
(%o15)　\(53\)

\(N=p^rq^s\)的質因數\(p,q\)
(%i17)
p: Q^(beta/gamma)*P^(b/gamma);
q: Q^(-alpha/gamma)*P^(-a/gamma);
(%o16)　\(61\)
(%o17)　\(53\)