/* eichler.gp */

/*
  This GP script implements Eichler orders over arbitrary central simple
  algebras, in the following sense.
  Let O be the order stored in an alginit structure A.
  Let pr be a prime ideal of the base field F that is unramified in O and k an integer.
  Let f: O/pr^k ---> M_d(ZF/pr^k) be an isomorphism.
  Then eichlerprimepower(A,pr,k) computes a basis for a preimage by some f of the
  set of matrices whose first column is 0 except the top entry:
  [* * * * ]
  [0 * * * ]
  [0 * * * ]
  [0 * * * ]
*/

algfromcenter(A,x) =
{
  my(y);
  y = vectorv(algdim(A,1));
  y = algbasistoalg(A,y);
  y[1] = nfbasistoalg(algcenter(A),x);
  algalgtobasis(A,y); \\We need an easier way to embed the center!!!
};

algmodpr(A,pr) =
{
  my(p,Ap,g,Q,pro,lif,map,mapi);
  p = pr.p;
  Ap = algtableinit(algmultable(A),p);
  g = algfromcenter(A,pr.gen[2]);
  g = algtomatrix(Ap,g);
  g = matimagemod(g,p); \\ideal generated by pr.gen[2]
  [Q,pro,lif] = algquotient(Ap,g,1);
  [map,mapi] = algsplit(Q);
  [map*pro,lif*mapi]
};

\\e in A such that e mod pr is an idempotent of maximal rank
eichleridempotent(A,pr) =
{
  my(map,mapi,d,n,e);
  [map, mapi] = algmodpr(A,pr);
  d = algdegree(A);
  n = pr.f;
  e = vectorv(n*d^2);
  for(i=1, d-1, e[1+n*(d+1)*i] = 1); \\We need an automatic conversion 
                                     \\  from matrices to vectors!!!
  mapi*e
};

eichlerprimepower(A,pr,k) =
{
  my(e,p,polidem = 3*x^2-2*x^3,B,N,Me,Mzk,nf,Mprk);
  p = pr.p;
  e = eichleridempotent(A,pr);
  i = 1;
  while(i<k,
    e = algpoleval(A,polidem,e); \\Newton iteration
    i *= 2;
    e %= p^ceil(i/pr.e);
  );
  N = algdim(A,1);
  B = Vec(matid(N));
  Me = Mat(vector(N,i,algmul(A,B[i],e))); \\We need access to the right multiplication table!!!
  nf = algcenter(A);
  Mzk = Mat([algfromcenter(A,z) | z <- nf.zk]);
  prk = idealtwoelt(nf,idealpow(nf,pr,k));
  Mprk = algtomatrix(A,algfromcenter(A,prk[2]),1);
  mathnfmodid(matconcat([Me,Mzk,Mprk]),prk[1]) 
};