Re: Difficulty with binary quadratic equation solver

[Bill Allombert <Bill.Allombert@math.u-bordeaux1.fr>]

On Fri, Nov 28, 2008 at 09:36:47AM -0700, Kurt Foster wrote:
> Bill provides the example
> 
> [a,b,c,d,e,f] = [1,0,-5,1,1,2]

(I'd like to note that http://www.alpertron.com.ar/QUAD.HTM
does not handle that equation either, so we cannot compare the results).

> for which the field discriminant N = 5, the index m is 2, U =  
> (B.fu)^2, the smallest value of k is 3, and plowing through the  
> calculations by hand shows that there are two "inequivalent" solutions
> 
> [x,y] = [1,1] and [-5, -2]
> 
> For a given positive fundamental discriminant N and index m > 1,  
> determining the exact smallest value of k for which U^k is in O has  
> proven to be rather exasperating, although certain multiples of it are  
> easy to find.  I've dug into it a fair bit, but I'm sure this is a  
> known calculation.  Can someone provide a reference?

Well, we can use the ray class number formula that give the ray class number
in term of the class number:

   h_f=h*Phi(f)/phi(f)/k

So the following should compute k:

  K=bnfinit(x^2-b*x+a*c);
  B=bnrinit(K,K.index);
  k=K.no*B.bid.no/(eulerphi(K.index)*B.no);

I join a new version of queq.gp that use this formula. 

Cheers,
Bill.

qeeval(Q,v)=local(a=Q[1],b=Q[2],c=Q[3],d=Q[4],e=Q[5],f=Q[6],x=v[1],y=v[2]);a*x^2+b*x*y+c*y^2+d*x+e*y+f

qered(Q)=
{
  local(a=Q[1],b=Q[2],c=Q[3],d=Q[4],e=Q[5],f=Q[6],F,u,v);
  F=[2*a,b;b,2*c]^-1*[-d,-e]~;u=F[1];v=F[2];
  D=denominator(content(F));
  [D^2*(u^2*a+v*u*b+v^2*c+u*d+v*e+f),D,D*u,D*v]
}

qesolve(Q)=
{
  local(a,b,c,f,V,K,B,d,L,w,U,Ud,S=List(),idx);
  local(R,D,M,u,v);
  a=Q[1];b=Q[2];c=Q[3];
  if (issquare(b^2-4*a*c),error("reducible equation"));
  R=qered(Q); f=R[1]; D=R[2];
  K=bnfinit(x^2-b*x+a*c);
  L=bnfisintnorm(K,-f*a);
  if (#K.fu, 
    idx=K.index; B=bnrinit(K,idx);
    d=K.no*B.bid.no/(eulerphi(idx)*B.no);
    U=K.fu[1]; if (norm(U)==-1,U=U^2);
   ,d=1; U=1);
  for(i=1,#L,
    for(j=0,K.tu[1]-1,
      for(k=0,d-1,
        w=U^k*L[i]*K.tu[2]^j;
        w=lift(w); u=polcoeff(w,0); v=polcoeff(w,1);
        if (denominator(u)!=1 || denominator(v)!=1,next);
        if (u%a!=0,next,u\=a);
        if ((R[3]-u)%D!=0,next,u=(R[3]-u)\D);
        if ((R[4]-v)%D!=0,next,v=(R[4]-v)\D);
        listput(S,[u,v]);
  )));Vec(S)
}