\\ This file contains PARI routines for working with (real) quadratic
\\ forms and computing the Shintani cycles associated to them. 


\\ Routine to set up some initial gloabal data - must be run on the 
\\ discrimainant of a real quadratic order before any other computation can
\\ be carried out. 
setup(d0,f,u,d)=\
 gc=f;\
 gd0=d0;\
 gd=d0*f^2;\
 gs=sqrt(gd);\
 gsp=factmod(x^2-gd,11);\
 gsp=compo(gsp,1);\
 gsp=compo(gsp,1);\
 gsp=compo(gsp,1);\
 gsp=compo(gsp,2);\
 u=unit(gd)^2;\
  print(u);\
 gu=compo(u,2);gv=compo(u,3);\
 if(gd%4==1,gu=gu+gv/2; gv=gv/2;gid=try(1,1),\
   gv=gv/2;gid=try(0,1));\
 print(gu^2-gd*gv^2);\

\\ List manipulation routines: last extracts the last element from a 
\\ list v, and rest returns the rest of the list. 
last(v)=compo(v,length(v));
rest(v,n,res)=\
  res=[];\
  n=length(v)-1;\
  for(j=1,n,res=concat(res,[compo(v,j)]));\
  res
  

\\ Modifies the quadratic form f by the transformation (1,k;0,1).
shft(f,k,a,b,c,bb)=\
 a=f[1]; b=f[2]; c=f[3];\
 bb = b+2*a*k; c = a*k^2 + b*k+c;\
 b=bb;\
 [a,b,c]

\\ Returns 1 if f is reduced, 0 otherwise. 
reduced(f,a,b,c,res)=\
 a=f[1];b=f[2];c=f[3];\
 res=1;\
 if(gs-b<0,res=0,);\
 if(gs-b>2*abs(a),res=0,);\
 if(gs+b<2*abs(a),res=0,);\
 res

\\ reduces the form f with positive discriminant. 
reduce(f)=\
 until(reduced(f),f=next(f));\
 f

\\ returns the form adjacent to f in the cycle of reduced forms.
next(f,a,b,c,k)=\
 a=f[3];b=-f[2];c=f[1];\
 k=trunc( (-b+gs)/(2*a) );\
 shft([a,b,c],k)

\\ Returns 1 if f1 is equivalent to f2 in the class group, 0 otherswise. 
eqf(f1,f2,res,ff)=\
 f1=reduce(f1);\
 ff=f2;\
 until(eqqf(ff,f1) || eqqf(ff,f2), ff=next(ff));\
 if(eqqf(ff,f1),res=1,res=0);\
 res

\\ Returns 1 if f1 and f2 are equal, 0 otherwise. 
eqqf(f1,f2,res)=\
 res=1;\
 if(f1[1]<>f2[1],res=0,);\
 if(f1[2]<>f2[2],res=0,);\
 if(f1[3]<>f2[3],res=0,);\
 res

\\ The folllowing three routines impelement the Guassian composition of
\\ binary quadrartic forms, as explained in Duncan Buell. 
pair(a,b,u,v,n)=\
 n=gcd(a,b); a=a/n;b=b/n;\
 u=(1/a)%b;\
 v=(1-a*u)/b;\
 [u,v]

triple(a,b,c,u,c1,c2)=\
  u=gcd(a,b);\
  c1=pair(u,c);\
  c2=pair(a,b);\
  [c1[1]*c2[1],c1[1]*c2[2],c1[2]]

cmp(f1,f2,\
   a1,b1,c1,a2,b2,c2,discr,beta,n,s,t,u,v,a,b,c)=\
   a1=f1[1];b1=f1[2];c1=f1[3];\
   a2=f2[1];b2=f2[2];c2=f2[3];\
   discr=b1^2-4*a1*c1;\
   beta=(b1+b2)/2;\
   n=gcd(beta,gcd(a1,a2));\
   s=triple(a1,a2,beta);\
   t=s[1];u=s[2];v=s[3];\
   a=a1*a2/n^2;\
   b=(a1*b2*t+a2*b1*u+v*(b1*b2+discr)/2)/n;\
   c=(b^2-discr)/(4*a);\
   reduce([a,b,c])

\\ Returns the form f1 to the power n, in the class group.
pwr(f1,n,res,res2)=\
 if(n==1,res=f1,\
   if(n%2==0,res2=pwr(f1,n/2);res=cmp(res2,res2),\
             res2=pwr(f1,(n-1)/2);res=cmp(cmp(res2,res2),f1)));\
  res

\\ utility routine to find the exact order of f in the class group. 
\\ Needs the form f, and a bound n on the order of f. Tries all divisors of
\\ n. 
exactorder(f,n,id,d,ff,res)=\
 res=0;\
 fordiv(n,d,if(res==0,ff=pwr(f,d);\
                print1(d);print1("  ");res=eqf(gid,ff);print(res);\
                if(res==0,,res=d),\
                ));\
  res

\\ Given forms f1 and f2 of exact order n, returns a value k such that
\\ f1=k f2. If no such exists, returns k=0. 
multipl(f1,f2,n,ff,res)=\
 ff=f2;\
 res=0;\
for(k=1,n,if(res==0,res=eqf(ff,f1)*k;ff=cmp(ff,f2),));\
 res

\\ The function divisor(n,j) returns the jth divisor of n. 
divisr(n,j,d,k)=\
k=0;\
fordiv(n,d,if(k<j,res=d;k=k+1,));\
res

\\ Generate some quadratic forms. 
\\ Input is b (middle coeff) and j.
try(b,j,a,c,ac)= ac=(b^2-gd)/4;\
    a=divisr(ac,j);\
    c=ac/a;\
 reduce([a,b,c]);

\\ Generate the form which corresponds to the idele (...,1,a mod f,1,... ).
\\   Of discriminant D. 
ram(a,res)=\
 res=[gc^2,trace(a)*gc,norm(a)] ;\
reduce(res)

\\ Generate the form which corresponds to the idele ( ..., p^j, ...).
unram(p,a,j,res)=\
 res=[p,trace(a)*gc,norm(a)*gc^2/p];\
 print(res);\
 res=reduce(res);\
 print(res);\
 pwr(res,j) 



\\ spanned(b,f): Given a basis for a space of forms of order 2, this
\\ function returns 1 if f is spanned by b, and 0 otherwise. 
spanned(b,f,res,ff,bb)=\
 ff=last(b);\
 bb=rest(b);\
 print(ff);print(bb);\
 if(eqf(f,gid),res=1,\
    if(length(b)==1,print(f);print(ff); res=eqf(f,ff),\
      if(spanned(bb,f),res=1,res=spanned(bb,cmp(f,ff)))));\
 res