1   word 80
    2   point 8
    3   Peye=#pi:Gamma=#euler
    4   Eeps=1000*#eps:Rteps=sqrt(#eps)
    5   clr time
  100   dim A(6),Aa(6),Bddat(10,5),Perd(2),Bb(8),U(4),Ibb(3),Torpt(16,2)
  110   dim F(20),B(8)
  150   a(1)=
  160   a(2)=
  170   a(3)=
  180   a(4)=
  190   a(6)=
  210   gosub *Printeq(A())
  300   gosub *Reduce(A(),&Aa(),&Bddat())
  310   print:print "Minimal equation":gosub *Printeq(Aa())
  312   gosub *Urst(A(),Aa(),&U())
  313   print U(1),U(2),U(3),U(4)
  320   Ncon=Bddat(0,1)
  325   print "Conductor=",Ncon
  340   gosub *Printred(Bddat())
  350   Tam=1
  360   for I=1 to Bddat(0,0):Tam=Tam*Bddat(I,2):next I
  365   print:print "Tamagawa number=",Tam
  370   print:print "Time=";time1000:print " "
  402   print "Input no of torsion points known"
  404   input Ntor
  406   if Ntor>0 then goto 417
  410   gosub *Torspts(Aa(),&Ntor)
  415   print "No of torsion points = ";Ntor
  417   print 1,"Infinity"
  420   if Ntor=1 then jump
  422   for I=2 to Ntor
  424   Ia1=Torpt(I,1):Ia2=Torpt(I,2)
  425   Ia3=U(1)*U(1)*Ia1+U(2)
  426   Ia4=U(1)^3*Ia2+U(3)*U(1)*U(1)*Ia1+U(4)
  428   print I,"(";Ia3;",";Ia4;")"
  429   Torpt(I,1)=Ia3:Torpt(I,2)=Ia4
  430   next I
  435   **:print:print "Time=";time1000:print " "
  440   gosub *Period(Aa(),&Perd())
  445   W1=Perd(1)
  452   print " w1 = ";W1
  454   print " w2 = ";Perd(2)
  456   print " Omega = ";Perd(0)
  457   Q=-exp(-2*Peye*im(Perd(2))/Perd(1))
  458   gosub *Form71(Aa(),&Bb())
  459   if Bb(1)>0 then Q=abs(Q)
  460   print " q = ";Q
  465   Tau=Perd(2)/W1
  468   print " tau = ";Tau
  482   print "how many a(i) values do you want for conductor = ";Ncon
  484   input Na
  486   if Na>0 then goto 500
  490   stop
  500   dim An%(Na)
  502   gosub *Anvals(Aa(),Bddat(),Na,&An%())
  515   Rb1=2*Peye/sqrt(Ncon)
  520   Rb2=0.90*Rb1:Rb3=Rb1/0.90:Rb4=1.10*Rb1:Rb5=Rb1/1.10
  525   Rb6=0:Rb7=0:Rb8=0:Rb9=0
  530   for I=1 to Na
  532   if An%(I)=0 then goto 545
  533   Rb10=exp(-I*Rb2):Rb11=exp(-I*Rb3):Rb12=exp(-I*Rb4):Rb13=exp(-I*Rb5)
  535   Rb6=Rb6+An%(I)*(Rb10+Rb11)/I:Rb7=Rb7+An%(I)*(Rb12+Rb13)/I
  540   Rb8=Rb8+An%(I)*(Rb10-Rb11)/I:Rb9=Rb9+An%(I)*(Rb12-Rb13)/I
  545   next I
  550   print "Sign =+1",Rb6,Rb7
  555   print "Sign =-1",Rb8,Rb9
  560   Rb10=0:Rb11=0
  565   if abs(Rb8)<0.001 then Rb10=1
  570   if abs(Rb9)<0.001 then Rb11=1
  575   if (Rb10+Rb11)=2 then goto 585
  580   print "Not enough evidence that sign =-1. Try with more ai values"
  582   stop
  585   print:print "program assumes sign =  -1"
  600   print:print "Computing L'(E,1)"
  610   Rc1=0.0:Ra1=2*Peye/sqrt(Ncon)
  615   for I=1 to Na
  617   if An%(I)=0 then goto 625
  620   Ra2=Ra1*I:Ra3=fnEi(Ra2)
  622   Rc1=Rc1+An%(I)*Ra3/I
  625   next I
  630   Rc1=Rc1+Rc1
  635   print " L'(E,1) = ";Rc1
  650   print:print " Height = ";
  655   Ht=Rc1*Ntor*Ntor/(2*Perd(0)*Tam)
  670   print Ht
  680   print:print "Time=";time1000
  700   Ic1=1:Ic2=0
  705   for Ic3=1 to Bddat(0,0)
  710   if Bddat(Ic3,3)=1 then goto 790
  715   if Bddat(Ic3,3)>2 then goto 720
  716   Ic4=Bddat(Ic3,4)\2+1
  717   if Ic4>Ic2 then Ic2=Ic4
  718   Ic1=Ic1*Ic4
  719   goto 790
  720   if Bddat(Ic3,3)>3 then goto 725
  721   Ic4=3
  722   if Ic4>Ic2 then Ic2=Ic4
  723   Ic1=Ic1*Ic4
  724   goto 790
  725   if Bddat(Ic3,3)>4 then goto 740
  726   Ic4=1
  727   if Ic4>Ic2 then Ic2=Ic4
  728   Ic1=Ic1*Ic4
  729   goto 790
  740   if Bddat(Ic3,3)>=10 then goto 750
  741   Ic4=2
  742   if Ic4>Ic2 then Ic2=Ic4
  743   Ic1=Ic1*Ic4
  744   goto 790
  750   Ic4=1
  752   if Ic4>Ic2 then Ic2=Ic4
  753   Ic1=Ic1*Ic4
  754   goto 790
  790   next Ic3
  792   print:print "Number of local heights = ";Ic1
  795   dim Lbadt(Ic2),Lbad(Ic1),Dlim(Ic1),Stp(Ic1)
  800   print:print "Computing Local heights"
  802   Nlb=1:Lbad(1)=0
  804   for I=1 to Bddat(0,0)
  806   if Bddat(I,3)=1 then goto *Lbend1
  810   Ra1=log(Bddat(I,1)):Ib1=Bddat(I,0):Ra2=Ib1/12.0
  812   if Bddat(I,3)>2 then jump
  814   Ib2=Ib1\2+1
  816   for Ib4=0 to Ib2-1:Lbadt(Ib4+1)=Ra1*(Ra2-Ib4/2.0+(Ib4*Ib4)/(Ib1+Ib1)):next Ib4
  818   goto *Lbend
  820   **:if Bddat(I,3)>3 then jump
  824   Ib2=3:Lbadt(1)=Ra2*Ra1
  826   Ic1=Bddat(I,4):Lbadt(2)=Ra1*(Ic1/12.0+(Ib1-Ic1-6)/12.0)
  828   Lbadt(3)=Ra1*(-Ic1/24.0+(Ib1-Ic1-6)/12.0)
  830   goto *Lbend
  834   **:if Bddat(I,3)>4 then jump
  836   Ib2=1:Lbadt(1)=Ra2*Ra1:goto *Lbend
  838   **:if Bddat(I,3)>5 then jump
  840   Lbadt(1)=Ra2*Ra1
  844   Ib2=2:Lbadt(2)=Ra1*(Ib1-3)/12.0:goto *Lbend
  846   **:if Bddat(I,3)>6 then jump
  847   Lbadt(1)=Ra2*Ra1
  849   Ib2=2:Lbadt(2)=Ra1*(Ib1-4)/12.0:goto *Lbend
  850   **:if Bddat(I,3)>7 then jump
  851   Lbadt(1)=Ra1*Ra2
  853   Ib2=2:Lbadt(2)=Ra1*(Ra2-0.5):goto *Lbend
  854   **:if Bddat(I,3)>8 then jump
  856   Lbadt(1)=Ra2*Ra1
  860   Ib2=2:Lbadt(2)=Ra1*(Ib1-8)/12:goto *Lbend
  862   **:if Bddat(I,3)>9 then jump
  864   Lbadt(1)=Ra2*Ra1
  868   Ib2=2:Lbadt(2)=(Ib1-9)*Ra1/12.0:goto *Lbend
  870   **:Ib2=1:Lbadt(1)=Ra2*Ra1
  872   *Lbend:Ib6=0
  874   for Ib5=1 to Nlb:Ra1=Lbad(Ib5)
  875   for Ib3=1 to Ib2:Ra2=Ra1+Lbadt(Ib3)
  876   Ib6=(Ib3-1)*Nlb+Ib5:Lbad(Ib6)=Ra2
  877   next Ib3
  878   next Ib5
  880   Nlb=Nlb*Ib2
  890   *Lbend1:next I
  900   if Nlb=1 then goto 1000
  910   for Ic1=1 to Nlb-1
  915   Ic3=Ic1:Rc1=Lbad(Ic1)
  920   for Ic2=Ic1+1 to Nlb
  923   Rc2=Lbad(Ic2)
  926   if Rc2>=Rc1 then goto 940
  930   Ic3=Ic2:Rc1=Rc2
  940   next Ic2
  945   if Ic3=Ic1 then goto 960
  950   Lbad(Ic3)=Lbad(Ic1)
  955   Lbad(Ic1)=Rc1
  960   next Ic1
 1000   Ie2=0
 1001   Id3=0:Id4=0:Id5=0
 1005   if Q<0 then goto 1060
 1010   if Ntor>1 then goto 1020
 1015   Ie2=2
 1016   goto 1060
 1020   Rb1=1.0:Rb2=A(1)*A(1)/4+A(2)
 1021   Rb3=A(1)*A(3)/2+A(4):Rb4=A(3)*A(3)/4+A(6)
 1030   gosub *Cubsol(Rb1,Rb2,Rb3,Rb4,&Ie3,&Ibb())
 1035   Ie2=2
 1040   for I=2 to Ntor
 1045   if Torpt(I,1)<=Ibb(2) then Ie2=1
 1050   next I
 1060   if ie2>0 then goto 1070
 1063   print " No second component"
 1066   goto 1100
 1070   if ie2=1 then goto 1080
 1073   print " Must do second component section"
 1076   goto 1100
 1080   print " Could skip second component"
 1082   print " Input 2 to force second component"
 1085   input ie3
 1090   if ie3=2 then ie2=2
 1100   Ie1=fnGetnq(Q)
 1105   Rc1=sqrt(Q)
 1110   Rc2=(abs(Q)^(1.0/12.0))
 1115   Rc3=sqrt(Rc2)
 1120   for I=1 to Nlb
 1122   Lbad(I)=exp(Lbad(I)-Ht)/Rc2
 1125   next I
 1127   Tam=0
 1128   if Ie2<2 then goto 1140
 1130   gosub *Gsdeval(0,Rc1,&Rc4,&Rc5)
 1135   gosub *Gsdeval(2,Rc1,&Tam,&Rc5)
 1140   gosub *Fsdeval(2,Q,&Rc5,&Rc6)
 1150   Rc8=-0.107653919
 1155   Rc9=(abs(Q))^(1.0/8.0)
 1160   Ie4=0
 1161   if Q>=Rc8 then Ie4=1
 1170   if Ie4=1 then goto 1200
 1175   Rb10=1.55
 1176   gosub *Fsdeval(Rb10,Q,&Rb12,&Rb11)
 1177   Rb12=1.75
 1178   gosub *Fsdeval(Rb12,Q,&Rb14,&Rb13)
 1180   Rb14=Rb12-Rb13*(Rb12-Rb10)/(Rb13-Rb11)
 1182   gosub *Fsdeval(Rb14,Q,&Rb16,&Rb15)
 1184   if abs(Rb14-Rb12)<10^(-10) then goto 1190
 1186   Rb10=Rb12:Rb11=Rb13
 1187   Rb12=Rb14:Rb13=Rb15
 1188   goto 1180
 1190   Rc6=Rb14
 1195   Rc7=Rb16
 1200   Maxf=Rc5
 1210   if Ie4=1 then goto 1250
 1220   Maxf=Rc7
 1250   Rb1=Maxf
 1260   Rb2=Tam/Rc9
 1265   if Rb2>Rb1 then Rb1=Rb2
 1270   for I=1 to Nlb
 1280   Dlim(I)=int(Rb1/Lbad(I))
 1290   next I
 1300   Rb11=1.0
 1305   Rb12=Q
 1310   for I=1 to Ie1
 1315   Rb11=Rb11*(1-Rb12)
 1320   Rb12=Rb12*Q
 1325   next I
 1330   Rb11=Rb11*Rb11
 1335   for I=1 to Nlb
 1340   Stp(I)=Lbad(I)/Rb11
 1345   next I
 1350   Id1=0
 1360   for I=Nlb to 1 step -1
 1365   if Dlim(I)=0 then Id1=I
 1370   next I
 1380   if Id1>0 then Nlb=Id1-1
 1390   Id1=Dlim(1)
 1400   gosub *Form71(Aa(),&Bb())
 2000   for D=1 to Id1
 2001   print d;
 2005   Ibb(1)=Bb(2)*D*D:Ibb(2)=2*Bb(4)*D^4:Ibb(3)=Bb(6)*D^6
 2010   for Id2=1 to Nlb
 2020   if D<=Dlim(Id2) then goto 2100
 2030   Nlb=Id2-1
 2035   cancel for:goto 3960
 2100   Rb19=D*Lbad(Id2)
 2105   if Rb19>Maxf then goto 3000
 2110   Rb11=Stp(Id2)
 2200   gosub *Fsdeval(Rb11,Q,&Rb12,&Rb20)
 2205   Rb12=Rb12-Rb19
 2210   Rb18=Rb11-Rb12/Rb20
 2220   if abs(Rb18-Rb11)<Eeps then goto 2300
 2230   Rb11=Rb18
 2240   goto 2200
 2300   Stp(Id2)=Rb18+Lbad(Id2)/Rb20
 2320   Rb11=asin(Rb18/2.0)/Peye
 2325   Rb12=D*D*(fnWp(Rb11,Q,W1)-Bb(2)/12)
 2330   Rb13=re(Rb12):Rb14=round(Rb13)
 2340   Rb15=abs(Rb13-Rb14)
 2355   if Rb15>0.1 then goto 2500
 2400   Ib1=((4*Rb14+Ibb(1))*Rb14+Ibb(2))*Rb14+Ibb(3)
 2405   if Ib1<=0 then goto 2500
 2410   Ib2=isqrt(Ib1)
 2415   if res>0 then 2500
 2420   print:print " Solution with f(s) first point"
 2430   Ia1=Rb14//D^2
 2435   Ia3=Ib2-Rb14*Aa(1)*D-Aa(3)*D*D*D
 2440   Ia2=Ia3//(2*D^3)
 2445   Ia3=U(1)*U(1)*Ia1+U(2)
 2450   Ia4=U(1)^3*Ia2+U(3)*U(1)*U(1)*Ia1+U(4)
 2455   print "x=";Ia3
 2460   print "y=";Ia4
 2470   print time1000
 2480   goto 3980
 2500   if Ie4=1 then goto 3000
 2510   if Rb19<Rc5 then goto 3000
 2520   Rb11=2*Rc6-Rb18
 2530   if Rb11>2 then Rb11=1.95
 2600   gosub *Fsdeval(Rb11,Q,&Rb12,&Rb13)
 2605   Rb12=Rb12-Rb19
 2610   Rb17=Rb11-Rb12/Rb13
 2620   if abs(Rb17-Rb11)<Eeps then goto 2700
 2630   Rb11=Rb17
 2640   goto 2600
 2700   Rb11=asin(Rb17/2.0)/Peye
 2705   Rb12=D*D*(fnWp(Rb11,Q,W1)-Bb(2)/12)
 2710   Rb13=re(Rb12):Rb14=round(Rb13)
 2720   Rb15=abs(Rb13-Rb14)
 2725   if Rb15>0.1 then goto 3000
 2800   Ib1=((4*Rb14+Ibb(1))*Rb14+Ibb(2))*Rb14+Ibb(3)
 2805   if Ib1<=0 then 3000
 2810   Ib2=isqrt(Ib1)
 2815   if res>0 then 3000
 2820   print:print " Solution with f(s) second point"
 2830   Ia1=Rb14//D^2
 2835   Ia3=Ib2-Rb14*Aa(1)*D-Aa(3)*D*D*D
 2840   Ia2=Ia3//(2*D^3)
 2845   Ia3=U(1)*U(1)*Ia1+U(2)
 2850   Ia4=U(1)^3*Ia2+U(3)*U(1)*U(1)*Ia1+U(4)
 2855   print "x=";Ia3
 2860   print "y=";Ia4
 2870   print time1000
 2880   goto 3980
 3000   if Ie2<2 then goto 3950
 3100   Rb19=Rc9*Rb19
 3105   if Rb19<=Rc4 then goto 3950
 3106   if Rb19>=Tam then goto 3950
 3110   Rb11=Stp(Id2)-Lbad(Id2)/Rb20
 3125   if Stp(Id2)>2.0 then Stp(Id2)=1.95
 3200   gosub *Gsdeval(Rb11,Rc1,&Rb12,&Rb13)
 3205   Rb12=Rb12-Rb19
 3210   Rb18=Rb11-Rb12/Rb13
 3220   if abs(Rb18-Rb11)<Eeps then goto 3300
 3230   Rb11=Rb18
 3235   if Rb11<0 then Rb11=-Rb11
 3240   goto 3200
 3300   Rb11=asin(Rb18/2.0)/Peye+Tau/2
 3305   Rb12=D*D*(fnWp(Rb11,Q,W1)-Bb(2)/12)
 3310   Rb13=re(Rb12):Rb14=round(Rb13)
 3320   Rb15=abs(Rb13-Rb14)
 3325   if Rb15>0.1 then goto 3950
 3400   Ib1=((4*Rb14+Ibb(1))*Rb14+Ibb(2))*Rb14+Ibb(3)
 3405   if Ib1<=0 then goto 3950
 3410   Ib2=isqrt(Ib1)
 3415   if res>0 then 3950
 3420   print:print " Solution with g(s)"
 3430   Ia1=Rb14//D^2
 3435   Ia3=Ib2-Rb14*Aa(1)*D-Aa(3)*D*D*D
 3440   Ia2=Ia3//(2*D^3)
 3445   Ia3=U(1)*U(1)*Ia1+U(2)
 3450   Ia4=U(1)^3*Ia2+U(3)*U(1)*U(1)*Ia1+U(4)
 3455   print "x=";Ia3
 3460   print "y=";Ia4
 3470   print time1000
 3480   goto 3980
 3950   next Id2
 3960   if Nlb>0 then goto 3975
 3970   print "no more solutions possible so stopping"
 3971   cancel for:goto 3999
 3975   next D
 3980   print Id3,Id4,Id5
 3995   print
 3999   end
 4000   *Fsdeval(S,Q,&Fs,&Fds)
 4001   inc Id3
 4010   Fs=1:Fds=0
 4020   Rb2=Q:Rb4=S*S:Rb3=Rb4-2
 4030   for Ib1=1 to Ie1
 4040   Rb1=Rb2*(Rb2+Rb3)+1
 4045   Fs=Fs*Rb1
 4050   Fds=Fds+Rb2/Rb1
 4055   Rb2=Rb2*Q
 4060   next Ib1
 4070   Fds=Fs*(1+2*Rb4*Fds)
 4080   Fs=Fs*S
 4099   return
 4100   *Gsdeval(S,Bq,&Gs,&Gsd)
 4101   inc Id4
 4105   Gs=1:Gsd=0
 4110   Rb1=Bq:Rb4=Bq*Bq
 4115   Rb2=S*S-2
 4120   for Ib1=1 to Ie1
 4130   Rb3=Rb1*(Rb1+Rb2)+1
 4140   Gs=Gs*Rb3
 4145   Gsd=Gsd+Rb1/Rb3
 4150   Rb1=Rb1*Rb4
 4160   next Ib1
 4170   Gsd=2*Gs*S*Gsd
 4199   return
 4200   fnGetnq(Q)
 4210   Rb1=int(log(Eeps)/log(abs(Q)))+2
 4299   return(Rb1)
 4500   fnWp(Z,Q,W1)
 4501   inc Id5
 4525   Rb3=0.0:Rb1=1.0:Rb2=exp(2*Peye*Z*#i)
 4530   for Ib1=1 to Ie1
 4535   Rb1=Q*Rb1:Rb4=(1-Rb1*Rb2)^2:Rb5=(Rb1-Rb2)^2
 4540   Rb4=1.0/Rb4+1.0/Rb5:Rb5=Rb2*Rb4
 4545   Rb4=(1-Rb1)^2:Rb4=2.0/Rb4
 4550   Rb3=Rb3+Rb1*(Rb5-Rb4)
 4555   next Ib1
 4560   Rb4=(1-Rb2)^2:Rb3=1.0/12.0+Rb2/Rb4+Rb3
 4565   Rb5=2*Peye*#i/W1:Rb4=Rb5*Rb5*Rb3
 4599   return(Rb4)
 5000   fnEi(X)
 5015   if X>4 then goto 5100
 5020   Rb10=1
 5022   Rb8=log(Eeps):Rb7=log(X):Rb5=-Rb8
 5024   while (Rb10>0.1)
 5026   Rb9=log(Rb5):Rb2=Rb5*(Rb7-1.0)-(Rb5+1.5)*Rb9-Rb8
 5028   Rb3=Rb7-2.0-3.0/(Rb5+Rb5)-Rb9
 5030   Rb4=Rb5-Rb2/Rb3:Rb10=abs(Rb5-Rb4):Rb5=Rb4
 5032   wend
 5034   Ib2=int(Rb5)+5
 5036   if odd(Ib2)=0 then Ib2=Ib2+1
 5038   Rb6=0.0:Rb1=1.0
 5040   for Ib1=Ib2 to 1 step -1
 5042   Rb6=(Rb6+Rb1/Ib1)*X/Ib1
 5044   Rb1=-Rb1
 5046   next Ib1
 5050   Rb19=-Gamma-Rb7+Rb6
 5052   goto *Eiend
 5100   Ib1=int((log(Eeps)+10)/(-0.4))+5
 5112   Rb12=X
 5114   for Ib2=Ib1 to 1 step -1
 5116   Rb10=1.0+Ib2/Rb12:Rb11=X+Ib2/Rb10:Rb12=Rb11
 5118   next Ib2
 5120   Rb19=exp(-X)/Rb12
 5149   *Eiend:return(Rb19)
 6000   *Anvals(A(),Bad(),Na,&An%())
 6020   gosub *Form71(A(),&B())
 6022   Ic4=B(2):Ic5=2*B(4):Ic6=B(6)
 6024   An%(1)=1
 6025   for I=2 to Na
 6028   Ind=0
 6030   for J=1 to Bad(0,0)
 6032   if I=Bad(J,1) then Ind=J
 6034   next J
 6036   if Ind=0 then jump
 6038   An%(I)=Bad(Ind,5):goto *Anj1
 6040   **:Ind=fnRabmil(I)
 6042   if Ind=0 then goto 6050
 6044   An%(I)=fnApval(I)
 6046   goto *Anj1
 6050   gosub *Factor(I,&F())
 6052   if F(0)>1 then 6060
 6053   Ind=1
 6054   for J=1 to Bad(0,0):if F(1)=Bad(J,1) then Ind=0
 6055   next J
 6056   P%=F(1):I1%=I\P%:I2%=I1%\P%
 6057   An%(I)=An%(P%)*An%(I1%)-Ind*P%*An%(I2%)
 6058   goto *Anj1
 6060   P%=F(1)^F(2):I1%=I\P%
 6062   An%(I)=An%(P%)*An%(I1%)
 6070   *Anj1:next I
 6099   return
 6200   fnApval(P)
 6224   Ic2=0
 6226   for Ic1=0 to P-1
 6228   Ic3=((4*Ic1+Ic4)*Ic1+Ic5)*Ic1+Ic6
 6230   Ic2=Ic2+kro(Ic3,P)
 6232   next Ic1
 6249   return(-Ic2)
 6500   *Period(A(),&Perd())
 6510   ' perd(1)=w1  perd(2)=w2  perd(0)=Omega
 6530   gosub *Form71(A(),&B())
 6532   Rb1=B(6)/4.0:Rb3=B(4)/2.0:Rb2=B(2)/4.0
 6534   Rb4=Rb3/3.0-Rb2*Rb2/9.0
 6536   Rc9=(Rb3*Rb2-3.0*Rb1)/6.0-(Rb2^3)/27.0
 6538   if B(1)<0 then goto 6600
 6540   Rb5=-(Rb4*Rb4*Rb4+Rc9*Rc9):Rb6=sqrt(Rb5):Rb7=sqrt(Rc9*Rc9+Rb5)
 6542   if Rc9=0.0 then Rb8=Peye/2.0
 6544   :else Rb8=atan(Rb6/abs(Rc9))
 6546   if Rc9<0.0 then Rb8=Peye-Rb8
 6550   Rb9=Rb7^(1.0/3.0)
 6551   Rb17=cos(Rb8/3.0):Rb18=sin(Rb8/3.0)
 6552   Rb10=2.0*Rb9*Rb17-Rb2/3.0
 6554   Rb11=-Rb9*Rb17-Rb2/3.0-sqrt(3.0)*Rb9*Rb18
 6556   Rb12=-Rb9*Rb17-Rb2/3.0+sqrt(3.0)*Rb9*Rb18
 6558   if Rb11>Rb10 then swap Rb10,Rb11
 6560   if Rb11>=Rb12 then goto 6570
 6562   if Rb10>Rb12 then jump
 6564   swap Rb11,Rb12:swap Rb10,Rb11
 6566   goto 6570
 6568   **:swap Rb11,Rb12
 6570   Rb1=sqrt(Rb10-Rb11):Rb3=sqrt(Rb10-Rb12)
 6572   Perd(1)=Peye/fnAgm(Rb1,Rb3)
 6574   Rb1=sqrt(Rb11-Rb12)
 6576   Perd(2)=0.0+#i*(Peye/fnAgm(Rb1,Rb3))
 6578   Perd(0)=2.0*Perd(1)
 6580   goto *Perend
 6600   Rb5=Rb4*Rb4*Rb4+Rc9*Rc9
 6601   Rb16=Rc9+sqrt(Rb5)
 6602   if Rb16>=0.0 then Rb14=Rb16^(1.0/3.0)
 6603   if Rb16<0.0 then Rb14=-((-Rb16)^(1.0/3.0))
 6604   Rb16=Rc9-sqrt(Rb5)
 6606   if Rb16>=0.0 then Rb15=Rb16^(1.0/3.0)
 6608   if Rb16<0.0 then Rb15=-((-Rb16)^(1.0/3.0))
 6610   Rb10=Rb14+Rb15-Rb2/3.0
 6612   Rb11=3.0*Rb10+B(2)/4.0
 6614   Rb12=sqrt(3.0*Rb10*Rb10+B(2)*Rb10/2.0+B(4)/2.0)
 6616   Rb1=2.0*sqrt(Rb12):Rb3=sqrt(Rb12+Rb12+Rb11)
 6618   Perd(1)=2.0*Peye/fnAgm(Rb1,Rb3)
 6620   Rb3=sqrt(Rb12+Rb12-Rb11)
 6622   Perd(2)=Perd(1)/2.0+#i*Peye/fnAgm(Rb1,Rb3)
 6624   Perd(0)=Perd(1)
 6649   *Perend:return
 6700   *Cubsol(A,B,C,D,&Nrl,&R())
 6704   Rc1=B/A:Rc2=C/A:Rc3=D/A
 6706   Rc4=Rc2/3-Rc1*Rc1/9
 6708   Rc5=(Rc1*Rc2-3*Rc3)/6-Rc1^3/27
 6710   Rc6=Rc4*Rc4*Rc4+Rc5*Rc5
 6712   if Rc6<=0.0 then goto 6740
 6714   Rc7=Rc5+sqrt(Rc6)
 6716   if Rc7>=0 then Rc8=Rc7^(1.0/3.0)
 6718   if Rc7<0 then Rc8=-((-Rc7)^(1.0/3.0))
 6720   Rc7=Rc5-sqrt(Rc6)
 6722   if Rc7>=0 then Rc9=Rc7^(1.0/3.0)
 6724   if Rc7<0 then Rc9=-((-Rc7)^(1.0/3.0))
 6726   Nrl=1
 6728   R(1)=Rc8+Rc9-Rc1/3
 6730   Rc6=-(Rc8+Rc9)/2-Rc1/3
 6732   Rc7=sqrt(3.0)*(Rc8-Rc9)/2
 6734   R(2)=Rc6+Rc7*#i
 6736   R(3)=Rc6-Rc7*#i
 6738   goto *Cusolend
 6740   if Rc6<0.0 then goto 6754
 6742   Nrl=3
 6744   if Rc5>=0 then Rc6=Rc5^(1.0/3.0)
 6746   if Rc5<0 then Rc6=-((-Rc5)^(1.0/3.0))
 6748   R(1)=2*Rc6-Rc1/3
 6750   R(2)=-Rc6-Rc1/3:R(3)=R(2)
 6752   goto *Cusort
 6754   Nrl=3
 6756   Rc7=sqrt(-Rc6)
 6758   Rc8=sqrt(Rc5*Rc5-Rc6)
 6760   if Rc5=0.0 then Rc9=Peye/2
 6762   if Rc5<>0.0 then Rc9=atan(Rc7/Rc5)
 6764   if Rc5<0 then Rc9=Rc9+Peye
 6766   Rc6=Rc8^(1.0/3.0):Rc7=Rc9/3
 6768   Rc8=2*Rc6*cos(Rc7):Rc9=2*Rc6*sin(Rc7)
 6770   R(1)=Rc8-Rc1/3
 6772   Rc6=-Rc9*sqrt(3.0)/2.0
 6774   Rc7=-Rc8/2.0-Rc1/3.0
 6776   R(2)=Rc7+Rc6:R(3)=Rc7-Rc6
 6778   *Cusort:for Ic1=1 to 2:Rc1=R(Ic1):Ic3=Ic1
 6780   for Ic2=Ic1+1 to 3
 6782   if R(Ic2)<=Rc1 then jump
 6784   Ic3=Ic2:Rc1=R(Ic2)
 6786   **:next Ic2
 6788   if Ic3=Ic1 then jump
 6790   Rc2=R(Ic3):R(Ic3)=R(Ic1):R(Ic1)=Rc2
 6792   **:next Ic1
 6799   *Cusolend:return
 6900   fnAgm(A,B)
 6920   Rc2=A:Rc4=B:Rc5=Eeps+Eeps
 6925   while (Rc5>Eeps)
 6930   Rc1=(Rc2+Rc4)/2.0:Rc3=sqrt(Rc2*Rc4)
 6935   Rc5=abs((Rc1-Rc3)/Rc1)
 6940   Rc2=Rc1:Rc4=Rc3
 6945   wend
 6950   return(Rc2)
 7000   *Torspts(A(),&Ntor)
 7010   ' ntor=no. of points found
 7020   dim Bet(10),R(3),R1(3),Xpts(15)
 7030   Ntor=1
 7032   gosub *Form71(A(),&B())
 7034   Ib2=B(2):Ib3=2*B(4):Ib4=B(6)
 7036   gosub *Ratroot3(Ib2,Ib3,Ib4,&R())
 7038   if R(0)=0 then goto 7050
 7040   for Ib1=1 to R(0)
 7042   Ntor=Ntor+1:Ib23=-(A(1)*R(Ib1)+A(3))//2
 7044   Torpt(Ntor,1)=R(Ib1):Torpt(Ntor,2)=Ib23
 7046   next Ib1
 7050   Ib1=abs(B(1))
 7052   gosub *Factor(Ib1,&F())
 7054   for Ib1=1 to F(0):F(Ib1+Ib1)=F(Ib1+Ib1)\2:next Ib1
 7056   Ib5=0
 7058   for Ib1=1 to F(0)
 7060   if F(Ib1+Ib1)>Ib5 then Ib5=F(Ib1+Ib1)
 7062   next Ib1
 7064   Ib5=Ib5+1:Ib7=1:Ib6=0
 7066   for Ib1=1 to F(0)
 7068   Ib6=Ib6+F(Ib1+Ib1)*Ib7:Ib7=Ib5*Ib7
 7070   next Ib1
 7072   for Ib1=0 to Ib6
 7074   for Ib14=1 to F(0):Bet(Ib14)=0:next Ib14
 7076   Ib11=Ib1:Ib20=1
 7078   if Ib11=0 then jump
 7080   Ib12=Ib11\Ib5:Bet(Ib20)=res:inc Ib20:Ib11=Ib12
 7082   goto 7078
 7084   **:Ib13=1
 7086   for Ib14=1 to F(0)
 7088   if Bet(Ib14)<=F(Ib14+Ib14) then jump
 7090   cancel for:goto *Torj1
 7092   **:Ib13=Ib13*(F(Ib14+Ib14-1)^Bet(Ib14))
 7094   next Ib14
 7100   for Ib14=0 to 2:Ib11=Ib13*(2^Ib14)
 7102   Ib12=Ib4-Ib11*Ib11
 7104   gosub *Ratroot3(Ib2,Ib3,Ib12,&R1())
 7106   if R1(0)=0 then goto 7240
 7107   for Ib8=1 to R1(0)
 7108   Ib17=R1(Ib8):Ib20=0
 7109   if Ntor=1 then goto 7114
 7110   for Ib10=2 to Ntor:if Ib17=Xpts(Ib10) then Ib20=1
 7111   next Ib10
 7112   if Ib20=0 then goto 7114
 7113   goto 7220
 7114   Ib23=(Ib11-A(1)*Ib17-A(3))//2
 7115   Ib16=(3*Ib17*Ib17+2*A(2)*Ib17+A(4)-A(1)*Ib23)//(2*Ib23+A(1)*Ib17+A(3))
 7116   Ib18=-Ib17-Ib17-A(2)+Ib16*(Ib16+A(1))
 7118   Ib21=-Ib23-A(3)-A(1)*Ib18+Ib16*(Ib17-Ib18)
 7120   if R(0)=0 then goto 7140
 7122   Ib20=0
 7124   for Ib10=1 to R(0)
 7126   if Ib18=R(Ib10) then Ib20=1
 7128   next Ib10
 7130   if Ib20=0 then goto 7140
 7132   Ntor=Ntor+1:Xpts(Ntor)=Ib17
 7134   Torpt(Ntor,1)=Ib17:Torpt(Ntor,2)=Ib23
 7136   Ntor=Ntor+1:Ib23=-(Ib11+A(1)*Ib17+A(3))//2:Xpts(Ntor)=Ib17
 7138   Torpt(Ntor,1)=Ib17:Torpt(Ntor,2)=Ib23
 7139   goto 7220
 7140   for Ib10=3 to 6
 7142   Ib24=abs(Ib17-Ib18)+abs(Ib23-Ib21)
 7144   if Ib24>0 then goto 7154
 7146   Ib15=2*Ib23+A(1)*Ib17+A(3)
 7148   if Ib15=0 then goto 7190
 7150   Ib16=(3*Ib17*Ib17+2*A(2)*Ib17+A(4)-A(1)*Ib23)//Ib15
 7152   goto 7160
 7154   if Ib17=Ib18 then goto 7190
 7156   Ib16=(Ib23-Ib21)//(Ib17-Ib18)
 7160   Ib19=-Ib17-Ib18-A(2)+Ib16*(Ib16+A(1))
 7162   Ib22=-Ib23-A(3)-A(1)*Ib19+Ib16*(Ib17-Ib19)
 7164   if R(0)=0 then goto 7180
 7166   Ib20=0
 7168   for Ib9=1 to R(0)
 7170   if Ib19=R(Ib9) then Ib20=1
 7172   next Ib9
 7174   if Ib20=1 then goto 7190
 7180   if Ib19=Ib18 then goto 7190
 7182   Ib18=Ib19:Ib21=Ib22
 7184   goto 7200
 7190   Ntor=Ntor+1:Xpts(Ntor)=Ib17
 7192   Torpt(Ntor,1)=Ib17:Torpt(Ntor,2)=Ib23
 7194   Ntor=Ntor+1:Ib23=-(Ib11+A(1)*Ib17+A(3))//2:Xpts(Ntor)=Ib17
 7196   Torpt(Ntor,1)=Ib17:Torpt(Ntor,2)=Ib23
 7198   cancel for:goto 7220
 7200   next Ib10
 7220   next Ib8
 7240   next Ib14
 7250   *Torj1:next Ib1
 7299   return
 7500   *Ratroot3(A,B,C,&Rts())
 7505   'rational roots of 4x^3+ax^2+bx+c=0
 7511   dim Rr(3)
 7512   if C<>0 then goto 7550
 7514   if B<>0 then goto 7530
 7516   Rts(0)=3:Rts(1)=0:Rts(2)=0:Rts(3)=-A//4
 7518   goto *Rt3end
 7530   Rts(0)=1:Rts(1)=0:Ic8=A*A-16*B
 7532   if Ic8<0 then goto *Rt3end
 7534   Ic6=isqrt(Ic8):Ic7=Ic8-Ic6*Ic6
 7536   if Ic7>0 then goto *Rt3end
 7538   Rts(0)=3:Rts(2)=(-A+Ic6)//8:Rts(3)=(-A-Ic6)//8
 7540   goto *Rt3end
 7550   gosub *Cubsol(4,A,B,C,&Ic1,&Rr())
 7552   Ic6=0
 7554   for Ic2=1 to Ic1
 7556   Rc1=4*Rr(Ic2):Ic3=round(Rc1):Rc2=abs(Rc1-Ic3)
 7558   if Rc2>Rteps then goto 7568
 7559   Ic3=Ic3//4
 7560   Ic4=((4*Ic3+A)*Ic3+B)*Ic3+C
 7562   if Ic4<>0 then goto 7568
 7564   inc Ic6
 7566   Rts(Ic6)=Ic3
 7568   next Ic2
 7570   Rts(0)=Ic6
 7599   *Rt3end:return
 7700   *Rt3ch(A,B,C,M,N,&Ind)
 7705   ' checks if (m/n) is a root of 4x^3+ax^2+bx+c=0
 7715   Ind=0
 7720   Id1=C\M:Id2=res
 7722   if Id2>0 then goto *Rt3chend
 7724   Id3=B+N*Id1:Id1=Id3\M:Id2=res
 7726   if Id2>0 then goto *Rt3chend
 7728   Id3=A+N*Id1:Id1=Id3\M:Id2=res
 7730   if Id2>0 then goto *Rt3chend
 7732   Id3=4+N*Id1
 7734   if Id3<>0 then goto *Rt3chend
 7736   Ind=1
 7749   *Rt3chend:return
 7800   *Reduce(Ain(),&Am(),&Baddat())
 7820   ' baddat(0,0)=no. of bad primes  baddat(0,1)=conductor
 7825   ' baddat(i,1)=prime  baddat(i,2)=values of cp
 7830   ' baddat(i,3)=type   baddat(i,4)=numerical value for some types
 7835   ' baddat(i,0)=power of prime in delta  baddat(i,5)=split/non-split/add
 7840   local Del,J,Nbadp,I
 7845   dim V(4),A(6),D(5)
 7850   for I=1 to 6:A(I)=Ain(I):next I
 7855   gosub *Minimal(A(),&Am())
 7860   for I=1 to 6:A(I)=Am(I):next I
 7865   gosub *Form71(A(),&B())
 7870   Del=B(1)
 7875   if Del<0 then Del=-Del
 7880   gosub *Factor(Del,&F())
 7885   Nbadp=F(0):Baddat(0,0)=Nbadp:Ncon=1
 7890   for I=1 to Nbadp:Baddat(I,1)=F(I+I-1):Baddat(I,0)=F(I+I):next I
 7895   for I=1 to Nbadp
 7900   if Baddat(I,1)>3 jump
 7905   gosub *Red23(A(),&D(),Baddat(I,1))
 7910   goto 7920
 7915   **:gosub *Redp(A(),&D(),Baddat(I,1))
 7920   Ncon=Ncon*Baddat(I,1)^D(3)
 7925   Baddat(I,2)=D(4):Baddat(I,3)=D(1):Baddat(I,4)=D(2)
 7930   Baddat(I,5)=D(5)
 7935   next I
 7940   Baddat(0,1)=Ncon
 7999   return
 8000   *Minimal(A(),&Am())
 8015   gosub *Form71(A(),&B())
 8020   C4=B(3):C6=B(5):Ic7=C6*C6:Ic8=B(1):Ic5=1
 8025   Ic9=gcd(Ic7,Ic8)
 8030   Ic8=6*Ic9:gosub *Factor(Ic8,&F())
 8035   for Ic1=1 to F(0)
 8037   Ic10=F(Ic1+Ic1-1):Ic6=F(Ic1+Ic1)
 8039   if Ic10=2 then Ic6=Ic6-1
 8040   if Ic10=3 then Ic6=Ic6-1
 8042   Ic6=Ic6\12:Ic5=Ic5*(Ic10^Ic6)
 8045   if Ic10>2 then jump
 8050   Ic3=C4\(2^(4*Ic6)):Ic4=C6\(2^(6*Ic6))
 8055   Ic7=fnMods(Ic3,2):Ic8=fnMods(Ic4,4)
 8060   if Ic7<>1 then goto 8075
 8065   if Ic8<>-1 then goto 8075
 8070   Ic2=-(Ic4+1)\4:goto 8135
 8075   Ic7=fnMods(Ic3,16):Ic8=fnMods(Ic4,16)
 8080   if Ic7<>0 then goto 8100
 8085   if Ic8=0 then goto 8095
 8088   if Ic8=8 then goto 8095
 8090   goto 8100
 8095   Ic2=Ic4\8:goto 8135
 8100   Ic5=Ic5\2:Ic2=0:goto 8135
 8105   **:if Ic10>3 then jump
 8110   Ic7=0:Ic8=C6:if Ic8=0 then goto 8125
 8115   Ic9=Ic8\3:if res>0 then goto 8125
 8120   inc Ic7:Ic8=Ic9:goto 8115
 8125   Ic8=Ic7-(6*Ic6+2)
 8130   if Ic8=0 then Ic5=Ic5\3
 8135   **:next Ic1
 8140   C4=C4\Ic5^4:C6=C6\Ic5^6
 8145   Am(1)=fnMods(C4,2)
 8150   Ic7=Am(1)+C6:Am(2)=-fnMods(Ic7,3)
 8155   Ic7=Ic2+Am(1)*Am(2):Am(3)=fnMods(Ic7,2)
 8160   Ic7=Am(1)+4*Am(2):Ic8=Ic7*Ic7-C4-24*Am(1)*Am(3)
 8165   Am(4)=Ic8\48
 8170   Ic8=36*Ic7*(Am(1)*Am(3)+2*Am(4))-Ic7^3-C6-216*Am(3)
 8175   Am(6)=Ic8\864:Am(5)=0
 8199   return
 8300   *Redp(A(),&D(),P)
 8325   gosub *Form71(A(),&B())
 8330   Del=B(1):C4=B(3):J=C4^3//Del:C6=B(5):D(5)=0
 8335   Ic2=abs(num(J)):Ic3=abs(den(J))
 8337   Ic4=0
 8338   if J=0 then jump
 8340   Ic4=fnOrdp(Ic2,P)
 8350   if Ic4>0 then jump
 8360   Ic4=-fnOrdp(Ic3,P)
 8365   **:Ic7=abs(Del)
 8370   Ic5=fnOrdp(Ic7,P)
 8380   Ic6=Ic5
 8385   if Ic4<0 then Ic6=Ic6+Ic4
 8390   ' Ic7=1:Ic8=0:Ic9=0:Ic10=0
 8395   ' if Ic6<12 then jump
 8400   Ic7=P^(Ic6\12)
 8405   Ic2=odd(A(1))
 8410   if Ic2=1 then Ic9=(Ic7-A(1))\2
 8415   if Ic2=0 then Ic9=-A(1)\2
 8420   Ic2=A(2)-Ic9*A(1)-Ic9*Ic9
 8425   Ic3=Ic2\3:Ic11=res
 8430   if Ic11=0 then Ic8=-Ic2\3
 8435   if Ic11=1 then Ic8=(Ic7*Ic7-Ic2)\3
 8440   if Ic11=2 then Ic8=(-Ic7*Ic7-Ic2)\3
 8445   Ic2=A(3)+Ic8*A(1)
 8450   Ic3=odd(Ic2)
 8455   if Ic3=1 then Ic10=(Ic7^3-Ic2)\2
 8460   if Ic3=0 then Ic10=-Ic2\2
 8465   Ic2=Ic6\12:Ic6=res:Del=Del\(Ic7^12)
 8470   C4=C4\Ic7^4:C6=C6\Ic7^6
 8475   **:if Ic4>=0 then goto 8545
 8480   Ic2=-Ic4
 8485   if Ic6=6 then goto 8515
 8490   D(3)=1:Ic3=kro(-C6,P):D(5)=Ic3
 8495   if Ic3=1 then D(4)=Ic5
 8500   if Ic3=-1 then D(4)=gcd(2,Ic5)
 8505   D(1)=2:D(2)=Ic2
 8510   goto *Redpend
 8515   D(3)=2:D(1)=3:D(2)=Ic2
 8520   Ic3=odd(Ic2)
 8525   if Ic3=1 then Ic5=Del*C6\P^(9+Ic2)
 8530   if Ic3=0 then Ic5=Del\P^(6+Ic2)
 8535   D(4)=3+kro(Ic5,P)
 8540   goto *Redpend
 8545   D(3)=2
 8550   if Ic6=0 then D(3)=0
 8555   if Ic6=0 then D(1)=1:D(4)=1
 8556   if Ic6=2 then D(1)=4:D(4)=1
 8557   if Ic6=3 then D(1)=5:D(4)=2
 8558   if Ic6=4 then D(1)=6
 8559   if Ic6=6 then D(1)=7
 8560   if Ic6=8 then D(1)=8
 8561   if Ic6=9 then D(1)=9:D(4)=2
 8562   if Ic6=10 then D(1)=10:D(4)=1
 8565   D(2)=0
 8570   if Ic6<>4 goto 8575
 8572   Ic2=-6*C6\(P^2):D(4)=2+kro(Ic2,P)
 8575   if Ic6<>8 goto 8580
 8577   Ic2=-6*C6\(P^4):D(4)=2+kro(Ic2,P)
 8580   if Ic6<>6 goto *Redpend
 8582   Ic2=3*C4\(P^2):Ic3=C6\(P^3)
 8583   Ic5=0
 8584   for Ic1=0 to P-1
 8585   Ic6=4*Ic1*Ic1*Ic1-Ic2*Ic1-Ic3
 8586   Ic4=Ic6\P
 8587   if res=0 then Ic5=Ic5+1
 8588   next Ic1
 8589   D(4)=1+Ic5
 8599   *Redpend:return
 8600   *Red23(A(),&D(),P)
 8610   local Bt,Nu,Fp,Cp,M
 8615   dim Ad(6),Ad1(6),V(4),Rt(3,2)
 8620   for Ic4=1 to 6:Ad(Ic4)=A(Ic4):next Ic4
 8624   gosub *Form71(Ad(),&B())
 8626   Ic5=abs(B(1)):Nu=-1:Ic11=0:Ic18=0:D(5)=0
 8628   while Ic11=0
 8630   Ic6=Ic5\P:Ic11=res:Nu=Nu+1:Ic5=Ic6
 8632   wend
 8634   if Nu>0 then jump
 8636   Bt=1:Fp=0:Cp=1:D(5)=0:goto *Red23end
 8638   **:if P=3 then goto 8646
 8640   Ic5=B(2)\P:if res>0 then jump
 8642   Ic1=fnMods(Ad(4),P):Ic3=fnMods(Ic1*(1+Ad(2)+Ad(4))+Ad(6),P):goto 8656
 8644   **:Ic1=fnMods(Ad(3),P):Ic3=fnMods(Ic1+Ad(4),P):goto 8656
 8646   Ic5=B(2)\P:if res>0 then jump
 8648   Ic1=-fnMods(B(6),P):goto 8654
 8652   **:Ic1=-fnMods((B(2)*B(4)),P)
 8654   Ic5=Ad(1)*Ic1+Ad(3):Ic3=fnMods(Ic5,P)
 8656   V(1)=1:V(2)=Ic1:V(3)=0:V(4)=Ic3
 8657   gosub *Form72(Ad(),V(),&Ad1())
 8658   for Ic4=1 to 6:Ad(Ic4)=Ad1(Ic4):next Ic4
 8660   gosub *Form71(Ad(),&B())
 8662   Ic5=B(3)\P:if res=0 then jump
 8664   V(1)=1:V(2)=Ad(1):V(3)=-Ad(2):gosub *Roots2(V(),P,&Rt())
 8666   if Rt(0,0)>0 then Cp=Nu:D(5)=1:goto 8670
 8668   Ic5=Nu\2:Cp=2-res:D(5)=-1
 8670   Bt=2:Fp=1:Ic18=Nu:goto *Red23end
 8672   **:Ic5=Ad(6)\(P*P):if res=0 then jump
 8674   Bt=4:Fp=Nu:Cp=1:D(5)=0:goto *Red23end
 8676   **:Ic5=B(8)\(P^3):if res=0 then jump
 8678   Bt=5:Fp=Nu-1:Cp=2:D(5)=0:goto *Red23end
 8680   **:Ic5=B(6)\(P*P*P):if res=0 then jump
 8682   Bt=6:Fp=Nu-2:V(1)=1:V(2)=Ad(3)\P:V(3)=-Ad(6)\P^2
 8684   gosub *Roots2(V(),P,&Rt())
 8686   Cp=1:if Rt(0,0)>0 then Cp=3
 8688   goto *Red23end
 8690   **:if P=3 then jump
 8692   Ic2=fnMods(Ad(2),2):Ic3=2*fnMods((Ad(6)\4),2):goto 8696
 8694   **:Ic2=-Ad(1)*modinv(2,P):Ic3=-Ad(3)*modinv(2,P)
 8696   V(1)=1:V(2)=0:V(3)=Ic2:V(4)=Ic3:gosub *Form72(Ad(),V(),&Ad1())
 8698   for Ic4=1 to 6:Ad(Ic4)=Ad1(Ic4):next Ic4
 8700   gosub *Form71(Ad(),&B())
 8702   Ic7=Ad(2)\P:Ic8=Ad(4)\(P*P):Ic9=Ad(6)\(P^3)
 8704   Ic10=27*Ic9*Ic9-Ic7^2*Ic8^2+4*Ic7^3*Ic9-18*Ic7*Ic8*Ic9+4*Ic8^3
 8706   Ic11=3*Ic8-Ic7^2
 8708   Ic5=Ic10\P:if res=0 then goto 8716
 8710   Bt=7:Fp=Nu-4:V(1)=1:V(2)=Ic7:V(3)=Ic8:V(4)=Ic9
 8712   gosub *Roots3(V(),P,&Rt())
 8714   Cp=1+Rt(0,1):D(5)=0:goto *Red23end
 8716   Ic5=Ic11\P:if res=0 then goto 8775
 8718   if P=2 then Ic1=Ic8
 8720   if P=3 then Ic1=Ic7*Ic8
 8722   Ic1=P*fnMods(Ic1,P)
 8724   V(1)=1:V(2)=Ic1:V(3)=0:V(4)=0:gosub *Form72(Ad(),V(),&Ad1())
 8726   for Ic4=1 to 6:Ad(Ic4)=Ad1(Ic4):next Ic4
 8728   gosub *Form71(Ad(),&B())
 8730   M=1:Ic12=P*P:Ic13=P*P:Cp=0
 8732   while Cp=0
 8734   Ic14=Ad(2)\P:Ic15=Ad(3)\Ic13:Ic16=Ad(4)\(P*Ic12):Ic17=Ad(6)\(Ic12*Ic13)
 8736   Ic5=Ic15*Ic15+4*Ic17:Ic6=Ic5\P
 8738   if res=0 then jump
 8740   V(1)=1:V(2)=Ic15:V(3)=-Ic17
 8742   gosub *Roots2(V(),P,&Rt())
 8744   Cp=2:if Rt(0,0)>0 then Cp=4
 8746   goto 8772
 8748   **:if P=2 then Ic3=Ic13*Ic17
 8750   if P=3 then Ic3=Ic13*fnMods((-Ic15*modinv(2,P)),P)
 8752   V(1)=1:V(2)=0:V(3)=0:V(4)=Ic3
 8754   gosub *Form72(Ad(),V(),&Ad1())
 8755   for Ic4=1 to 6:Ad(Ic4)=Ad1(Ic4):next Ic4
 8756   gosub *Form71(Ad(),&B())
 8758   Ic13=Ic13*P:inc M:Ic14=Ad(2)\P:Ic15=Ad(3)\Ic13
 8759   Ic16=Ad(4)\(P*Ic12):Ic17=Ad(6)\(Ic12*Ic13)
 8760   Ic5=Ic16*Ic16-4*Ic14*Ic17:Ic6=Ic5\P
 8761   if res=0 then jump
 8762   V(1)=Ic14:V(2)=Ic16:V(3)=Ic17
 8763   gosub *Roots2(V(),P,&Rt())
 8764   Cp=2:if Rt(0,0)>0 then Cp=4
 8765   goto 8772
 8766   **:if P=2 then Ic1=Ic12*fnMods((Ic17*Ic14),2)
 8767   if P=3 then Ic1=Ic12*fnMods((-Ic16*modinv(2*Ic14,P)),P)
 8768   V(1)=1:V(2)=Ic1:V(3)=0:V(4)=0:gosub *Form72(Ad(),V(),&Ad1())
 8769   for Ic4=1 to 6:Ad(Ic4)=Ad1(Ic4):next Ic4
 8770   gosub *Form71(Ad(),&B())
 8771   Ic12=Ic12*P:inc M
 8772   wend
 8773   Bt=3:Fp=Nu-M-4:Ic18=M:goto *Red23end
 8775   if P=2 then Ic5=-Ic7
 8776   if P=3 then Ic5=-Ic9
 8777   Ic1=P*fnMods(Ic5,P):V(1)=1:V(2)=Ic1:V(3)=0:V(4)=0
 8778   gosub *Form72(Ad(),V(),&Ad1())
 8779   for Ic4=1 to 6:Ad(Ic4)=Ad1(Ic4):next Ic4
 8780   gosub *Form71(Ad(),&B())
 8781   Ic15=Ad(3)\(P*P):Ic17=Ad(6)\(P^4)
 8782   Ic5=Ic15*Ic15+4*Ic17:Ic6=Ic5\P
 8783   if res=0 then jump
 8784   V(1)=1:V(2)=Ic15:V(3)=-Ic17
 8785   gosub *Roots2(V(),P,&Rt())
 8786   Cp=1:if Rt(0,0)>0 then Cp=3
 8787   Bt=8:Fp=Nu-6:goto *Red23end
 8788   **:if P=2 then Ic3=Ic17
 8789   if P=3 then Ic3=2*Ic15
 8790   Ic3=-P*P*fnMods(Ic3,P)
 8791   V(1)=1:V(2)=0:V(3)=0:V(4)=Ic3:gosub *Form72(Ad(),V(),&Ad1())
 8792   for Ic4=1 to 6:Ad(Ic4)=Ad1(Ic4):next Ic4
 8793   Ic5=Ad(4)\P^4:if res=0 then jump
 8794   Bt=9:Fp=Nu-7:Cp=2:goto *Red23end
 8795   **:Bt=10:Fp=Nu-8:Cp=1
 8798   *Red23end:D(1)=Bt:D(2)=Ic18:D(3)=Fp:D(4)=Cp
 8799   return
 8800   *Form71(A(),&B())
 8810   ' b(2)=b2  b(4)=b4  b(6)=b6  b(8)=b8
 8815   ' b(1)=delta  b(3)=c4    b(5)=c6
 8820   B(2)=A(1)*A(1)+4*A(2):B(4)=A(1)*A(3)+2*A(4)
 8825   B(6)=A(3)*A(3)+4*A(6)
 8830   B(8)=A(1)*A(1)*A(6)+4*A(2)*A(6)-A(1)*A(3)*A(4)
 8835   B(8)=B(8)+A(2)*A(3)*A(3)-A(4)*A(4)
 8840   B(3)=B(2)*B(2)-24*B(4)
 8845   B(5)=-B(2)^3+36*B(2)*B(4)-216*B(6)
 8850   B(1)=-B(2)*B(2)*B(8)-8*B(4)^3-27*B(6)*B(6)+9*B(2)*B(4)*B(6)
 8899   return
 8900   *Form72(A(),V(),&Ad())
 8905   ' implements formulas 7.2 from Cohen
 8930   Ad(1)=(A(1)+V(3)+V(3))\V(1)
 8935   Ad(2)=(A(2)-V(3)*A(1)+3*V(2)-V(3)*V(3))\(V(1)*V(1))
 8940   Ad(3)=(A(3)+V(2)*A(1)+2*V(4))\V(1)^3
 8945   Ad(4)=A(4)-V(3)*A(3)+2*V(2)*A(2)-(V(4)+V(2)*V(3))*A(1)
 8947   Ad(4)=(Ad(4)+3*V(2)*V(2)-2*V(3)*V(4))\V(1)^4
 8950   Ad(6)=A(6)+V(2)*A(4)+V(2)*V(2)*A(2)+V(2)^3-V(4)*A(3)
 8955   Ad(6)=(Ad(6)-V(4)*V(4)-V(2)*V(4)*A(1))\V(1)^6
 8999   return
 9000   fnRabmil(N)
 9010   ' Result=0 Composite     Result=1 Prime
 9020   Ie2=even(N):Ie1=1
 9025   if Ie2=0 jump
 9030   Ie1=0:if N=2 then Ie1=1
 9032   goto *Finrm
 9035   **:if N=1 goto *Finrm
 9040   Ie3=N-1:Ie4=0:Ie2=0
 9045   while Ie2=0
 9050   Ie6=Ie3\2:Ie2=res
 9055   Ie3=Ie6:inc Ie4
 9060   wend
 9065   Ie3=Ie3+Ie3+1:dec Ie4
 9070   randomize
 9075   for Ie7=1 to 20
 9080   Ie6=int((N-2)*rnd+2):Ie2=0
 9085   Ie5=modpow(Ie6,Ie3,N)
 9090   if Ie5=1 jump
 9095   Ie8=modpow(Ie5,1,N)
 9100   while and{Ie8>1,Ie8<N-1,Ie2<=Ie4-2}
 9105   Ie5=modpow(Ie5,2,N):inc Ie2
 9107   Ie8=modpow(Ie5,1,N)
 9110   wend
 9115   if Ie5=N-1 jump
 9117   Ie1=0
 9120   cancel for
 9125   goto *Finrm
 9130   **:next Ie7
 9140   Ie1=1
 9150   *Finrm:return(Ie1)
 9200   *Factor(N,&Fac())
 9210   ' fac(0)=no. of factors
 9215   '  n = fac(1)^fac(2)*fac(3)^fac(4)*.....
 9220   if N>1 goto 9240
 9225   Fac(0)=1:Fac(1)=1:Fac(2)=1
 9230   goto *Finfac
 9240   Id6=fnRabmil(N)
 9245   if Id6=0 goto 9260
 9250   Fac(0)=1:Fac(1)=N:Fac(2)=1
 9255   goto *Finfac
 9260   Id5=0
 9265   for Id1=1 to 10000000
 9267   if Id1>12251 goto 9273
 9270   Id2=prm(Id1)
 9271   goto 9275
 9273   Id2=Id2+2
 9275   Id3=N\Id2:Id4=res
 9280   if Id4>0 goto *Nxfac2
 9285   Id5=Id5+1:Fac(0)=Id5
 9290   Fac(Id5+Id5-1)=Id2:Fac(Id5+Id5)=1
 9295   N=Id3
 9300   Id3=N\Id2:Id4=res
 9305   if Id4>0 goto *Nxfac1
 9310   Fac(Id5+Id5)=Fac(Id5+Id5)+1:N=Id3
 9315   if N>1 goto 9300
 9320   cancel for
 9325   goto *Finfac
 9330   *Nxfac1:Id6=fnRabmil(N)
 9332   if Id6=0 goto *Nxfac2
 9334   Id5=Id5+1:Fac(0)=Id5
 9336   Fac(Id5+Id5-1)=N:Fac(Id5+Id5)=1
 9337   cancel for
 9338   goto *Finfac
 9340   *Nxfac2:next Id1
 9349   *Finfac:return
 9350   fnOrdp(N,P)
 9360   Id1=N:Id2=0
 9365   Id3=Id1\P
 9370   if res>0 then jump
 9375   Id2=Id2+1:Id1=Id3
 9380   goto 9365
 9399   **:return(Id2)
 9400   *Roots2(C(),P,&R())
 9402   ' finds roots of c(1)x^2+c(2)x+c(3) mod p
 9404   ' r(0,0)=no of roots  r(0,1)=no. of distinct roots
 9406   ' r(i,0)=root no i  r(i,1)=degree of root
 9413   Id5=0:Id6=0
 9415   for Id3=0 to P-1
 9420   Id1=C(1)*Id3*Id3+C(2)*Id3+C(3)
 9422   Id4=Id1\P:Id2=res
 9425   if Id2<>0 then jump
 9430   Id5=Id5+1:Id6=Id6+1:R(Id6,0)=Id3:R(Id6,1)=1
 9435   Id1=2*C(1)*Id3+C(2)
 9440   Id4=Id1\P:Id2=res
 9450   if Id2<>0 then jump
 9455   R(Id6,1)=2:Id5=Id5+1
 9460   **:next Id3
 9465   R(0,0)=Id5:R(0,1)=Id6
 9470   return
 9500   *Roots3(C(),P,&R())
 9502   ' finds roots of c(1)x^3+c(2)x^2+c(3)x+c(4)  mod p
 9504   ' r(0,0)=no of roots  r(0,1)=no. of distinct roots
 9506   ' r(i,0)=root no i  r(i,1)=degree of root
 9513   Id5=0:Id6=0
 9515   for Id3=0 to P-1
 9520   Id1=C(1)*Id3*Id3*Id3+C(2)*Id3*Id3+C(3)*Id3+C(4)
 9522   Id4=Id1\P:Id2=res
 9525   if Id2<>0 then jump
 9530   Id5=Id5+1:Id6=Id6+1:R(Id6,0)=Id3:R(Id6,1)=1
 9535   Id1=3*C(1)*Id3*Id3+2*C(2)*Id3+C(3)
 9540   Id4=Id1\P:Id2=res
 9550   if Id2<>0 then jump
 9555   R(Id6,1)=2:Id5=Id5+1
 9560   Id1=6*C(1)*Id3+2*C(2)
 9565   Id4=Id1\P:Id2=res
 9570   if Id2<>0 then jump
 9575   R(Id6,1)=3:Id5=Id5+1
 9580   **:next Id3
 9585   R(0,0)=Id5:R(0,1)=Id6
 9590   return
 9600   *Printeq(A())
 9605   ' prints out equation
 9610   print "  y^2 ";
 9612   if A(1)=0 then jump
 9614   Ib1=abs(A(1))
 9616   if A(1)>0 then print " + ";Ib1;" xy ";
 9618   if A(1)<0 then print " - ";Ib1;" xy ";
 9620   **:if A(3)=0 then jump
 9622   Ib1=abs(A(3))
 9624   if A(3)>0 then print " + ";Ib1;" y ";
 9626   if A(3)<0 then print " - ";Ib1;" y ";
 9630   **:print " = x^3 ";
 9632   if A(2)=0 then jump
 9634   Ib1=abs(A(2))
 9636   if A(2)>0 then print " + ";Ib1;" x^2 ";
 9638   if A(2)<0 then print " - ";Ib1;" x^2 ";
 9640   **:if A(4)=0 then jump
 9642   Ib1=abs(A(4))
 9644   if A(4)>0 then print " + ";Ib1;" x ";
 9646   if A(4)<0 then print " - ";Ib1;" x ";
 9650   **:if A(6)=0 then jump
 9652   Ib1=abs(A(6))
 9654   if A(6)>0 then print " + ";Ib1;
 9656   if A(6)<0 then print " - ";Ib1;
 9660   **:print " "
 9662   print " "
 9699   return
 9700   *Printred(Bad())
 9720   print 
 9725   print "Prime","Kodaira","Reduction","Cp"
 9728   for Ib1=1 to Bad(0,0)
 9730   print Bad(Ib1,1),
 9735   Ib2=Bad(Ib1,3)
 9740   if Ib2=1 then print "I0",
 9741   if Ib2=2 then print "I";Bad(Ib1,4),
 9742   if Ib2=3 then print "I*";Bad(Ib1,4),
 9743   if Ib2=4 then print "II",
 9744   if Ib2=5 then print "III",
 9745   if Ib2=6 then print "IV",
 9746   if Ib2=7 then print "I0*",
 9747   if Ib2=8 then print "IV*",
 9748   if Ib2=9 then print "III*",
 9749   if Ib2=10 then print "II*",
 9750   if Bad(Ib1,5)=-1 then print "Non-split",
 9751   if Bad(Ib1,5)=0 then print "Additive",
 9752   if Bad(Ib1,5)=1 then print "Split",
 9760   print Bad(Ib1,2),
 9769   print
 9770   next Ib1
 9775   print
 9799   return
 9800   fnMods(N,P)
 9805   ' finds symmetric residue mod p
 9815   Id1=N\P:Id2=res
 9820   Id1=P\2
 9825   if Id2>Id1 then Id2=Id2-P
 9830   return(Id2)
 9900   *Urst(A(),Am(),&U())
 9910   dim B(8)
 9915   gosub *Form71(A(),&B())
 9920   Rb1=B(3):gosub *Form71(Am(),&B())
 9925   Rb2=B(3)
 9927   if Rb2<>0 then goto 9930
 9928   U(1)=1:goto 9935
 9930   Rb3=Rb1\Rb2:Rb4=isqrt(Rb3):U(1)=isqrt(Rb4)
 9935   U(3)=(U(1)*Am(1)-A(1))\2
 9940   U(2)=(U(1)*U(1)*Am(2)-A(2)+U(3)*A(1)+U(3)*U(3))\3
 9945   U(4)=(U(1)^3*Am(3)-A(3)-U(2)*A(1))\2
 9949   return