10 ' r_sievex.ub 20 ' y2=x3+ax+b , 4a3+27b2<>0 30 ' y2=x3+az4x+bz6, (X/Z2, Y/Z3), h=log max{|x|,Z2} 40 Va%=9*7*11*13:Vb%=128*5*17:Vc%=720 50 dim Ta%(Va%-1),Tb%(Vb%-1),Tc%(Vc%),Td%(Vc%),E%(1000),Te%(256,256) 60 print "making sieve tables" 70 for I%=0 to Va%-1:Ta%(I%)=1:next I% 80 for I%=0 to Va%\2+1:Ta%(modpow(I%,2,Va%))=0:next I% 90 for I%=0 to Vb%-1:Tb%(I%)=1:next I% 100 for I%=0 to Vb%\2+1:Tb%(modpow(I%,2,Vb%))=0:next I% 110 for I%=0 to Vc%-1:Tc%(I%)=1:next I% 120 for I%=0 to Vc%\2:Tc%(modpow(I%,2,Vc%))=0:next I% 130 for I%=2 to 256 140 for J%=0 to I%-1:Te%(I%,J%)=gcd(I%,J%):next J% 150 next I% 160 ' 170 input "search level (1..4 ; 0:help) = ";L% 180 if L%<>0 then 240 190 print "level 1 : max(x)=10000, max(z)=100, H=";log(10000) 200 print "level 2 : max(x)=22026, max(z)=148, H=10" 210 print "level 3 : max(x)=65536, max(z)=256, H=";log(65536) 220 print "level 4 : max(x)=1000000, max(z)=1000, H=";log(1000000) 230 goto 170 240 if L%=1 then Ex=10000:Ez%=100:H=log(Ex) 250 if L%=2 then Ex=22026:Ez%=148:H=10 260 if L%=3 then Ex=65536:Ez%=256:H=log(Ex) 270 if L%=4 then Ex=1000000:Ez%=1000:H=log(Ex) 280 if or{L%<0,L%>4} then 230 290 print "search range :" 300 print "max(x) =";Ex;", max(z) =";Ez%;", naive height =";H 310 ' 320 print "y2 = x3 + ax + b : input a, b" 330 input "a=";A:input "b=";B 340 D=4*A^3+27*B^2:if D=0 then print "not irreducible":goto 320 350 print "a=";A;", b=";B;", discriminant=";D 360 clr time 370 ' case : z=1 (integer) 380 print "integer point : x <= 1000000" 390 X=-Ex:H%=0 400 U=X^3+A*X+B:V=3*X^2+A:inc H%:if H%>30 then Sx=-Ex:goto 440 410 if V=0 then Sx=1:goto 440 else W=U/V 420 if abs(W)>0.5 then X=X-W:goto 400 430 Sx=int(X):if Sx<-Ex then Sx=-Ex 440 J%=0 450 for I%=0 to Vc%-1:W%=(I%^3+A*I%+B)@Vc% 460 if Tc%(W%)=0 then inc J%:Td%(J%)=I% 470 next I% 480 ' 490 Base=Vc%*(Sx\Vc%):if Base<0 then Base=Base-Vc% 500 for I%=1 to J%:X=Base+Td%(I%) 510 Y2=X^3+A*X+B:if Y2<0 then 560 520 if Tb%(Y2@Vb%)=1 then 560 530 if Ta%(Y2@Va%)=1 then 560 540 Y=isqrt(Y2):if res>0 then 560 550 inc F%:R%=fnSub(X,Y) 560 next I% 570 Base=Base+Vc%:if Base30 then X=-Ex:goto 650 630 if V=0 then Sx=1:goto 660 else W=U/V 640 if abs(W)>0.5 then X=X-W:goto 620 650 Sx=int(X):if Sx<-Ex then Sx=-Ex 660 J%=0 670 for I%=0 to Z%-1:E%(I%)=0:next I% 680 for I%=0 to Z%\2:E%(modpow(I%,2,Z%))=1:next I% 690 for I%=0 to Z%-1:if Te%(Z%,I%)>1 then 720 700 W%=modpow(I%,3,Z%) 710 if E%(W%)=1 then inc J%:Td%(J%)=I% 720 next I% 730 Base=Z%*(Sx\Z%):if Base<0 then Base=Base-Z% 740 for I%=1 to J%:X=Base+Td%(I%) 750 Y2=X^3+La*X+Lb:if Y2<0 then 800 760 if Tb%(Y2@Vb%)=1 then 800 770 if Ta%(Y2@Va%)=1 then 800 780 Y=isqrt(Y2):if res>0 then 800 790 inc F%:print F%;": (";X;",";Y;",";Z%;")" 800 next I% 810 Base=Base+Z%:if Base1 then print:return(0) 920 X1=M^2-2*X:Y1=M*(X-X1)-Y:print:print spc(4);"2P = (";X1;",";Y1;")"; 930 if and{X=X1,Y=(-Y1)} then print " order";C%+1:return(0) 940 M=(Y1-Y)//(X1-X):if den(M)<>1 then print:return(0) 950 Xn=M^2-X-X1:Yn=M*(X-Xn)-Y:inc C%:X1=Xn:Y1=Yn 960 print:print spc(4);cutspc(str(C%));"P = (";X1;",";Y1;")"; 970 goto 930