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