\\ This file contains the PARI calculations that were performed 
\\ in section  5 of the paper 
\\ "Integration on H_p x H and arithmetic applications"
\\  
\\
\\ The calculation can be divided into three parts:
\\ 
\\ Part I: Initialization of the data concerning X_0(11) 
\\ Part 2: Initialization of the data on the real quadratic order
\\ Part 3: Computing the points

\\ Global variables: ggp and ggn. 
\\ Calculations on the curve will be done modulo ggp^ggn.
ggp = 11;
ggn =  3;

\\ *********************************************************************
\\ ** PART I : The curve X_0(11) 
\\ *********************************************************************
\\ ggmsymbol = ([C, 0, A, B, B-A, -A, -A, -B, A-B, A, 0, -C]);\
\\
gge = initell([0,-1,1,-10,-20]);\
 ggmsymbol = ([[0,0,1], [0,0,0], [1,0,0], [0,1,0], [-1,1,0], [-1,0,0], [-1,0,0], [0,-1,0], [1,-1,0], [1,0,0], [0,0,0], [0,0,-1]]);\

\\ *********************************************************************
\\ ** Modular symbol routines 
\\ *********************************************************************
\\
\\ path(a,b,c,d): computes the symbol { a/b, c/d } as in Cremona, i.e.
\\ the path from (a/b) to (c/d)
path(a,b,c,d)= pathinfty(c/d) -  pathinfty(a/b)

\\ pathinfty(a/b): computes the symbol :  { infty, c/d }
pathinfty(x,\
   c,ll,res,qleft,qmiddle,qright,sgn,j) = \
 c = cf(x) ;\
 ll=length(c);\
 res=0;\
 qleft = 1; qmiddle = 0;\
 sgn = 1;\
 for(j=1,ll,\
  qright = c[j]*qmiddle + qleft;\
  if(sgn==1,res = res - smb(qmiddle,qright),res = res + smb(qright,qmiddle));\
  sgn=-sgn;\
  qleft=qmiddle; qmiddle = qright);\
 res;


smb(c,d)= \
  if( (d% ggp)==0, ggmsymbol[ggp+1], ggmsymbol[((c/d)% ggp)+1] );



\\ **********************************************************************
\\ ** Initialization routine: InitForms(D), where D is a (not necessarily
\\ **  fundamental) discriminant. This sets up certain global
\\ **  variables: 
\\ ** ggD: The discriminant D of the order.
\\ ** ggw: The element such that [1,gw] is a basis for the order.
\\ ** ggUnit: The fundamental unit of norm 1 in the order of
\\ **              discriminant D. 
\\ ** ggMUnit: The matrix giving the multiplication by the unit, in the 
\\ **       basis  [1, sqrt{D}]
\\ ** ggH: The class number of the order.
\\ ** ggClassGroup: The class group of the order, as an array of 
\\ **   representative forms. 
\\ ** ggi: an index which divides 12 = 11+1, and is used to make the unit 
\\ **      congruent to a scalar in O^*. 
\\ ** ggsd: the square root of d
\\
InitForms(d,\
            temp,structure,generators,j,u,v)=\
 ggD = d;\
 ggw = quadgen(d);\
 ggi = 1;\
 ggunit = unit(d)^ggi;\
 u = real(ggunit); v= imag(ggunit);\
 if( (ggD % 4) ==1, u = u+ v/2; v=v/2, v=v/2);\
 ggMunit = [ u, ggd * v;\
	     v, u ];\
 temp=buchreal(d,1,0.5,0.5,20);\
 ggh = temp[1];\
 if(norm(ggunit)==1,ggh=ggh*2,);\
 structure=temp[2];\
 generators=length(structure);\
 if(generators==0,gggen=pf(ggd,2),gggen=temp[3][1]);\
 gggen=reduce(stripd(gggen));\
 ggclassgroup=temp[3];\
  if((d%4)==0,ggsd=2*w,ggsd=2*w-1);\
 for(j=1,generators,ggclassgroup[j]=stripd(ggclassgroup[j]));

