/* Deduction of minimal polynomials for coordinates of Heegner points */

realtofrac(x) =
{
   local(L,M);
   L = contfrac(x);
   M = contfracpnqn(L);
   return(M[1,1]/M[2,1]);
}

heegpoly_fieldelt(x,D) =
{
   /* this program assumes D is a negative fundamental discriminant and gives output in the form (a + b omega) where omega is canonical integral basis of the ring of integers */
   local(p,q,a,b,t);
   p = realtofrac(real(x));
   q = realtofrac(imag(x)^2);
   t = -4*q/D;
   b = sign(imag(x))*sqrtint(numerator(t))/sqrtint(denominator(t));
   if((b^2 != t),
      print([p,q,t,b]);
      error("something went wrong when finding element of imaginary quadratic extension")
   );
   if((D % 4) == 0, a = p, a = p - b/2);
   return(a+b*quadgen(D));
}

heegpoly(P,D) =
{
   local(d,j,PP);
   d = poldegree(P);
   PP = sum(j=0,d,heegpoly_fieldelt(polcoeff(P,j),D) * x^j);
   return(PP);
}

ellHeegnerPolys(E,D,s=0) =
{
   local(pts,Fx,Fy,j);
   if(s==0,pts = ellHeegnerPoints(E,D),pts = ellHeegnerPoints(E,D,s));
   Fx = prod(j=1,length(pts),(X-pts[j][1]));
   Fy = prod(j=1,length(pts),(X-pts[j][2]));
   Fx = heegpoly(Fx,coredisc(D));
   Fy = heegpoly(Fy,coredisc(D));
   return([Fx,Fy]);
}