Subject: Congruent numbers
Date: Fri, 02 Apr 1999 11:39:04 +0000
From: Allan MacLeod
To: kc2h-msm@asahi-net.or.jp
Hisanori,
I finally got organised enough to get together this mail.
There are 3 separate components:
(a) A further list of solutions. These give N, and the numerator
and denominator of ONE side of the right-angled triangle of area N.
Thus the only values under 500 left are 367 and 373. Computations
with the L-series of the elliptic curve suggest they both will
have about 45-55 digits in the solution, but I'm working on them.
(b) The Latex2e file for a description of the 4-descent method.
I have been extending this method to an 8-descent, but I have not
yet written this up. The 8-descent generated most of the larger
solutions I have sent.
(c) A UBASIC file to do the 4-descent. It is set up for the
congruent number elliptic curve, but can be used for any curve
of the form y^2 = x ( x^2 + A x + B ). All you do is change
lines 110 and 120, and define a different set of torsion points
in 200-220. It is possible for the descent method to find a known
torsion point so the program rejects these - assuming the correct
values have been given. NTOR is the number of DISTINCT x-values
giving torsion points, not the total number of torsion points.
Please note that the code outputs results in blocks of 4-digits
but UBASIC does not include leading zeros, so 123 should be read
as 0123, and 45 as 0045, and so on.
If you have any problems with the code please contact me.
Allan
Results with leading zeros inserted manually
N = 302
NUM= 3222 2699 3797 7270 9043
DEN= 651 2555 7288 7396 9390
N = 358
NUM= 1 1495 8231 9978 8309
DEN= 448 0287 0612 8675
N = 446
NUM= 2926 1625 4502 0512
DEN= 188 2156 1327 6049
N = 461
NUM= 35 4245 8984 9428 7468 7926 3101
DEN= 1 9986 6684 4073 6116 9045 4150
N = 478
NUM= 38 9304 0879 5252 2553 1965 4000
DEN= 7 4978 2399 6743 8022 9625 4569
N = 502
NUM= 302 9541 8177 4340
DEN= 21 8061 8794 2363
N = 542
NUM= 202 5303 8190 7827 1177 0377 6560
DEN= 8 1428 1712 5003 0879 6745 0503
N = 566
NUM= 6 3472 4216 3764
DEN= 1 1567 9795 0661
N = 573
NUM= 73 1041 7244 9217 3487
DEN= 3 4633 3902 3546 6360
N = 597
NUM= 2 2346 2941 9695 0994 9337
DEN= 929 4420 8391 2610 3590
N = 613
NUM= 84 4763 3466 8325 7932 2518 6681 0020
DEN= 13 4241 0851 0489 8238 0780 2413 9997
N = 614
NUM= 11 2157 3825 6258 2251
DEN= 3429 2682 8998 1509
N = 623
NUM= 582 6343 8089 5689 6576
DEN= 36 8607 9972 4308 8375
N = 631
NUM= 2 2846 9608 8751 5232 2597 6894 1083
DEN= 7092 1720 7605 9524 0258 4250 2860
N = 661
NUM= 37 1925 6188 2781 4761 0879 7023 9363
DEN= 7 0314 5782 5766 9676 2277 1438 8190
N = 662
NUM= 477 2458 5785 9600 0046 5937
DEN= 24 2705 4669 7268 1246 2015
N = 685
NUM= 1981 9227 6722 9340
DEN= 100 3896 8550 8003
N = 717
NUM= 169 4663 6219 4960
DEN= 349 4898 5679 4871
N = 766
NUM= 36 0680 4383 4916 0272
DEN= 3 5335 6948 0956 2929
N = 773
NUM= 3478 0904 7047 5122 8519 0874 5980 3540
DEN= 90 6544 1731 4969 3914 3615 2596 9763
N = 878
NUM= 32 1083 2481 3901
DEN= 3 3714 4829 3270
N = 893
NUM= 15 5732 2789 2649 5436
DEN= 1 4429 4881 2221 8085
N = 917
NUM= 6 5979 5110 7197 5972
DEN= 2305 6284 1217 1265
N = 926
NUM= 5 4778 4812 8494 0868 1737 5197 5728
DEN= 4014 2816 3974 8811 4352 4809 4809
N = 933
NUM= 14 9556 8797 6751 4191 8679 6773
DEN= 5684 2573 4010 4853 4094 7100
N = 941
NUM= 5639 4866 5024 8149 9426 5529 7021
DEN= 1475 7347 8114 3394 3670 3517 4410
N = 965
NUM= 1953 2964 9309 0138 2740
DEN= 136 6639 8835 9937 4483
N = 974
NUM= 1 7412 1182 7075 3276 1341 4072
DEN= 1325 2392 3110 3802 0930 3837
N = 991
NUM= 1 8494 5993 7818 8935 0214 4063
DEN= 1324 4708 4668 7755 9474 8680
----------------------------- Latex2e file ----------------------------
\documentclass[12pt]{article}
\begin{document}
\begin{center}
{\bf DESCENT FORMULAE FOR ELLIPTIC CURVES}
\end{center}
\bigskip
\section{Four-descent}
The following set of formulae describe the algebra necessary to perform
a general four-descent procedure for the elliptic curve
%eqn 1
\begin{equation}
y^2 = x^3 + a x^2 + b x
\end{equation}
so the curve is assumed to have at least one point of order 2.
We can assume without loss of generality that $x = du^2/v^2$ and $y=duw/v^3$,
with $d$ squarefree and $(d,v)=1$, giving
%eqn 2
\begin{equation}
d^2 w^2 = d^3 u^4 + d^2 a u^2 v^2 + d b v^4
\end{equation}
which implies that $d|b$, so that $b=de$. There are thus a finite set of
possible values for d, which we can work out easily from the factors of
$b$, unless $b$ is very large and difficult to factor. Remember that $d$
can be negative.
We consider first the simpler quadratic
%eqn 3
\begin{equation}
h^2 = d f^2 + a f g + e g^2
\end{equation}
and look for a solution $(h_0,f_0,g_0)$ either by searching or by more
advanced methods. We assume $g_0 \ne 0$.
If we find a solution, we can parameterise as follows. Define $x=f/g$ and
$y=h/g$, so that the equation is
\begin{displaymath}
y^2=d x^2+a x+e
\end{displaymath}
then the line through $(x_0,y_0)=(f_0/g_0,h_0/g_0)$, with gradient $m$,
meets the curve again where
%eqn 4
\begin{equation}
x=\frac{a + d x_0 + m (m x_0 - 2 y_0)}{m^2-d}
\end{equation}
Assuming $m=p/q$ is rational, we simplify this to
%eqn 5
\begin{equation}
\frac{a g_0 q^2 + d f_0 q^2 + p (f_0 p - 2 h_0 q)}{g_0 (p^2 - d q^2)}
\end{equation}
To solve our original problem we look for values of the parameters
giving $f/g=u^2/v^2$, which leads to the two quadratics
%eqn 6
\begin{equation}
k u^2 = g_0 p^2 - g_0 d q^2
\end{equation}
and
%eqn 7
\begin{equation}
k v^2 = f_0 p^2 - 2 h_0 p q + (a g_0 + d f_0)q^2
\end{equation}
with $k$ squarefree.
The possible values of $k$ are those which divide the resultant of the
two quadratics, which means that $k$ must divide the determinant
\begin{center}
\begin{array}{|cccc|}
g_0&0&-d g_0&0 \\
0&g_0&0&-d g_0 \\
f_0&-2 h_0&a g_0+d f_0&0 \\
0&f_0&-2 h_0&a g_0 + d f_0
\end{array}
\end{center}
which reduces to $g_0^4(a^2-4b)$. Let $k_0$ be a squarefree divisor (possibly
negative).
We search the first of these quadratics to find a solution $(u_0,p_0,q_0)$.
If we find one, we can perform another simple line-quadratic intersection
analysis to characterise the solutions to the first quadratic as
%eqn 8
\begin{equation}
p = 2 u_0 k_0 r s - p_0 (g_0 s^2 + k_0 r^2)
\end{equation}
and
%eqn 9
\begin{equation}
q = q_0 ( g_0 s^2 - k_0 r^2)
\end{equation}
Substitute these into $k_0 * (7)$ gives the quartic
%eqn 10
\begin{equation}
z^2 = z_1 r^4 + z_2 r^3 s + z_3 r^2 s^2 + z_4 r s^3 + z_5 s^4
\end{equation}
with
%eqn 11
\begin{equation}
z_1 = k_0^3 (a g_0 q_0^2 + d f_0 q_0^2 + f_0 p_0^2 - 2 h_0 p_0 q_0)
\end{equation}
%eqn12
\begin{equation}
z_2 = 4 u_0 k_0^3 (h_0 q_0 - f_0 p_0)
\end{equation}
%eqn 13
\begin{equation}
z_3 = 2 k_0^2 (f_0 (2 u_0^2 k_0 - d g_0 q_0^2 + g_0 p_0^2) -
a g_0^2 q_0^2)
\end{equation}
%eqn 14
\begin{equation}
z_4 = -(4 u_0 g_0 k_0^2 (f_0 p_0 + h_0 q_0))
\end{equation}
%eqn 15
\begin{equation}
z_5 = g_0^2 k_0 (a g_0 q_0^2 + d f_0 q_0^2 + f_0 p_0^2 + 2 h_0 p_0 q_0)
\end{equation}
\section{COMPUTING}
The above formulae form the basis of a simple code written by the author
in UBASIC. This system is fast, simple, free, and runs on a multitude of PC's,
old and new. The major advantage, however, is that large numbers can be dealt
with by one single statement at the start of the program - the rest of the
code is standard Basic.
The basic structure of the algorithm is
\bigskip
\begin{flushleft}
find divisors d of b
for s1 = s1a to s1b
\hspace{2em}test if $h^2=df^2+afg+eg^2$ soluble with $|f|+|g|=s1$
\hspace{2em}if $(d,h_0,f_0,g_0)$ a solution then
\hspace{4em}find divisors k of $g_0(a^2-4b)$
\hspace{4em}for s2 =2 to s2b
\hspace{6em}test if $k u^2 = g_0 p^2 - g_0 d q^2$ with $|p|+|q|=s2$
\hspace{6em}if $(k,u_0,p_0,q_0)$ a solution, then
\hspace{8em}form quartic equation from equations (10) to (15)
\hspace{8em}for s3=2 to s3b
\hspace{10em}test if quartic is square for $|r|+|s|=s3$
\hspace{10em}if it is, determine point on curve, and stop
\hspace{8em}next s3
\hspace{4em}next s2
next s1
\end{flushleft}
We now make some comments on this method, culled from extensive practical
experience.
If the curve has a rational point of moderate height, the above search-based
procedure works very well. We have found the method to be effective for points
with height of up to about 15 (in the Silverman height normalisation), if we
select s2b=99 and s3b=199. It is impossible to predict in advance the best
choice of s1a and s1b.
For larger heights, the quartic search needs a larger value of s3b and
can take a considerable time. In many
cases, the quartic is insoluble so searching is futile. Cremona describes
how to test this, and we have implemented this test as a precursor to
searching. This method is probably the most advanced mathematically
in the whole procedure.
Cremona also describes various methods for searching the quartic
which reduce the time taken. We have chosen NOT to implement these, partly
to keep the code reasonably understandable by amateurs, and partly because
UBASIC has restrictions on the number of variables.
For investigating families of curves where $a$ and $b$ are functions of $n$,
we found the simple search for $(h_0,f_0,g_0)$ could fail to find any
solutions in the specified range. In such a case we can adapt the search as
follows. If
\begin{displaymath}
h^2 = d f^2 +a f g + b g^2 / d
\end{displaymath}
then
\begin{displaymath}
d(2h)^2 = ( 2df+ ag)^2 - (a^2-4b)g^2
\end{displaymath}
and if we factor $a^2-4b = \alpha \beta^2$ with $\alpha$ squarefree,
%eqn 16
\begin{equation}
dH^2 = F^2 - \alpha G^2
\end{equation}
with $F = 2df+ag , \; G=\beta g, \; H = 2h$. We can recover integer solutions
to the original equations as $f_0=\beta F_0 - a G_0 , \; g_0 = 2 d G_0, \;
h_0 = d \beta H_0$.
In equation (6), it is clear that $(\pm u_0,\pm p_0, \pm q_0)$ satisfy the
equation. We do not, however, need to consider all 8 possible combinations of
sign. From (11) to (15), the sign of $u_0$ is irrelevant if we allow
negative values of r or s. Similarly, $-p_0$ and $-q_0$ lead to the same
quartic as $p_0$ and $q_0$, and the other two possible combinations give
the same quartic. Thus there are only two essentially different quartics,
and it is easy to show the relationship between them implies that they are
both soluble or both insoluble. It is worthwhile to search both (if soluble),
as the coefficients are different so the fixed search range might give a
solution for one but not the other.
At several points in our code, we need to factor numbers to find divisors.
This is done by an extremely simple-minded search procedure, without recourse
to any modern factorisation techniques. So far, this has not had a severe
effect on the performance of the code, but one could easily devise an
elliptic curve where the factorisation time was dominant. The current code
has the advantage that is is short.
Finally, users should note that the code could easily find a torsion point
of the curve, and not a point of infinite order. Since the latter is
usually wanted, the code checks that the point found is not a torsion
point from a list of x-values provided by the user.
\end{document}
------------------------- UBASIC CODE ------------------------------
1 word 40
10 data 2,9999
11 data 199
12 data 399
20 Facsize=128
50 dim Dval(Facsize),Kval(Facsize)
55 dim Fac1(20),Fac2(20),Wait(20,10)
60 dim Z0(5),Tor(9)
70 dim Badp(10),Ibin(20),Pro%(50)
100 input N
110 A=0
120 B=-N*N
130 Ab=A*A-4*B
200 Ntor=3
210 Tor(1)=0
211 Tor(2)=-N
212 Tor(3)=N
500 read Ser1a:read Ser1b:read Ser2b:read Ser3b
600 Numd=2:Dval(1)=1:Dval(2)=-1
610 if abs(B)=1 then goto 1000
615 E=abs(B)
620 gosub *Factor(E,&Fac1())
630 gosub *Divisors(Fac1(),&Dval())
640 Numd=Dval(0)
1000 for Ser1=Ser1a to Ser1b
1002 print Ser1;
1005 for Af=1 to Ser1-1
1010 G0=Ser1-Af
1015 if gcd(Af,G0)>1 then goto 9890
1020 G02=G0*G0
1025 for F0=-Af to Af step (2*Af)
1030 F02=F0*F0:Fg0=F0*G0
1035 for Id=1 to Numd
1040 D=Dval(Id):E=B\D
1045 Ib1=D*F02+A*Fg0+E*G02
1050 if Ib1<=0 then goto 9870
1055 H0=isqrt(Ib1)
1060 if res>0 then goto 9870
1070 print:print "a:";D,H0,F0,G0
1100 Resl=abs(G0*Ab)
1105 gosub *Factor(Resl,&Fac1())
1110 gosub *Divisors(Fac1(),&Kval())
1120 Numk=Kval(0)
1200 for Ser2=2 to Ser2b
1210 for P0=1 to Ser2-1
1220 Q0=Ser2-P0
1225 if gcd(P0,Q0)>1 then goto 9780
1230 for Ik=1 to Numk
1235 Ib1=G0*(P0*P0-D*Q0*Q0)
1240 K0=Kval(Ik)
1245 Ib2=Ib1\K0
1250 if res>0 then goto 9770
1255 if Ib2<=0 then goto 9770
1260 U0=isqrt(Ib2)
1270 if res>0 then goto 9770
1280 print " b:";K0,U0,P0,Q0
1300 Z0(1)=K0^3*(A*G0*Q0^2+D*F0*Q0^2+F0*P0^2-2*H0*P0*Q0)
1305 Z0(2)=4*U0*K0^3*(H0*Q0-F0*P0)
1310 Z0(3)=2*K0^2*(F0*(2*U0^2*K0-D*G0*Q0^2+G0*P0^2)-A*G0^2*Q0^2)
1315 Z0(4)=-4*U0*G0*K0^2*(F0*P0+H0*Q0)
1320 Z0(5)=G0^2*K0*(A*G0*Q0^2+D*F0*Q0^2+F0*P0^2+2*H0*P0*Q0)
1325 Ib1=gcd(Z0(1),Z0(2)):Ib2=gcd(Ib1,Z0(3))
1330 Ib3=gcd(Ib2,Z0(4)):Ib4=gcd(Ib3,Z0(5))
1335 if Ib4=1 then goto 1400
1340 Ib1=1
1350 gosub *Factor(Ib4,&Fac1())
1355 for I=1 to Fac1(0)
1360 Ib2=Fac1(I+I)\2
1365 Ib1=Ib1*(Fac1(I+I-1)^Ib2)
1370 next I
1380 for I=1 to 5
1385 Z0(I)=Z0(I)\(Ib1*Ib1)
1390 next I
1400 gosub *Testsol(Z0(),&Is1)
1410 if Is1=1 then goto 9770
1420 for Iq0=1 to 2
1430 if Iq0=1 then goto 1500
1435 Q0=-Q0
1450 Z0(1)=K0^3*(A*G0*Q0^2+D*F0*Q0^2+F0*P0^2-2*H0*P0*Q0)
1455 Z0(2)=4*U0*K0^3*(H0*Q0-F0*P0)
1460 Z0(3)=2*K0^2*(F0*(2*U0^2*K0-D*G0*Q0^2+G0*P0^2)-A*G0^2*Q0^2)
1465 Z0(4)=-4*U0*G0*K0^2*(F0*P0+H0*Q0)
1470 Z0(5)=G0^2*K0*(A*G0*Q0^2+D*F0*Q0^2+F0*P0^2+2*H0*P0*Q0)
1500 for Ser3=2 to Ser3b
1510 for R=1 to Ser3-1
1515 S=Ser3-R
1520 if gcd(R,S)>1 then goto 9670
1530 Ib1=R*R:Ib2=R*S:Ib3=S*S
1540 Ic1=Ib1*(Z0(1)*Ib1+Z0(3)*Ib3)+Z0(5)*Ib3*Ib3
1545 Ic2=Ib2*(Z0(2)*Ib1+Z0(4)*Ib3)
1550 Id1=Ic1+Ic2
1555 if Id1<=0 then goto 1600
1560 Id2=isqrt(Id1)
1570 if res>0 then goto 1600
1575 print " c:";R;S
1580 gosub *Testpt(&Indpt)
1590 if Indpt=1 then goto 6000
1600 Id1=Ic1-Ic2
1610 if Id1<=0 then goto 9670
1620 Id2=isqrt(Id1)
1630 if res>0 then goto 9670
1640 S=-S
1645 print " c:";R;S
1650 gosub *Testpt(&Indpt)
1660 if Indpt=0 then goto 9670
6000 print:print
6010 print "SOLUTION FOUND"
6015 print:print "N = ";N
6020 print:print "d = ";D
6025 print
6030 print "u = ";:gosub *Prout(U)
6035 print:print
6040 print "v = ";:gosub *Prout(V)
6045 print:print
6050 print "w = ";:gosub *Prout(W)
6055 print:print
6060 X=D*U^2//V^2:Y=D*U*W//V^3
6065 Ib1=num(X):Ib2=den(X)
6070 print "x = ";:gosub *Prout(Ib1)
6075 print:print "over"
6080 print " ";:gosub *Prout(Ib2)
6085 print
6090 Ib1=num(Y):Ib2=den(Y)
6092 print
6095 print "y = ";:gosub *Prout(Ib1)
6100 print:print "over"
6105 print " ";:gosub *Prout(Ib2)
6110 print
6115 print:print
6199 stop
9670 next R
9680 next Ser3
9690 next Iq0
9770 next Ik
9780 next P0
9790 next Ser2
9870 next Id
9880 next F0
9890 next Af
9900 next Ser1
9999 end
11000 *Testpt(&Indp)
11005 Indp=0
11010 Ie1=Q0*(G0*S^2-K0*R^2)
11015 if Ie1=0 then goto 11199
11020 Ie2=2*U0*K0*R*S-P0*(G0*S^2+K0*R^2)
11025 Ie3=Ie2//Ie1
11030 Ie1=num(Ie3):Ie2=den(Ie3)
11050 Ie3=A*G0*Ie2^2+D*F0*Ie2^2+Ie1*(F0*Ie1-2*H0*Ie2)
11055 Ie4=G0*(Ie1^2-D*Ie2^2)
11060 if Ie4=0 then goto 11199
11070 Ie5=Ie3//Ie4
11075 Ie1=num(Ie5):Ie2=den(Ie5)
11080 if Ie1<=0 then goto 11199
11085 Ie3=isqrt(Ie1)
11090 if res>0 then goto 11199
11100 Ie4=isqrt(Ie2)
11105 if res>0 then goto 11199
11110 Ie5=D*Ie1^2+A*Ie1*Ie2+E*Ie2^2
11115 if Ie5<=0 then goto 11199
11120 Ie6=isqrt(Ie5)
11125 if res>0 then goto 11199
11130 Ie1=D*Ie3^2//(Ie4^2):Ie2=0
11135 for Id1=1 to Ntor
11140 if Ie1=Tor(Id1) then Ie2=1
11145 next Id1
11150 if Ie2=1 then goto 11199
11160 Ie2=D*Ie3*Ie6//(Ie4^3)
11170 Id1=Ie2*Ie2-Ie1*(Ie1*Ie1+A*Ie1+B)
11175 if Id1<>0 then goto 11199
11180 Indp=1
11190 U=Ie3:V=Ie4:W=Ie6
11199 return
12000 *Prout(Nn)
12010 N9=abs(Nn)
12015 if Nn<0 then print " - ";
12020 Nd9%=0
12030 N91=N9\10000:N92=res
12040 inc Nd9%
12050 Pro%(Nd9%)=N92
12060 if N91=0 then goto 12100
12070 N9=N91
12080 goto 12030
12100 for I9%=Nd9% to 1 step -1
12110 print Pro%(I9%);
12120 next I9%
12199 return
13000 *Divisors(F1(),&F2())
13005 Ib1=F1(0)
13010 F2(0)=2:F2(1)=1:F2(2)=-1
13015 Ib2=2^(Ib1+1)
13020 if Ib1>1 then goto 13100
13025 if F1(1)=1 then goto 13199
13030 F2(0)=4:F2(3)=F1(1):F2(4)=-F2(3)
13035 goto 13199
13100 for Ic1=1 to Ib1:Ibin(Ic1)=0:next Ic1
13105 Ibin(1)=1
13110 if Ib2>Facsize then Ib2=Facsize
13115 for Ic2=3 to Ib2 step 2
13120 Id1=1
13125 for Id2=1 to Ib1
13130 Id1=Id1*(F1(Id2+Id2-1)^Ibin(Id2))
13135 next Id2
13140 F2(Ic2)=Id1:F2(Ic2+1)=-Id1
13145 Id1=1
13150 for Id2=1 to Ib1
13155 Id3=Ibin(Id2)+Id1
13160 if Id3=2 then goto 13175
13165 Ibin(Id2)=1
13170 cancel for:goto 13190
13175 Ibin(Id2)=0
13180 next Id2
13190 next Ic2
13195 F2(0)=Ib2
13199 return
14000 *Testsol(Z(),&Inn)
14005 Inn=0
14010 Bi=12*Z(1)*Z(5)-3*Z(2)*Z(4)+Z(3)*Z(3)
14015 Bj=72*Z(1)*Z(3)*Z(5)+9*Z(2)*Z(3)*Z(4)
14020 Bj=Bj-27*Z(1)*Z(4)*Z(4)-27*Z(5)*Z(2)*Z(2)-2*Z(3)*Z(3)*Z(3)
14025 Delta=4*Bi*Bi*Bi-Bj*Bj
14030 Ia1=abs(Delta):gosub *Factor(Ia1,&Fac1())
14035 Nbadp=Fac1(0)
14040 for I=1 to Nbadp
14045 Badp(I)=Fac1(I+I-1)
14050 next I
14060 for Ip=1 to Nbadp
14065 Ip1=Badp(Ip)
14070 gosub *Qpsol(Z(1),Z(2),Z(3),Z(4),Z(5),Ip1,&Qeq)
14075 if Qeq=1 then goto 14099
14080 Inn=1
14090 cancel for:goto 14199
14099 next Ip
14199 return
14300 *Qpsol(A,B,C,D,E,P,&Resq)
14310 gosub *Zpsol(A,B,C,D,E,P,&Resq)
14320 if Resq=1 then goto 14349
14330 gosub *Zpsol(E,D,C,B,A,P,&Resq)
14349 return
14400 *Zpsol(A,B,C,D,E,P,&Resz)
14402 Ib2=0
14404 Ib3=0:Ib4=1:Ib5=0:Ib6=1
14410 Xx=Ib3+Ib5*Ib6
14420 if P=2 then Ib1=fnLemma7(A,B,C,D,E,Ib4,Xx)
14422 if P>2 then Ib1=fnLemma6(A,B,C,D,E,P,Ib4,Xx)
14430 if Ib1<1 then goto 14440
14432 Resz=1:goto *Zpend
14440 if Ib1=0 then goto 14460
14442 Ib5=Ib5+1
14444 if Ib5<P then goto 14410
14445 if Ib2>0 then goto 14450
14447 Resz=0:goto *Zpend
14450 Ib3=Wait(Ib2,1):Ib4=Wait(Ib2,2):Ib5=Wait(Ib2,3):Ib6=Wait(Ib2,4)
14452 dec Ib2
14454 goto 14410
14460 if Ib5<(P-1) then goto 14470
14462 goto 14480
14470 inc Ib2
14472 Wait(Ib2,1)=Ib3:Wait(Ib2,2)=Ib4:Wait(Ib2,3)=Ib5+1:Wait(Ib2,4)=Ib6
14480 inc Ib4:Ib5=0:Ib6=P*Ib6
14482 goto 14410
14499 *Zpend:return
14700 fnLemma6(A,B,C,D,E,P,N,X)
14702 Ic1=(((A*X+B)*X+C)*X+D)*X+E
14704 Ic2=fnPadsq(Ic1,P)
14706 if Ic2=0 then goto 14710
14708 Resl=1:goto *Lem6end
14710 Ic3=((4*A*X+3*B)*X+2*C)*X+D
14712 Ic4=fnOrdp(Ic1,P)
14714 Ic5=fnOrdp(Ic3,P)
14720 if Ic4<(Ic5+N) then goto 14730
14722 if N<=Ic5 then goto 14730
14724 Resl=1:goto *Lem6end
14730 if Ic4<(2*N) then goto 14740
14732 if Ic5<N then goto 14740
14734 Resl=0:goto *Lem6end
14740 Resl=-1
14749 *Lem6end:return(Resl)
14800 fnLemma7(A,B,C,D,E,N,X)
14801 P2=2
14802 Ic1=(((A*X+B)*X+C)*X+D)*X+E
14804 Ic2=fnPadsq(Ic1,P2)
14806 if Ic2=0 then goto 14810
14808 Resl=1:goto *Lem7end
14810 Ic3=((4*A*X+3*B)*X+2*C)*X+D
14812 Ic4=fnOrdp(Ic1,P2)
14814 Ic5=fnOrdp(Ic3,P2)
14816 Ic6=Ic1
14820 if odd(Ic6)=1 then goto 14830
14822 Ic6=Ic6\2
14824 goto 14820
14830 Ic7=Ic6\4:Ic6=res
14832 if Ic4<(Ic5+N) then goto 14840
14834 if N<=Ic5 then goto 14840
14836 Resl=1:goto *Lem7end
14840 if N<=Ic5 then goto 14850
14842 if Ic4<>(Ic5+N-1) then goto 14850
14844 if odd(Ic4)=1 then goto 14850
14846 Resl=1:goto *Lem7end
14850 if N<=Ic5 then goto 14860
14852 if Ic4<>(Ic5+N-2) then goto 14860
14854 if Ic6<>1 then goto 14860
14856 Resl=1:goto *Lem7end
14860 if Ic5<N then goto 14870
14862 if Ic4<(N+N) then goto 14870
14864 Resl=0:goto *Lem7end
14870 if Ic5<N then goto 14880
14872 if Ic4<>(N+N-2) then goto 14880
14874 Resl=0:goto *Lem7end
14880 Resl=-1
14899 *Lem7end:return(Resl)
14900 fnPadsq(A,P)
14902 if A<>0 then jump
14904 Resl=1:goto *Padend
14910 **:Id1=fnOrdp(A,P)
14912 if odd(Id1)=0 then jump
14914 Resl=0:goto *Padend
14920 **:Id2=A\(P^Id1)
14922 if P>2 then goto 14940
14924 Id1=Id2\8:Id2=res
14926 Resl=0
14928 if Id2=1 then Resl=1
14930 goto *Padend
14940 Resl=0
14942 if kro(Id2,P)=1 then Resl=1
14949 *Padend:return(Resl)
19000 fnRabmil(N)
19005 ' rabin-miller test for n prime
19010 ' Result=0 Composite Result=1 Prime
19020 Ie2=even(N):Ie1=1
19025 if Ie2=0 jump
19030 Ie1=0:if N=2 then Ie1=1
19032 goto *Finrm
19035 **:if N=1 goto *Finrm
19040 Ie3=N-1:Ie4=0:Ie2=0
19045 while Ie2=0
19050 Ie6=Ie3\2:Ie2=res
19055 Ie3=Ie6:inc Ie4
19060 wend
19065 Ie3=Ie3+Ie3+1:dec Ie4
19070 randomize
19075 for Ie7=1 to 20
19080 Ie6=int((N-2)*rnd+2):Ie2=0
19085 Ie5=modpow(Ie6,Ie3,N)
19090 if Ie5=1 jump
19095 Ie8=modpow(Ie5,1,N)
19100 while and{Ie8>1,Ie8<N-1,Ie2<=Ie4-2}
19105 Ie5=modpow(Ie5,2,N):inc Ie2
19107 Ie8=modpow(Ie5,1,N)
19110 wend
19115 if Ie5=N-1 jump
19117 Ie1=0
19120 cancel for
19125 goto *Finrm
19130 **:next Ie7
19140 Ie1=1
19150 *Finrm:return(Ie1)
19200 *Factor(N,&Fac())
19205 ' factors n - very simple search procedure
19210 ' fac(0)=no. of factors
19215 ' n = fac(1)^fac(2)*fac(3)^fac(4)*.....
19220 if N>1 goto 19240
19225 Fac(0)=1:Fac(1)=1:Fac(2)=1
19230 goto *Finfac
19240 Id6=fnRabmil(N)
19245 if Id6=0 goto 19260
19250 Fac(0)=1:Fac(1)=N:Fac(2)=1
19255 goto *Finfac
19260 Id5=0
19265 for Id1=1 to 1000000
19267 if Id1>12251 goto 19273
19270 Id2=prm(Id1)
19271 goto 19275
19273 Id2=Id2+2
19275 Id3=N\Id2:Id4=res
19280 if Id4>0 goto *Nxfac2
19285 Id5=Id5+1:Fac(0)=Id5
19290 Fac(Id5+Id5-1)=Id2:Fac(Id5+Id5)=1
19295 N=Id3
19300 Id3=N\Id2:Id4=res
19305 if Id4>0 goto *Nxfac1
19310 Fac(Id5+Id5)=Fac(Id5+Id5)+1:N=Id3
19315 if N>1 goto 19300
19320 cancel for
19325 goto *Finfac
19330 *Nxfac1:Id6=fnRabmil(N)
19332 if Id6=0 goto *Nxfac2
19334 Id5=Id5+1:Fac(0)=Id5
19336 Fac(Id5+Id5-1)=N:Fac(Id5+Id5)=1
19337 cancel for
19338 goto *Finfac
19340 *Nxfac2:next Id1
19349 *Finfac:return
19350 fnOrdp(N,P)
19352 ' finds order of p in N
19354 if abs(N)>1 then goto 19360
19356 Ie2=0
19358 if N=0 then Ie2=123456789
19359 goto 19399
19360 Ie1=abs(N):Ie2=0
19365 Ie3=Ie1\P
19370 if res>0 then jump
19375 Ie2=Ie2+1:Ie1=Ie3
19380 goto 19365
19399 **:return(Ie2)
19800 fnMods(N,P)
19805 ' finds symmetric residue mod p
19815 Id1=N\P:Id2=res
19820 Id1=P\2
19825 if Id2>Id1 then Id2=Id2-P
19830 return(Id2)
E-mail : kc2h-msm@asahi-net.or.jpHisanori Mishima