\\ This file contains the PARI calculations that were performed 
\\ in section 5.1 of the paper
\\ " Stark-Heegner Points over real quadratic fields " 
\\
\\ 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 =  2;

\\ *********************************************************************
\\ ** 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^*. 
\\
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);\
 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];\
 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.
\\ ************************************************************************
\\
\\ geodesic(f,alpha) : takes as input a form f of discriminant D, an 
\\   element alpha in O/p^nO^*. Returns the geodesic cycle associated
\\   to the lattice of discriminant Dp^{2n} associated to F and alpha
\\ (For details, see the paper...)
geodesic(f,alpha,\
   u,v,a,b,lambda,lambda0,m) =\
 u = real(alpha);\
 v = imag(alpha);\
 if( (ggD % 4)==1, u=u+v/2; v=v/2 , v=v/2);\
 a=f[1]; b=f[2];\
 if( (v% ggp)==0,\
       lambda  = (2*a*v)/(u+b*v) % (ggp^ggn);\
       lambda0 = lambda - ggp^(ggn-1);\
       m = [ 2*a*ggp - b* lambda0 * ggp , 2*a- b*lambda ;\
	     -ggp*lambda0                , -lambda ];\
       m = m^(-1) * ggmunit * m,\
       lambda  = (u+b*v)/(2*a*v) % (ggp^ggn);\
       lambda0 = lambda + ggp^(ggn-1);\
       m = [ 2*a*ggp* lambda0 - b* ggp   , 2*a*lambda - b ;\
	     -ggp                        , -1 ];\
       m = m^(-1) * ggmunit * m);\
     pathinfty(m[1,1]/m[2,1])  
       
geodesics(f,\
   n,p,j) = \
   n=ggn;\
   p=ggp;\
   j=1;\
   while(ordr(j+ ggw)< (p+1)*p^(n-1), j=j+1);\
   ggalpha = j+ggw;\
   print("Primitive Element: ",ggalpha);\
   vector((p+1)*p^(n-1),j,geodesic(f,ggalpha^j));

ordr(x,\
  j) =\
 j=1;\
 while(imag(x^j % (ggp^ggn)),j=j+1);\
 j








check(v,\
 j,res)=\
  res=0;\
 for(j=1,(ggp+1)*ggp^(ggn-1),res=res+v[j]);\
 res;





koly(v,\
  j,res)=\
  res=0;\
  for(j=1,(ggp+1)*ggp^(ggn-1),res=res+j*v[j]);\
  res % ((ggp+1)*ggp^(ggn-1))