H = alginit(1.,1/2);
algdim(H)
algisdivision(H)
algiscommutative(H)

w = [Pi,2,sqrt(3),-7]~
algmul(,w,w)
alginv(,w)
alginvol(,w)
algtrace(,w)
algnorm(,w)

W = algtomatrix(,w)
trace(W) - algtrace(,w)
matdet(W) - algnorm(,w)
W2 = algtomatrix(,w,1)
trace(W2) - 2*algtrace(,w)
matdet(W2) - algnorm(,w)^2

u = w/sqrt(algnorm(,w))
U = algtomatrix(,u)
exponent(conj(U)~ * U - matid(2))

C = alginvol(H);
rot(w) = ((algtomatrix(,w,1)*C)^2)[^1,^1];
R = rot(u)
exponent(R~*R - matid(3))

nf = nfinit(y^4-y-1);
al = alginit(nf,[-7,y]);
algdim(al)
algdim(al,1)
algiscommutative(al)
algissimple(al)

[a,b] = algisquatalg(al)

z = algquattobasis(al, [-2,y,1+y,y/2]~)
lift(algbasistoquat(al, alginvol(al, z)))
algtrace(al, z)
algpoleval(al, algcharpoly(al,z), z)

Z = algtomatrix(al,z); liftall(Z)
trace(Z) - algtrace(al,z)
matdet(Z) - algnorm(al,z)
Z2 = algtomatrix(al,z,1); matsize(Z2)
trace(Z2) - 2*nfelttrace(nf,algtrace(al,z))
matdet(Z2) - nfeltnorm(nf,algnorm(al,z))^2

algisramified(al, 1)
pr = idealprimedec(nf,2)[1];
algisramified(al, pr)

nfhilbert(nf, a, b, pr)

algramifiedplaces(al)
algisdivision(al)

al2 = alginit(nf, [[pr],[0,1]]);
#algramifiedplaces(al2)
algisramified(al2,pr)
algisquatalg(al2)

mt = algmultable(al);
denominator(mt)

algdisc(al)
D = idealnorm(nf, algramifiedplaces(al)[2]);
2^algdim(al,1) * (nf.disc^2 * D)^2

pr3 = idealprimedec(nf,3)[1];
pr3.f
modP3 = algmodprinit(al, pr3);
pr7 = algramifiedplaces(al)[2];
pr7.f
modP7 = algmodprinit(al, pr7);

algmodpr(al, z, modP3)
algmodpr(al, z, modP7)
t = algquattobasis(al,[0,1,2,1/7*y^3+1/7*y-2/7]~);
algmodpr(al, t, modP7)

li1 = algmodprlift(al, [1,x;0,1], modP3)
algmodpr(al, li1, modP3)
li2 = algmodprlift(al, Mat(x), modP7)
algmodpr(al, li2, modP7)

quatembed(al,x,pl) =
{
  my([a,b]=algisquatalg(al),nf=algcenter(al),x1,x2,x3,x4,sa,sb);
  x = algbasistoquat(al,x);
  [a,b,x1,x2,x3,x4] = apply(y->nfeltembed(nf,y,pl), concat([a,b],x~));
  if(pl>nf.r1 || a>0,
    sa=sqrt(a);
    return([x1+x2*sa,b*(x3+x4*sa);x3-x4*sa,x1-x2*sa]));
  if(b>0,
    sb=sqrt(b);
    return([x1+x3*sb,a*(x2-x4*sb);x2+x4*sb,x1-x3*sb]));
  sa = sqrt(-a);
  sb = sqrt(-b);
  [x1,x2*sa,x3*sb,x4*sa*sb]~
};

quatembed(al,z,1)
quatembed(al,z,2)
quatembed(al,z,3)

eich = algeichlerbasis(al, Mat([pr3,3]))

lat1 = alglathnf(al,z);
lat1[2]
lat2 = alglathnf(al,eich);
lat2[2]

alglatsubset(al,lat1,lat2)
alglatsubset(al,lat2,lat1)
alglatindex(al,lat2,lat1)

alglatlefttransporter(al,lat2,lat2) == lat2
triv = alglathnf(al,matid(16));
lat1inv = alglatlefttransporter(al,lat1,triv);
alglatmul(al,lat1inv,lat1) == triv