\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ ss.gp -- pari-gp v2 program                           \\
\\                                                       \\
\\ Date: June 9, 1999                                    \\
\\                                                       \\
\\ Author: William A. Stein (was@math.berkeley.edu)      \\
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

len(v)=matsize(v)[2];
nrows(v)=matsize(v)[1];

{composite_finitefield(p, f, g,
    v)=
\\print("comp p=",p,", f=", f,", g=",g);
   if(poldegree(f)<=1,return(g));
   if(poldegree(g)<=1,return(f));
   v=polcompositum(f*Mod(1,p),g*Mod(1,p));
\\   if(len(v)>1, print("v = ",v));
   return(lift(v[1]));
}

\\ a root of g(x) in Fp[x]/(f(x)). 
{rootin_finitefield(p, f, g,
    fac)=
    fac=factorff(g, p, subst(f,x,t));
\\    print("fac = ",fac);
    if(poldegree(fac[1,1])!=1,
       print("ERROR"));
    -subst(fac[1,1],x,0);
}

\\ sum n_E/(j-j_E) \in Fpbar(j).
{fdiv(p,i,sq,
   ans, f, alp, bet, ratfun1, ratfun2, fac, m, r, n_E, j_E)=

   fac = lift(factormod(sseigen[p,i][1], p));
   ans = vector(nrows(fac),m,0);
   for(r=1,nrows(fac),
     if(fac[r,2]!=1,print("WARNING, fac = ",fac));
     f = composite_finitefield(p, fac[r,1], ssbasis[p][1]);
     alp = rootin_finitefield(p, f, fac[r,1]);
     bet = rootin_finitefield(p, f, ssbasis[p][1]);
     ratfun1 = 0;   \\ this will be sum n_E/(j-j_E).
     ratfun2 = 0;
     for(m=1,len(ssbasis[p][2]),
       n_E = subst(sseigen[p,i][2][m],x,alp)*Mod(1,p);
       j_E = subst(ssbasis[p][2][m],x,bet)*Mod(1,p);
       ratfun1 += n_E/(x-j_E);
       if(sq==1, ratfun2 += n_E/((x-j_E)^2));
     );
     if(sq==1,
       ans[r] = [subst(f,x,t), lift(lift(ratfun1)),
                             lift(lift(ratfun2))],
       ans[r] = [subst(f,x,t), lift(lift(ratfun1))];
     );
   );
   return(ans);
}

\\ Compute the Fp-rational zeros of the rational function
\\ g(x) in Fp[t]/f(t). 
{zeros(p, f, g,
   fac,c,v,i)=
   ans=[];
   g=numerator(g);
   for(j=0,p-1,
     if(Mod(subst(g,x,j),f)*Mod(1,p)==0,
        ans=setunion(ans,[j]);
     )
   );
   return(vecsort(eval(ans)));
}

\\ fdivzeros
{fdivzeros(p,i,
   v)=
   v=fdiv(p,i);
   for(i=1,len(v),
     v[i]=[v[i][1],zeros(p,v[i][1],v[i][2])];
   );
   return(v);
}

\\ fdivzeros_table
{fdivzeros_table(pstart, pstop,
   p, i, v)=
   forprime(p=pstart, pstop,
     i=1;
     while(type(sseigen[p,i])=="t_VEC",
        print("rf[",p,",",i,"] = ",fdivzeros(p,i),";");
        i++;
     )
   );
}

\\ fdiv2zeros
{fdiv2zeros(p,i,
   v)=
   v=fdiv(p,i);
   for(i=1,len(v),
     v[i]=[v[i][1],zeros(p,v[i][1],v[i][2])];
   );
   return(v);
}

{fdiv2zeros_table(pstart, pstop,
   p, i, v)=
   forprime(p=pstart, pstop,
     i=1;
     while(type(sseigen[p,i])=="t_VEC",
        print("rf[",p,",",i,"] = [",fdiv2zeros(p,i));
        i++;
     )
   );
}


\\ apEj: returns a list of the set of _all_ pairs [j, a_p(E_j)] where 
\\ j lies in F_p and E_j is some elliptic curve of j-invariant
\\ j. If several curves has the same j-invariant but different
\\ a_p they will appear on the list.  Note, a_p := p+1 - #E(Fp). 
\\ The algorithm is as stupid as imaginable, but who cares 
\\ since p < 1000. 
{jap(p,   
  ans, a, b, E, cn)=
\\  print1 ("ja, p=",p,": [");
  ans=[];
  cn=floor(p/40)+1;
  for(a=0,p-1,
\\    if(a%cn==0,print1("."));
    for(b=0,p-1,
      if(poldisc((x^3+a*x+b)*Mod(1,p))!=0,
\\         print(poldisc((x^3+a*x+b)*Mod(1,p)));
         E=ellinit([0,0,0,Mod(a,p),Mod(b,p)]);
         ans = setunion(ans, [[lift(E.j), lift(ellap(E,p))]]);
      );
    );
  );
  return(vecsort(eval(ans),1));
}

\\ Returns all j's for which there exists SOME elliptic curve
\\ over E/Fp of j-invariant j for which ap=1, i.e.,
\\ |E(Fp)| = p+1-ap = p..
\\ v = output of jap.
{jwithapone(v,
  i, ans)=
  ans=[];
  for(i=1,len(v),if(v[i][2]==1,ans=setunion(ans,[v[i][1]])));
  return(vecsort(eval(ans)));
}

{jwithapone_table(pstart, pstop,
  p)=
  japs=vector(pstop);
  forprime(p=pstart,pstop,
    japs[p]=jap(p);
    print("jap1[",p,"] = ",jwithapone(japs[p]),";")
  );
}



{rflist(p,
   v, i, j)=
   i=1;
   v=[];
   while(type(rf[p,i])=="t_VEC",
      for(j=1,len(rf[p,i]),
         v=setunion(v,rf[p,i][j][2]));
      i++;
   );
   return(eval(v));
}

{rf_0jap1(pstart,pstop,
   p, v, w, z)=
   forprime(p=pstart,pstop,
     v=Set(rflist(p));  
     w=Set(jap1[p]);
     z=eval(setintersect(v,w));
     if(len(z)!=0,
       print(p, ":  ",v,", ", w,", ", z)
     );
   );