UBASICからの移植願い

 投稿者：しばっち  投稿日：2016年 3月29日(火)20時53分23秒
  下記プログラム(UBASIC)の移植をどなたかお願いします

   10   'ecm
   20   '
   30   '  Prime Decomposition by ECM
   40   '        1987,88,89 by Yuji KIDA
   50   '
   60   word -120:point -2
   70   dim SS(20),PD(20):'中間約数用&素因数用
   80   print "素因数分解（楕円曲線法）"
   90   '
  100   '  各数ごとの分解
  110   '
  120   print:input "整数=";N
  130   if N=0 then end
  140   gosub *Factoring_ECM(N)
  150   goto 120
  160   '
  170   '  メイン・ルーチン
  180   '
  190   *Factoring_ECM(N)
  200   NN=N:S1%=1:SP%=0:DP%=0
  210   '
  220   '  小さい素数で割る
  230   '
  240   D=prmdiv(N):if D=0 goto 330
  250   if N=D goto 280
  260   print D;"*";:N=N\D:S1%=0
  270   goto 240
  280   if S1% then print N;"は素数です":goto 820
  290   print N:goto 820
  300   '
  310   '  大きい素因数を持つ場合
  320   '
  330   print
  340   '
  350   '  ADLEMAN TEST
  360   '
  370   if fnAdleman(N)=0 then clr S1%:goto 470
  380   '
  390   '  素数ならば PD に登録
  400   '
  410   PD(DP%)=N:inc DP%
  420   if SP% then dec SP%:N=SS(SP%):goto 370
  430   goto 600
  440   '
  450   '  ある数のべきであるかチェックする
  460   '
  470   N1=fnPowerCheck(N):N2=N\N1
  480   if N1=1 then 540
  490   if N1>N2 then swap N1,N2
  500   SS(SP%)=N2:inc SP%:N=N1:goto 370
  510   '
  520   '  EC 法で約数を見つける
  530   '
  540   N1=fnECM(N):N2=N\N1
  550   if N1>N2 then swap N1,N2
  560   SS(SP%)=N2:inc SP%:N=N1:goto 370
  570   '
  580   '  結果の表示
  590   '
  600   if S1% then print N;"は素数です":goto 820
  610   '
  620   '  新めて小さい方から表示
  630   '
  640   N=NN:print:print N;"="
  650   D=prmdiv(N):if D=0 goto 680
  660   print D;"*";:N=N\D
  670   goto 650
  680   if DP%=1 then print PD(0):goto 820
  690   '
  700   '  大きい方はSORT して表示
  710   '
  720   repeat
  730      clr SW%
  740      for I=0 to DP%-2
  750         if PD(I)>PD(I+1) then swap PD(I),PD(I+1):SW%=1
  760      next
  770   until SW%=0
  780   for I=0 to DP%-2
  790      print PD(I);"*";
  800   next
  810   print PD(DP%-1)
  820   return
  830   '
  840   '
  850   fnADLEMAN(N)
  860   local E%,ET%,G%,H%,I%,J%,MW%,NP%,NQ%,OK%,P%,Q%,SW%,W%,WK%,X%
  870   local E,IN,M1,P1,S,T,U,U1,V,V1,W,U#,W#
  880   dim Index%(9282):'index mod Q with base G
  890   dim Phi%(9282),Psi%(9282)
  900   dim Tau(16),Tau1(16),Tau2(16),Tau_e(16),Sum(16),Tau_p(16):'16=max of P-1
  910   dim P%(6),Q%(27),G%(27):'6=no of P -1,27=no of Q -1
  920   dim WK%(6,27),Eta%(6,27)
  930   '
  940   print "ADLEMAN 素数判定 for";N
  950   restore *PrimeData1
  960   for I%=0 to 6:read P%(I%):next
  970   restore *PrimeData2
  980   for I%=0 to 27:read Q%(I%),G%(I%):next
  990   '
