\\ File tools: contains the software tools for dealing with real
\\ quadratic forms.
\\ 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 
\\ discriminant of a real quadratic order before any other computation can
\\ be carried out. Syntax: setupdisc(D0,f), D0=fund. disc.
setupdisc(d0,f,u,d)=\
 gc=f;\
 gd0=d0;\
 w=quadgen(d0);\
 gd=d0*f^2;\
 gs=sqrt(gd);\
 gsp=0;\
 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);


\\ Routine which, given an element alpha in k_p*/Fp*,
\\ returns the form associated to alpha
\\ of discriminant Dp^2 associated to it; it is assumed that
\\ D is a fundamental discriminant of narrow class number 1, so that
\\ k_p^{*}/Fp* is the narrow class group.
formof(alpha,res)=\
 if(imag(alpha)%37==0,res=[1,trace(w)*gc,norm(w)*gc^2],\
 alpha = alpha/(imag(alpha)) % 37;\
 res=[37^2,37*trace(alpha),norm(alpha)]);\
 res=reduce(res);
 

\\ 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 otherwise. 
eqf(f1,f2,res,ff)=\
 f1=reduce(f1);\
 f2=reduce(f2);\
 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.  syntax: exactorder(f,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

\\ utility routine to find the exact order of alpha in k_p*/Fp*, p=37.
\\ Needs alpha, and a bound n on the order of alpha. Tries all divisors of
\\ n.  syntax: eo(alpha,n)
eo(alpha,n,id,d,ff,res)=\
 res=0;\
 fordiv(n,d,if(res==0,ff= (alpha^d) % 37;\
                print1(d);print1("  ");res=(imag(ff)==0);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

\\ Routine to find a generator for discriminant of the real quadratic
\\ order.
generator(test,x)=\
 x=test;\
 until(eo(x+w,38)==38,x=x+1);\
 (x+w)% 37;

\\ 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




\\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ File setup37



smb(a,b,q,s1,s2,q1,q2,p1,p2,n,lambda,c1,c2,res) =\
 print1("Computing symbol for ");print(a);\
 b= (-1/a)% level;\
 q=exp(-2*pi/level);\
 s1=0;s2=0;\
 q1 = exp(2*pi*i*a/level)*q;\
 q2 = exp(2*pi*i*b/level)*q;\
 p1=q1;p2=q2;\
 for(n=1,200,\
  s1 = s1 + p1*f[n]/n;\
  s2 = s2 + p2*f[n]/n;\
  p1 = p1 * q1;\
  p2 = p2 * q2);\
 lambda=s1-s2 ;\
 c2 = imag(lambda)/imag(omega2);\
 print(c2);\
 c2=round(c2);\
 lambda=lambda-c2*omega2;\
 c1=real(lambda)/real(omega1);\
 print(c1);\
 c1=round(c1);\
 if(delta>0,res=c1+c2*i,\
   res=c1 + c2*(1/2+1/2*i));\
 res

\\ The routine setupcurve() sets up a global table  of modular symbols
\\  to compute the theta-element. Must be invoked...
setupcurve(n,u) =\
   table = vector(level,dummy,x);\
   put(1,0);\
   for(n=2,level-1,\
    u = (1-1/n)% level;\
    if(table[u]<>x,put(n,-table[u]+table[n-1]),\
      if(table[n]<>x,,put(n,smb(n)));\
       put(u,-table[n]+table[n-1])))
 
put(n,s,u) =\
 print1(" Entering the symbol ");print(n);\
 print1(n); print1(" ");print1(table[n]); print1(" -->  ");\
 table[n]=s;\
 print(table[n]);\
 print1(level-n);print1(" ");print1(table[level-n]);print1(" -->  ");\
 table[level-n]= conj(s);\
 print(table[level-n]);\
 u = (1/n) % level;\
 print1(u);print1(" ");\
 print1(table[u]);print1(" --> ");\
 table[u]= -conj(s);\
 print(table[u]);\
 print1(level-u);print1(" ");\
 print1(table[level-u]);print1(" --> ");\
 table[level-u]= -s;\
 print(table[level-u]); 

\\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ 
\\ File: 37 (Computation of the modular symbols) 

md(a,b,res)=\
 res=a%b;\
 if(2*res>abs(b),res=res-abs(b),);\
 res


symbol(a,b,w,ww,aa,bb,res) =\
 until( ((2*abs(a) <= abs(b)*level) && (2*abs(b)<= abs(a))) || (b==0),\
   a=md(a,b*level);\
   b=md(b,a) );\
 if(b<0,b=-b;a=-a,);\
 if(b==0,res=0,\
  if(b==1,res=table[a%level],\
     w=round(-a/b);\
     ww = (1/w)% level;\
     bb= b*w + a;\
     aa= (a* ww + b*(w*ww-1));\
     res=symbol(aa,bb)-table[ww] ));\
 res
  

adjust(f,a,b,c,res)=\
  a=f[1];b=f[2];c=f[3];\
  if (a%level==0,if( (b+gsp) % level ==0,res=f,res=adjust([c,-b,a])),\
    k=( ( -b+gsp)/(2*a)) % level;\
    res=shft(f,k);\
    res=[res[3],-res[2],res[1]]);\
   res

geodesic(f,ff,a,b)=\
 ff=adjust(f);\
  a=ff[1];b=ff[2];\
  symbol(gu+gv*b,-2*a*gv/level)

computetheta(f,n, th,dummy,ff,j)=\
 th = vector(n,dummy,0);\
 ff=f;\
 for(j=1,n,\
  th[j]=geodesic(ff);\
  print1(j);print1("  ");print(th[j]);\
  ff=cmp(f,ff));\
 th



    
main(d,test)=\
 e=initell([0,0,1,-1,0]);\
 level=37;\
 f = anell(e,200);\
 omega1=e[15];\
 omega2=e[16];\
 delta = e[12];\
 setupcurve();\
 setupdisc(d,37);\
 alpha=generator(test);\
 f=formof(alpha);\
 th=computetheta(f,38);\
 print("th= ",th);\
 s=0;\
 for(j=1,38,s=s+th[j]);\
 print("Check Sum: ",s);\
 for(j=1,38,s=s+j*th[j]);\
 print(" s= ",s);\
 s = real(s);\
 s = s % 38;\
 print("s reduced =",s);\
 koly=(alpha^s) % 37 ;\
 print("Kolyvagin derivative: ",koly);\
 koly = koly/conj(koly) % 37;\
 print(" projected: ", koly);

\\ p=phitate(koly)