/* subroutines */


/* a tuples of Fq  */

Range(a)={
local(F0,i0,ia0,ia1,RA0,RA1,ra0);
RA0=vector(q,i0,[i0-1]);
for(i0=1,a-1,ra0=length(RA0);
RA1=vector(q*ra0,i0,
ia0=(i0-1) % q;
ia1=(i0-1) \ q;
concat(RA0[1+ia1],ia0)
);
RA0=RA1
);
return(RA0);
}

/* deg <=a polys over Fq  */

PolRange(a)={
local(RA0,RA1,i0);
RA0=Range(a+1);RA1=[];
for(i0=1,length(RA0),
RA1=concat(RA1,Polrev(RA0[i0],t));
);
return(RA1);
}




/* Polynomial Legendre Symbol */

LegPol(PA,QA)={
local(cA,iA);
cA=poldegree(QA);
return(lift(Mod(PA,QA)^((q^cA-1)/2)));
}






/* List of deg=a irred poly neq P&Q*/


IrredNIN(a)={
local(RA0,RA1);
RA0=Range(a);RA1=[];
for(i=1,q^a,
g=Mod(t,q)^a+Mod(Polrev(RA0[i],t),q);
if((polisirreducible(g)==1 && g<>P && g<>Q),RA1=concat(RA1,g))
);
return(RA1);
}


/* valuation inf=0*/
ValPol(PA,QA)={
local(i0);
i0=0;
if(PA==Mod(0,q),return(0));
while(PA%QA==Mod(0,q),i0=i0+1;PA=PA/QA);
return(i0);
}


/*  main   */

