/* test-eichler.gp */

\r eichler.gp

testalgfromcenter(A) =
{
  my(nf,n,x,a,pol,y,b,c,d);
  nf = algcenter(A);
  n = poldegree(nf.pol);
  x = vectorv(n,i,random(-10));
  a = algfromcenter(A,x);
  pol = minpoly(nfbasistoalg(nf,x));
  y = vectorv(n,i,random(-10));
  b = algfromcenter(A,y);
  c = algfromcenter(A,nfeltmul(nf,x,y));
  d = algfromcenter(A,x+y);
  [algpoleval(A,pol,a) == 0, \\image vanishes on minpoly
   c == algmul(A,a,b), \\map is multiplicative
   d == algadd(A,a,b)] \\map is additive
};

vec2mat(d,n,V,g) =
{
  my(M = matrix(d,d), c);
  for(i=0,d^2*n-1,
    c = V[i+1];
    M[i\(n*d) + 1, (i\n)%d + 1] += c*g^(i%n)
  );
  M
};

testalgmodpr(A) =
{
  my(map,mapi,nf,dec,pr,N,x,y,xy,fx,fy,fxy,d,n,L,M1,M2,g);
  L = List();
  nf = algcenter(A);
  forprime(p=1,10,
    dec = idealprimedec(nf,p);
    pr = dec[#dec];
    if(algtype(A)==2,
      if(algdisc(A)%p==0, next);
    ,\\else
      if(algisramified(A,pr), next);
    );
    d = algdegree(A);
    n = pr.f;
    [map,mapi] = algmodpr(A,pr);
    listput(~L, matsize(mapi)[2] == n*d^2);
    N = algdim(A,1);
    x = vectorv(N,i,random(-10));
    y = vectorv(N,i,random(-10));
    fx = map*x;
    fy = map*y;
    xy = algmul(A,x,y);
    fxy = map*xy;
    listput(~L, fxy == fx*fy); \\map is multiplicative
    x = vectorv(d^2*n,i,random(p));
    g = ffgen(map[1][1,1]);
    M1 = map*(mapi*x);
    M2 = vec2mat(d,n,x,g);
    listput(~L, M1 == M2); \\maps are inverse of each other
    x = idealhnf(nf,pr)*vectorv(poldegree(nf.pol),i,random(p));
    listput(~L,map*algfromcenter(A,x) == 0); \\pr maps to 0
  );
  Vec(L);
};

testeichleridempotent(A) =
{
  my(L,nf,dec,pr,e,map,mapi,M);
  L = List();
  nf = algcenter(A);
  forprime(p=1,10,
    dec = idealprimedec(nf,p);
    pr = dec[#dec];
    if(algtype(A)==2,
      if(algdisc(A)%p==0, next);
    ,\\else
      if(algisramified(A,pr), next);
    );
    e = eichleridempotent(A,pr);
    [map,mapi] = algmodpr(A,pr);
    M = map*e;
    listput(~L, M^2 == M); \\idempotent mod pr
    listput(~L, matrank(1-M) == 1); \\1-e mod pr has rank 1
  );
  Vec(L);
};

diag(M) = vector(#M~,i,M[i,i]);

testeichlerprimepower(A) =
{
  my(L,nf,dec,pr,n,ord,lat,N,d,Ap,mt,B,Bi,P,AP,comm,Pk);
  L = List();
  nf = algcenter(A);
  N = algdim(A,1);
  d = algdegree(A);
  forprime(p=1,10,
  \\forprime(p=1,30,
    dec = idealprimedec(nf,p);
    pr = dec[#dec];
    n = pr.f;
    if(algtype(A)==2,
      if(algdisc(A)%p==0, next);
    ,\\else
      if(algisramified(A,pr), next);
    );
    for(k=1,3,
    \\for(k=1,20,
      ord = eichlerprimepower(A,pr,k);
      listput(~L, ord[,1]==1); \\contains 1
      listput(~L, denominator(ord)==1); \\contained in maxord
      lat = [ord,1];
      listput(~L, vecprod(diag(ord)) == idealnorm(nf,pr)^(k*(d-1))); \\index
      B = matimagemod(ord,p);
      Bi = lift(Mod(B,p)^(-1));
      mt = [Bi*algtomatrix(A,b,1)*B | b <- Vec(B)];
      listput(~L, denominator(mt)==1); \\stable under mult
      Ap = algtableinit(mt,p);
      P = matimagemod(algtomatrix(A,algfromcenter(A,pr.gen[2]),1),p); \\pr*maxord
      P = lift(Bi*matintersect(Mod(B,p),Mod(P,p)));
      AP = algquotient(Ap,P);
      listput(~L, algdim(AP) == n*(d^2-(d-1))); \\dim ord/pr*maxord
      my([J,ssdec] = algsimpledec(AP));
      listput(~L, (d==1 && J==0) || matsize(J)[2] == (d-1)*n); \\dim radical
      listput(~L, d==1 || (#ssdec == 2 \\2 semisimple factors
        && algdim(ssdec[1])==n \\first is residue field
        && #algcenter(ssdec[2])==n)); \\second is M_{d-1}(residue field)
      Pk = algfromcenter(A,idealtwoelt(nf,idealpow(nf,pr,k))[2]);
      Pk = mathnfmodid(algtomatrix(A,Pk,1),p^k);
      listput(~L, alglatsubset(A,alglathnf(A,Pk),lat)); \\contains pr^k*maxord
      comm = matimagemod(Mat([algmul(A,b1,b2)-algmul(A,b2,b1) | b1 <- Vec(ord); b2 <- Vec(ord)]),p^k);
      comm = mathnfmodid(concat(comm,Pk),p^k);
      listput(~L, vecprod(diag(comm)) == idealnorm(nf,pr)^(k*if(d==1,1,d+1))); \\index commutator
    );
  );
  Vec(L);
};

alltests(A) =
{
  print(testalgfromcenter(A));
  print(testalgmodpr(A));
  print(testeichleridempotent(A));
  print(testeichlerprimepower(A));
};

nf = nfinit(y);
A = alginit(nf,1);
print("trivial algebra");
alltests(A);

A = alginit(nf,[-1,-1]);
print("quatalg/Q");
alltests(A);

A = alginit(nf,algmultable(A));
print("table quatalg/Q");
alltests(A);

pr1 = idealprimedec(nf,3)[1];
pr2 = idealprimedec(nf,5)[1];
A = alginit(nf, [3,[[pr1,pr2],[1/3,-1/3]],[0]]);
print("degree 3 alg/Q");
alltests(A);

A = alginit(nf, [5,[[pr1,pr2],[1/5,-1/5]],[0]]);
print("degree 5 alg/Q");
alltests(A);

nf = nfinit(y^2+3);
A = alginit(nf,1);
print("degree 1 alg/IQF");
alltests(A);

A = alginit(nf,[y,7+9*y]);
print("quatalg/IQF");
alltests(A);

\\TODO table CSA/IQF

pr1 = idealprimedec(nf,3)[1];
pr2 = idealprimedec(nf,7)[1];
A = alginit(nf, [3,[[pr1,pr2],[1/3,-1/3]],[]]);
print("degree 3 alg/IQF");
alltests(A);

A = alginit(nf, [5,[[pr1,pr2],[1/5,-1/5]],[]],,1);
print("degree 5 alg/IQF");
alltests(A);