\\ stripd: A utility routine to convert a real quadratic form
\\ qfr(a,b,c,d) to the usual representation [a,b,c]
stripd(form)=\
 [compo(form,1),compo(form,2),compo(form,3)];

\\ convertqfr: Performs the converse conversion,i.e,. turns a 
\\ quadratic form of the type [a,b,c] into qfr(a,b,c,0).
convertqfr(form)=\
 qfr(form[1],form[2],form[3],0.0);


\\ reduces the form f with positive discriminant. 
reduce(f)=\
 stripd(redreal(convertqfr(f)) );

\\ next(f):
\\ returns the form adjacent to f in the cycle of reduced forms.
next(f)=\
 stripd(rhoreal(convertqfr(f)) );
  
\\ eqf(f1,f2):
\\ Returns 1 if f1 is equivalent to f2 in the class group, 0 otherwise. 
eqf(f1,f2,res,ff)=\
 f1=reduce(f1);\
 ff=f2;\
 until(ff==f1 || ff==f2, ff=next(ff));\
 if(ff==f1,res=1,res=0);\
 res

\\ The folllowing routine implements the Guassian composition of
\\ binary quadrartic forms. 
compose(f1,f2)=\
  stripd(redreal(comprealraw( convertqfr(f1),convertqfr(f2) )));

\\ pwr(f1,n):
\\ Returns the form f1 to the power n, in the class group.
pwr(f1,n)=\
   stripd(redreal(powrealraw(convertqfr(f1),n)));

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

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

\\ Generate some quadratic forms. 
\\ Input is a prime p (first coeff) and j.
try(p)= reduce(stripd(pf(ggd,p)))



\\ ************************************************************************
\\ **  Routines for computing the geodesic cycles and the local points.
\\ ************************************************************************
\\
\\
\\ IPsi(F): takes as input a binary quadratic form F=(a,b,c) of discriminant
\\              D, and returns the period attached to the embdding Psi sending
\\              sqrt(D) to the matrix    | b      2c |, and hence 
\\                                                  | -2a  -b  |.
\\             w to   | (b+1)/2          c     |
\\                      | -a          (-b+1)/2  |
IPsi(F,\
  a,b,c,matw,gammaPsi,zPsi) =\
  a=F[1]; b=F[2]; c=F[3];\
  d= (b^2-4*a*c);\
  if( (d%4)==0, matw=[b,2*c;-2*a,-b],\
 matw = [(b+1)/2,c; -a, (-b+1)/2 ]);\
 gammaPsi = real(ggUnit)*[1,0;0,1] + imag(ggunit)*matw;\
  zPsi = (- b + ggsd )/(2*a);\
  dmint(zPsi, gammaPsi[1,1]/gammaPsi[2,1]));\

\\ The function dmint(z,q) implements the double integral introduced in the 
\\  paper, in an important special case: that where z is defined over the quadratic
\\ unramified extension of Q_p, and is of the form  a+bw, where a is a rational
\\ number in Zp and b is in Zp*. 
\\ The integral computed is dmint(z-bar,z; infinity,q).
dmint(z,q,\
  res,p,n,u,alpha,lhs,rhs,Nalpha,M,v) =\
  res=1;\
  p=ggp; n=ggn;\
  for(u=0,p^n-1,\
   alpha = u+w;\
    lhs=alpha-conj(alpha);\
    rhs = (alpha*z)-conj(alpha*z);\
    Nalpha = (rhs/lhs) % (p^n);\
     M = [1,-Nalpha; 0,p^n];\
     res=(res* (alpha/conj(alpha))^pathinfty( (q-Nalpha)/p^n) ) % p^n);\
  for(v=0,p^(n-1)-1,\
    alpha = 1+p*v*w;\
     lhs = (alpha-conj(alpha));\
     rhs = (alpha*z)-conj(alpha*z);\
     Nalpha = ((-lhs/rhs) % (p^n));\
     res=(res*(alpha/conj(alpha))^path(-Nalpha/p^n,(-Nalpha*q-1)/(p^n*q))) % (p^n) );\
  res)