1000   *PrimeData1
1010   data 2,3,5,7,11,13,17
1020   '
1030   *PrimeData2
1040   data 3,2,7,3,11,2,31,3
1050   data 43,3,71,7,211,2
1060   data 23,5,67,2,331,3,463,3,2311,3
1070   data 79,3,131,2,547,2,859,2,911,17,2003,5,2731,3,6007,3
1080   data 103,5,239,7,443,2,1123,2,1327,3,1871,14,3571,2,9283,2
1090   '
1100   if N=0 then return(N)
1110   '
1120   S=1
1130   for NQ%=0 to 27
1140      S=S*Q%(NQ%):if S^2>N then cancel for:jump
1150   next
1160   print:print "数値が大き過ぎます":stop:end
1170   **
1180   NP%=2+(NQ%>3)+(NQ%>6)+(NQ%>11)+(NQ%>19)
1190   '
1200   W=prmdiv(N)
1210   if W goto *DIVISOR?
1220   '
1230   '☆☆ TEST FOR P=2
1240   '
1250   gosub *Test4_2
1260   if H%=1 then SW%=1 else gosub *Test3_2:SW%=OK%
1270   MW%=1
1280   for J%=0 to NQ%
1290     Q%=Q%(J%):G%=G%(J%)
1300     gosub *Test2_2
1310     if not OK% then cancel for:goto *FAILED
1320     if SW% then WK%=H%
1330     :else gosub *Test5_2
1340       :if not OK% then cancel for:goto *DIVISOR
1350     WK%(0,J%)=WK%:Eta%(0,J%)=ET%:if WK%>MW% then MW%=WK%
1360   next
1370   '☆ W の最大値と η(χ) の計算
1380   if MW%<>1 then
1390     :for J%=0 to NQ%
1400       :if MW%>WK%(0,J%) then Eta%(0,J%)=0 endif
1410     :next
1420   '
1430   '☆☆ TEST for P>2
1440   '
1450   for I%=1 to NP%
1460     P%=P%(I%)
1470     gosub *Test4_P
1480     if H%=1 then SW%=1 else gosub *Test3_P:SW%=OK%
1490     '
1500     MW%=1
1510     for J%=0 to NQ%
1520       Q%=Q%(J%):G%=G%(J%)
1530       if (Q%-1)@P% then clr WK%:clr ET%:jump
1540       gosub *Test2_P
1550       if not OK% then cancel for,for:goto *FAILED
1560       if SW% then WK%=H%:ET%=Index%(N@Q%)@P%:jump
1570       gosub *Test5_P
1580       if not OK% then cancel for,for:goto *FAILED
1590       **
1600       WK%(I%,J%)=WK%:Eta%(I%,J%)=ET%
1610       if WK%>MW% then MW%=WK%
1620     next J%
1630      '☆ W の最大値と η(χ) の計算
1640     if MW%<>1 then
1650       :for J%=0 to NQ%
1660          :WK%=WK%(I%,J%)
1670          :if WK% then if MW%>WK% then Eta%(I%,J%)=0 endif endif
1680       :next J%
1690   next I%
1700   '
1710   '☆☆ 最終テスト
1720   '
1730   gosub *GET_V
1740   '
1750   W=1
1760   T=1:for I%=0 to NP%:T=T*P%(I%):next
1770   while T
1780     W=W*V@S:if W=1 goto *ALLOK
1790     if N@W=0 goto *GETDIVISOR
1800     dec T
1810   wend
1820   '
1830   *ALLOK
1840   return(1)
1850   '
1860   *FAILED
1870   return(0)
1880   '
1890   *DIVISOR?
1900   if W=N goto *ALLOK
1910   *GETDIVISOR
1920   return(0)
1930   '
1940   '☆☆ サブルーチン for P=2
1950   '
1960   '☆ 条件 3 のチェック for P=2
1970   '
1980   *Test3_2
1990   local J%
2000   for J%=0 to NQ%
2010     Q%=Q%(J%)
2020     if kro(N,Q%)<>1 then cancel for:OK%=1:jump
2030   next:OK%=0
2040   **
2050   return
2060   '
2070   '☆ 条件 4 のチェック for P=2
2080   '
2090   *Test4_2
2100   W#=N-1:H%=-1
2110   repeat U#=W#:W#=W#\2:inc H% until res
2120   return
2130   '
2140   '☆ 条件 2 のチェック for P=2
2150   '
2160   *Test2_2
2170   W=Q%:if (Q%-1)@4 then W=N-W
2180   U=modpow(W,(N-1)\2,N)
2190   if kro(N,Q%)<0 then ET%=1:OK%=(U=N-1)
2200   :else ET%=0:OK%=(U=1)
2210   return
2220   '
2230   '☆条件 5 のチェック for P=2
2240   '
2250   *Test5_2
2260   T=kro(-1,Q%)*Q%:'τ^2
2270   T=modpow(T,U#,N):WK%=1:OK%=1:ET%=0
2280   if T=1 then jump
2290   if T=N-1 then ET%=1:jump
2300   repeat
2310     inc WK%:W=T:T=T^2@N
2320     if T=N-1 then ET%=1:jump
2330   until T=1
2340   W=gcd(W-1,N):OK%=(W=1)
2350   **
2360   return
2370   '
2380   '☆☆ サブルーチン for P>2
2390   '
2400   '☆ 条件 3 のチェック for P>2
2410   '
2420   *Test3_P
2430   local J%
2440   for J%=0 to NQ%
2450     Q%=Q%(J%):E%=(Q%-1)\P%
2460     if not res then
2470       :if modpow(N,E%,Q%)<>1 then cancel for:OK%=1:jump
2480   next:OK%=0
2490   **
2500   return
2510   '
2520   '☆ 条件 4 のチェック for P>2 (N^(P%-1)-1 を P%^H%*U と分解)
2530   '
2540   *Test4_P
2550   W#=N^(P%-1)-1:H%=-1
2560   repeat U#=W#:W#=W#\P%:inc H% until res
2570   return
2580   '
2590   '☆ 条件 2 のチェック for P>2
2600   '
2610   ' τ(χ)^N/τ(χ^N) の計算
2620   *Test2_P
2630   local I%
2640   W%=1:for I%=1 to Q%-1:W%=W%*G%@Q%:Index%(W%)=I%@P%:next:'ind の計算
2650   gosub *TAU_P
2660   if not SW% then block Tau_p(0..P%-2)=block Tau(0..P%-2)
2670   '
2680   E=N\P%:gosub *T_E
2690   block Tau_e(0..P%-2)=block Tau(0..P%-2)
2700   E%=N@P%:gosub *JACOBI
2710   gosub *T_TM
2720   'MAIN OF TEST2
2730   E%=Index%(N@Q%)*(-N)@P%
2740   if E%=P%-1 then W=-(N-1)*(P%-1)
2750   :else if Tau(E%)<>1 then clr OK%:jump endif
2760      :dec Tau(E%):clr W
2770   for I%=0 to P%-2:W=W+Tau(I%):next
2780   OK%=(W=0)
2790   **
2800   return
2810   '
2820   '☆条件 5 のチェック for P>2
2830   '
2840   ' (τ(χ)^P)^U の計算
2850   *Test5_P
2860   local I%,J%
2870   block Tau1(0..P%-2)=block Tau_p(0..P%-2)
2880   Tau(0)=1:block Tau(1..P%-2)=0
2890   for I%=0 to len(U#)-2
2900     if bit(I%,U#) gosub *T_T1
2910     gosub *T1_T1
2920   next
2930   gosub *T_T1
2940   ' (τ(χ)^P)^U を P 乗しながら条件 5 をチェック
2950   WK%=1:gosub *W_KAI:if OK% then jump
2960   repeat
2970     block Tau1(0..P%-2)=block Tau(0..P%-2)
2980     block Tau_e(0..P%-2)=block Tau(0..P%-2)
2990     for I%=1 to len(P%)-2
3000       gosub *T1_T1
3010       if bit(I%,P%) gosub *T_T1
3020     next
3030     gosub *T1_T1
3040     gosub *T_T1
3050     inc WK%:gosub *W_KAI
3060   until OK%
3070   if ET%>0 then OK%=1:jump
3080   ' MAIN OF TEST 5
3090   I%=0
3100   loop
3110     T=Tau_e(I%)
3120     if gcd(T,N)=1 goto 3170
3130     if T then W=gcd(T,N):clr OK%:jump
3140     inc I%
3150   endloop
3160   '
3170   if T<>1 then
3180     :if gcd(T-1,N)=1 goto 3270 endif
3190     :W=gcd(T-1,N):clr OK%:jump
3200   :else J%=I%+1
3210      :loop
3220         :T=Tau_e(J%):if gcd(T,N)=1 goto 3270 endif
3230         :if T then W=gcd(T,N):clr OK%:jump endif
3240         :inc J%
3250      :endloop
3260   '
3270   I%=0
3280   loop
3290     T=Tau_e(I%)
3300     if gcd(T+1,N)=1 then OK%=1:jump
3310     if T<N-1 then W=gcd(T+1,N):clr OK%:jump
3320     inc I%
3330   endloop
3340   **
3350   return
3360   '
3370   '* W(χ) を決定するためのチェック
3380   '
3390   *W_KAI
3400     local I%
3410     T=Tau(0)
3420     if T<=1 then
3430       :W=T:clr ET%
3440       :for I%=1 to P%-2
3450         :T=Tau(I%):if T then W=W+T:ET%=I% endif
3460       :next
3470       :OK%=(W=1)
3480     :elseif T<N-1 then
3490       :clr OK%
3500     :else for I%=1 to P%-2
3510             :if Tau(I%)<T then cancel for:clr OK%:jump endif
3520          :next
3530          :OK%=1:ET%=P%-1
3540     :endif
3550     **
3560   return
3570   '
3580   '* τ(χ)^P の計算
3590   '
3600   *TAU_P
3610     local I%,J%
3620     block Psi%(1..Q%-1)=block Index%(1..Q%-1)
3630     block Phi%(1..Q%-1)=block Index%(1..Q%-1)
3640     Tau(0)=1:block Tau(1..P%-2)=0:'SET 1
3650     Tau1(0)=1:block Tau1(1..P%-2)=0:'SET 1
3660     for I%=1 to len(P%)-2
3670       gosub *PSI_PSI
3680       if bit(I%,P%) gosub *PHY_PSI
3690     next
3700     gosub *PSI_PSI
3710     gosub *T_T1
3720     block Tau1(0..P%-2)=block Tau(0..P%-2)
3730     clr Tau1(P%-1)
3740     J%=Phi%(Q%-1)
3750     for I%=0 to P%-1:W%=(I%+J%)@P%:Tau(W%)=Tau1(I%):next
3760     W=Tau(P%-1)
3770     for I%=0 to P%-2:Tau(I%)=(Tau(I%)-W)*Q%@N:next
3780   return
3790   '
3800   '* T^E の計算
3810   '
3820   *T_E
3830     local I%,J%
3840     block Tau1(0..P%-2)=block Tau(0..P%-2)
3850     J%=0
3860     if even(E) then
3870       :while bit(J%,E)=0:gosub *T1_T1:inc J%:wend
3880       :block Tau(0..P%-2)=block Tau1(0..P%-2)
3890     for I%=J%+1 to len(E)-1
3900       gosub *T1_T1
3910       if bit(I%,E) gosub *T_T1
3920     next
3930   return
3940   '
3950   '* τ(χ)^E/τ(χ^E) の計算
3960   '
3970   *JACOBI
3980     local I%,J%
3990     block Psi%(0..Q%-1)=block Index%(0..Q%-1)
4000     Tau1(0)=1:block Tau1(1..P%-2)=0
4010     J%=0:while bit(J%,E%)=0:gosub *PSI_PSI:inc J%:wend
4020     block Tau(0..P%-2)=block Tau1(0..P%-2)
4030     block Phi%(0..Q%-1)=block Psi%(0..Q%-1)
4040     for I%=J%+1 to len(E%)-1
4050       gosub *PSI_PSI
4060       if bit(I%,E%) gosub *PHY_PSI
4070     next
4080   return
4090   '
4100   '* T=T*TM の計算
4110   '
4120   *T_TM
4130     local I%,J%
4140     for J%=0 to P%-2:Tau2(J%)=Tau(0)*Tau_e(J%)@N:next:clr Tau2(P%-1)
4150     for I%=1 to P%-2
4160       for J%=0 to P%-2
4170         W%=(I%+J%)@P%:Tau2(W%)=Tau2(W%)+Tau(I%)*Tau_e(J%)@N
4180       next
4190     next
4200     W=Tau2(P%-1)
4210     for I%=0 to P%-2:Tau(I%)=(Tau2(I%)-W)@N:next
4220     block Tau2(0..P%-2)=block Tau(0..P%-2)
4230   return
4240   '
4250   '* PSI(A)*PSI(1-A) の和の計算
4260   '
4270   *PSI_PSI
4280     local I%,J%
4290     block Sum(0..P%-1)=0
4300     for I%=2 to Q%-1:inc Sum((Psi%(I%)+Psi%(Q%+1-I%))@P%):next
4310     ' T1=T1*T1*U の計算
4320     for I%=0 to (P%-1)\2:Tau2(2*I%)=Tau1(I%)^2@N:next
4330     for I%=(P%+1)\2 to P%-2:Tau2(2*I%-P%)=Tau1(I%)^2@N:next:clr Tau2(P%-2)
4340     for I%=0 to P%-3
4350       for J%=I%+1 to P%-2
4360         W%=(I%+J%)@P%:Tau2(W%)=Tau2(W%)+2*Tau1(I%)*Tau1(J%)@N
4370       next
4380     next
4390     for J%=0 to P%-1:Tau1(J%)=Tau2(0)*Sum(J%)@N:next
4400     for I%=1 to P%-1
4410       for J%=0 to P%-1
4420         W%=(I%+J%)@P%:Tau1(W%)=Tau1(W%)+Tau2(I%)*Sum(J%)@N
4430       next
4440     next
4450     W=Tau1(P%-1)
4460     for I%=0 to P%-2:Tau1(I%)=(Tau1(I%)-W)@N:next
4470     ' PSI^2 の計算
4480     for I%=0 to Q%-1:Psi%(I%)=2*Psi%(I%)@P%:next
4490   return
4500   '
4510   '* PHY(A)*PSI(1-A) の和の計算
4520   '
4530   *PHY_PSI
4540     local I%,J%
4550     block Sum(0..P%-1)=0
4560     for I%=2 to Q%-1:inc Sum((Phi%(I%)+Psi%(Q%+1-I%))@P%):next
4570     ' T=T*T1*SU の計算
4580     for J%=0 to P%-2:Tau2(J%)=Tau(0)*Tau1(J%)@N:next:clr Tau2(P%-1)
4590     for I%=1 to P%-2
4600       for J%=0 to P%-2
4610         W%=(I%+J%)@P%:Tau2(W%)=Tau2(W%)+Tau(I%)*Tau1(J%)@N
4620       next
4630     next
4640     for J%=0 to P%-1:Tau(J%)=Tau2(0)*Sum(J%):next
4650     for I%=1 to P%-1
4660       for J%=0 to P%-1
4670         W%=(I%+J%)@P%:Tau(W%)=Tau(W%)+Tau2(I%)*Sum(J%)@N
4680       next
4690     next
4700     W=Tau(P%-1)
4710     for I%=0 to P%-2:Tau(I%)=(Tau(I%)-W)@N:next
4720     ' PHY=PHY*PSY の計算
4730     for I%=0 to Q%-1:Phi%(I%)=(Phi%(I%)+Psi%(I%))@P%:next
4740   return
4750   '
4760   '* T=T*T1 の計算
4770   '
4780   *T_T1
4790     local I%,J%
4800     for J%=0 to P%-2:Tau2(J%)=Tau(0)*Tau1(J%)@N:next:clr Tau2(P%-1)
4810     for I%=1 to P%-2
4820       for J%=0 to P%-2
4830         W%=(I%+J%)@P%:Tau2(W%)=Tau2(W%)+Tau(I%)*Tau1(J%)@N
4840       next
4850     next
4860     W=Tau2(P%-1)
4870     for I%=0 to P%-2:Tau(I%)=(Tau2(I%)-W)@N:next
4880   return
 		 

Re: UBASICからの移植願い

 投稿者：しばっち  投稿日：2016年 3月29日(火)20時54分12秒
  > No.4019[元記事へ]

続き

4890   '
4900   '* T1=T1*T1 の計算
4910   '
4920   *T1_T1
4930     local I%,J%
4940     for I%=0 to (P%-1)\2:Tau2(2*I%)=Tau1(I%)^2@N:next
4950     for I%=(P%+1)\2 to P%-2:Tau2(2*I%-P%)=Tau1(I%)^2@N:next:clr Tau2(P%-2)
4960     for I%=0 to P%-3
4970       for J%=I%+1 to P%-2
4980          W%=(I%+J%)@P%:Tau2(W%)=Tau2(W%)+2*Tau1(I%)*Tau1(J%)@N
4990       next
5000     next
5010     W=Tau2(P%-1)
5020     for I%=0 to P%-2:Tau1(I%)=(Tau2(I%)-W)@N:next
5030   return
5040   '
5050   '☆☆ サブルーチン for 最終テスト
5060   '
5070   ' ☆ V mod S の計算
5080   '
5090   *GET_V
5100     local I%,J%
5110     for J%=0 to NQ%
5120       IN=Eta%(0,J%):P1=2
5130       for I%=1 to NP%
5140         if WK%(I%,J%) then
5150           :P%=P%(I%):U1=P1*modinv(P1,P%):V1=1-U1
5160           :P1=P1*P%:IN=(IN*V1+Eta%(I%,J%)*U1)@P1
5170       next
5180       Q%=Q%(J%):X%=modpow(G%(J%),IN,Q%)
5190       if J% then
5200         :U1=M1*modinv(M1,Q%):V1=1-U1
5210         :M1=M1*Q%:V=(V*V1+X%*U1)@M1
5220       :else M1=Q%:V=X%
5230     next
5240   return
5250   '
5260   'POWERCHK V1.1
5270   '  N がある数の巾かチェック
5280   '  1987,89 by Yuji KIDA
5290   '
5300   fnPowerCheck(N)
5310      local E%,D#,EP#,W#
5320   print "巾乗判定 ";
5330   for E%=(len(N)-1)\17 to 2 step -1:print E%;
5340   '
5350   ' NEWTON法で近似値を求める
5360   '
5370   D#=exp(log(N)/E%):W#=D#^E%:'NEWTON法の初期値
5380   EP#=(D#+1)^E%-W#:'許容誤差
5390   '
5400   repeat
5410   D#=((E%-1)*W#+N)*D#/(E%*W#):W#=D#^E%:'NEWTON法メイン
5420   until abs(W#-N)<EP#
5430   '
5440   ' 整数値を一つずつ調べる
5450   '
5460   D#=int(D#):W#=D#^E%
5470   if W#>N goto 5550
5480   if W#=N then cancel for:goto 5600
5490   '小さい場合
5500   inc D#:W#=D#^E%
5510   if W#=N then cancel for:goto 5600
5520   if W#<N goto 5500
5530   goto 5580:'次の E% へ
5540   '大きい場合
5550   dec D#:W#=D#^E%
5560   if W#=N then cancel for:goto 5600
5570   if W#>N goto 5550
5580   next
5590   D#=1:'巾でない
5600   print
5610   return(D#)
10000   '
10010   'ELLIPTIC CURVE METHOD ver 4.2
10020   ' originally   by H.W.LENSTRA
10030   ' fast version by P.L.Montgomery
10040   ' fast version by H.Suyama
10050   ' inplemented  by Y.Kida
10060   '
10070   fnECM(N)
10080   local EC%,I%,IP%,J%,Q%
10090   local A,A0,AA,DX,DZ,GD,L1,L2,LG,LS,N1,N2,P
10100   local TX,TZ,UX,UZ,W1,W2,W3,W4,WX,WZ,X,Z
10110   dim X(47),V0(11),V1(11),V2(11),V3(11)
10120   '
10130   print "楕円曲線法 with パラメーター ";
10140   L1=int(log(N)^2.65/10):L1=min(130000,max(L1,500))
10150   L2=40*L1
10160   print "( B1=";L1;"B2=";L2;")"
10170   LG=log(L1):LS=isqrt(L1)
10180   EC%=0
10190   print "曲線";
10200   '
10210   *SET_CURVE
10220   '
10230   inc EC%:print EC%;
10240   A0=2*(EC%+1)*modinv(3*(EC%+1)^2-1,N)@N
10250   if A0*(A0^2-1)*(9*A0^2-1)@N=0 then 10230
10260   A=(-3*A0^4-6*A0^2+1)*modinv(4*A0^3,N)@N
10270   AA=(A+2)*modinv(4,N)@N
10280   X=(3*A0^2+1)@N:Z=4*A0@N
10290   '
10300   '* 1st step
10310   '
10320   'for power of 2
10330   '
10340   for I%=1 to len(L1)
10350     W1=(X+Z)^2@N:W2=(X-Z)^2@N
10360     X=W1*W2@N:Z=(W1-W2)*(W2+AA*(W1-W2)@N)@N
10370   next
10380   GD=gcd(Z,N):if GD>1 then *MITUKETA
10390   '
10400   'for powers of odd primes
10410   '
10420   IP%=2
10430   repeat
10440     P=prm(IP%):M#=P^int(LG/log(P)):inc IP%
10450     gosub *ECM_SUB:GD=gcd(Z,N):if GD>1 then *MITUKETA
10460   until P>LS
10470   '
10480   repeat
10490     M#=1:for IP%=IP% to IP%+9:M#=M#*prm(IP%):next
10500     gosub *ECM_SUB:GD=gcd(Z,N):if GD>1 then *MITUKETA
10510   until prm(IP%+9)>L1
10520   '
10530   '* 2nd STEP
10540   '
10550   UX=X:UZ=Z:'                                =Q
10560   X(0)=X*modinv(Z,N)@N:J%=0
10570   W1=(X+Z)^2@N:W2=(X-Z)^2@N
10580   TX=W1*W2@N:TZ=(W1-W2)*(W2+AA*(W1-W2)@N)@N:'=2Q
10590   W1=(X-Z)*(TX+TZ)@N:W2=(X+Z)*(TX-TZ)@N
10600   X=(W1+W2)^2@N*UZ@N:Z=(W1-W2)^2@N*UX@N:'    =3Q
10610   for I%=5 to 209 step 2
10620     WX=X:WZ=Z
10630     W1=(X-Z)*(TX+TZ)@N:W2=(X+Z)*(TX-TZ)@N
10640     X=(W1+W2)^2@N*UZ@N:Z=(W1-W2)^2@N*UX@N:' =7Q,9Q,...
10650     if I%=105 then DX=X:DZ=Z:'              =105Q
10660     if gcd(I%,210)=1 then inc J%:X(J%)=X*modinv(Z,N)@N
10670     UX=WX:UZ=WZ
10680   next
10690   '
10700   W1=(DX+DZ)^2@N:W2=(DX-DZ)^2@N
10710   X=W1*W2@N:Z=(W1-W2)*(W2+AA*(W1-W2)@N)@N:'  =210Q
10720   '
10730   for J%=0 to 11:I%=J%*4
10740     W1=X(I%)+X(I%+1):W2=X(I%)*X(I%+1)@N
10750     W3=X(I%+2)+X(I%+3):W4=X(I%+2)*X(I%+3)@N
10760     V3(J%)=W1+W3:V2(J%)=W1*W3@N+W2+W4
10770     V1(J%)=(W2*W3+W1*W4)@N:V0(J%)=W2*W4@N
10780   next
10790   '
10800   UX=X:UZ=Z:'                                =210Q
10810   W1=(X+Z)^2@N:W2=(X-Z)^2@N
10820   TX=W1*W2@N:TZ=(W1-W2)*(W2+AA*(W1-W2)@N)@N:'=420Q
10830   W1=(X-Z)*(TX+TZ)@N:W2=(X+Z)*(TX-TZ)@N
10840   X=(W1+W2)^2@N*UZ@N:Z=(W1-W2)^2@N*UX@N:'    =630Q
10850   '
10860   for Q%=2 to L1\420
10870     WX=X:WZ=Z
10880     W1=(X-Z)*(TX+TZ)@N:W2=(X+Z)*(TX-TZ)@N
10890     X=(W1+W2)^2@N*UZ@N:Z=(W1-W2)^2@N*UX@N
10900     UX=WX:UZ=WZ
10910   next
10920   '
10930   for Q%=Q% to L2\420
10940     W1=X*modinv(Z,N)@N:W2=W1^2@N:W3=W1*W2@N:W4=W2^2@N
10950     M#=(W4-V3(0)*W3+V2(0)*W2-V1(0)*W1+V0(0))@N
10960     for J%=1 to 11
10970       M#=((W4-V3(J%)*W3+V2(J%)*W2-V1(J%)*W1+V0(J%))@N)*M#@N
10980     next
10990     GD=gcd(M#,N):if GD>1 then cancel for:goto *MITUKETA
11000     WX=X:WZ=Z
11010     W1=(X-Z)*(TX+TZ)@N:W2=(X+Z)*(TX-TZ)@N
11020     X=(W1+W2)^2@N*UZ@N:Z=(W1-W2)^2@N*UX@N
11030     UX=WX:UZ=WZ
11040   next:Q%=0
11050   '
11060   goto *SET_CURVE:'次の curve へ
11070   '
11080   *MITUKETA
11090   if and{GD=N,len(N)<50} then *SET_CURVE:'次の curve へ
11100   print:return(GD)
11110   '
11120   '* main subroutine
11130   '
11140   *ECM_SUB
11150   TX=X:TZ=Z:UX=X:UZ=Z
11160   W3=TX+TZ:W4=TX-TZ
11170   for I%=1 to len(M#)-1
11180     W1=W3^2@N:W2=W4^2@N
11190     TX=W1*W2@N:TZ=(W1-W2)*(W2+AA*(W1-W2)@N)@N
11200     W3=TX+TZ:W4=TX-TZ
11210     if bit(I%,M#) then
11220       :W1=W4*(X+Z)@N:W2=W3*(X-Z)@N
11230       :X=(W1+W2)^2@N*UZ@N:Z=(W1-W2)^2@N*UX@N
11240     :else W1=W4*(UX+UZ)@N:W2=W3*(UX-UZ)@N
11250       :UX=(W1+W2)^2@N*Z@N:UZ=(W1-W2)^2@N*X@N
11260   next
11270   return
 		 

戻る