\\ ---------------  GP code  ---------------------------------------
\\
\\ Time-stamp: <Thu, Dec 19, 2002 - 10:57:58 - apacetti>
\\
\\ Description: Routines for operating with binary positive definite
\\ quadratic forms.
\\
\\ File: /home/villegas/gp/bforms.gp
\\
\\ Original Author:     Fernando Rodriguez-Villegas 
\\                      villegas@math.utexas.edu
\\                      University of Texas at Austin
\\
\\
\\ With the assistance of:  Ariel Pacetti
\\			    apacetti@math.utexas.edu
\\			    University of Texas at Austin
\\
\\ Created:             Wed Feb  19  2000
\\
\\-----------------------------------------------------------------
\\ Elementary routine to determine if a rational number is integer

isint(n)=(denominator(n)==1)

\\====================================================================
\\computes the power of the prime in the number

primepower(p,a)=
	{local(cont,aux);
	cont=0;
	aux=a;
	while(aux%p==0,aux=aux/p;cont=cont+1);
	cont
}

\\====================================================================
\\ Computes the square root of a number module another number

sqrtmod(a,b)=
	{local(aux,ans,n,fac,aux2,pow,a2);
	ans=Mod(1,1);
	if(b<0,b=-b);
	fac=factor(b);
	n=length(fac[,1]);
	for(k=1,n,if(fac[k,2]==1,aux=truncate(sqrt(a+O(fac[k,1]))),
		if(fac[k,1]<>2,
			pow=primepower(fac[k,1],a);
			if(pow>=fac[k,2],aux=0,
			if(pow%2==0,,error("no solution"));
			a2=a/(fac[k,1]^pow);
			aux=truncate(sqrt(a2+O(fac[k,1])));
			for(j=1,fac[k,2]-1-pow,
			aux2=(bezout(2*aux,fac[k,1])[1]*(a2-sqr(aux))/(fac[k,1])^j)%fac[k,1];
			aux=(aux+aux2*fac[k,1]^j)%(fac[k,1]^fac[k,2]));
			aux=aux*fac[k,1]^(pow/2)),
		if(fac[k,2]<=3,aux=sqrtmodbrute(a,2^(fac[k,2])),
			aux=sqrtmodbrute(a,8);
			for(j=1,fac[k,2]-3,
			if(aux%2==0,if((sqr(aux)-a)/2^(j+2)%2==1,error("not a
			square")),
				aux2=-(sqr(aux)-a)/2^(j+2);
			aux=aux+aux2*2^(j+1))))));
	ans=chinese(ans,Mod(aux,fac[k,1]^(fac[k,2]))));
	lift(ans)
}

\\======================================================================
\\ Brute force

sqrtmodbrute(a,b)=
	{local(aux1,aux2);
	aux1=abs(b);
	aux2=a%aux1;
	for(k=0,aux1-1,if((k^2%aux1)==aux2,return(k)));
	error("not a square")
}



\\====================================================================
\\ evalfm(q,v)
\\ Evaluate quadratic form q at vector v

evalfm(q,v)=q[1]*v[1]^2+q[2]*v[1]*v[2]+q[3]*v[2]^2


\\====================================================================
\\ primefm(a,d) returns a form [a,b,c] of discriminant d.  If d is
\\ not a square mod a or a is not prime or 1 then gives error.

primefm(a,d)=
	{local(b,aux);
	if(a==1,
	  if(sqr(d%4)-(d%4)==0,[1,d%2,(d%2^2-d)\4],error("no solution")),
            if(isprime(a)&&(sqr(d%4)-(d%4))%4==0,
		if(d%a==0,if(d%4==0,return([a,0,-d/4/a]),
			return([a,a,(a^2-d)/4/a])),);
		if(kronecker(d,a)==1,
			aux=sqrt(d+O(a)); 
			b=truncate(aux);
			if(b%2!=d%2,b=b+a,);
			[a,b,(sqr(b)-d)\4\a],
        error("Error in primefm: no such form or not a prime input"))))
}
\\======================================================================
\\ Computes the same but for any a if exists

fm(a,d)=
	{local(b,c);
	if(d%4==2,error("error no solution"));
	if(d%4==3,error("error no solution"));
	b=sqrtmod(d,4*a);
	if(a<>0,c=(b^2-d)/4/a;
	[a,b,c],if(d<0,error("No
	solution"),if(d==sqrtint(d)^2,[0,sqrtint(d),1],error("no solution"))))

}


\\=====================================================================
\\ discrf(q)  discriminant of the binary form q 
\\ q=[a,b,c] and d= 4*a*c-b*b.

discrf(q)=sqr(q[2])-4*q[1]*q[3]

\\=======================================================================
\\ compf(q1,q2)= compositum of q1 and q2.
\\ discriminants must be the same.

compf(q1,q2)=
      {local(aux,mumu,a,b,c,dd,nunu);
      dd=-discrf(q1);
  if(discrf(q2)!=-dd, error(" *** Error in compf: not the same discr"),
      aux=bezout(q1[1],q2[1]);
      mumu=aux[2];
      aux=bezout(aux[3],(q1[2]+q2[2])\2);
      mumu=mumu*aux[1];nunu=aux[2];
      a=q2[1]\aux[3];
      b=q2[2]+(mumu*(q1[2]-q2[2])\2-nunu*q2[3])*a*2;
      a=a*q1[1]\aux[3];
      b=b%(2*a);
      c=(sqr(b)+dd)\(4*a);
      [a,b,c];
     )
}

\\======================================================================
\\ redf(q) returns a vector with [1]= q_0 = reduced q and [2]= a vector
\\ [u,v] so that q_0(u,v)=q[1]

redf(q)=
	{local(aa,a,b,c,dd,top,bot,n,aux,aux1,s);
	dd=-discrf(q);
	a=q[1];b=q[2];c=q[3];
	aa= 2*a;
	top=[-b,1];
	bot=[aa,0];
	n=(a+b)\aa;
	if(n,
	    aux=a*n;
	    c=c+(aux-b)*n;
	    b=b-2*aux;
	    top=top+n*bot);
	while(c<a,
               aux=c;c=a;
               a=aux;aa=2*a;
               b=-b;
               aux1=top;
               top=-bot;
               bot=aux1;
               n=(a+b)\aa;
               if(n,
                    aux=a*n;
                    c=c+(aux-b)*n;
                    b=b-2*aux;
                   top=top+n*bot;
                 ,);
           );
     if(a==c && b>0,
               b=-b;
               aux1=top;
               top=-bot;
               bot=aux1;
        ,);
            q=[a,b,c];
            s=kronecker(bot[1],dd);
            if(s==-1,bot=-bot,);
	    bot[1]=(bot[1]-q[2]*bot[2])/2/q[1];
            [q,bot]
}


\\=====================================================================
\\ compredf(q1,q2) gives the reduced compositum of the forms.
\\ error if discr are not the same.

compredf(q1,q2)=
        {local(q);
              q= compf(q1,q2);
                 red1f(q)
}

\\======================================================================
\\ powredf(q,n) gives the reduced from equiv to q^n.
\\ uses the  binary decomposition of n and reduces at every step.

powredf(q,n)=
     {local(aux,qq,k,nn);
     aux=binary(n);
     nn=2; qq=q;
     k=1; if(n==1,return(redf(q)[1]));
       while(nn<n,nn=2*nn;
                  k=k+1;
                  qq=compredf(qq,qq);
                    if(aux[k]==1,qq=compredf(q,qq);
                       ,);
             );
        qq
 }
        
\\=======================================================================
\\ ordf(q,n)= order of q in its class group. q^n must be principal,
\\ searches over divisors of n.

ordf(q,n)=
       {local(k,a,aux);
                  aux=divisors(n);
                  k=1;
                  a= powredf(q,aux[k])[1];
              while(a!=1 ,k=k+1;
                          a=powredf(q,aux[k])[1];
                    );
                  aux[k]
}

\\======================================================================
\\ powerf(q,n) gives  q^n.
\\ uses the  binary decomposition of n. does not reduce at any time.

powerf(q,n)=
     {local(aux,qq,k,nn);
     if(n,
     aux=binary(n);
     nn=2; qq=q;
     k=1;
       while(nn<=n,nn=nn+nn;
                  k=k+1;
                  qq=compf(qq,qq);
                    if(aux[k]==1,qq=compf(q,qq);
                       ,);
             );
        qq,primefm(1,discrf(q)))
}

\\=======================================================================
\\ powlistf(q,n) 
\\ list all reduced forms q^k for k=1 to n.

powlistf(q,n)=
          {local(qq);
              qq=red1f(q);
              for(k=1,n,
                    print(k,"   ", qq);
                    qq=compredf(q,qq);
                  );

}

\\=======================================================================
\\ red1f(q) returns  reduced q.

red1f(q)=
         {local(aa,a,b,c,n,aux);
         a=q[1];b=q[2];c=q[3];
         aa= a+a;
         n=(a+b)\aa;
          if(n,
               aux=a*n;
               c=c+(aux-b)*n;
               b=b-2*aux;
             ,);
     while(c<a,
               aux=c;c=a;
               a=aux;aa=2*a;
               b=-b;
               n=(a+b)\aa;
               if(n,
                    aux=a*n;
                    c=c+(aux-b)*n;
                    b=b-2*aux;
                 ,);
           );
       if(a==c && b>0,
            b=-b;
          ,);
          [a,b,c]
}

\\=====================================================================
\\ genf(dd,hh) returns a generator for the class group of discr dd and
\\ order hh when such group is cyclic

genf(dd,hh)=
    {local(k,q,sq,aux,a,nn,b,n,p,aux1);
	  dd=-dd;
          aux=divisors(hh);
          sq=sqrt(dd);
          n=0;k=1;
          q=[1,-1,(1+dd)\4];
            while(k!=hh,
                         n=n+1;
                         p=prime(n);
                         while(kronecker(-dd,p)!=1,
                                            n=n+1;
                                            p=prime(n);
                               );
                          b=-dd+O(p^(hh+2));
                          aux1=sqrt(b);
                          b=truncate(aux1);
                            if(b%2==0,b=truncate(-aux1),);
                          aux1=(b*b+dd)\4;
                          nn=2;k=aux[nn];
                          a=p^k;
                          q=[a,b,aux1\a];
                          while(red1f(q)[1]!=1,
                                             nn=nn+1;
                                             k=aux[nn];
                                             a=p^k;
                                             q=[a,b,aux1\a];
                                );
                    );
                red1f([p,b,aux1\p])
}

\\=======================================================================
\\ powvecf(q,n)
\\ vector with all reduced forms q^k for k=1 to n.

powvecf(q,n)=
          {local(qq,vecf);
              vecf=[];
              qq=red1f(q);
              for(k=1,n,
                    vecf= concat(vecf,[qq]);
                    qq=compredf(q,qq);
                  );
              vecf
}

\\==========================================================================
\\ indf(q,vecf) returns k such that vecf[k]=red1f(q) if it exists, otherwise 
\\ returns zero. 

indf(q,vecf)=
      {local(testf,qq,k);
          k=0;
          testf=0;
          qq=red1f(q);
          until(testf!=0 || k==length(vecf),
                     k=k+1;
                       if(vecf[k]==qq,testf=k,);
                );
           testf
}
    

\\========================================================================
\\ hvecf(dd)
\\ Vector of reduced forms.

hvecf(dd)=
    {local(f,be,dv);
	dd=-dd;
	if(dd%4==0||dd%4==3,
      f=[];
      forstep(x=dd%2,sqrtint(dd\3),2,
        be=(sqr(x)+dd)\4;
                for(a=max(x,1),sqrtint(be),
                          dv=divrem(be,a);
                            if(dv[2]==0,
                                  if(gcd(gcd(a,x),dv[1])==1,
                                           f=concat(f,[[a,-x,dv[1]]]);
                                             if(x!=0 && a!=x && a*a!=be,
                                           f=concat(f,[[a,x,dv[1]]]))))));
     f,
	print("  ***   Error in hvecf: not a discriminant"))
}

\\========================================================================
\\ clnof(dd)
\\ Class number of primitive forms of dicsr d.

clnof(dd)=
    {local(be,dv,ct);
	dd=-dd;
	if(dd%4==0||dd%4==3,
      forstep(x=dd%2,sqrtint(dd\3),2,
        be=(sqr(x)+dd)\4;
                for(a=max(x,1),sqrtint(be),
                          dv=divrem(be,a);
                            if(dv[2]==0,
                                  if(gcd(gcd(a,x),dv[1])==1,
                                           ct=ct+1;
                                             if(x!=0 && a!=x && a*a!=be,
                                           ct=ct+1)))));
     ct,
	print("  ***   Error in clnof: not a discriminant"))
}


\\========================================================================
\\ cprf(q,plist)
\\ Compute a form equivalent to q whose norm a is relatively
\\ prime to the primes in plist; if plist is not given then
\\ plist=[2,3,p_1,...,p_r] is assumed, where p_1,...,p_r are the prime
\\ factors of the discriminant of q. This implementation actually
\\ gives forms with prime norm with this property.


cprf(q,plist)=
	{local(d,p,qq,qqred,a,b);
	if(q[1]==1,q,
		d=discrf(q); q=red1f(q);
		if(plist,,
			plist=rad(d);
			if(d%2,plist=concat([2],plist),);
			if(d%3,plist=concat([3],plist),));
			plist=vecsort(plist);
		p=nextprime(q[1]+1);
		until(kronecker(d,p)==1&&isnotin(p,plist),
			p=nextprime(p+1));
			qq=fm(p,d);
			qqred=red1f(qq);
			a=qqred[1];b=qqred[2];
			while(a!=q[1]||abs(b)!=abs(q[2]),
				p=nextprime(p+1);
				until(kronecker(d,p)==1&&isnotin(p,plist),
					p=nextprime(p+1));
					qq=fm(p,d);
					qqred=red1f(qq);
					a=qqred[1];b=qqred[2]);
			[qq[1],qq[2]*if(b==q[2],1,-1),qq[3]])
}

\\======================================================================
\\ Given a quadratic form q, and an integer N, computes a form
\\ equivalent to q, but such that the first coordinate is prime to N

cprf2(q,N)=
	{local(ans);
	if(gcd(q[1],N)==1,return(q));
	ans=q;while(gcd(ans[3],N)>1,ans=[ans[1],ans[2]+2*ans[1],ans[1]+ans[2]+ans[3]]);
	[ans[3],-ans[2],ans[1]]}


\\========================================================================
\\ Test for not-membership of p in plist; plist is a list of distinct
\\ numbers.

\\isnotin(p,plist)=plist=vecsort(plist);1-length(setintersect([p],plist))

isnotin(p,plist)=
	{local(ans);
	ans=1;
	for(k=1,length(plist),ans=ans&&(plist[k]!=p));
	ans
}

\\========================================================================
\\ rad(m)
\\ Get radical of an integer as a plist


rad(m)=
	{local(f);
	f=factor(m);
	vector(length(f~),r,f[r,1])
}

\\========================================================================
\\ cprvecf(d,plist)
\\ Compute a list of inequivalent forms [a,b,c] of discriminant d
\\ with a relatively prime to all primes in the list plist.


cprvecf(d,plist)=
	{local(v);
	v=hvecf(d);
	for(k=1,length(v),v[k]=cprf(v[k],plist));
	v
}

\\========================================================================
\\ orderf(q) gives the order in the equivalence relation of the form q.


orderf(q)=
	{local(cont,qq,aux);
	cont=1;
	qq=red1f(q);
	aux=compredf(qq,qq);
	while (qq != aux, aux=compredf(qq,aux); cont=cont+1);
	cont
}



addhelp(sqrtmod,"Computes the square root of a number module another number. For more references see the bforms.txt file");

addhelp(evalfm,"Evaluate quadratic form q at vector v. For more references see the bforms.txt file");

addhelp(primefm,"primefm(a,d) returns a form [a,b,c] of discriminant d.  If d is not a square mod a or a is not prime or 1 then gives error. For more references see the bforms.txt file");

addhelp(fm,"fm(a,d) returns a form [a,b,c] of discriminant d for any d, if it exists. For more references see the bforms.txt file")

addhelp(discrf,"discrf(q) discriminant of the binary form q if q=[a,b,c]  then b^2-4*a*c. For more references see the bforms.txt file")

addhelp(compf,"compf(q1,q2) computes the compositum of the forms q1 and q2. For doing this, the discriminants must be the same.")

addhelp(redf,"redf(q) returns a vector with [1]= q_0 = reduced q and [2]= a vector [u,v] so that q_0(u,v)=q[1]")

addhelp(compredf,"compredf(q1,q2) gives the reduced compositum of the forms q1 and q2. Answers error if discr are not the same.")

addhelp(powredf,"powredf(q,n) gives the reduced from equiv to q^n. It uses the  binary decomposition of n and reduces at every step.")

addhelp(orf,"ordf(q,n) gives the  order of q in its class group. q^n must be principal, it just searches over the divisors of n.")

addhelp(powerf,"powerf(q,n) gives  q^n. It uses the  binary decomposition of n. Note that it does not reduce at any time.")

addhelp(powlistf,"powlistf(q,n) list all reduced forms q^k for k=1 to n.")

addhelp(red1f,"red1f(q) returns  just the reduced form of q.")

addhelp(genf,"genf(dd,hh) returns a generator for the class group of discr dd and order hh when such group is cyclic.")

addhelp(powvecf,"powvecf(q,n)gives  vector with all reduced forms q^k for k=1 to n.")

addhelp(indf,"indf(q,vecf) returns k such that vecf[k]=red1f(q) if it exists, otherwise returns zero.")

addhelp(hvecf,"hvecf(d) gives a vector of all reduced forms of discriminant d")

addhelp(clnof,"clnof(d) gives the class number of primitive forms of dicsr d. Remember that d must be negative")

addhelp(cprf,"cprf(q,plist) compute a form equivalent to q whose norm a is relatively prime to the primes in plist, if plist is not given then plist=[2,3,p_1,...,p_r] is assumed, where p_1,...,p_r are the prime factors of the discriminant of q. This implementation actually gives forms with prime norm with this property.")

addhelp(cprvecf,"cprvecf(d,plist) computes a list of inequivalent forms [a,b,c] of discriminant d with a relatively prime to all primes in the list plist.");

addhelp(orderf,"orderf(q) gives the order in the equivalence relation of the form q.")

{helpbforms()=
print("-------------------------------------------------------------------");
print("      		Computational number Theory                        "); 
print("-------------------------------------------------------------------");
print("");
print("Description: Routines for operating with binary positive definite");
print("quadratic forms.");
print("");
print("Original Author:     Fernando Rodriguez-Villegas");
print("                     villegas@math.utexas.edu");
print("                     University of Texas at Austin");
print("");
print("With the assistance of:  Ariel Pacetti");
print("                         apacetti@math.utexas.edu");
print("                         University of Texas at Austin");
print("-------------------------------------------------------------------");
print("");
print("List of routines'" );
print("");
print("sqrtmod   evalfm    primefm    fm    discrf");
print("compf    redf    compredf    powredf    ordf");
print("powerf    powlistf    red1f    genf    powvecf");
print("indf    hvecf    clnof    cprf    cprvecf");
print("orderf");
print("")}