{



Y=t;
q=3;P=t^2+1;Q=t^2+t+2;
q=3;P=t^3+2*t+1;Q=t+2;
q=3;P=t^4+t+2;Q=t^4+2*t^3+1*t^2+2*t+1;
q=3;P=t^5+2*t+1;Q=t+2;

q=5;P=t^3+t+1;Q=t+2;
q=5;P=t^4+2;Q=t^2+t+1;

q=3;P=t^2+1;Q=t^2+t+2;
q=3;P=t^3+2*t+1;Q=t+2;
q=3;P=t^5+2*t+1;Q=t+2;
q=5;P=t^3+t+1;Q=t+2;
q=5;P=t^4+2;Q=t^2+t+1;
q=3;P=t^4 + t^3 + 2*t + 1;Q=t^2 + 1;

q=3;P=t^3 + t^2 + t + 2;Q=t + 1;
q=5;P=t^3 + t^2 + 4*t + 1;Q=t+2;
q=5;P=t^3 + 3*t + 2;Q=t+4;

q=3;P=t^4+t+2;Q=t^4+2*t^3+1*t^2+2*t+1;
q=3;P=t^3 + t^2 + t + 2;Q=t + 1;
q=3;P=t^4 + t^3 + 2*t + 1;Q=t^2 + 1;
q=3;P=t^5+2*t+1;Q=t+2;
q=5;P=t^3 + t^2 + 4*t + 1;Q=t+2;
q=5;P=t^4+2;Q=t^2+t+1;
q=7;P=t^3 + 2;Q=t+3;

q=3;P=t^6 + t^2 + t + 1;Q=t^2 + 1;
q=3;P=t^6 + t^2 + t + 1;Q=t^2 + t + 2;
q=3;P=t^6 + t^2 + t + 1;Q=t^4 + t^3 + t^2 + 1;

q=5;P=t^4+2;Q=t^4 + 3*t^2 + t + 1;

q=7;P=t^4 + t + 1;Q=t^2 + 1;

q=11;P=t^3 + t + 4;Q=t + 2;
q=11;P=t^3 + t + 4;Q=t^3 + t^2 + 2;

q=3;P=t^3 + t^2 + t + 2;Q=t + 1;
q=3;P=t^4 + t^3 + 2*t + 1;Q=t^2 + 1;

q=3;P=t^3 + t^2 + t + 2;Q=t + 1;
q=3;P=t^4 + t^3 + 2*t + 1;Q=t^2 + 1;
q=3;P=t^5+2*t+1;Q=t+2;
q=5;P=t^3 + t^2 + 4*t + 1;Q=t+2;
q=5;P=t^4+2;Q=t^2+t+1;
q=7;P=t^3 + 2;Q=t+3;



r=poldegree(P);s=poldegree(Q);
L=2*(r+s)-2;

e=znprimroot(q);

P=Mod(P,q);Q=Mod(Q,q);Y=Mod(Y,q);


print("(q,P,Q,Y)=("q","lift(P)","lift(Q)","lift(Y)")");


/* validity check  */

if(polisirreducible(P)<>1 || polisirreducible(Q)<>1 || polisirreducible(Y)<>1, print("invalid (P,Q)");
return); 

fr=0;
if((LegPol(Mod(Y,q),Q)<>Mod(1,q)) && (LegPol(Mod(Y,q),P)<>Mod(1,q)), fr=1);
if((LegPol(e*Mod(Y,q),Q)<>Mod(1,q)) && (LegPol(e*Mod(Y,q),P)<>Mod(1,q)), fr=1);

if(((LegPol(P,Q)==Mod(1,q)) || ((poldegree(P)%2)==0)) && (LegPol(e*P,Q)==Mod(1,q)), fr=1);
if(((LegPol(Q,P)==Mod(1,q)) || ((poldegree(Q)%2)==0)) && (LegPol(e*Q,P)==Mod(1,q)), fr=1);
if(fr==1,print("(P,Q,Y) doesn't satisfy conditions");return);


/* easy check  */


M=r+s-1; 
G=PolRange(M);  /* elts of A/(PQ)  */

for(m=0,M-1,
F=PolRange(m);


k0=0; /* k0 is k at which G[k] fails the test */
for(k=1,length(G),
if(Mod(Mod(G[k],q),P)==0, next());
if(Mod(Mod(G[k],q),Q)==0, next());

j0=0; /* j0 is j such that a=F[j] is good for G[k] */
for(j=1,length(F),
g=Mod(F[j],q)^2-Mod(G[k],q);

LgP=LegPol(g,P);
LgQ=LegPol(g,Q);

if(LgP==Mod(-1,q) && LgQ==Mod(-1,q),j0=j;break(1));
);
if(j0==0,k0=k;break(1));  /* the case where G[k] fails the test   */
);
if(k0==0,L=2*m;break(1)); /* the case where m passes the test  */
);

print("degree >= "L+1" checked"); 



RCR=[];RRS=[];
for(l=1,L,
m=if(l%2==0,l/2,(l-1)/2);
F=PolRange(m);
R=IrredNIN(l);
print("irred poly of degree "l" computed");


/* check the conditions */

for(k=1,length(R),


j0=0;
for(j=1,length(F),
g=Mod(F[j],q)^2-4*R[k];
LgP=LegPol(g,P);
LgQ=LegPol(g,Q);
h=Mod(F[j],q)^2-4*e*R[k];
LhP=LegPol(h,P);
LhQ=LegPol(h,Q);

if((poldegree(g)%2==1 || pollead(g)^((q-1)/2)<>Mod(1,q)) &&
(LgP==Mod(-1,q) || (LgP==Mod(0,q) && ValPol(g,P)%2==1)) && 
(LgQ==Mod(-1,q) || (LgQ==Mod(0,q) && ValPol(g,Q)%2==1)),j0=j;break(1));
if((poldegree(h)%2==1 || pollead(h)^((q-1)/2)<>Mod(1,q)) &&
(LhP==Mod(-1,q) || (LhP==Mod(0,q) && ValPol(h,P)%2==1)) && 
(LhQ==Mod(-1,q) || (LhQ==Mod(0,q) && ValPol(h,Q)%2==1)),j0=j;break(1));


);
if(j0==0,RCR=concat(RCR,R[k]);print(R[k]);break(2))
);
print("degree "l" checked")
);

}