\\ General purpose routines for working in the ring of Hamilton
\\ quaternions.

hadd(p,q,x,y,z,w)= x=p[1]+q[1];\
           y=p[2]+q[2];\
           z=p[3]+q[3];\
           w=p[4]+q[4];\
           [x,y,z,w];

hmult(p,q,x,y,z,w)= x=p[1]*q[1] - p[2]*q[2] - p[3]*q[3] - p[4]*q[4];\
           y=p[1]*q[2] + p[2]*q[1] +p[3]*q[4] - p[4]*q[3];\
           z=p[1]*q[3]+p[3]*q[1] -p[2]*q[4] + p[4]*q[2];\
           w=p[1]*q[4] + q[1]*p[4] +p[2]*q[3] - p[3]*q[2];\
           [x,y,z,w];

hnorm(p) = p[1]^2 + p[2]^2 + p[3]^2 + p[4]^2;

hconj(p) = [p[1],-p[2],-p[3],-p[4]];

hinv(p) = hconj(p)/hnorm(p);

hpow(p,n,q) = if(n==0,[1,0,0,0],\
               q=hpow(p,n\2);\
               if(n%2==0,hmult(q,q),\
                   hmult(p,hmult(q,q))))

hround(p,pp) = pp=round(2*p);\
            if((pp%2 ==[0,0,0,0]) || (pp%2==[1,1,1,1]),\
                pp/2,round(p));

hgcd(a,b,q,r) = q = hround(hmult(a,hinv(b)));\
                r = a-hmult(q,b);\
                if(r==[0,0,0,0],b,hgcd(b,r));

\\ The function "adjust" modifies an "adelic" quaternion concentrated at
\\ 7 by multiplying it on the right by global elements of norm 7 until
\\ it becomes a unit at 7. 
adjust(q) = if(hnorm(q)%7==0,adjust(hmult(q,hinv(hgcd(q,[7,0,0,0])))),q);


\\ The function "class" returns the class in R-hat*\B-hat*/B* that an idele
\\  which is concentrated at 7 belongs to. There are two classes in this
\\ space (i.e., two components on Gross' Shimura curve) and these classes
\\ are denoted symbolically by -1 and 1. We first modify the idele by global 
\\ elements so that it is a unit at 7, using the function "adjust". 
\\ Then, we compute the image of the CONJUGATE of the
\\ resulting matrix in GL_2(F7), using
\\ the  local embedding
\\ i --> 0  1       j --> 2  3       k --> 3 -2
\\      -1  0             3 -2            -2 -3
\\ If this matrix sends infinity to 0,2,3 or infinity, then it is the class 1.
\\ If it maps infinity to 1,4,5,6, then it is the class -1. 

class(q, qq, num,den,res,tst) = qq=adjust(q);\
           print(" Adjusted: ",qq);\
            tst= hmult(hinv(qq),q);\
           print ("        Ratio: ", tst, "   has norm: ",hnorm(tst));\
           num = qq[1] - 2*qq[3] - 3*qq[4];\
           den = qq[2] - 3*qq[3] + 2*qq[4];\
           print("num: ",num, "      den: ",den);\
           if(den%7==0,1,res=num/den % 7;\
                    if(res==0 || res==2 || res==3,1,-1));


\\ The function symbol(a,b,n) returns the component on which the special 
\\ point of level 7^n, corresponding to the idele a+b.w concentrated at
\\ 7, lands.  Here w is the element (1+sqrt(-11))/2 in the quadratic field
\\ of class number 1. We start with an initial embedding sending w to 
\\ [1,1,1,3]/2, and we conjugate this by [2,1,1,1]^n to obtain an
\\ optimal embedding of the order of conductor 7^(n+1).

symbol(a,b,n,q,qq) =  q=[a,0,0,0]+ b*[1/2,1/2,1/2,3/2];\
                  print("q= ",q);\
                  qq=hmult(hpow([2,1,1,1],n),q);\
                  print("qq= ",qq,"   norm: ", hnorm(qq));\
                  class(qq)
 

\\ The program go() takes as  input parameter the level n of 7-adic accuracy
\\ that is desired. It assumes that there is 
\\ a global array gross() containing 8*7^n entries, and it computes
\\ the special points  corresponding to the 8*7^n optimal embeddings of
\\ the order of conductor 7^(n+1).
go(n, alpha, w, lim1, lim2, j, t, a, b) = w=quadgen(-11);\
        lim1 = 8*7^n;\
        lim2 = 7^(n+1);\
        for(j=1,lim1,\
            t = alpha^j % lim2;\
            a=real(t); b=imag(t);\
            gross[j]=symbol(a,b,n) );