1al = alginit(nfinit(y), [Mat(1)]); \\reported by Bill
2J = varhigher("J",y);
3default(realprecision,38);
4default(parisize,"100M");
5Q=nfinit(y); T=polcyclo(5, 'x); F=rnfinit(Q, T);
6A = alginit(F, [Mod(x^2,T), 3]);
7print("contains nfabs: ", algsplittingfield(A)[12][1] != 0);
8A[1][12][1] = 0;
9A
10
11
12do(name, test) = {
13 setrand(1);
14 print(name,": ", iferr(test(), E, E));
15}
16gusuite(name, tests) = print("Suite: ", name); tests();
17
18searchin(hf,pr,h) =
19{
20  my(i,n);
21  for(i=1,#hf[1],
22    if(hf[1][i]==pr, return(hf[2][i]==h))
23  );
24  return(h==0);
25};
26
27samehasse(hf1,hi1,hf2,hi2) =
28{
29  my(i,n);
30  if(hi1 != hi2, return(0));
31
32  n = #hf1[1];
33  for(i=1,n,
34    if(!searchin(hf2,hf1[1][i],hf1[2][i]), return(0));
35  );
36
37  n = #hf2[1];
38  for(i=1,n,
39    if(!searchin(hf1,hf2[1][i],hf2[2][i]), return(0));
40  );
41
42  return(1);
43};
44
45hassetriv(hf,hi) = samehasse(hf,hi,[[],Vecsmall([])],Vecsmall(vector(#hi,i,0)));
46altriv(al) = hassetriv(alghassef(al),alghassei(al));
47alsame(al,hf,hi) = samehasse(alghassef(al),alghassei(al),hf,hi);
48
49testcharpol(al,a)=
50{ my (T = algcharpoly(al,a));
51  print(T);
52  !algpoleval(al,T,a);
53}
54get() = gusuite("get", ()->{
55  my(s='s,i='i,poli,nf,ii,pols,rnf,ss,al);
56  poli = i^2+1;
57  nf = nfinit(poli);
58  ii = Mod(i,poli);
59  pols = s^2+2;
60  rnf = rnfinit(nf, pols);
61  ss = Mod(s,pols);
62  al = alginit(rnf,[-ss,ii]);
63
64  do("degree", ()->algdegree(al)==2);
65  do("center", ()->algcenter(al)==nf);
66  do("splitting", ()->algsplittingfield(al)==rnf);
67  do("automorphism", ()->algaut(al)==-ss);
68  do("b", ()->algb(al)==ii);
69  do("trivial hasse invariants", ()->altriv(al));
70  do("charac", ()->algchar(al)==0);
71  do("dim", ()->algdim(al)==4);
72  do("absdim", ()->algdim(al,1)==8);
73  do("basis", ()->#algbasis(al)==8);
74  do("invbasis", ()->#alginvbasis(al)==8);
75  do("basis*invbasis", ()->algbasis(al)*alginvbasis(al)==matid(8));
76  do("iscyclic", ()->algtype(al)==3);
77  do("radical", ()->algradical(al)==0); \\cyclic => simple
78});
79
80operations() = gusuite("operations", ()->{
81  my(s='s,i='i,poli,nf,ii,pols,rnf,ss,al,n,un,u,j);
82  poli = i^2+1;
83  nf = nfinit(poli);
84  ii = Mod(i,poli);
85  pols = s^2+2;
86  rnf = rnfinit(nf, pols);
87  ss = Mod(s,pols);
88  ss = ss*Mod(1,poli);
89  al = alginit(rnf,[-ss,ii]);
90
91  do("radical", ()->algradical(al)==0); \\cyclic => simple
92
93  nfii = ii;
94  ii = ii*Mod(1,pols);
95
96  n = [-ss,ii-1]~;
97  un = [1,0]~;
98  u = [1-ss,ii-1]~;
99  j = [0,1]~;
100  do("addition", ()->algadd(al,n,un)==u);
101  do("negation", ()->algneg(al,u)==[ss-1,1-ii]~);
102  do("soustraction", ()->algsub(al,u,n)==un);
103  do("multiplication", ()->algmul(al,n,un)==n);
104  do("non-commutativity", ()->algmul(al,n,j)==algmul(al,j,n));
105  do("left division", ()->algdivl(al,u,n)==n);
106  do("right division", ()->algdivr(al,n,u)==n);
107  do("noncommutative left division", ()->algdivl(al,u,j)==[ii+1,1-ss]~);
108  do("noncommutative right division", ()->algdivr(al,j,u)==[ii+1,1+ss]~);
109  do("division by non-invertible", ()->algdivl(al,n,u));
110  do("nilpotent", ()->algsqr(al,n)==[0,0]~);
111  do("square", ()->algsqr(al,u)==algadd(al,u,n));
112  do("square j", ()->algsqr(al,j)==[ii,0]~);
113  do("inverse", ()->alginv(al,u)==algsub(al,un,n));
114  do("powers", ()->algpow(al,u,124)==algadd(al,un,algmul(al,[124,0]~,n)));
115  do("negative powers", ()->algpow(al,j,-56)==un);
116  do("multiplication table j", ()->algtomatrix(al,j)==[0,ii;1,0]);
117  do("multiplication table", ()->algtomatrix(al,u)==[1-ss,-1-ii;ii-1,1+ss]);
118  do("characteristic polynomial", ()->algcharpoly(al,u,'y)==y^2-2*y+1);
119  do("characteristic polynomial j", ()->algcharpoly(al,j,'y)==y^2-nfii);
120  do("trace zero", ()->algtrace(al,n)==0);
121  do("trace commutator", ()->algtrace(al,algsub(al,algmul(al,j,u),algmul(al,u,j)))==0);
122  do("trace", ()->algtrace(al,algmul(al,u,j))==-2-2*nfii);
123  do("norm zero", ()->algnorm(al,n)==0);
124  do("norm one", ()->algnorm(al,u)==1);
125  do("norm j", ()->algnorm(al,j)==-nfii);
126  a = algadd(al,u,j); b = algadd(al, n, [12-4*ii+ss*(4-ii),7+3*ii+ss*(1+7*ii)]~);
127  do("norm is multiplicative a*b", ()->nfalgtobasis(nf,algnorm(al,algmul(al,a,b)))==nfalgtobasis(nf,nfeltmul(nf,algnorm(al,a),algnorm(al,b))));
128  do("norm is multiplicative b*a", ()->nfalgtobasis(nf,algnorm(al,algmul(al,b,a)))==nfalgtobasis(nf,nfeltmul(nf,algnorm(al,b),algnorm(al,a))));
129  do("poleval", ()->algbasistoalg(al,algpoleval(al,x^3-2,u)) ==\
130algsub(al,algpow(al,u,3),[2,0]~));
131  do("poleval b", ()->algbasistoalg(al,algpoleval(al,x^3+i*x-3*i,u)) ==\
132algsub(al,algadd(al,algpow(al,u,3), algmul(al,[i,0]~,u)), [3*i,0]~));
133});
134
135tensor() = gusuite("tensor product of cyclic algebras", ()->{
136  my(nf,pol,jj,al1,al2,al3,hf12,hf23,p7,p7b,p3,p5);
137  pol = y^2+y+1;
138  nf = nfinit(pol);
139  jj = Mod(y,pol);
140  al1 = alginit(rnfinit(nf,x^2-(1+jj)), [-x, 4*jj+5]);
141  al2 = alginit(rnfinit(nf,x^3-2), [jj*x, 7]);
142  al3 = alginit(rnfinit(nf,x^4+x^3+x^2+x+1), [x^2, 7]);
143  do("radical 1", ()->algradical(al1)==0); \\cyclic => simple
144  do("radical 2", ()->algradical(al2)==0); \\cyclic => simple
145  do("radical 3", ()->algradical(al3)==0); \\cyclic => simple
146  p3 = idealprimedec(nf,3)[1];
147  p5 = idealprimedec(nf,5)[1];
148  p7 = idealprimedec(nf,7)[1];
149  p7b = idealprimedec(nf,7)[2];
150
151  hf12 = [[p3,p7,p7b],Vecsmall([3,4,5])];
152  hf23 = [[p5,p7,p7b],Vecsmall([6,11,7])];
153
154  do("tensor of degree 2 and 3 no mo", ()->alsame(algtensor(al1,al2,0),hf12,Vecsmall([])));
155});
156
157rnfprimedec2(rnf,pp,nf2) = {
158  my(nf = rnf.nf,pp2,pp3);
159  pp2 = rnfidealup(rnf,pp);
160  pp3 = idealadd(nf2,pp2[1],pp2[2]);
161  for(i=3, #pp2,
162    pp3 = idealadd(nf2,pp3,pp2[i]));
163  return(idealfactor(nf2,pp3));
164};
165
166rnfprimedec(rnf,pp) = rnfprimedec2(rnf,pp,nfinit(rnf.polabs));
167
168testgwa(nf,Lpr,Ld,pl) = ()->{
169  my(rnf, d, n, nf2);
170  rnf = rnfinit(nf,nfgrunwaldwang(nf,Lpr,Ld,pl,x));
171  d = rnf.disc[1];
172  n = poldegree(rnf.pol);
173  return(1);
174  nf2 = nfinit(rnf.polabs);
175  for(i=1,#Lpr,
176    if(#mattranspose(rnfprimedec2(rnf,Lpr[i],nf2))*Ld[i] != n, return(0)));
177  return(1);
178};
179
180gw() = gusuite("Grunwald-Wang", ()->{
181  my(p2,p3,p5,p7,p11,p13,p17,pp,ppp,nf);
182  nf = nfinit(y);
183  p2 = idealprimedec(nf,2)[1];
184  p3 = idealprimedec(nf,3)[1];
185  p5 = idealprimedec(nf,5)[1];
186  pp = idealprimedec(nf,nextprime(1234))[1];
187  ppp = idealprimedec(nf,nextprime(4321))[1];
188
189  do("A quadratic over Q, 2 large inert, imaginary", testgwa(nf,[pp,ppp],Vecsmall([2,2]),Vecsmall([-1])));
190  do("A quartic over Q, 2 large inert, imaginary", testgwa(bnfinit(nf),[pp,ppp],Vecsmall([4,4]),Vecsmall([-1])));
191
192  nf = nfinit(y^2+1);
193  p2 = idealprimedec(nf,2)[1];
194  p3 = idealprimedec(nf,3)[1];
195  p5 = idealprimedec(nf,5)[1];
196  do("A : degree 4 over Q(i), local degrees [4,1,1]", testgwa(nf,[p2,p3,p5],Vecsmall([4,1,1]),Vecsmall([])));
197
198  nf = nfinit(y^2+y+1);
199  p2 = idealprimedec(nf,2)[1];
200  p3 = idealprimedec(nf,3)[1];
201  p5 = idealprimedec(nf,5)[1];
202  p7 = idealprimedec(nf,7)[1];
203  pp = idealprimedec(nf,nextprime(1248))[1];
204  ppp = idealprimedec(nf,nextprime(7531))[1];
205
206  do("A degree 3 over Q(j), local degrees [3,3] larger primes", testgwa(nf,[pp,ppp],[3,3],[]));
207
208  nf = nfinit(y^2-5);
209  p2 = idealprimedec(nf,2)[1];
210  p3 = idealprimedec(nf,3)[1];
211  p5 = idealprimedec(nf,5)[1];
212  p7 = idealprimedec(nf,7)[1];
213  p11 = idealprimedec(nf,11)[1];
214  p13 = idealprimedec(nf,13)[1];
215  p17 = idealprimedec(nf,17)[1];
216  pp = idealprimedec(nf,nextprime(1248))[1];
217  ppp = idealprimedec(nf,nextprime(4897))[1];
218
219  do("A : degree 3 over Q(sqrt(5)), local degrees [3,3] [0,0], larger primes",
220testgwa(nf,[pp,ppp],[3,3],[0,0]));
221  /* TODO check what happens with [-1,-1] */
222
223  nf = nfinit(y^2-7);
224  p2 = idealprimedec(nf,2)[1];
225  p3 = idealprimedec(nf,3)[1];
226  p5 = idealprimedec(nf,5)[1];
227  p7 = idealprimedec(nf,7)[1];
228  p11 = idealprimedec(nf,11)[1];
229  p13 = idealprimedec(nf,13)[1];
230  p17 = idealprimedec(nf,17)[1];
231
232  do("A : degree 5 over Q(sqrt(7)), local degrees [5,5,5,5,5,5,5] [0,0]",
233testgwa(nf,[p2,p3,p5,p7,p11,p13,p17],Vecsmall([5,5,5,5,5,5,5]),Vecsmall([0,0])));
234  /* TODO check what happens with [-1,-1] */
235
236  nf = nfinit(polcyclo(9,y));
237  p2 = idealprimedec(nf,2)[1];
238  p3 = idealprimedec(nf,3)[1];
239  p5 = idealprimedec(nf,5)[1];
240  p7 = idealprimedec(nf,7)[1];
241  do("A : degree 9 over Q(zeta_9), local degrees [9,9,9,9]", testgwa(nf,[p2,p3,p5,p7],Vecsmall([9,9,9,9]),Vecsmall([])));
242
243  nf = nfinit(y^6 -y^5 -7*y^4 +2*y^3 +7*y^2 -2*y -1);
244  p2 = idealprimedec(nf,2)[1];
245  p3 = idealprimedec(nf,3)[1];
246  p5 = idealprimedec(nf,5)[1];
247  p7 = idealprimedec(nf,7)[1];
248  p11 = idealprimedec(nf,11)[1];
249  p13 = idealprimedec(nf,13)[1];
250  p17 = idealprimedec(nf,17)[1];
251  pp = idealprimedec(nf,nextprime(1357))[1];
252  ppp = idealprimedec(nf,nextprime(853))[1];
253
254  do("A degree 2 over totally real sextic, local degrees [2,2] [2,2,2,2,2,2], larger primes", testgwa(nf,[pp,ppp],Vecsmall([2,2]),Vecsmall([-1,-1,-1,-1,-1,-1])));
255
256  do("A degree 2 over totally real sextic, local degrees [] [2,2,2,2,2,2]", testgwa(nf,[],Vecsmall([]),Vecsmall([-1,-1,-1,-1,-1,-1])));
257});
258
259moreoperations() = gusuite("more operations", ()->{
260  my(x='x, y='y, nf, p1, p2, iinv, finv, al, u, t, b, w, un, pol1, pol2, cp,
261    mul, rnf, aut, ww, tt, ord, invord, Y = varhigher("Y",y));
262
263  pol1 = y;
264  nf = nfinit(pol1);
265  p1 = idealprimedec(nf,2)[1];
266  p2 = idealprimedec(nf,3)[1];
267  iinv = [0];
268  finv = [[p1,p2],[-2/3,-1/3]];
269  setrand(3);
270  al = alginit(nf,[3,finv,iinv],x);
271
272  do("construct algebra", ()->al);
273
274  pol2 = algsplittingfield(al).pol;
275  mul = Mod(1,pol1)*Mod(1,pol2);
276  u = [0,1,0]~*mul;
277  t = [x,0,0]~*mul;
278  b = [algb(al),0,0]~*mul;
279  w = [x+x^2-7,1/2+x-3/7*x^2,12*x-1]~*mul;
280  un = [1,0,0]~*mul;
281
282  do("norm(u)", ()->lift(algnorm(al,u))==lift(algb(al)));
283  do("norm(t)", ()->algnorm(al,t)==rnfeltnorm(algsplittingfield(al),x));
284  do("trace(u)", ()->algtrace(al,u)==0);
285  do("trace(t)", ()->algtrace(al,t)==rnfelttrace(algsplittingfield(al),x));
286  do("u+t", ()->algadd(al,u,t)==algadd(al,t,u));
287  do("u*t", ()->algmul(al,u,t)!=algmul(al,t,u));
288  do("u^3", ()->algpow(al,u,3)==b);
289  do("w^-1 L", ()->algmul(al,w,alginv(al,w))==un);
290  do("w^-1 R", ()->algmul(al,alginv(al,w),w)==un);
291  do("w^-1*u", ()->algmul(al,w,algdivl(al,w,u)));
292  do("u*w^-1", ()->algmul(al,algdivr(al,u,w),w));
293  cp = algcharpoly(al,w,Y);
294  do("charpol(w)", ()->cp);
295  do("eval charpol", ()->algadd(al,algadd(al,algpow(al,w,3),algmul(al,[polcoeff(cp,2),0,0]~*mul,algsqr(al,w))),algadd(al,algmul(al,[polcoeff(cp,1),0,0]~*mul,w),[polcoeff(cp,0),0,0]~*mul))==0);
296  do("trace(w)", ()->algtrace(al,w)==-polcoeff(cp,2));
297  do("norm(w)", ()->algnorm(al,w)==-polcoeff(cp,0));
298  do("dim", ()->algdim(al)==9);
299  do("absdim", ()->algdim(al,1)==9);
300  do("iscommutative", ()->algiscommutative(al)==0);
301  do("issemisimple", ()->algissemisimple(al)==1);
302  do("issimple", ()->algissimple(al)==1);
303
304  pol1 = y^2+1;
305  nf = nfinit(pol1);
306  pol2 = x^3 + x^2 - 2*x - 1;
307  rnf = rnfinit(nf,pol2);
308  aut = x^2-2;
309  al = alginit(rnf,[aut,Mod(2,pol1)]);
310  mul = Mod(1,pol1)*Mod(1,pol2);
311  u = [0,1,0]~*mul;
312  t = [x,0,0]~*mul;
313  tt = [y*x,0,0]~*mul;
314  b = [2,0,0]~*mul;
315  un = [1,0,0]~*mul;
316  w = [y*x-1/3*x^2, 2+y/7-x, -12*x-3*y*x^2]~*mul;
317  ww = [-x^2*y, 1/13+y+x+4*x^2, (-2+y)*x^2+(y/3+1/5)*x]~*mul;
318  ord = algbasis(al);
319  invord = alginvbasis(al);
320  do("algleftmultable w+ww", ()->algtomatrix(al,algadd(al,w,ww),1)==(algtomatrix(al,w,1)+algtomatrix(al,ww,1)));
321  do("algleftmultable w*ww", ()->algtomatrix(al,algmul(al,w,ww),1)==(algtomatrix(al,w,1)*algtomatrix(al,ww,1)));
322  do("alg(basis(w))", ()->algbasistoalg(al,algalgtobasis(al,w))==w);
323  do("alg(basis(ww))", ()->algbasistoalg(al,algalgtobasis(al,ww))==ww);
324  do("basis(w)+ww", ()->algadd(al,algalgtobasis(al,w),ww)==algalgtobasis(al,algadd(al,w,ww)));
325  do("basis(w)-ww", ()->algsub(al,algalgtobasis(al,w),ww)==algalgtobasis(al,algsub(al,w,ww)));
326  do("w+basis(ww)", ()->algadd(al,w,algalgtobasis(al,ww))==algalgtobasis(al,algadd(al,w,ww)));
327  do("w-basis(ww)", ()->algsub(al,w,algalgtobasis(al,ww))==algalgtobasis(al,algsub(al,w,ww)));
328  do("basis(w)*ww", ()->algmul(al,algalgtobasis(al,w),ww)==algalgtobasis(al,algmul(al,w,ww)));
329  do("w*basis(ww)", ()->algmul(al,w,algalgtobasis(al,ww))==algalgtobasis(al,algmul(al,w,ww)));
330  do("basis(w)^2", ()->algsqr(al,algalgtobasis(al,w))==algalgtobasis(al,algsqr(al,w)));
331  do("basis(ww)^2", ()->algsqr(al,algalgtobasis(al,ww))==algalgtobasis(al,algsqr(al,ww)));
332  do("basis(w)\\ww", ()->algdivl(al,algalgtobasis(al,w),ww)==algalgtobasis(al,algdivl(al,w,ww)));
333  do("w\\basis(ww)", ()->algdivl(al,w,algalgtobasis(al,ww))==algalgtobasis(al,algdivl(al,w,ww)));
334  do("basis(ww)\\w", ()->algdivl(al,algalgtobasis(al,ww),w)==algalgtobasis(al,algdivl(al,ww,w)));
335  do("ww\basis(w)", ()->algdivl(al,ww,algalgtobasis(al,w))==algalgtobasis(al,algdivl(al,ww,w)));
336  do("basis(w)^-1", ()->alginv(al,algalgtobasis(al,w))==algalgtobasis(al,alginv(al,w)));
337  do("basis(ww)^-1", ()->alginv(al,algalgtobasis(al,ww))==algalgtobasis(al,alginv(al,ww)));
338  do("basis(w)/ww", ()->algdivr(al,algalgtobasis(al,w),ww)==algalgtobasis(al,algdivr(al,w,ww)));
339  do("w/basis(ww)", ()->algdivr(al,w,algalgtobasis(al,ww))==algalgtobasis(al,algdivr(al,w,ww)));
340  do("basis(ww)/w", ()->algdivr(al,algalgtobasis(al,ww),w)==algalgtobasis(al,algdivr(al,ww,w)));
341  do("ww/basis(w)", ()->algdivr(al,ww,algalgtobasis(al,w))==algalgtobasis(al,algdivr(al,ww,w)));
342  do("trace(basis(w))", ()->algtrace(al,w)==algtrace(al,algalgtobasis(al,w)));
343  do("trace(basis(ww))", ()->algtrace(al,ww)==algtrace(al,algalgtobasis(al,ww)));
344  w = [0,0,x*y]~*mul;
345  ww = [1+y,1+x,1+x^2]~*mul;
346  do("alg(basis(w)) 2", ()->algbasistoalg(al,algalgtobasis(al,w))==w);
347  do("alg(basis(ww)) 2", ()->algbasistoalg(al,algalgtobasis(al,ww))==ww);
348  do("basis(w)+ww 2", ()->algadd(al,algalgtobasis(al,w),ww)==algalgtobasis(al,algadd(al,w,ww)));
349  do("basis(w)-ww 2", ()->algsub(al,algalgtobasis(al,w),ww)==algalgtobasis(al,algsub(al,w,ww)));
350  do("w+basis(ww) 2", ()->algadd(al,w,algalgtobasis(al,ww))==algalgtobasis(al,algadd(al,w,ww)));
351  do("w-basis(ww) 2", ()->algsub(al,w,algalgtobasis(al,ww))==algalgtobasis(al,algsub(al,w,ww)));
352  do("basis(w)*ww 2", ()->algmul(al,algalgtobasis(al,w),ww)==algalgtobasis(al,algmul(al,w,ww)));
353  do("w*basis(ww) 2", ()->algmul(al,w,algalgtobasis(al,ww))==algalgtobasis(al,algmul(al,w,ww)));
354  do("basis(w)^2 2", ()->algsqr(al,algalgtobasis(al,w))==algalgtobasis(al,algsqr(al,w)));
355  do("basis(ww)^2 2", ()->algsqr(al,algalgtobasis(al,ww))==algalgtobasis(al,algsqr(al,ww)));
356  do("basis(w)\ww 2", ()->algdivl(al,algalgtobasis(al,w),ww)==algalgtobasis(al,algdivl(al,w,ww)));
357  do("w\basis(ww) 2", ()->algdivl(al,w,algalgtobasis(al,ww))==algalgtobasis(al,algdivl(al,w,ww)));
358  do("basis(ww)\w 2", ()->algdivl(al,algalgtobasis(al,ww),w)==algalgtobasis(al,algdivl(al,ww,w)));
359  do("ww\basis(w) 2", ()->algdivl(al,ww,algalgtobasis(al,w))==algalgtobasis(al,algdivl(al,ww,w)));
360  do("basis(w)^-1 2", ()->alginv(al,algalgtobasis(al,w))==algalgtobasis(al,alginv(al,w)));
361  do("basis(ww)^-1 2", ()->alginv(al,algalgtobasis(al,ww))==algalgtobasis(al,alginv(al,ww)));
362  do("basis(w)/ww 2", ()->algdivr(al,algalgtobasis(al,w),ww)==algalgtobasis(al,algdivr(al,w,ww)));
363  do("w/basis(ww) 2", ()->algdivr(al,w,algalgtobasis(al,ww))==algalgtobasis(al,algdivr(al,w,ww)));
364  do("basis(ww)/w 2", ()->algdivr(al,algalgtobasis(al,ww),w)==algalgtobasis(al,algdivr(al,ww,w)));
365  do("ww/basis(w) 2", ()->algdivr(al,ww,algalgtobasis(al,w))==algalgtobasis(al,algdivr(al,ww,w)));
366  do("trace(basis(w)) 2", ()->algtrace(al,w)==algtrace(al,algalgtobasis(al,w)));
367  do("trace(basis(ww)) 2", ()->algtrace(al,ww)==algtrace(al,algalgtobasis(al,ww)));
368  w = [1/2,1/3*x,1/5]~*mul;
369  ww = [1+y,1+x,1+x^2]~*mul;
370  do("alg(basis(w)) 3", ()->algbasistoalg(al,algalgtobasis(al,w))==w);
371  do("alg(basis(ww)) 3", ()->algbasistoalg(al,algalgtobasis(al,ww))==ww);
372  do("basis(w)+ww 3", ()->algadd(al,algalgtobasis(al,w),ww)==algalgtobasis(al,algadd(al,w,ww)));
373  do("basis(w)-ww 3", ()->algsub(al,algalgtobasis(al,w),ww)==algalgtobasis(al,algsub(al,w,ww)));
374  do("w+basis(ww) 3", ()->algadd(al,w,algalgtobasis(al,ww))==algalgtobasis(al,algadd(al,w,ww)));
375  do("w-basis(ww) 3", ()->algsub(al,w,algalgtobasis(al,ww))==algalgtobasis(al,algsub(al,w,ww)));
376  do("basis(w)*ww 3", ()->algmul(al,algalgtobasis(al,w),ww)==algalgtobasis(al,algmul(al,w,ww)));
377  do("w*basis(ww) 3", ()->algmul(al,w,algalgtobasis(al,ww))==algalgtobasis(al,algmul(al,w,ww)));
378  do("basis(w)^2 3", ()->algsqr(al,algalgtobasis(al,w))==algalgtobasis(al,algsqr(al,w)));
379  do("basis(ww)^2 3", ()->algsqr(al,algalgtobasis(al,ww))==algalgtobasis(al,algsqr(al,ww)));
380  do("basis(w)\ww 3", ()->algdivl(al,algalgtobasis(al,w),ww)==algalgtobasis(al,algdivl(al,w,ww)));
381  do("w\basis(ww) 3", ()->algdivl(al,w,algalgtobasis(al,ww))==algalgtobasis(al,algdivl(al,w,ww)));
382  do("basis(ww)\w 3", ()->algdivl(al,algalgtobasis(al,ww),w)==algalgtobasis(al,algdivl(al,ww,w)));
383  do("ww\basis(w) 3", ()->algdivl(al,ww,algalgtobasis(al,w))==algalgtobasis(al,algdivl(al,ww,w)));
384  do("basis(w)^-1 3", ()->alginv(al,algalgtobasis(al,w))==algalgtobasis(al,alginv(al,w)));
385  do("basis(ww)^-1 3", ()->alginv(al,algalgtobasis(al,ww))==algalgtobasis(al,alginv(al,ww)));
386  do("basis(w)/ww 3", ()->algdivr(al,algalgtobasis(al,w),ww)==algalgtobasis(al,algdivr(al,w,ww)));
387  do("w/basis(ww) 3", ()->algdivr(al,w,algalgtobasis(al,ww))==algalgtobasis(al,algdivr(al,w,ww)));
388  do("basis(ww)/w 3", ()->algdivr(al,algalgtobasis(al,ww),w)==algalgtobasis(al,algdivr(al,ww,w)));
389  do("ww/basis(w) 3", ()->algdivr(al,ww,algalgtobasis(al,w))==algalgtobasis(al,algdivr(al,ww,w)));
390  do("trace(basis(w)) 3", ()->algtrace(al,w)==algtrace(al,algalgtobasis(al,w)));
391  do("trace(basis(ww)) 3", ()->algtrace(al,ww)==algtrace(al,algalgtobasis(al,ww)));
392  do("radical", ()->algradical(al)==0); \\cyclic => simple
393  do("iscommutative cyc 3", ()->algiscommutative(al)==0);
394  do("issemisimple cyc 3", ()->algissemisimple(al)==1);
395  do("issimple cyc 3", ()->algissimple(al)==1);
396
397  pol1 = y;
398  nf = nfinit(pol1);
399  pol2 = x^2-2;
400  rnf = rnfinit(nf,pol2);
401  aut = -x;
402  al = alginit(rnf,[aut,Mod(5,pol1)]);
403  mul = Mod(1,pol1)*Mod(1,pol2);
404  u = [0,1]~*mul;
405  t = [x,0]~*mul;
406  b = [5,0]~*mul;
407  un = [1,0]~*mul;
408  w = [-1/3*x^2, 2/7-x]~*mul;
409  ww = [-x^2*4, 1/13+x+4*x^2]~*mul;
410  ord = algbasis(al);
411  invord = alginvbasis(al);
412  do("algleftmultable/Q w+ww", ()->algtomatrix(al,algadd(al,w,ww),1)==(algtomatrix(al,w,1)+algtomatrix(al,ww,1)));
413  do("algleftmultable/Q w*ww", ()->algtomatrix(al,algmul(al,w,ww),1)==(algtomatrix(al,w,1)*algtomatrix(al,ww,1)));
414  do("alg(basis(w))/Q", ()->algbasistoalg(al,algalgtobasis(al,w))==w);
415  do("alg(basis(ww))/Q", ()->algbasistoalg(al,algalgtobasis(al,ww))==ww);
416  do("basis(w)+ww/Q", ()->algadd(al,algalgtobasis(al,w),ww)==algalgtobasis(al,algadd(al,w,ww)));
417  do("basis(w)-ww/Q", ()->algsub(al,algalgtobasis(al,w),ww)==algalgtobasis(al,algsub(al,w,ww)));
418  do("w+basis(ww)/Q", ()->algadd(al,w,algalgtobasis(al,ww))==algalgtobasis(al,algadd(al,w,ww)));
419  do("w-basis(ww)/Q", ()->algsub(al,w,algalgtobasis(al,ww))==algalgtobasis(al,algsub(al,w,ww)));
420  do("basis(w)*ww/Q", ()->algmul(al,algalgtobasis(al,w),ww)==algalgtobasis(al,algmul(al,w,ww)));
421  do("w*basis(ww)/Q", ()->algmul(al,w,algalgtobasis(al,ww))==algalgtobasis(al,algmul(al,w,ww)));
422  do("basis(w)^2/Q", ()->algsqr(al,algalgtobasis(al,w))==algalgtobasis(al,algsqr(al,w)));
423  do("basis(ww)^2/Q", ()->algsqr(al,algalgtobasis(al,ww))==algalgtobasis(al,algsqr(al,ww)));
424  do("basis(w)\ww/Q", ()->algdivl(al,algalgtobasis(al,w),ww)==algalgtobasis(al,algdivl(al,w,ww)));
425  do("w\basis(ww)/Q", ()->algdivl(al,w,algalgtobasis(al,ww))==algalgtobasis(al,algdivl(al,w,ww)));
426  do("basis(ww)\w/Q", ()->algdivl(al,algalgtobasis(al,ww),w)==algalgtobasis(al,algdivl(al,ww,w)));
427  do("ww\basis(w)/Q", ()->algdivl(al,ww,algalgtobasis(al,w))==algalgtobasis(al,algdivl(al,ww,w)));
428  do("basis(w)^-1/Q", ()->alginv(al,algalgtobasis(al,w))==algalgtobasis(al,alginv(al,w)));
429  do("basis(ww)^-1/Q", ()->alginv(al,algalgtobasis(al,ww))==algalgtobasis(al,alginv(al,ww)));
430  do("basis(w)/ww/Q", ()->algdivr(al,algalgtobasis(al,w),ww)==algalgtobasis(al,algdivr(al,w,ww)));
431  do("w/basis(ww)/Q", ()->algdivr(al,w,algalgtobasis(al,ww))==algalgtobasis(al,algdivr(al,w,ww)));
432  do("basis(ww)/w/Q", ()->algdivr(al,algalgtobasis(al,ww),w)==algalgtobasis(al,algdivr(al,ww,w)));
433  do("ww/basis(w)/Q", ()->algdivr(al,ww,algalgtobasis(al,w))==algalgtobasis(al,algdivr(al,ww,w)));
434  do("trace(basis(w))/Q", ()->algtrace(al,w)==algtrace(al,algalgtobasis(al,w)));
435  do("trace(basis(ww))/Q", ()->algtrace(al,ww)==algtrace(al,algalgtobasis(al,ww)));
436  do("radical/Q", ()->algradical(al)==0); \\cyclic => simple
437  do("iscommutative /Q", ()->algiscommutative(al)==0);
438  do("issemisimple /Q", ()->algissemisimple(al)==1);
439  do("issimple /Q", ()->algissimple(al)==1);
440});
441
442tablealg() = gusuite("table algebra", ()->{
443  my(x='x, al, mt, p, un, a, b, c, d, e, ss, projm, liftm, pa, sc);
444  mt = [[1,0,0;0,1,0;0,0,1],[0,0,0;1,0,1;0,0,0],[0,0,0;0,0,0;1,0,1]]; \\Matrices [*,*;0,*]
445  un = [1,0,0]~;
446  a = [1,0,-1]~;
447  b = [0,-1,1]~;
448  do("algisassociative 0.0", ()->algisassociative(mt));
449  do("algisassociative 0.1", ()->algisassociative(mt[2..3]));
450  my (mt0 = mt);
451  mt0[3][3,1] = 1;
452  do("algisassociative 0.2", ()->algisassociative(mt0));
453  do("algisassociative 0.3", ()->algisassociative('x));
454  al = algtableinit(mt,0);
455  do("construction 0", ()->al);
456  do("iscyclic 0", ()->algtype(al)==1);
457  do("dim 0", ()->algdim(al,1)==3);
458  do("dim 0b", ()->algdim(al)==3);
459  do("char 0", ()->algchar(al)==0);
460  do("a+b 0", ()->algadd(al,a,b)==[1,-1,0]~);
461  do("a-b 0", ()->algsub(al,a,b)==[1,1,-2]~);
462  do("a*b 0", ()->algmul(al,a,b)==[0,-1,0]~);
463  do("b*a 0", ()->algmul(al,b,a)==[0,0,0]~);
464  do("a^2 0", ()->algsqr(al,a)==a);
465  do("b^2 0", ()->algsqr(al,b)==b);
466  e = [1,1,0]~;
467  do("e^691691 0", ()->algpow(al,e,691691)==[1,691691,0]~);
468  d = [1,0,1]~;
469  do("d^101 0", ()->algpow(al,d,101)==[1,0,2^101-1]~);
470  do("multable(a) 0", ()->algtomatrix(al,a,1)==[1,0,0;0,1,0;-1,0,0]);
471  do("multable(b) 0", ()->algtomatrix(al,b,1)==[0,0,0;-1,0,-1;1,0,1]);
472  do("divl(d,a) 0", ()->algdivl(al,d,a)==a);
473  do("divl(d,b) 0", ()->algdivl(al,d,b)==[0,-1,1/2]~);
474  do("d^-1 0", ()->alginv(al,d)==[1,0,-1/2]~);
475  do("divr(a,d) 0", ()->algdivr(al,a,d)==a);
476  do("divr(b,d) 0", ()->algdivr(al,b,d)==[0,-1/2,1/2]~);
477  c = [0,7,0]~;
478  do("rad(al) 0", ()->#algradical(al)==1); \\matrices [0,*;0,0]
479  do("ss(al) 0", ()->#algradical(algquotient(al,algradical(al)))==0);
480  [ss,projm,liftm] = algquotient(al,algradical(al),1);
481  pa = projm*a;
482  do("proj(a) idem 0", ()->algsqr(ss,pa)==pa);
483  do("idemproj 0", ()->algcentralproj(ss,[pa,algsub(ss,projm*un,pa)]));
484  sc = algcentralproj(ss,[pa,algsub(ss,projm*un,pa)]);
485  do("simple components 0", ()->algmultable(sc[1])==[Mat(1)] && algmultable(sc[2])==[Mat(1)]);
486  do("center al 0", ()->#algcenter(al)==1);
487  do("center ss 0", ()->#algcenter(ss)==2);
488  do("primesubalg ss 0", ()->#algprimesubalg(ss)==-1);
489  do("charpol annihil(a) 0", ()->testcharpol(al,a));
490  do("charpol annihil(b) 0", ()->testcharpol(al,b));
491  do("charpol annihil(c) 0", ()->testcharpol(al,c));
492  do("charpol annihil(d) 0", ()->testcharpol(al,d));
493  do("charpol annihil(e) 0", ()->testcharpol(al,e));
494  do("random 0", ()->algrandom(al,1));
495  do("algsimpledec 0", ()->#algsimpledec(ss)[2]==2);
496  do("alg_decomposition 0", ()->dec=algsimpledec(al); #dec[1]==1 && #dec[2]==2);
497  do("iscommutative 0", ()->algiscommutative(al)==0);
498  do("issemisimple 0", ()->algissemisimple(al)==0);
499  do("issimple 0", ()->algissimple(al)==0);
500  do("issimple ss 0", ()->algissimple(ss)==0);
501  do("isdivision 0", ()->algisdivision(al)==0);
502
503  p = 2;
504  al = algtableinit(mt,p);
505  do("algisassociative 2", ()->algisassociative(mt,p));
506  do("construction 2", ()->al);
507  do("iscyclic 2", ()->algtype(al)==1);
508  do("dim 2", ()->algdim(al,1)==3);
509  do("char 2", ()->algchar(al)==p);
510  do("a+b 2", ()->algadd(al,a,b)==[1,1,0]~);
511  do("a-b 2", ()->algsub(al,a,b)==algadd(al,a,b));
512  do("a*b 2", ()->algmul(al,a,b)==[0,p-1,0]~);
513  do("b*a 2", ()->algmul(al,b,a)==[0,0,0]~);
514  do("a^2 2", ()->algsqr(al,a)==a*Mod(1,p));
515  do("b^2 2", ()->algsqr(al,b)==b*Mod(1,p));
516  do("multable(a) 2", ()->algtomatrix(al,a,1)==[1,0,0;0,1,0;-1,0,0]*Mod(1,p));
517  do("multable(b) 2", ()->algtomatrix(al,b,1)==[0,0,0;-1,0,-1;1,0,1]*Mod(1,p));
518  do("divl(un,a) 2", ()->algdivl(al,un,a)==a*Mod(1,p));
519  do("divl(un,b) 2", ()->algdivl(al,un,b)==b*Mod(1,p));
520  do("un^-1 2", ()->alginv(al,un)==un);
521  do("divr(a,un) 2", ()->algdivr(al,a,un)==a*Mod(1,p));
522  do("divr(b,un) 2", ()->algdivr(al,b,un)==b*Mod(1,p));
523  do("rad(al) 2", ()->#algradical(al)==1); \\matrices [0,*;0,0]
524  do("ss(al) 2", ()->#algradical(algquotient(al,algradical(al)))==0);
525  [ss,projm,liftm] = algquotient(al,algradical(al),1);
526  pa = projm*a;
527  do("proj(a) idem 2", ()->algsqr(ss,pa)==pa*Mod(1,p));
528  do("idemproj 2", ()->algcentralproj(ss,[pa,algsub(ss,projm*un,pa)]));
529  sc = algcentralproj(ss,[pa,algsub(ss,projm*un,pa)]);
530  do("simple components 2", ()->algmultable(sc[1])==[Mat(Mod(1,p))] && algmultable(sc[2])==[Mat(Mod(1,p))]);
531  do("center al 2", ()->#algcenter(al)==1);
532  do("center ss 2", ()->#algcenter(ss)==2);
533  do("primesubalg ss 2", ()->#algprimesubalg(ss)==2);
534  do("charpol annihil(a) 2", ()->testcharpol(al,a));
535  do("charpol annihil(b) 2", ()->testcharpol(al,b));
536  do("charpol annihil(c) 2", ()->testcharpol(al,c));
537  do("random 2", ()->algrandom(al,1));
538  do("algsimpledec 2", ()->#algsimpledec(ss)[2]==2);
539  do("alg_decomposition 2", ()->dec=algsimpledec(al); #dec[1]==1 && #dec[2]==2);
540  do("iscommutative 2", ()->algiscommutative(al)==0);
541  do("issemisimple 2", ()->algissemisimple(al)==0);
542  do("issimple 2", ()->algissimple(al)==0);
543  do("issimple ss 2", ()->algissimple(ss)==0);
544  do("matrix trace 2", ()->algtrace(al,[un,vector(3)~;vector(3)~,un])==0);
545  do("matrix norm 2", ()->algnorm(al,[un,vector(3)~;vector(3)~,un])==1);
546  do("norm 2", ()->algnorm(al,un)==1);
547
548  p = 3;
549  mt = [[1,0;0,1],[0,0;1,0]];\\F3[x]/(x^2)
550  un = [1,0]~;
551  a = [1,-1]~;
552  b = [0,1]~;
553  al = algtableinit(mt,p);
554  do("construction 3", ()->al);
555  do("iscyclic 3", ()->algtype(al)==1);
556  do("dim 3", ()->algdim(al,1)==2);
557  do("char 3", ()->algchar(al)==p);
558  do("a+b 3", ()->algadd(al,a,b)==un);
559  do("a-b 3", ()->algsub(al,a,b)==[1,1]~);
560  do("a*b 3", ()->algmul(al,a,b)==[0,1]~);
561  do("b*a 3", ()->algmul(al,b,a)==[0,1]~);
562  do("a^2 3", ()->algsqr(al,a)==[1,1]~);
563  do("b^2 3", ()->algsqr(al,b)==[0,0]~);
564  do("a^691691 3", ()->algpow(al,a,691691)==[1,-691691]~*Mod(1,p));
565  do("multable(a) 3", ()->algtomatrix(al,a,1)==[1,0;-1,1]*Mod(1,p));
566  do("multable(b) 3", ()->algtomatrix(al,b,1)==[0,0;1,0]);
567  do("algdivl(a,b) 3", ()->algdivl(al,a,b)==b);
568  do("a^-1 3", ()->alginv(al,a)==[1,1]~);
569  do("algdivr(b,a) 3", ()->algdivr(al,b,a)==b);
570  do("rad(al) 3", ()->#algradical(al)==1); \\ideal (x)
571  do("ss(al) 3", ()->#algradical(algquotient(al,algradical(al)))==0);
572  [ss,projm,liftm] = algquotient(al,algradical(al),1);
573  do("center al 3", ()->#algcenter(al)==2);
574  do("center ss 3", ()->#algcenter(ss)==1);
575  do("primesubalg ss 3", ()->#algprimesubalg(ss)==1);
576  do("charpol annihil(a) 3", ()->testcharpol(al,a));
577  do("charpol annihil(b) 3", ()->testcharpol(al,b));
578  do("random 3", ()->algrandom(al,1));
579  do("algsimpledec 3", ()->#algsimpledec(ss)[2]==1);
580  do("alg_decomposition 3", ()->dec=algsimpledec(al); #dec[1]==1 && #dec[2]==1);
581  do("iscommutative 3", ()->algiscommutative(al)==1);
582  do("issemisimple 3", ()->algissemisimple(al)==0);
583  do("issemisimple ss 3", ()->algissemisimple(ss)==1);
584  do("issimple 3", ()->algissimple(al)==0);
585  do("issimple ss 3", ()->algissimple(ss)==1);
586
587  p = 3;
588  mt = [[1,0,0;0,1,0;0,0,1],[0,0,0;1,0,0;0,1,0],[0,0,0;0,0,0;1,0,0]];\\F3[x]/(x^3)
589  un = [1,0,0]~;
590  a = [1,-1,0]~;
591  b = [0,1,0]~;
592  al = algtableinit(mt,p);
593  do("construction 3c", ()->al);
594  do("iscyclic 3c", ()->algtype(al)==1);
595  do("dim 3c", ()->algdim(al,1)==3);
596  do("char 3c", ()->algchar(al)==p);
597  do("a+b 3c", ()->algadd(al,a,b)==un);
598  do("a-b 3c", ()->algsub(al,a,b)==[1,1,0]~);
599  do("a*b 3c", ()->algmul(al,a,b)==[0,1,p-1]~);
600  do("b*a 3c", ()->algmul(al,b,a)==[0,1,p-1]~);
601  do("a^2 3c", ()->algsqr(al,a)==[1,1,1]~);
602  do("b^2 3c", ()->algsqr(al,b)==[0,0,1]~);
603  do("a^691691 3c", ()->algpow(al,a,691691)==[1,-691691,(691691*691690)\2]~*Mod(1,p));
604  do("multable(a) 3c", ()->algtomatrix(al,a,1)==[1,0,0;-1,1,0;0,-1,1]*Mod(1,p));
605  do("multable(b) 3c", ()->algtomatrix(al,b,1)==[0,0,0;1,0,0;0,1,0]);
606  do("algdivl(a,b) 3c", ()->algdivl(al,a,b)==[0,1,1]~);
607  do("a^-1 3c", ()->alginv(al,a)==[1,1,1]~);
608  do("algdivr(b,a) 3c", ()->algdivr(al,b,a)==[0,1,1]~);
609  do("rad(al) 3c", ()->#algradical(al)==2); \\ideal (x), basis (x,x^2)
610  do("ss(al) 3c", ()->#algradical(algquotient(al,algradical(al)))==0);
611  [ss,projm,liftm] = algquotient(al,algradical(al),1);
612  do("center al 3c", ()->#algcenter(al)==3);
613  do("center ss 3c", ()->#algcenter(ss)==1);
614  do("primesubalg ss 3c", ()->#algprimesubalg(ss)==1);
615  do("charpol annihil(a) 3c", ()->testcharpol(al,a));
616  do("charpol annihil(b) 3c", ()->testcharpol(al,b));
617  do("random 3c", ()->algrandom(al,1));
618  do("algsimpledec 3c", ()->#algsimpledec(ss)[2]==1);
619  do("alg_decomposition 3c", ()->dec=algsimpledec(al); #dec[1]==2 && #dec[2]==1);
620  do("iscommutative 3c", ()->algiscommutative(al)==1);
621  do("issemisimple 3c", ()->algissemisimple(al)==0);
622  do("issemisimple ss 3c", ()->algissemisimple(ss)==1);
623  do("issimple 3c", ()->algissimple(al)==0);
624  do("issimple ss 3c", ()->algissimple(ss)==1);
625
626  p = 2;
627  mt = [[1,0;0,1],[0,1;1,1]]; \\F2[x]/(x^2+x+1)
628  un = [1,0]~;
629  a = [0,1]~;
630  b = [1,1]~;
631  al = algtableinit(mt,p);
632  do("construction 2b", ()->al);
633  do("iscyclic 2b", ()->algtype(al)==1);
634  do("dim 2b", ()->algdim(al,1)==2);
635  do("char 2b", ()->algchar(al)==p);
636  do("a+b 2b", ()->algadd(al,a,b)==un);
637  do("a-b 2b", ()->algsub(al,a,b)==un);
638  do("a*b 2b", ()->algmul(al,a,b)==un);
639  do("b*a 2b", ()->algmul(al,b,a)==un);
640  do("a^2 2b", ()->algsqr(al,a)==b);
641  do("b^2 2b", ()->algsqr(al,b)==a);
642  do("a^691691 2b", ()->algpow(al,a,691691)==b);
643  do("multable(a) 2b", ()->algtomatrix(al,a,1)==[0,1;1,1]);
644  do("multable(b) 2b", ()->algtomatrix(al,b,1)==[1,1;1,0]);
645  do("divl(a,b) 2b", ()->algdivl(al,a,b)==a);
646  do("a^-1 2b", ()->alginv(al,a)==b);
647  do("divr(b,a) 2b", ()->algdivr(al,b,a)==a);
648  do("rad(al) 2b", ()->#algradical(al)==0); \\separable extension of F2
649  do("center al 2b", ()->#algcenter(al)==2);
650  do("primesubalg al 2b", ()->#algprimesubalg(al)==1);
651  do("charpol annihil(a) 2b", ()->testcharpol(al,a));
652  do("charpol annihil(b) 2b", ()->testcharpol(al,b));
653  do("random 2b", ()->algrandom(al,1));
654  do("algsimpledec 2b", ()->#algsimpledec(al)[2]==1);
655  do("alg_decomposition 2b", ()->dec=algsimpledec(al); dec[1]==0 && #dec[2]==1 && algdim(dec[2][1],1)==2);
656  do("iscommutative 2b", ()->algiscommutative(al)==1);
657  do("issemisimple 2b", ()->algissemisimple(al)==1);
658  do("issimple 2b", ()->algissimple(al)==1);
659  do("issimple,1 2b", ()->algissimple(al,1)==1);
660
661
662  p = 3;
663  mt = [matid(4),
664
665         [0,1,0,0;
666          1,0,0,0;
667          0,0,1,0;
668          0,0,0,-1],
669
670         [0,0,0,1/2;
671          0,0,0,1/2;
672          1,-1,0,0;
673          0,0,0,0],
674
675         [0,0,1/2,0;
676          0,0,-1/2,0;
677          0,0,0,0;
678          1,1,0,0]]*Mod(1,p); \\M_2(F3)
679  un = [1,0,0,0]~;
680  a = [0,1,-1,0]~;
681  b = [1,1,0,1]~;
682  al = algtableinit(mt,p);
683  do("construction 3b", ()->al);
684  do("iscyclic 3b", ()->algtype(al)==1);
685  do("dim 3b", ()->algdim(al,1)==4);
686  do("char 3b", ()->algchar(al)==p);
687  do("a+b 3b", ()->algadd(al,a,b)==[1,-1,2,1]~*Mod(1,p));
688  do("a-b 3b", ()->algsub(al,a,b)==[2,0,2,2]~);
689  do("a*b 3b", ()->algmul(al,a,b)==[2,2,0,2]~);
690  do("b*a 3b", ()->algmul(al,b,a)==[2,0,1,1]~);
691  do("a^2 3b", ()->algsqr(al,a)==un);
692  do("b^2 3b", ()->algsqr(al,b)==-b*Mod(1,p));
693  do("a^691691 3b", ()->algpow(al,a,691691)==a*Mod(1,p));
694  do("b^691691 3b", ()->algpow(al,b,691691)==b);
695  do("multable(a) 3b", ()->algtomatrix(al,a,1)==[0,1,0,1;1,0,0,1;2,1,1,0;0,0,0,2]);
696  do("multable(b) 3b", ()->algtomatrix(al,b,1)==[1,1,2,0;1,1,1,0;0,0,2,0;1,1,0,0]);
697  do("divl(a,b) 3b", ()->algdivl(al,a,b)==[2,2,0,2]~);
698  do("a^-1 3b", ()->alginv(al,a)==[0,1,2,0]~);
699  do("divr(b,a) 3b", ()->algdivr(al,b,a)==[2,0,1,1]~);
700  c = [0,0,1,0]~;
701  do("rad(al) 3b", ()->#algradical(al)==0); \\matrix ring is semisimple
702  do("center al 3b", ()->#algcenter(al)==1);
703  do("primesubalg al 3b", ()->#algprimesubalg(al)==1);
704  do("charpol annihil(a) 3b", ()->testcharpol(al,a));
705  do("charpol annihil(b) 3b", ()->testcharpol(al,b));
706  do("charpol annihil(c) 3b", ()->testcharpol(al,c));
707  do("random 3b", ()->algrandom(al,1));
708  do("algsimpledec 3b", ()->#algsimpledec(al)[2]==1);
709  do("alg_decomposition 3b", ()->dec=algsimpledec(al); dec[1]==0 && #dec[2]==1 && #algcenter(dec[2][1])==1);
710  do("iscommutative 3b", ()->algiscommutative(al)==0);
711  do("issemisimple 3b", ()->algissemisimple(al)==1);
712  do("issimple 3b", ()->algissimple(al)==1);
713
714  p = 2;
715  mt = [matid(4),
716
717     [0,0,1,0;
718      1,0,0,1;
719      0,0,0,0;
720      0,0,-1,0],
721
722     [0,0,0,0;
723      0,0,0,0;
724      1,0,0,0;
725      0,1,0,0],
726
727     [0,0,0,0;
728      0,0,0,0;
729      0,0,1,0;
730      1,0,0,1]]*Mod(1,p); \\M_2(F2)
731  un = [1,0,0,0]~;
732  a = [0,1,0,0]~;
733  b = [1,0,0,1]~;
734  al = algtableinit(mt,p);
735  do("construction 2c", ()->al);
736  do("iscyclic 2c", ()->algtype(al)==1);
737  do("dim 2c", ()->algdim(al,1)==4);
738  do("char 2c", ()->algchar(al)==p);
739  do("a+b 2c", ()->algadd(al,a,b)==[1,1,0,1]~);
740  do("a-b 2c", ()->algsub(al,a,b)==[1,1,0,1]~);
741  do("a*b 2c", ()->algmul(al,a,b)==[0,0,0,0]~);
742  do("b*a 2c", ()->algmul(al,b,a)==a);
743  do("a^2 2c", ()->algsqr(al,a)==[0,0,0,0]~);
744  do("b^2 2c", ()->algsqr(al,b)==b);
745  c = [1,1,1,1]~;
746  do("a^691691 2c", ()->algpow(al,a,691691)==[0,0,0,0]~);
747  do("b^691691 2c", ()->algpow(al,b,691691)==b);
748  do("c^691691 2c", ()->algpow(al,c,691691)==[0,1,1,1]~);
749  do("multable(a) 2c", ()->algtomatrix(al,a,1)==[0,0,1,0;1,0,0,1;0,0,0,0;0,0,1,0]);
750  do("multable(b) 2c", ()->algtomatrix(al,b,1)==[1,0,0,0;0,1,0,0;0,0,0,0;1,0,0,0]);
751  do("divl(c,a) 2c", ()->algdivl(al,c,a)==[0,0,0,1]~);
752  do("divl(c,b) 2c", ()->algdivl(al,c,b)==[0,0,1,0]~);
753  do("c^-1 2c", ()->alginv(al,c)==[0,1,1,1]~);
754  do("divr(a,c) 2c", ()->algdivr(al,a,c)==[1,1,0,1]~);
755  do("divr(b,c) 2c", ()->algdivr(al,b,c)==[0,1,0,0]~);
756  do("rad(al) 2c", ()->#algradical(al)==0); \\matrix ring is semisimple
757  do("center al 2c", ()->#algcenter(al)==1);
758  do("primesubalg al 2c", ()->#algprimesubalg(al)==1);
759  do("charpol annihil(a) 2c", ()->testcharpol(al,a));
760  do("charpol annihil(b) 2c", ()->testcharpol(al,b));
761  do("charpol annihil(c) 2c", ()->testcharpol(al,c));
762  do("random 2c", ()->algrandom(al,1));
763  do("algsimpledec 2c", ()->#algsimpledec(al)[2]==1);
764  do("alg_decomposition 2c", ()->dec=algsimpledec(al); dec[1]==0 && #dec[2]==1 && #algcenter(dec[2][1])==1);
765  do("iscommutative 2c", ()->algiscommutative(al)==0);
766  do("issemisimple 2c", ()->algissemisimple(al)==1);
767  do("issimple 2c", ()->algissimple(al)==1);
768
769  p = 5;
770  mt = [Mat(Mod(1,p))];
771  un = [1]~;
772  a = [2]~;
773  b = [3]~;
774  al = algtableinit(mt,p);
775  do("construction 5", ()->al);
776  do("iscyclic 5", ()->algtype(al)==1);
777  do("dim 5", ()->algdim(al,1)==1);
778  do("char 5", ()->algchar(al)==p);
779  do("a+b 5", ()->algadd(al,a,b)==[Mod(0,p)]~);
780  do("a-b 5", ()->algsub(al,a,b)==[Mod(4,p)]~);
781  do("a*b 5", ()->algmul(al,a,b)==[Mod(1,p)]~);
782  do("b*a 5", ()->algmul(al,b,a)==[Mod(1,p)]~);
783  do("a^2 5", ()->algsqr(al,a)==[Mod(4,p)]~);
784  do("b^2 5", ()->algsqr(al,b)==[Mod(-1,p)]~);
785  do("a^691691 5", ()->algpow(al,a,691691)==b);
786  do("multable(a) 5", ()->algtomatrix(al,a,1)==Mat(Mod(2,p)));
787  do("multable(b) 5", ()->algtomatrix(al,b,1)==Mat(Mod(3,p)));
788  do("divl(a,b) 5", ()->algdivl(al,a,b)==[Mod(4,p)]~);
789  do("a^-1 5", ()->alginv(al,a)==[Mod(3,p)]~);
790  do("divr(a,b) 5", ()->algdivr(al,a,b)==[Mod(4,p)]~);
791  do("rad(al) 5", ()->#algradical(al)==0); \\F5, dim 1
792  do("center al 5", ()->#algcenter(al)==1);
793  do("primesubalg al 5", ()->#algprimesubalg(al)==1);
794  do("charpol annihil(a) 5", ()->testcharpol(al,a));
795  do("charpol annihil(b) 5", ()->testcharpol(al,b));
796  do("random 5", ()->algrandom(al,1));
797  do("algsimpledec 5", ()->#algsimpledec(al)[2]==1);
798  do("alg_decomposition 5", ()->dec=algsimpledec(al); dec[1]==0 && #dec[2]==1 && algdim(dec[2][1],1)==1);
799  do("iscommutative 5", ()->algiscommutative(al)==1);
800  do("issemisimple 5", ()->algissemisimple(al)==1);
801  do("issimple 5", ()->algissimple(al)==1);
802
803  p = 0; \\M_2(Q)+Q
804  mt = [matid(5),
805
806     [0,0,1,0,0;
807      1,0,0,1,0;
808      0,0,0,0,0;
809      0,0,-1,0,0;
810      0,1,-1,-1,1],
811
812     [0,0,0,0,0;
813      0,0,0,0,0;
814      1,0,0,0,0;
815      0,1,0,0,0;
816      0,0,0,0,0],
817
818     [0,0,0,0,0;
819      0,0,0,0,0;
820      0,0,1,0,0;
821      1,0,0,1,0;
822      0,0,0,0,0],
823
824     [0,0,0,0,0;
825      0,0,0,0,0;
826      0,0,0,0,0;
827      0,0,0,0,0;
828      1,1,0,0,1]
829  ];
830  un = [1,0,0,0,0]~;
831  a = [1,0,0,0,-1]~;
832  b = [1,1,0,0,0]~;
833  al = algtableinit(mt,p);
834  do("construction 0b", ()->al);
835  do("iscyclic 0b", ()->algtype(al)==1);
836  do("dim 0b", ()->algdim(al,1)==5);
837  do("char 0b", ()->algchar(al)==p);
838  do("a+b 0b", ()->algadd(al,a,b)==[2,1,0,0,-1]~);
839  do("a-b 0b", ()->algsub(al,a,b)==[0,-1,0,0,-1]~);
840  do("a*b 0b", ()->algmul(al,a,b)==[1,1,0,0,-2]~);
841  do("b*a 0b", ()->algmul(al,b,a)==algmul(al,a,b));\\a central
842  do("a^2 0b", ()->algsqr(al,a)==a);
843  do("b^2 0b", ()->algsqr(al,b)==[1,2,0,0,1]~);
844  do("a^691691 0b", ()->algpow(al,a,691691)==a);
845  do("b^691 0b", ()->algpow(al,b,691)==[1,691,0,0,2^691-1-691]~);
846  do("multable(a) 0b", ()->algtomatrix(al,a,1)==
847      [1,0,0,0,0;
848       0,1,0,0,0;
849       0,0,1,0,0;
850       0,0,0,1,0;
851       -1,-1,0,0,0]);
852  do("multable(b) 0b", ()->algtomatrix(al,b,1)==
853      [1,0,1,0,0;
854       1,1,0,1,0;
855       0,0,1,0,0;
856       0,0,-1,1,0;
857       0,1,-1,-1,2]);
858  do("divl(b,a) 0b", ()->algdivl(al,b,a)==[1,-1,0,0,0]~);
859  do("b^-1 0b", ()->alginv(al,b)==[1,-1,0,0,1/2]~);
860  do("divr(a,b) 0b", ()->algdivr(al,a,b)==algdivl(al,b,a));
861  do("rad(al) 0b", ()->#algradical(al)==0);
862  do("idemproj 0b", ()->algcentralproj(al,[a,algsub(al,un,a)]));
863  sc = algcentralproj(al,[a,algsub(al,un,a)]);
864  do("simple components 0b", ()->algdim(sc[1],1)==4 && algdim(sc[2],1)==1);
865  do("mt M2 component 0b", ()->algmultable(sc[1])[1]==matid(4));
866  do("center al 0b", ()->#algcenter(al)==2);
867  do("primesubalg al 0b", ()->#algprimesubalg(al)==-1);
868  do("charpol annihil(a) 0b", ()->testcharpol(al,a));
869  do("charpol annihil(b) 0b", ()->testcharpol(al,b));
870  do("random 0b", ()->algrandom(al,1));
871  do("algsimpledec 0b", ()->#algsimpledec(al)[2]==2);
872  do("alg_decomposition 0b", ()->dec=algsimpledec(al); dec[1]==0 && #dec[2]==2
873    && #algcenter(dec[2][1])==1 && #algcenter(dec[2][2])==1 &&
874    (algdim(dec[2][1],1)==4 || algdim(dec[2][2],1)==4));
875  do("subalg M2(Q)", ()->sal=algsubalg(al,[1,0,0,0; 0,1,0,0; 0,0,1,0; 0,0,0,1;\
876    0,0,0,0])[1]; algisassociative(sal) && algradical(sal)==0 &&\
877    #algsimpledec(sal)[2]==1);
878  do("iscommutative 0b", ()->algiscommutative(al)==0);
879  do("issemisimple 0b", ()->algissemisimple(al)==1);
880  do("issimple 0b", ()->algissimple(al)==0);
881
882  p = 3;
883  al = algtableinit(mt,p);
884  do("construction 3d", ()->al);
885  do("iscyclic 3d", ()->algtype(al)==1);
886  do("dim 3d", ()->algdim(al,1)==5);
887  do("char 3d", ()->algchar(al)==p);
888  do("a+b 3d", ()->algadd(al,a,b)==[2,1,0,0,-1]~*Mod(1,p));
889  do("a-b 3d", ()->algsub(al,a,b)==[0,-1,0,0,-1]~*Mod(1,p));
890  do("a*b 3d", ()->algmul(al,a,b)==[1,1,0,0,-2]~*Mod(1,p));
891  do("b*a 3d", ()->algmul(al,b,a)==algmul(al,a,b));\\a central
892  do("a^2 3d", ()->algsqr(al,a)==a*Mod(1,p));
893  do("b^2 3d", ()->algsqr(al,b)==[1,2,0,0,1]~);
894  do("a^691691 3d", ()->algpow(al,a,691691)==a*Mod(1,p));
895  do("b^691 3d", ()->algpow(al,b,691)==[1,691,0,0,2^691-1-691]~*Mod(1,p));
896  do("multable(a) 3d", ()->algtomatrix(al,a,1)==
897      [1,0,0,0,0;
898       0,1,0,0,0;
899       0,0,1,0,0;
900       0,0,0,1,0;
901       -1,-1,0,0,0]*Mod(1,p));
902  do("multable(b) 3d", ()->algtomatrix(al,b,1)==
903      [1,0,1,0,0;
904       1,1,0,1,0;
905       0,0,1,0,0;
906       0,0,-1,1,0;
907       0,1,-1,-1,2]*Mod(1,p));
908  do("divl(b,a) 3d", ()->algdivl(al,b,a)==[1,-1,0,0,0]~*Mod(1,p));
909  do("b^-1 3d", ()->alginv(al,b)==[1,-1,0,0,1/2]~*Mod(1,p));
910  do("divr(a,b) 3d", ()->algdivr(al,a,b)==algdivl(al,b,a));
911  do("rad(al) 3d", ()->#algradical(al)==0);
912  do("idemproj 3d", ()->algcentralproj(al,[a,algsub(al,un,a)]));
913  sc = algcentralproj(al,[a,algsub(al,un,a)]);
914  do("simple components 3d", ()->algdim(sc[1],1)==4 && algdim(sc[2],1)==1);
915  do("mt M2 component 3d", ()->algmultable(sc[1])[1]==matid(4));
916  do("center al 3d", ()->#algcenter(al)==2);
917  do("primesubalg al 3d", ()->#algprimesubalg(al)==2);
918  do("charpol annihil(a) 3d", ()->testcharpol(al,a));
919  do("charpol annihil(b) 3d", ()->testcharpol(al,b));
920  do("random 3d", ()->algrandom(al,1));
921  do("algsimpledec 3d", ()->#algsimpledec(al)[2]==2);
922  do("alg_decomposition 3d", ()->dec=algsimpledec(al); dec[1]==0 && #dec[2]==2
923    && #algcenter(dec[2][1])==1 && #algcenter(dec[2][2])==1 &&
924    (algdim(dec[2][1],1)==4 || algdim(dec[2][2],1)==4));
925  do("subalg M2(F3)", ()->sal=algsubalg(al,[1,0,0,0; 0,1,0,0; 0,0,1,0; 0,0,0,1;\
926    0,0,0,0]); algisassociative(sal[1]) && algradical(sal[1])==0 &&\
927    #algsimpledec(sal[1])[2]==1);
928  do("iscommutative 3d", ()->algiscommutative(al)==0);
929  do("issemisimple 3d", ()->algissemisimple(al)==1);
930  do("issimple 3d", ()->algissimple(al)==0);
931  do("issimple,1 3d", ()->algissimple(al,1)==0);
932});
933
934all() = gusuite("all", ()->{
935  get();
936  operations();
937  tensor();
938  gw();
939  moreoperations();
940  tablealg();
941});
942
943all();
944
945nf = nfinit(y^2-2);
946al = alginit(nf, [-3,-5], x,1);
947print("maxorder assoc: ", algisassociative(al[9]));
948al0 = alginit(nf, [-3,-5], x,0);
949print("natorder assoc: ", algisassociative(al0[9]));
950un = [1,0]~;
951ii = [x,0]~;
952jj = [0,1]~;
953kk = algmul(al,ii,jj);
954print("spl(1): ", algtomatrix(al,un)==matid(2));
955print("spl(i): ", algtomatrix(al,ii)==[x,0;0,-x]);
956print("spl(j): ", algtomatrix(al,jj)==[0,-5;1,0]);
957print("spl(k): ", algtomatrix(al,kk)==[0,-5*x;-x,0]);
958print("spl(basis(1)): ", algtomatrix(al,algalgtobasis(al,un))==matid(2));
959print("spl(basis(i)): ", algtomatrix(al,algalgtobasis(al,ii))==[x,0;0,-x]);
960print("spl(basis(j)): ", algtomatrix(al,algalgtobasis(al,jj))==[0,-5;1,0]);
961print("spl(basis(k)): ", algtomatrix(al,algalgtobasis(al,kk))==[0,-5*x;-x,0]);
962a = y+1;
963b = 1/3;
964c = -y/5+1/2;
965print("spl(a*1): ", algtomatrix(al,a*un)==a*matid(2));
966print("spl(a*i): ", algtomatrix(al,a*ii)==a*[x,0;0,-x]);
967print("spl(a*j): ", algtomatrix(al,a*jj)==a*[0,-5;1,0]);
968print("spl(a*k): ", algtomatrix(al,a*kk)==a*[0,-5*x;-x,0]);
969print("spl(b*1): ", algtomatrix(al,b*un)==b*matid(2));
970print("spl(b*i): ", algtomatrix(al,b*ii)==b*[x,0;0,-x]);
971print("spl(b*j): ", algtomatrix(al,b*jj)==b*[0,-5;1,0]);
972print("spl(b*k): ", algtomatrix(al,b*kk)==b*[0,-5*x;-x,0]);
973
974
975ord = algbasis(al);
976invord = alginvbasis(al);
977setrand(1);
978x1 = algrandom(al,1);
979ax1 = algbasistoalg(al,x1);
980nx1 = algalgtobasis(al0,ax1);
981print("nattomax 1: ", nx1==ord*x1);
982setrand(2);
983x2 = algrandom(al,1);
984ax2 = algbasistoalg(al,x2);
985nx2 = algalgtobasis(al0,ax2);
986print("nattomax 2: ", nx2==ord*x2);
987print("ord*invord=id: ", ord*invord == matid(8));
988print("spl additive: ", algtomatrix(al,x1) + algtomatrix(al,x2) == algtomatrix(al, algadd(al,x1,x2)));
989print("spl multiplicative: ", algtomatrix(al,x1) * algtomatrix(al,x2) == algtomatrix(al, algmul(al,x1,x2)));
990print("changebasis bug 1: ", algalgtobasis(al,algbasistoalg(al,algmul(al,x1,x2)))==algmul(al,x1,x2));
991print("changebasis bug 2: ", algalgtobasis(al0,algmul(al0,ax1,ax2)) == algmul(al0,nx1,nx2));
992print("changebasis bug 3: ", invord*algmul(al0,nx1,nx2) == algmul(al,x1,x2));
993print("changebasis bug 4: ", algtomatrix(al,x1,1) == invord*algtomatrix(al0,nx1,1)*ord);
994
995print("algtableinit segfault bug: ");
996alt = algtableinit(al[9]);
997print(alt != 'alt);
998print("center of CSA: ", #algcenter(alt)==2);
999print("radical of CSA: ", algradical(alt)==0);
1000print("decomposition of CSA: ", #algsimpledec(alt)[2]==1);
1001dec = algsimpledec(alt);
1002{print("alg_decomposition of CSA: ", #dec==2 && dec[1]==0 && #dec[2]==1 &&
1003  #algcenter(dec[2][1])==2 && algdim(dec[2][1],1)==8);}
1004
1005print("alsimple bug");
1006mt = [matid(3), [0,0,0;1,1,0;0,0,0], [0,0,0;0,0,0;1,0,1]];
1007A = algtableinit(mt);
1008algissimple(A)
1009
1010print("tests for al_CSA: ");
1011T = y^3-y+1;
1012nf = nfinit(T);
1013   m_i = [0,-1,0, 0;\
1014          1, 0,0, 0;\
1015          0, 0,0,-1;\
1016          0, 0,1, 0];
1017   m_j = [0, 0,-1,0;\
1018          0, 0, 0,1;\
1019          1, 0, 0,0;\
1020          0,-1, 0,0];
1021   m_k = [0, 0, 0, -1;\
1022          0, 0,-1, 0;\
1023          0, 1, 0, 0;\
1024          1, 0, 0, 0];
1025mt = [matid(4), m_i, m_j, m_k];
1026print(algisassociative(mt));
1027al = alginit(nf, mt, 'x);
1028print(al != 0);
1029print("algebra:");
1030print("csa getcenter: ", algcenter(al) == nf);
1031print("csa getsplitting: ", algsplittingfield(al) != 0);
1032print("getrelmultable: ", algrelmultable(al) == mt);
1033print("getsplittingdata:");
1034print(#algsplittingdata(al) == 3);
1035print(#algsplittingdata(al)[1] == 12);
1036print(#algsplittingdata(al)[2] == 2);
1037print(#algsplittingdata(al)[3][1,] == 12);
1038print(#algsplittingdata(al)[3][,1] == 2);
1039print(al[3][3]*al[3][2][,1] == [1,0]~);
1040print(al[3][3]*al[3][2][,2] == [0,1]~);
1041polabs = al[1][12][1].pol;
1042for(i=1,10,\
1043print(al[3][3]*algmul(al, al[3][2][,1], algpow(al,al[3][1],i)) == [Mod(x^i,polabs),0]~);\
1044print(al[3][3]*algmul(al, al[3][2][,2], algpow(al,al[3][1],i)) == [0,Mod(x,polabs)^i]~)\
1045);
1046print("hasse invariants:");
1047do("hassei csa", () -> alghassei(al) == 0);
1048do("hassef cas", () -> alghassef(al) == 0);
1049do("hasse csa", () -> alghasse(al,1) == 0);
1050print("csa splitting pol: ", poldegree(al[1][12][1].pol) == 6);
1051print("csa basis: ", matsize(algbasis(al)) == [12,12]);
1052print("csa invbasis: ", matsize(alginvbasis(al)) == [12,12]);
1053print("csa absdim: ", #algmultable(al) == algdim(al,1));
1054print("csa char: ", algchar(al) == 0);
1055print("csa deg: ", algdegree(al) == 2);
1056print("csa dim: ", algdim(al) == 4);
1057print("csa absdim: ", algdim(al,1) == 12);
1058print("csa type: ", algtype(al) == 2); \\2==al_CSA
1059print("csa iscommutative: ", algiscommutative(al)==0);
1060print("csa issemisimple: ", algissemisimple(al)==1);
1061
1062print("elements:");
1063a = [0, Mod(y,T), 0, 0]~;
1064b = [0, -1, Mod(2*y^2,T), 0]~;
1065c = [Mod(1-y+2*y^2,T), 3, 0, Mod(-3*y,T)]~;
1066mynorm(aa) = sum(i=1,4,aa[i]^2);
1067algbasistoalg(al,a)
1068algalgtobasis(al,[1,2,3,4,5,6,7,8,9,10,11,12]~)
1069
1070print("csa add: ", algadd(al,a,b) == a+b);
1071print("csa neg: ", algneg(al,a) == -a);
1072print("csa neg 2: ", algneg(al,b) == -b);
1073print("csa sub: ", algsub(al,a,b) == a-b);
1074print("csa mul: ", algmul(al,a,b) == [Mod(y,T), 0, 0, Mod(2*y-2,T)]~);
1075print("csa mul 2: ", algmul(al,b,a) == [Mod(y,T), 0, 0, Mod(2-2*y,T)]~);
1076print("csa sqr: ", algsqr(al,a) == [Mod(-y^2,T),0,0,0]~);
1077print("csa sqr 2: ", algsqr(al,b) == [Mod(-1-4*y^4,T),0,0,0]~);
1078print("csa mt: ", algrelmultable(al)*b == -m_i + Mod(2*y^2,T)*m_j);
1079print("csa inv: ", alginv(al,a) == -1/y^2*a);
1080print("csa inv 2: ", alginv(al,b) == -1/Mod(1+4*y^4,T)*b);
1081print("csa divl: ", algdivl(al,1+a+b,b) == algmul(al, alginv(al, 1+a+b), b));
1082print("csa pow: ", algpow(al, a, 5) == Mod(y^4,T)*a);
1083aa = algalgtobasis(al, a);
1084bb = algalgtobasis(al, b);
1085cc = algalgtobasis(al, c);
1086print("csa mul 3: ", algmul(al,aa,b) == algalgtobasis(al,algmul(al,a,b)));
1087print("csa mul 4: ", algmul(al,a,bb) == algalgtobasis(al,algmul(al,a,b)));
1088print("csa pow 2: ", algpow(al,aa,13) == algalgtobasis(al,algpow(al,a,13)));
1089print("csa sub 2: ", algsub(al,aa,b) == algalgtobasis(al,algsub(al,a,b)));
1090print("csa sub 3: ", algsub(al,bb,a) == algalgtobasis(al,algsub(al,b,a)));
1091print("csa inv 3: ", alginv(al,aa) == algalgtobasis(al,alginv(al,a)));
1092print("csa inv 4: ", alginv(al,bb) == algalgtobasis(al,alginv(al,b)));
1093print("csa inv 5: ", alginv(al,algadd(al,a,bb)) == algalgtobasis(al,alginv(al,algadd(al,a,b))));
1094print("csa trace: ", algtrace(al,cc) == 2*c[1]);
1095print("csa trace 2: ", algtrace(al,c) == 2*c[1]);
1096
1097D = 12;
1098flag = 1;
1099for(i=1, D,\
1100  for(j=i, D,\
1101    ei = matid(D)[,i];\
1102    ej = matid(D)[,j];\
1103    flag = flag && (algtomatrix(al,ei)*algtomatrix(al,ej) == algtomatrix(al, algmul(al, ei, ej)));\
1104    flag = flag && (algtomatrix(al,ei)+algtomatrix(al,ej) == algtomatrix(al, algadd(al, ei, ej)))\
1105    ));
1106print(flag);
1107
1108print("testcharpol");
1109testcharpol(al,elt) = print(algcharpoly(al,elt)==x^2-algtrace(al,elt)*x+mynorm(elt));
1110testcharpol(al,a);
1111testcharpol(al,b);
1112testcharpol(al,c);
1113
1114print("testcharpol2");
1115testcharpol2(al,elt) = print(algcharpoly(al,elt)==x^2-algtrace(al,elt)*x+mynorm(algbasistoalg(al,elt)));
1116testcharpol2(al,aa);
1117testcharpol2(al,bb);
1118testcharpol2(al,cc);
1119
1120print("testnorm");
1121testnorm(al,elt) = print(algnorm(al,elt) == mynorm(elt));
1122testnorm(al,a);
1123testnorm(al,b);
1124testnorm(al,c);
1125
1126print("testnorm2");
1127testnorm2(al,elt) = print(algnorm(al,elt) == mynorm(algbasistoalg(al,elt)));
1128testnorm2(al,aa);
1129testnorm2(al,bb);
1130testnorm2(al,cc);
1131
1132doubleindex(N,i,j) = (i-1)*N+j;
1133matrixringmt(N) =
1134{
1135  my(mt = [Mat([Col(0,N^2) | i<-[1..N^2]]) | j<-[1..N^2]], B, Bi, mt2=Vec(0,N^2));
1136  for(i=1,N,
1137    for(j=1,N,
1138      for(k=1,N,
1139        mt[doubleindex(N,i,j)][doubleindex(N,i,k),doubleindex(N,j,k)] = 1
1140      )
1141    )
1142  );
1143  B = matid(N^2);
1144  for(i=1, N, B[doubleindex(N,i,i),1]=1);
1145  Bi = B^(-1);
1146  mt2[1] = matid(N^2);
1147  for(i=2,N^2, mt2[i] = Bi*mt[i]*B);
1148  mt2;
1149}
1150
1151print("examples from docu");
1152setrand(1);
1153algtype([])
1154A = alginit(nfinit(y),[-1,1]);
1155algadd(A,[1,0]~,[1,2]~)
1156algisinv(A,[-1,1]~)
1157algisinv(A,[1,2]~,&ix)
1158ix
1159algisdivl(A,[x+2,-x-2]~,[x,1]~)
1160algisdivl(A,[x+2,-x-2]~,[-x,x]~)
1161algisdivl(A,[x+2,-x-2]~,[-x,x]~,&z)
1162z
1163
1164A = alginit(nfinit(y^2-5),[2,y]);
1165algalgtobasis(A,[y,1]~)
1166algbasistoalg(A,algalgtobasis(A,[y,1]~))
1167
1168algbasistoalg(A,[0,1,0,0,2,-3,0,0]~)
1169algalgtobasis(A,algbasistoalg(A,[0,1,0,0,2,-3,0,0]~))
1170
1171mt = [matid(3), [0,0,0; 1,1,0; 0,0,0], [0,0,1; 0,0,0; 1,0,1]];
1172A = algtableinit(mt,2);
1173e = [0,1,0]~;
1174one = [1,0,0]~;
1175e2 = algsub(A,one,e);
1176algcentralproj(A,[e,e2])
1177algprimesubalg(A)
1178algquotient(A,[0;1;0])
1179algsubalg(A,[1,0; 0,0; 0,1])
1180algiscommutative(A)
1181algissimple(A)
1182
1183mt = [matid(3),[0,0,0;1,0,1;0,0,0],[0,0,0;0,0,0;1,0,1]];
1184A = algtableinit(mt);
1185algiscommutative(A)
1186algissemisimple(A)
1187algissimple(A)
1188algissimple(A,1)
1189
1190nf = nfinit(y^2-5);
1191A = alginit(nf, [-3,1-y]);
1192alghassef(A)
1193algdegree(A)^algdim(A,1)*nf.disc^algdim(A)*idealnorm(nf,alghassef(A)[1][2])^algdegree(A)
1194algdisc(A)
1195
1196nf = nfinit(y^3-y+1);
1197A = alginit(nf, [-1,-1]);
1198algdim(A,1)
1199algcenter(A).pol
1200algdegree(A)
1201algdim(A)
1202
1203nf = nfinit(y);
1204p = idealprimedec(nf,7)[1];
1205p2 = idealprimedec(nf,11)[1];
1206A = alginit(nf,[3,[[p,p2],[1/3,2/3]],[0]]);
1207algaut(A)
1208algb(A)
1209
1210mt = [matid(3), [0,0,0; 1,1,0; 0,0,0], [0,0,1; 0,0,0; 1,0,1]];
1211A = algtableinit(mt,13);
1212algchar(A)
1213algtype(A)
1214
1215nf = nfinit(y^2-5);
1216A = alginit(nf, [-1,2*y-1]);
1217alghassef(A)
1218A = alginit(nf, [-1,y]);
1219alghassei(A)
1220alghasse(A, 1)
1221alghasse(A, 2)
1222alghasse(A, idealprimedec(nf,2)[1])
1223alghasse(A, idealprimedec(nf,5)[1])
1224algindex(A, 1)
1225algindex(A, 2)
1226algindex(A, idealprimedec(nf,2)[1])
1227algindex(A, idealprimedec(nf,5)[1])
1228algindex(A)
1229algisdivision(A, 1)
1230algisdivision(A, 2)
1231algisdivision(A, idealprimedec(nf,2)[1])
1232algisdivision(A, idealprimedec(nf,5)[1])
1233algisdivision(A)
1234algissplit(A, 1)
1235algissplit(A, 2)
1236algissplit(A, idealprimedec(nf,2)[1])
1237algissplit(A, idealprimedec(nf,5)[1])
1238algissplit(A)
1239algisramified(A, 1)
1240algisramified(A, 2)
1241algisramified(A, idealprimedec(nf,2)[1])
1242algisramified(A, idealprimedec(nf,5)[1])
1243algisramified(A)
1244algramifiedplaces(A)
1245
1246nf = nfinit(y^2-5);
1247pr = idealprimedec(nf,13)[1];
1248pol = nfgrunwaldwang(nf, [pr], [2], [0,-1], 'x)
1249
1250
1251A = alginit(nfinit(y), [-1,-1]);
1252alginvbasis(A)
1253algbasis(A)
1254algmultable(A)
1255alginv(A,[1,1,0,0]~)
1256algmul(A,[1,1,0,0]~,[0,0,2,1]~)
1257algtomatrix(A,[0,1,0,0]~,1)
1258algtomatrix(A,[0,x]~,1)
1259algneg(A,[1,1,0,0]~)
1260algtomatrix(A,[0,0,0,2]~)
1261algpow(A,[1,1,0,0]~,7)
1262algsub(A,[1,1,0,0]~,[1,0,1,0]~)
1263algtrace(A,[5,0,0,1]~)
1264algtype(A)
1265
1266nf = nfinit(y^3-5);
1267a = y; b = y^2;
1268{m_i = [0,a,0,0;
1269       1,0,0,0;
1270       0,0,0,a;
1271       0,0,1,0];}
1272{m_j = [0, 0,b, 0;
1273       0, 0,0,-b;
1274       1, 0,0, 0;
1275       0,-1,0, 0];}
1276{m_k = [0, 0,0,-a*b;
1277       0, 0,b,   0;
1278       0,-a,0,   0;
1279       1, 0,0,   0];}
1280mt = [matid(4), m_i, m_j, m_k];
1281A = alginit(nf,mt,'x);
1282algrelmultable(A)
1283algsplittingfield(A).pol
1284algsplittingdata(A)
1285algtype(A)
1286
1287mt = [matid(3), [0,0,0; 1,1,0; 0,0,0], [0,0,1; 0,0,0; 1,0,1]];
1288A = algtableinit(mt,19);
1289algnorm(A,[0,-2,3]~)
1290A = algtableinit(mt);
1291algnorm(A,[0,-2,3]~)
1292
1293m_i=[0,-1,0,0;1,0,0,0;0,0,0,-1;0,0,1,0];
1294m_j=[0,0,-1,0;0,0,0,1;1,0,0,0;0,-1,0,0];
1295m_k=[0,0,0,-1;0,0,-1,0;0,1,0,0;1,0,0,0];
1296mt = [matid(4), m_i, m_j, m_k];
1297A = algtableinit(mt);
1298algissemisimple(A)
1299algissimple(A)
1300algissimple(A,1)
1301\\end examples
1302
1303
1304
1305
1306
1307
1308print("matrices over algebras");
1309scal8(a) = vector(8,i,if(i==1,a,0))~;
1310setrand(1);
1311nf = nfinit(y^2-5);
1312al = alginit(nf, [-1,-7]);
1313setrand(1);
1314M1 = [algrandom(al,2),algrandom(al,2);algrandom(al,2),algrandom(al,2)]
1315M2 = [algrandom(al,2),algrandom(al,2);algrandom(al,2),algrandom(al,2)]
1316print("mul alM: ", algmul(al,M1,M2));
1317a = [1,2,3,4,5,6,7,8]~; M = Mat([0,0]); M[1,1] = a; M[1,2] = a;
1318algmul(al, M, [a,a,a;a,a,a]);
1319print("sqr alM: ", algsqr(al,M1) == algmul(al,M1,M1));
1320print("divl alM: ", algmul(al,M1,algdivl(al,M1,M2)) == M2);
1321print("divr alM: ", algmul(al,algdivr(al,M1,M2),M2) == M1);
1322print("isinv alM: ", algisinv(al, M1));
1323print("isinv alM 2: ", algisinv(al, M2));
1324un = [1,0,0,0,0,0,0,0]~; zero = [0,0,0,0,0,0,0,0]~; id = [un,zero;zero,un];
1325{print("inv alM: ", algmul(al,M1,alginv(al,M1)) == id &&
1326 algmul(al,alginv(al,M1),M1) == id);}
1327{print("inv alM 2: ", algmul(al,M2,alginv(al,M2)) == id &&
1328 algmul(al,alginv(al,M2),M2) == id);}
1329print("neg alM: ", algneg(al,M1) == -M1);
1330print("sub alM: ", algsub(al,M1,M2) == M1-M2);
1331print("add alM: ", algadd(al,M1,M2) == M1+M2);
1332print("algtobasis basistoalg alM 1: ", algalgtobasis(al, algbasistoalg(al, M1)) ==
1333 M1);
1334print("algtobasis basistoalg alM 2: ", algalgtobasis(al, algbasistoalg(al, M2)) ==
1335 M2);
1336print("algleftmultable add alM: ", algtomatrix(al, M1,1)+algtomatrix(al,M2,1) == algtomatrix(al, algadd(al, M1, M2),1));
1337print("algleftmultable mul alM: ", algtomatrix(al, M1,1)*algtomatrix(al,M2,1) == algtomatrix(al, algmul(al, M1, M2),1));
1338{print("algleftmultable sqr alM: ", algtomatrix(al, M1,1)^2 == algtomatrix(al, algsqr(al, M1),1));}
1339print("algsplitm add alM: ", algtomatrix(al, M1)+algtomatrix(al,M2) ==
1340 algtomatrix(al, algadd(al, M1, M2)));
1341print("algsplitm mul alM: ", algtomatrix(al, M1)*algtomatrix(al,M2) ==
1342 algtomatrix(al, algmul(al, M1, M2)));
1343print("algsplitm sqr alM: ", algtomatrix(al, M1)^2 ==
1344 algtomatrix(al, algsqr(al, M1)));
1345print("algsplitm sqr alM 2: ", algtomatrix(al, M2)^2 ==
1346 algtomatrix(al, algsqr(al, M2)));
1347{print("algtrace alM: ", algtrace(al,M1) == algtrace(al,M1[1,1]) +
1348 algtrace(al,M1[2,2]));}
1349{print("algtrace alM 2: ", algtrace(al,M2) == algtrace(al,M2[1,1]) +
1350 algtrace(al,M2[2,2]));}
1351{print("algtrace prod alM: ", algtrace(al, algmul(al,M1,M2)) == algtrace(al,
1352algmul(al,M2,M1)));}
1353{print("algnorm alM: ", algnorm(al,algmul(al,M1,M2)) == algnorm(al,M1) *
1354 algnorm(al,M2));}
1355{print("algnorm alM 2: ", algnorm(al,algmul(al,M2,M1)) == algnorm(al,M1) *
1356 algnorm(al,M2));}
1357{print("algcharpoly alM: ", poldegree(algcharpoly(al,M1))==4 &&
1358 polcoeff(algcharpoly(al,M1),3) == -algtrace(al,M1) &&
1359 polcoeff(algcharpoly(al,M1),0) == algnorm(al,M1));}
1360{print("algcharpoly alM 2: ", poldegree(algcharpoly(al,M2))==4 &&
1361 polcoeff(algcharpoly(al,M2),3) == -algtrace(al,M2) &&
1362 polcoeff(algcharpoly(al,M2),0) == algnorm(al,M2));}
1363m = 15; n = -8;
1364{print("pow alM: ", algmul(al, algpow(al,M1,m), algpow(al,M1,n)) ==
1365 algpow(al,M1,m+n));}
1366{print("pow alM 2: ", algmul(al, algpow(al,M2,m), algpow(al,M2,n)) ==
1367 algpow(al,M2,m+n));}
1368print("pow 0 alM: ", algpow(al,M1,0) == id);
1369algbasistoalg(al,M1)
1370algbasistoalg(al,M2)
1371m1 = [1,2;3,4];
1372m2 = [5,6;7,8];
1373M1 = apply(scal8,m1);
1374M2 = apply(scal8,m2);
1375print("mul scalar alM: ", algmul(al,M1,M2) == apply(scal8,m1*m2));
1376
1377m_i=[0,-1,0,0;1,0,0,0;0,0,0,-1;0,0,1,0];
1378m_j=[0,0,-1,0;0,0,0,1;1,0,0,0;0,-1,0,0];
1379m_k=[0,0,0,-1;0,0,-1,0;0,1,0,0;1,0,0,0];
1380mt = [matid(4), m_i, m_j, m_k];
1381al = algtableinit(mt);
1382setrand(10);
1383M1 = [algrandom(al,2),algrandom(al,2);algrandom(al,2),algrandom(al,2)]
1384M2 = [algrandom(al,2),algrandom(al,2);algrandom(al,2),algrandom(al,2)]
1385print("mul alM t: ", algmul(al,M1,M2));
1386print("sqr alM t: ", algsqr(al,M1) == algmul(al,M1,M1));
1387print("divl alM t: ", algmul(al,M1,algdivl(al,M1,M2)) == M2);
1388print("divr alM t: ", algmul(al,algdivr(al,M1,M2),M2) == M1);
1389print("isinv alM t: ", algisinv(al, M1));
1390print("isinv alM t 2: ", algisinv(al, M2));
1391un = [1,0,0,0]~; zero = [0,0,0,0]~; id = [un,zero;zero,un];
1392{print("inv alM t: ", algmul(al,M1,alginv(al,M1)) == id &&
1393 algmul(al,alginv(al,M1),M1) == id);}
1394{print("inv alM t 2: ", algmul(al,M2,alginv(al,M2)) == id &&
1395 algmul(al,alginv(al,M2),M2) == id);}
1396print("neg alM t: ", algneg(al,M1) == -M1);
1397print("sub alM t: ", algsub(al,M1,M2) == M1-M2);
1398print("add alM t: ", algadd(al,M1,M2) == M1+M2);
1399print("algleftmultable add alM t: ", algtomatrix(al,M1,1) + algtomatrix(al,M2,1) == algtomatrix(al, algadd(al,M1,M2),1));
1400print("algleftmultable mul alM t: ", algtomatrix(al,M1,1) * algtomatrix(al,M2,1) == algtomatrix(al, algmul(al,M1,M2),1));
1401print("algleftmultable sqr alM t: ", algtomatrix(al,M1,1)^2 == algtomatrix(al,algsqr(al,M1),1));
1402{print("algtrace alM t: ", algtrace(al,M1) == 2*(algtrace(al,M1[1,1]) +
1403 algtrace(al,M1[2,2])));}
1404{print("algtrace alM t 2: ", algtrace(al,M2) == 2*(algtrace(al,M2[1,1]) +
1405 algtrace(al,M2[2,2])));}
1406{print("algtrace prod alM t: ", algtrace(al, algmul(al,M1,M2)) == algtrace(al,
1407algmul(al,M2,M1)));}
1408{print("algnorm alM t: ", algnorm(al,algmul(al,M1,M2)) == algnorm(al,M1) *
1409 algnorm(al,M2));}
1410{print("algnorm alM t 2: ", algnorm(al,algmul(al,M2,M1)) == algnorm(al,M1) *
1411 algnorm(al,M2));}
1412{print("algcharpoly alM t: ", poldegree(algcharpoly(al,M1))==16 &&
1413 polcoeff(algcharpoly(al,M1),15) == -algtrace(al,M1) &&
1414 polcoeff(algcharpoly(al,M1),0) == algnorm(al,M1));}
1415{print("algcharpoly alM t 2: ", poldegree(algcharpoly(al,M2))==16 &&
1416 polcoeff(algcharpoly(al,M2),15) == -algtrace(al,M2) &&
1417 polcoeff(algcharpoly(al,M2),0) == algnorm(al,M2));}
1418m = 32; n = -63;
1419{print("pow alM t: ", algmul(al, algpow(al,M1,m), algpow(al,M1,n)) ==
1420 algpow(al,M1,m+n));}
1421{print("pow alM 2 t: ", algmul(al, algpow(al,M2,m), algpow(al,M2,n)) ==
1422 algpow(al,M2,m+n));}
1423print("pow 0 alM t: ", algpow(al,M2,0) == id);
1424
1425
1426
1427T = y^3-y+1;
1428nf = nfinit(T);
1429print("csa al2");
1430setrand(1);
1431al2 = alginit(nf, matrixringmt(2), 'x);
1432print("al2 contains nfabs: ", algsplittingfield(al2)[12][1] != 0);
1433al2b = al2; al2b[1][12][1] = 0; \\ depends on 32/64bit
1434print(al2b);
1435print("csa al3");
1436al3 = alginit(nf, matrixringmt(3), 'x);
1437print("al3 contains nfabs: ", algsplittingfield(al3)[12][1] != 0);
1438
1439\\limit cases
1440print("trivial algebra over a quadratic field");
1441al = alginit(rnfinit(nfinit(y^2+1),x),[y,1])
1442a = [y]~
1443b = [1-2*y]~
1444c = [-3,1]~
1445algadd(al,a,b)
1446algsub(al,c,a)
1447algmul(al,a,b)
1448algdivl(al,b,c)
1449algdivr(al,c,b)
1450alginv(al,a)
1451algalgtobasis(al,b)
1452algtomatrix(al,a)
1453algtomatrix(al,a,1)
1454algcharpoly(al,b)
1455algtrace(al,c)
1456algnorm(al,c)
1457algiscommutative(al)
1458algissemisimple(al)
1459algissimple(al)
1460alghasse(al,1)
1461alghasse(al,idealprimedec(nfinit(y^2+1),2)[1])
1462algindex(al)
1463algisdivision(al)
1464algissplit(al)
1465algisramified(al)
1466algramifiedplaces(al)
1467
1468
1469print("trivial algebra over Q");
1470al = alginit(rnfinit(nfinit(y),x),[y,1])
1471a = [-2]~
1472b = [1/3]~
1473c = [4/5]~
1474algadd(al,a,b)
1475algsub(al,c,a)
1476algmul(al,a,b)
1477algdivl(al,b,c)
1478algdivr(al,c,b)
1479alginv(al,a)
1480algalgtobasis(al,b)
1481algtomatrix(al,a,1)
1482algtomatrix(al,b)
1483algcharpoly(al,b)
1484algtrace(al,c)
1485algnorm(al,c)
1486algiscommutative(al)
1487algissemisimple(al)
1488algissimple(al)
1489alghasse(al,1)
1490alghasse(al,idealprimedec(nfinit(y),2)[1])
1491algindex(al)
1492algisdivision(al)
1493algissplit(al)
1494algisramified(al)
1495algramifiedplaces(al)
1496
1497print("trivial CSA over Q");
1498al = alginit(nfinit(y), [Mat(1)]);
1499algsqr(al,[Mod(3,y)]~)
1500algsqr(al,[2]~)
1501print("nontrivial CSA over Q");
1502{m_i = [0,-1,0, 0;
1503       1, 0,0, 0;
1504       0, 0,0,-1;
1505       0, 0,1, 0];
1506m_j = [0, 0,-1,0;
1507       0, 0, 0,1;
1508       1, 0, 0,0;
1509       0,-1, 0,0];
1510m_k = [0, 0, 0, -1;
1511       0, 0,-1, 0;
1512       0, 1, 0, 0;
1513       1, 0, 0, 0];
1514mt = [matid(4), m_i, m_j, m_k];}
1515al = alginit(nfinit(y), mt);
1516algsqr(al,[Mod(3,y),Mod(2,y),Mod(1,y),Mod(2,y)]~)
1517algsqr(al,[2,3,4,5]~)
1518
1519print("empty matrices");
1520al = alginit(nfinit(y), [-1,-1]);
1521print("-v: ", algneg(al,[;]) == [;]);
1522print("v^(-1): ", alginv(al,[;]) == [;]);
1523print("v^n: ", algpow(al, [;], 13) == [;]);
1524print("v^0: ", algpow(al, [;], 0) == [;]);
1525print("mt(v)", algtomatrix(al, [;], 1) == [;]);
1526print("spl(v)", algtomatrix(al, [;]) == [;]);
1527print("trace(v): ", algtrace(al, [;]) == 0);
1528print("norm(v): ", algnorm(al, [;]) == 1);
1529print("charpoly(v): ", algcharpoly(al, [;]) == 1 && type(algcharpoly(al,[;])) == "t_POL");
1530print("v+v: ", algadd(al,[;],[;]) == [;]);
1531print("v-v: ", algsub(al,[;],[;]) == [;]);
1532print("v*v: ", algmul(al,[;],[;]) == [;]);
1533print("v/v: ", algdivr(al,[;],[;]) == [;]);
1534print("v\\v: ", algdivl(al,[;],[;]) == [;]);
1535v1 = matrix(0,1);
1536print("v*nv: ", algmul(al,v1,matid(1))==v1);
1537print("v*v 2: ", algmul(al,[;],matrix(0,1))==matrix(0,1));
1538print("trace(v) char 2: ", algtrace(algtableinit([matid(1)],2), [;]) == 0);
1539
1540mt0 = [Mat([1])];
1541almt0 = algtableinit(mt0,0)
1542a = [12]~
1543b = [-1/7]~
1544algadd(almt0,a,b)
1545algsub(almt0,a,b)
1546algmul(almt0,a,b)
1547algneg(almt0,a)
1548alginv(almt0,a)
1549algsqr(almt0,b)
1550algdivl(almt0,a,b)
1551algtrace(almt0,a)
1552algnorm(almt0,b)
1553algcharpoly(almt0,a)
1554algtomatrix(almt0,b,1)
1555algpow(almt0,a,0)
1556algiscommutative(almt0)
1557algissemisimple(almt0)
1558algissimple(almt0)
1559algisdivision(almt0)
1560
1561
1562print("trivial tensor product");
1563al1 = alginit(nfinit(y),1);
1564al2 = alginit(nfinit(y),2);
1565print(algtensor(al1,al2)==al2);
1566print(algtensor(al2,al1)==al2);
1567
1568print("splitting a nasty commutative algebra");
1569{mt = [matid(8),
1570[0,2,0,0,0,0,0,0;
1571 1,0,0,0,0,0,0,0;
1572 0,0,0,2,0,0,0,0;
1573 0,0,1,0,0,0,0,0;
1574 0,0,0,0,0,2,0,0;
1575 0,0,0,0,1,0,0,0;
1576 0,0,0,0,0,0,0,2;
1577 0,0,0,0,0,0,1,0],
1578
1579[0,0,3,0,0,0,0,0;
1580 0,0,0,3,0,0,0,0;
1581 1,0,0,0,0,0,0,0;
1582 0,1,0,0,0,0,0,0;
1583 0,0,0,0,0,0,3,0;
1584 0,0,0,0,0,0,0,3;
1585 0,0,0,0,1,0,0,0;
1586 0,0,0,0,0,1,0,0],
1587
1588[0,0,0,6,0,0,0,0;
1589 0,0,3,0,0,0,0,0;
1590 0,2,0,0,0,0,0,0;
1591 1,0,0,0,0,0,0,0;
1592 0,0,0,0,0,0,0,6;
1593 0,0,0,0,0,0,3,0;
1594 0,0,0,0,0,2,0,0;
1595 0,0,0,0,1,0,0,0],
1596
1597[0,0,0,0,5,0,0,0;
1598 0,0,0,0,0,5,0,0;
1599 0,0,0,0,0,0,5,0;
1600 0,0,0,0,0,0,0,5;
1601 1,0,0,0,0,0,0,0;
1602 0,1,0,0,0,0,0,0;
1603 0,0,1,0,0,0,0,0;
1604 0,0,0,1,0,0,0,0],
1605
1606[0,0,0,0,0,10,0, 0;
1607 0,0,0,0,5, 0,0, 0;
1608 0,0,0,0,0, 0,0,10;
1609 0,0,0,0,0, 0,5, 0;
1610 0,2,0,0,0, 0,0, 0;
1611 1,0,0,0,0, 0,0, 0;
1612 0,0,0,2,0, 0,0, 0;
1613 0,0,1,0,0, 0,0, 0],
1614
1615[0,0,0,0,0,0,15, 0;
1616 0,0,0,0,0,0, 0,15;
1617 0,0,0,0,5,0, 0, 0;
1618 0,0,0,0,0,5, 0, 0;
1619 0,0,3,0,0,0, 0, 0;
1620 0,0,0,3,0,0, 0, 0;
1621 1,0,0,0,0,0, 0, 0;
1622 0,1,0,0,0,0, 0, 0],
1623
1624[0,0,0,0,0, 0, 0,30;
1625 0,0,0,0,0, 0,15, 0;
1626 0,0,0,0,0,10, 0, 0;
1627 0,0,0,0,5, 0, 0, 0;
1628 0,0,0,6,0, 0, 0, 0;
1629 0,0,3,0,0, 0, 0, 0;
1630 0,2,0,0,0, 0, 0, 0;
1631 1,0,0,0,0, 0, 0, 0]
1632];}
1633{chg =
1634[1, 0,0,0,0,0, 0, 0;
1635 0, 1,0,0,0,0, 0, 0;
1636 0,-2,1,1,0,0, 0, 0;
1637 0,-1,0,1,0,0, 0, 0;
1638 0, 0,0,0,1,1,-2, 0;
1639 0, 0,0,0,0,1, 0,-1;
1640 0, 0,0,0,0,0, 1, 0;
1641 0, 0,0,0,0,0, 0, 1];}
1642chgi = chg^(-1);
1643mt2 = vector(8,j,chgi*sum(i=1,8,chg[i,j]*mt[i])*chg);
1644algisassociative(mt2)
1645al = algtableinit(mt2);
1646algiscommutative(al)
1647algissemisimple(al)
1648setrand(9991);
1649algissimple(al,1)
1650
1651print("non associative algebra");
1652mt = [matid(3), [0,2,3;1,4,5;0,6,7], [0,8,9;0,10,11;1,12,13]];
1653algisassociative(mt)
1654
1655print("csa without maximal order");
1656alginit(nfinit(y), [matid(1)], 'x, 0);
1657
1658print("simplify bug #1671");
1659test(str,al)={
1660 my(sal = simplify(al), x, y);
1661 setrand(1);
1662 print("testing simplify: ", str);
1663 print(algtype(al) == algtype(sal));
1664 print(algdim(al) == algdim(sal));
1665 setrand(11);
1666 x = algrandom(al,3);
1667 y = algrandom(al,10);
1668 print(algsqr(al,x) == algsqr(sal,x));
1669 print(algtomatrix(al,x) == algtomatrix(sal,x));
1670 print(algcharpoly(al,x) == algcharpoly(sal,x));
1671 print(algnorm(al,x) == algnorm(sal,x));
1672 print(algtomatrix(al,x,1) == algtomatrix(sal,x,1));
1673 print(algtrace(al,x) == algtrace(sal,x));
1674 print(alginv(al,x) == alginv(sal,x));
1675 print(algpow(al,x,42) == algpow(sal,x,42));
1676 print(algmul(al,x,y) == algmul(sal,x,y));
1677 print(algdivl(al,x,y) == algdivl(sal,x,y));
1678
1679 print(algbasistoalg(al,x) == algbasistoalg(sal,x));
1680 x = algbasistoalg(al,x);
1681 print(algbasistoalg(al,y) == algbasistoalg(sal,y));
1682 y = algbasistoalg(al,y);
1683 print(algsqr(al,x) == algsqr(sal,x));
1684 print(algtomatrix(al,x) == algtomatrix(sal,x));
1685 print(algcharpoly(al,x) == algcharpoly(sal,x));
1686 print(algnorm(al,x) == algnorm(sal,x));
1687 print(algtomatrix(al,x,1) == algtomatrix(sal,x,1));
1688 print(algtrace(al,x) == algtrace(sal,x));
1689 print(alginv(al,x) == alginv(sal,x));
1690 print(algpow(al,x,42) == algpow(sal,x,42));
1691 print(algmul(al,x,y) == algmul(sal,x,y));
1692 print(algdivl(al,x,y) == algdivl(sal,x,y));
1693
1694 print(algadd(al,x,y) == algadd(al,x,simplify(y)));
1695 print(algmul(al,x,y) == algmul(al,x,simplify(y)));
1696 print(algmul(al,x,y) == algmul(al,simplify(x),simplify(y)));
1697};
1698test("degree 1 cyclic over Q", alginit(nfinit(y),1));
1699test("degree 1 cyclic over Q(i)", alginit(nfinit(y^2+1),1));
1700test("degree 1 csa over Q", alginit(nfinit(y), [matid(1)]));
1701test("degree 1 csa over Q(i)", alginit(nfinit(y^2+1), [matid(1)]));
1702test("quatalg over Q(s5)", alginit(nfinit(y^2-5), [-2,-3]));
1703{m_i = [0,-1,0, 0;
1704       1, 0,0, 0;
1705       0, 0,0,-1;
1706       0, 0,1, 0];
1707m_j = [0, 0,-1,0;
1708       0, 0, 0,1;
1709       1, 0, 0,0;
1710       0,-1, 0,0];
1711m_k = [0, 0, 0, -1;
1712       0, 0,-1, 0;
1713       0, 1, 0, 0;
1714       1, 0, 0, 0];
1715mt = [matid(4), m_i, m_j, m_k];}
1716test("quatalg csa over Q", alginit(nfinit(y), mt));
1717
1718m = matcompanion(x^4+1);
1719mt = [m^i | i <- [0..3]];
1720al = algtableinit(mt);
1721B = [1,0;0,0;0,1/2;0,0]
1722al2 = algsubalg(al,B);
1723
1724
1725\\ bad inputs
1726m_i = [0,-1,0, 0;\
1727       1, 0,0, 0;\
1728       0, 0,0,-1;\
1729       0, 0,1, 0];
1730m_j = [0, 0,-1,0;\
1731       0, 0, 0,1;\
1732       1, 0, 0,0;\
1733       0,-1, 0,0];
1734m_k = [0, 0, 0, -1;\
1735       0, 0,-1, 0;\
1736       0, 1, 0, 0;\
1737       1, 0, 0, 0];
1738mt = [matid(4), m_i, m_j, m_k];
1739almt = algtableinit(mt,0);
1740algsplittingfield(almt);
1741algdegree(almt);
1742alghassei(almt);
1743alghassef(almt);
1744algrandom(1,1)
1745algrandom(1,I)
1746algtype(1)
1747algdim([1,[1],0,0,0,0,0,0,0,0])
1748algdim([1,[1],0,0,0,0,0,0,0,0],1)
1749algtensor(al,al2)
1750algtensor(al2,al)
1751algtensor(1,z,1)
1752algisassociative([1],0)
1753algisassociative([[1,0;0,2],[0,0;0,0]]) \\valid input
1754algmul(almt,a,b)
1755algtomatrix(almt,a,1)
1756alginv(almt,a)
1757algalgtobasis(almt,a)
1758algbasistoalg(almt,[0,0,0,0]~)
1759algpoleval(almt,1,a)
1760zero = [0,0,0,0]~; m = Mat([1,1]); m[1,1]=zero; m[1,2]=zero;
1761algadd(almt, [zero;zero], m)
1762algadd(almt, [zero;zero;zero], [zero;zero]);
1763algsub(almt, [zero;zero], m)
1764algsub(almt, [zero;zero;zero], [zero;zero]);
1765algmul(almt, m, [zero;zero;zero]);
1766algsqr(almt, [zero;zero]);
1767algdivl(almt, m, zero);
1768algdivl(almt, m, [zero,zero;zero,zero]);
1769algdivl(almt, m, m);
1770alginv(almt, m);
1771algtomatrix(almt,m,1);
1772algpow(almt, m, 3);
1773algtrace(almt, m);
1774algcharpoly(almt, m);
1775algcharpoly(alginit(nfinit(y),[-1,-1]), m);
1776algnorm(almt, m);
1777algnorm(alginit(nfinit(y),[-1,-1]), m);
1778alginit(nfinit(y),[2,[[],[]],[x]])
1779alginit(nfinit(y),[2,[],[1,1]])
1780alginit(nfinit(y),[2,[[],[]],Vecsmall([1])])
1781alginit(y,[2,[[],[]],[1]])
1782alginit(nfinit(y), y)
1783alginit(nfinit(y), [1,2,3,4])
1784algtableinit(mt,y);
1785alginit(nfinit(y^2+1),-3);
1786alginit(nfinit(x^2+1),3);
1787highvar = varhigher("highvar");
1788alginit(nfinit(highvar^2+1),3);
1789al = alginit(nfinit(y^2-2),[-1,-1]); algrandom(al,-10)
1790al = alginit(nfinit(y^2-2),[-1,-1]);
1791algrelmultable(al);
1792algsplittingdata(al);
1793alghasse(almt, 1);
1794algindex(almt, 1);
1795algisdivision(almt);
1796algissplit(almt);
1797algisramified(almt);
1798algramifiedplaces(almt);
1799alghasse(al, -1);
1800alghasse(al, 3);
1801alghasse(al, 2^100);
1802alghasse(al, []);
1803alghasse(al, 1/3);
1804algtableinit([matid(2), [0,1/2;1,0]]);
1805Q = nfinit(y);
1806alginit(Q, [matid(2), [0,1/2;1,0]]);
1807alginit(Q, [-1/2, -1]);
1808alginit(Q, [-1, -1/2]);
1809alginit(rnfinit(Q, x^2+1), [-x,-1/2]);
1810algsqr([0,0,0,0,0,0,0,0,0,0,0],[]~);
1811algsqr([0,0,0,0,0,0,0,0,[],0,0],[]~);
1812algsqr([0,0,0,0,0,0,0,0,[0],0,0],[]~);
1813algsqr([0,0,0,0,0,0,0,0,[[;]],0,0],[]~);
1814algsqr([[],0,0,0,0,0,0,0,[[;]],0,0],[]~);
1815algsqr([[],[0],0,0,0,0,0,0,[[;]],0,0],[]~);
1816algdim([[],[0],0,0,0,0,0,0,[[;]],0,0]);
1817algdegree([[],[0],0,0,0,0,0,0,[[;]],0,0]);
1818algdegree([rnfinit(nfinit(y),x),[[]],0,0,0,0,0,0,[[;]],0,0]);
1819algcenter([rnfinit(nfinit(y),x),[[]],0,0,0,0,0,0,[[;]],0,0]);
1820algcentralproj(almt,0);
1821algcentralproj(almt,[zero,zero]);
1822algsubalg(almt,0);
1823algisassociative([]);
1824algisassociative([matid(2),Mat([1,1])]);
1825algisassociative([[1,2;3,4],matid(2)]) \\valid input
1826algisassociative([matid(1)],[]);
1827algsqr(algtableinit([matid(1)]),[1,2]~);
1828algsqr(al,vector(691)~);
1829algsqr(al,[1,2,3,4,5,6,7,f^2]~);
1830algsqr(al,[f^3,[]]~);
1831algmul(al,[;],[1,2]~);
1832algdivl(al,[;],matid(1));
1833algdivl(al,matid(1),matrix(1,2));
1834alginv(al,[0,0]~);
1835al0mt = algtableinit([matid(1)]);
1836algalgtobasis(al0mt,[1]~);
1837algbasistoalg(al0mt,[1]~);
1838nfgrunwaldwang(nfinit(y),0,[],[],'x);
1839nfgrunwaldwang(nfinit(y),[2],'x-'x,[1]);
1840alginit(rnfinit(nfinit(y),x),0);
1841alginit(rnfinit(nfinit(y),x),[1,2,3,4]);
1842alginit(nfinit(y), [matid(2),matid(2)]);
1843alginit(nfinit(y), [matid(2),[0,1;1,0]]);
1844
1845nfgrunwaldwang(nfinit(y), 0, [], [0]);
1846nfgrunwaldwang(nfinit(y), [2], [], [0]);
1847nfgrunwaldwang(nfinit(y), [2], [2], []);
1848nfgrunwaldwang(nfinit(y), [2], [6], [0]);
1849nfgrunwaldwang(nfinit(y), [2,3], [2,3], [0]);
1850nfgrunwaldwang(nfinit(y), [2], [3], [-1]);
1851nfgrunwaldwang(nfinit(y), [[]~], [3], [-1]);
1852nfgrunwaldwang(nfinit(y), [2], [9], [0]);
1853
1854mt=[matid(3), [0,0,0; 1,1,0; 0,0,0], [0,0,1; 0,0,0; 1,0,1]];
1855A=algtableinit(mt,2);
1856algdegree(A)
1857algsub(A,1,1)
1858algadd(A,1,1)
1859algneg(A,1)
1860algmul(A,1,1)
1861algsqr(A,1)
1862algdivl(A,1,1)
1863algdivr(A,1,1)
1864alginv(A,1)
1865a='a;
1866K=nfinit(a);PR=idealprimedec(K,2);A=alginit(K,[3,[PR,[1]],[-1]]);
1867K=nfinit(a);P2=idealprimedec(K,2);P3=idealprimedec(K,3);A=alginit(K,[3,[concat(P2,P3),[1/3,-2/3]],[1/3]]);
1868algtensor(alginit(nfinit(y),2),alginit(nfinit(y^2+1),3));
1869algtensor(alginit(nfinit(y),2),alginit(nfinit(y),2));
1870nf = nfinit(y); p2 = idealprimedec(nf,2)[1]; p3 = idealprimedec(nf,3)[1];
1871\\alginit(nf, [2, [[p2,p3],[1/2,1/2]], [0]]);
1872alginit(nf, [2, [[p2,p2],[1/2,1/2]], [0]]);
1873alginit(nf, [2, [[p2,p3],[1/2,1/2]], [0,0]]);
1874alginit(nf, [2, [[p2,p3],[1/2,1/2],0], [0]]);
1875alginit(nf, [2, [0,[1/2,1/2]], [0]]);
1876alginit(nf, [2, [[p2,p3],0], [0]]);
1877alginit(nf, [2, [[p2,p3],[1/2,1/2,0]], [0]]);
1878alginit(nf, [2, [[p2,p3],[1/2,1/2]], [1/3]]);
1879
1880al = alginit(nfinit(y),[-1,-1]);
1881setrand(23);
1882a = algrandom(al,2);
1883algcharpoly(al,a,'z);
1884
1885al = alginit(nfinit(y^2+7), [-1,-1]);
1886algcharpoly(al,[1,2,3]~);
1887
1888algindex(1, 1)
1889
1890al = alginit(nfinit(y), [Mat(1)]);
1891algsqr(al,[Mod(1,y),Mod(2,y)]~)
1892{m_i = [0,-1,0, 0;
1893       1, 0,0, 0;
1894       0, 0,0,-1;
1895       0, 0,1, 0];
1896m_j = [0, 0,-1,0;
1897       0, 0, 0,1;
1898       1, 0, 0,0;
1899       0,-1, 0,0];
1900m_k = [0, 0, 0, -1;
1901       0, 0,-1, 0;
1902       0, 1, 0, 0;
1903       1, 0, 0, 0];
1904mt = [matid(4), m_i, m_j, m_k];}
1905al = alginit(nfinit(y), mt);
1906algsqr(al,[Mod(1,y),Mod(2,y)]~)
1907
1908alfail = alginit(nf, [0,0], 'x);
1909algb(al);
1910algaut(al);
1911algtableinit([Mat(1)],1)
1912algtableinit([Mat(1)],4)
1913al = alginit(nfinit(y),[-1,-1]);
1914algpoleval(al,x+1,"toto")
1915algpoleval(al,x+1,[1,2,3])
1916algpoleval(al,x+1,[1,2])
1917a = [1..4]~;
1918b = [5..8]~;
1919mb = algtomatrix(al,b,1);
1920algpoleval(al,x+1,[a,mb]);
1921algpoleval(al,x+1,[1,mb]);
1922alginit(nfinit(y),["a",[[],[]],[]]);
1923alginit(nfinit(y),[1,[[],[]],[]]);
1924alginit(nfinit(y),[0,[[],[]],[]]);
1925\\end bad inputs
1926
1927
1928print("new algsimpledec");
1929\\ K^3
1930mt = [matid(3),[0,0,0;1,1,0;0,0,0],[0,0,0;0,0,0;1,0,1]];
1931al = algtableinit(mt);
1932algissimple(al)
1933algsimpledec(al,1)
1934al2 = algtableinit(mt,5);
1935algissimple(al2)
1936algsimpledec(al2,1)
1937\\ upper-tri in M2(K)
1938mt = [matid(3),[0,0,0;1,1,0;0,0,1],[0,0,0;0,0,0;1,0,0]];
1939al = algtableinit(mt);
1940algsimpledec(al,1)
1941al2 = algtableinit(mt,5);
1942algsimpledec(al2,1)
1943
1944print("norm(,1)");
1945nf = nfinit(y^2-5);
1946B = alginit(nf,[-1,y]);
1947b = [x,1]~;
1948algnorm(B,b,1)
1949n = algnorm(B,b)
1950nfeltnorm(nf,n)^2
1951algnorm(B,b/3,1)
1952algnorm(B,[3,5,1,6,7,-2,1/7,0]~,1)
1953m = [[1/3,y]~,[0,0]~;[y+1,y-1]~,[1..8]~];
1954algnorm(B,m,1)==matdet(algtomatrix(B,m,1))
1955a = algrandom(B,2);
1956algnorm(B,a,1)*algnorm(B,b,1)==algnorm(B,algmul(B,a,b),1)
1957
1958print("trace(,1)");
1959nf = nfinit(y^2-5);
1960A = alginit(nf,[-1,y]);
1961a = [1+x+y,2*y]~*Mod(1,y^2-5)*Mod(1,x^2+1);
1962t = algtrace(A,a)
1963algtrace(A,a,1)
1964algdegree(A)*nfelttrace(nf,t)
1965b = [1+x/3+y,y-x/7]~*Mod(1,y^2-5)*Mod(1,x^2+1);
1966algdegree(A)*nfelttrace(nf,algtrace(A,b))==algtrace(A,b,1)
1967c = [4/3,2,1,6,8,9,-1,3/2]~;
1968algdegree(A)*nfelttrace(nf,algtrace(A,c))==algtrace(A,c,1)
1969m = [[y,x]~*Mod(1,y^2-5)*Mod(1,x^2+1),[1,0]~;[1..8]~,[1/7,5,6,7,8,1/4,1/3,1/2]~];
19702*algdegree(A)*nfelttrace(nf,algtrace(A,m))==algtrace(A,m,1)
1971algtrace(A,[[1,0]~,[0,0]~;[0,0]~,[1,0]~],1)==4*algdim(A,1)
1972
1973print("charpoly(,1)");
1974nf = nfinit(y^2-5);
1975al = alginit(nf,[-3-y,y]);
1976pol = nf.pol;
1977polrel = algsplittingfield(al).pol;
1978a = [y,1+x]~*Mod(1,pol)*Mod(1,polrel);
1979P = lift(algcharpoly(al,a))
1980algcharpoly(al,a,,1)
1981lift(P*subst(P,y,-y)*Mod(1,pol))^2
1982b = [1+x/3+y/2,3*y-x/7]~*Mod(1,pol)*Mod(1,polrel);
1983P = lift(algcharpoly(al,b));
1984Q = algcharpoly(al,b,,1);
1985Q == lift(P*subst(P,y,-y)*Mod(1,pol))^2
1986c = [4/3,2,1/7,6,4,9,-1,3/2]~;
1987P = lift(algcharpoly(al,c));
1988Q = algcharpoly(al,c,,1);
1989Q == lift(P*subst(P,y,-y)*Mod(1,pol))^2
1990{m = [[y/3,x+1]~*Mod(1,pol)*Mod(1,polrel),[1,7/11]~;
1991    [1..8]~,[1/7,5,6,7/2,8,1/4,1,-1/2]~];}
1992P = lift(algcharpoly(al,m));
1993Q = algcharpoly(al,m,,1);
1994Q==lift(P*subst(P,y,-y)*Mod(1,pol))^4
1995
1996print("more al_MAT tests");
1997setrand(1);
1998nf=nfinit(y^2-5);
1999pol = nf.pol;
2000al=alginit(nf,[y,-1]);
2001polrel = algsplittingfield(al).pol;
2002m1 = matrix(2,2,i,j,algrandom(al,1)/(i+2*j));
2003{m2 = matrix(2,2,i,j,vector(2,k,random(2)+random(2)*x+random(2)*y+
2004  random(2)*x*y)~*Mod(1,polrel)*Mod(1,pol)/prime(i+2*j));}
2005{m3 = [[1,0]~,[1..8]~;[y+x,y-x]~*Mod(1,polrel)*Mod(1,pol),
2006  [1,3,5,0,-1,-2,0,-2]~/3];}
2007print("add");
2008algadd(al,m1,m2)==algadd(al,m2,m1)
2009algadd(al,m1,m3)==algadd(al,m3,m1)
2010algadd(al,m1,m2)==algadd(al,m2,m1)
2011algadd(al,algadd(al,m1,m2),m3) == algadd(al,m1,algadd(al,m2,m3))
2012print("alg/basis");
2013algalgtobasis(al,m1) == m1
2014algbasistoalg(al,m2) == m2
2015algalgtobasis(al,algbasistoalg(al,algalgtobasis(al,m1))) == algalgtobasis(al,m1)
2016algalgtobasis(al,algbasistoalg(al,algalgtobasis(al,m2))) == algalgtobasis(al,m2)
2017algalgtobasis(al,algbasistoalg(al,algalgtobasis(al,m3))) == algalgtobasis(al,m3)
2018algbasistoalg(al,algalgtobasis(al,algbasistoalg(al,m1))) == algbasistoalg(al,m1)
2019algbasistoalg(al,algalgtobasis(al,algbasistoalg(al,m2))) == algbasistoalg(al,m2)
2020algbasistoalg(al,algalgtobasis(al,algbasistoalg(al,m3))) == algbasistoalg(al,m3)
2021print("charpoly");
2022mid = [[1,0]~,[0,0]~;[0,0]~,[1,0]~];
2023t = 123/7;
2024testcp(m)=
2025{
2026  my(P,Q);
2027  P = algcharpoly(al,m);
2028  print(algnorm(al,m) == polcoeff(P,0));
2029  print(algtrace(al,m) == -polcoeff(P,3));
2030  print(algnorm(al,algsub(al,m,t*mid))==subst(P,x,t));
2031  Q = algcharpoly(al,m,,1);
2032  print(algnorm(al,m,1) == polcoeff(Q,0));
2033  print(algtrace(al,m,1) == -polcoeff(Q,31));
2034  print(algnorm(al,algsub(al,m,t*mid),1)==subst(Q,x,t));
2035}
2036testcp(m1);
2037testcp(m2);
2038testcp(m3);
2039print("inv/div");
2040algisinv(al,m1,&m1i)
2041m1i == alginv(al,m1)
2042algdivl(al,m1,m2) == algmul(al,m1i,m2)
2043algdivr(al,m2,m1) == algmul(al,m2,m1i)
2044algisdivl(al,m1,m3,&d13)
2045algdivl(al,m1,m3) == d13
2046algisinv(al,m2,&m2i)
2047m2i == alginv(al,m2)
2048algdivl(al,m2,m1) == algmul(al,m2i,m1)
2049algdivr(al,m1,m2) == algmul(al,m1,m2i)
2050algisdivl(al,m2,m3,&d23)
2051algdivl(al,m2,m3) == d23
2052algisinv(al,m3,&m3i)
2053m3i == alginv(al,m3)
2054algdivl(al,m3,m1) == algmul(al,m3i,m1)
2055algdivr(al,m1,m3) == algmul(al,m1,m3i)
2056algisdivl(al,m3,m2,&d32)
2057algdivl(al,m3,m2) == d32
2058midbasis = [[1,0,0,0,0,0,0,0]~,vector(8)~;vector(8)~,[1,0,0,0,0,0,0,0]~];
2059algmul(al,m1i,m1) == midbasis
2060algmul(al,m1,m1i) == midbasis
2061algmul(al,m2i,m2) == midbasis
2062algmul(al,m2,m2i) == midbasis
2063algmul(al,m3i,m3) == midbasis
2064algmul(al,m3,m3i) == midbasis
2065print("mul");
2066algmul(al,algmul(al,m1,m2),m3) == algmul(al,m1,algmul(al,m2,m3))
2067algmul(al,m1,algadd(al,m2,m3)) == algadd(al, algmul(al,m1,m2), algmul(al,m1,m3))
2068algmul(al,algadd(al,m1,m2),m3) == algadd(al,algmul(al,m1,m3),algmul(al,m2,m3))
2069print("neg");
2070algadd(al,m1,algneg(al,m1)) == 0
2071algadd(al,m2,algneg(al,m2)) == 0
2072algadd(al,m3,algneg(al,m3)) == 0
2073print("norm");
2074algnorm(al,m1)*algnorm(al,m2) == algnorm(al,algmul(al,m1,m2))
2075algnorm(al,m1)*algnorm(al,m3) == algnorm(al,algmul(al,m1,m3))
2076algnorm(al,m3)*algnorm(al,m2) == algnorm(al,algmul(al,m3,m2))
2077print("pow");
2078algpow(al,m1,5) == algmul(al,algpow(al,m1,2),algpow(al,m1,3))
2079algalgtobasis(al,algpow(al,m2,3)) == algmul(al,algpow(al,m2,5),algpow(al,m2,-2))
2080algpow(al,m3,3) == algmul(al,algpow(al,m3,-1),algpow(al,m3,4))
2081print("sqr");
2082algsqr(al,m1) == algpow(al,m1,2)
2083algsqr(al,m2) == algpow(al,m2,2)
2084algsqr(al,m3) == algpow(al,m3,2)
2085print("sub");
2086algsub(al,m2,m3) == algadd(al,m2,algneg(al,m3))
2087algsub(al,m1,m3) == algadd(al,m1,algneg(al,m3))
2088algsub(al,m1,m2) == algadd(al,m1,algneg(al,m2))
2089print("trace");
2090algtrace(al,m2)+algtrace(al,m3) == algtrace(al,algadd(al,m2,m3))
2091algtrace(al,m1)+algtrace(al,m3) == algtrace(al,algadd(al,m1,m3))
2092algtrace(al,m1)+algtrace(al,m2) == algtrace(al,algadd(al,m1,m2))
2093print("algtomatrix");
2094sm1 = algtomatrix(al,m1);
2095sm2 = algtomatrix(al,m2);
2096sm3 = algtomatrix(al,m3);
2097#sm1==4
2098#sm1[,1]==4
2099sm2+sm3 == algtomatrix(al,algadd(al,m2,m3))
2100sm2*sm3 == algtomatrix(al,algmul(al,m2,m3))
2101sm1+sm3 == algtomatrix(al,algadd(al,m1,m3))
2102sm1*sm3 == algtomatrix(al,algmul(al,m1,m3))
2103sm1+sm2 == algtomatrix(al,algadd(al,m1,m2))
2104sm1*sm2 == algtomatrix(al,algmul(al,m1,m2))
2105print("algleftmultable");
2106M1 = algtomatrix(al,m1,1);
2107M2 = algtomatrix(al,m2,1);
2108M3 = algtomatrix(al,m3,1);
2109#M1==32
2110#M1[,1]==32
2111#M2==32
2112#M2[,1]==32
2113#M3==32
2114#M3[,1]==32
2115M2+M3 == algtomatrix(al,algadd(al,m2,m3),1)
2116M2*M3 == algtomatrix(al,algmul(al,m2,m3),1)
2117M1+M3 == algtomatrix(al,algadd(al,m1,m3),1)
2118M1*M3 == algtomatrix(al,algmul(al,m1,m3),1)
2119M1+M2 == algtomatrix(al,algadd(al,m1,m2),1)
2120M1*M2 == algtomatrix(al,algmul(al,m1,m2),1)
2121
2122
2123print("more al_CSA tests");
2124setrand(1);
2125nf = nfinit(y^2-5);
2126mti = [0,-1,0,0;[1,0]~,0,0,0;0,0,0,-1;0,0,1,0];
2127mtj = [0,0,Mod(y,y^2-5),0;0,0,0,-y;1,0,0,0;0,-1,0,0];
2128mtk = [0,0,0,y;0,0,[1,2]~,0;0,1,0,0;1,0,0,0];
2129mt = [matid(4),mti,mtj,mtk];
2130al = alginit(nf,mt);
2131a1 = [1/2,-3,2,4/3,-1,0,1,2/5]~;
2132a2 = [1+2*y,y/3,5/2,-3*y/7]~*Mod(1,y^2-5);
2133a3 = algrandom(al,2)/13;
2134aid = [1,0,0,0]~;
2135algadd(al,algadd(al,a1,a2),a3) == algadd(al,a1,algadd(al,a2,a3))
2136algadd(al,a1,a2) == algadd(al,a2,a1)
2137algbasistoalg(al,algalgtobasis(al,algbasistoalg(al,a1))) == algbasistoalg(al,a1)
2138algalgtobasis(al,algbasistoalg(al,algalgtobasis(al,a2))) == algalgtobasis(al,a2)
2139t = 7/3;
2140testcp(a)=
2141{
2142  my(P,Q);
2143  P = algcharpoly(al,a);
2144  print(algnorm(al,a) == polcoeff(P,0));
2145  print(algtrace(al,a) == -polcoeff(P,1));
2146  print(algnorm(al,algsub(al,a,t*aid))==subst(P,x,t));
2147  Q = algcharpoly(al,a,,1);
2148  print(algnorm(al,a,1) == polcoeff(Q,0));
2149  print(algtrace(al,a,1) == -polcoeff(Q,7));
2150  print(algnorm(al,algsub(al,a,t*aid),1)==subst(Q,x,t));
2151}
2152print("charpoly");
2153testcp(a1);
2154testcp(a2);
2155testcp(a3);
2156print("inv/div");
2157algisinv(al,a1,&a1i)
2158a1i == alginv(al,a1)
2159algdivl(al,a1,a2) == algmul(al,a1i,a2)
2160algdivr(al,a2,a1) == algmul(al,a2,a1i)
2161algisdivl(al,a1,a3,&d13)
2162algdivl(al,a1,a3) == d13
2163aidbasis = vector(8,i,i==1)~;
2164algmul(al,a1i,a1) == aidbasis
2165algmul(al,a1,a1i) == aidbasis
2166a2i = alginv(al,a2);
2167print("mul");
2168algmul(al,a2i,a2) == aid
2169algmul(al,a2,a2i) == aid
2170algmul(al,algmul(al,a1,a2),a3) == algmul(al,a1,algmul(al,a2,a3))
2171print("neg");
2172algadd(al,a1,algneg(al,a1)) == 0
2173algadd(al,a2,algneg(al,a2)) == 0
2174algadd(al,a3,algneg(al,a3)) == 0
2175print("norm");
2176algnorm(al,a1)*algnorm(al,a2) == algnorm(al,algmul(al,a1,a2))
2177algnorm(al,a1)*algnorm(al,a3) == algnorm(al,algmul(al,a1,a3))
2178algnorm(al,a3)*algnorm(al,a2) == algnorm(al,algmul(al,a3,a2))
2179print("pow");
2180algpow(al,a1,5) == algmul(al,algpow(al,a1,2),algpow(al,a1,3))
2181algpow(al,a2,3) == algmul(al,algpow(al,a2,5),algpow(al,a2,-2))
2182algpow(al,a3,3) == algmul(al,algpow(al,a3,-1),algpow(al,a3,4))
2183print("sqr");
2184algsqr(al,a1) == algpow(al,a1,2)
2185algsqr(al,a2) == algpow(al,a2,2)
2186algsqr(al,a3) == algpow(al,a3,2)
2187print("sub");
2188algsub(al,a2,a3) == algadd(al,a2,algneg(al,a3))
2189algsub(al,a1,a3) == algadd(al,a1,algneg(al,a3))
2190algsub(al,a1,a2) == algadd(al,a1,algneg(al,a2))
2191print("trace");
2192algtrace(al,a2)+algtrace(al,a3) == algtrace(al,algadd(al,a2,a3))
2193algtrace(al,a1)+algtrace(al,a3) == algtrace(al,algadd(al,a1,a3))
2194algtrace(al,a1)+algtrace(al,a2) == algtrace(al,algadd(al,a1,a2))
2195print("algtomatrix");
2196sa1 = algtomatrix(al,a1);
2197sa2 = algtomatrix(al,a2);
2198sa3 = algtomatrix(al,a3);
2199#sa1==2
2200#sa1[,1]==2
2201sa2+sa3 == algtomatrix(al,algadd(al,a2,a3))
2202sa2*sa3 == algtomatrix(al,algmul(al,a2,a3))
2203sa1+sa3 == algtomatrix(al,algadd(al,a1,a3))
2204sa1*sa3 == algtomatrix(al,algmul(al,a1,a3))
2205sa1+sa2 == algtomatrix(al,algadd(al,a1,a2))
2206sa1*sa2 == algtomatrix(al,algmul(al,a1,a2))
2207print("algleftmultable");
2208ma1 = algtomatrix(al,a1,1);
2209ma2 = algtomatrix(al,a2,1);
2210ma3 = algtomatrix(al,a3,1);
2211#ma1==8
2212#ma1[,1]==8
2213#ma2==8
2214#ma2[,1]==8
2215ma2+ma3 == algtomatrix(al,algadd(al,a2,a3),1)
2216ma2*ma3 == algtomatrix(al,algmul(al,a2,a3),1)
2217ma1+ma3 == algtomatrix(al,algadd(al,a1,a3),1)
2218ma1*ma3 == algtomatrix(al,algmul(al,a1,a3),1)
2219ma1+ma2 == algtomatrix(al,algadd(al,a1,a2),1)
2220ma1*ma2 == algtomatrix(al,algmul(al,a1,a2),1)
2221
2222
2223print("csa pol/polmod bugs");
2224setrand(1);
2225nf = nfinit(y^2-5);
2226mti = [0,-1,0,0;1,0,0,0;0,0,0,-1;0,0,1,0];
2227mtj = [0,0,y,0;0,0,0,-y;1,0,0,0;0,-1,0,0];
2228mtk = [0,0,0,y;0,0,y,0;0,1,0,0;1,0,0,0];
2229mt = [matid(4),mti,mtj,mtk];
2230al = alginit(nf,mt*Mod(1,nf.pol));
2231algrelmultable(al)
2232a = [y,y,y,y/3]~;
2233algpow(al,a,4)
2234algpow(al,a,-2)
2235algmul(al,a,a)
2236algnorm(al,a) == algnorm(al,algalgtobasis(al,a))
2237algnorm(al,a,1) == algnorm(al,algalgtobasis(al,a),1)
2238algtrace(al,a) == algtrace(al,algalgtobasis(al,a))
2239algtrace(al,a,1) == algtrace(al,algalgtobasis(al,a),1)
2240algcharpoly(al,a) == algcharpoly(al,algalgtobasis(al,a))
2241algcharpoly(al,a,,1) == algcharpoly(al,algalgtobasis(al,a),,1)
2242algtomatrix(al,a) == algtomatrix(al,algalgtobasis(al,a))
2243
2244print("csa: denom over Z[y] but not over ZK");
2245setrand(1);
2246nf = nfinit(y^2-5);
2247mti = [0,-1,0,0;1,0,0,0;0,0,0,-1;0,0,1,0];
2248mtj = [0,0,(y-1)/2,0;0,0,0,-(y-1)/2;1,0,0,0;0,-1,0,0];
2249mtk = [0,0,0,(y-1)/2;0,0,(y-1)/2,0;0,1,0,0;1,0,0,0];
2250mt = [matid(4),mti,mtj,mtk];
2251al = alginit(nf,mt*Mod(1,nf.pol));
2252algrelmultable(al)
2253mti = [0,-1,0,0;1,0,0,0;0,0,0,-1;0,0,1,0];
2254mtj = [0,0,(y-1)/3,0;0,0,0,-(y-1)/3;1,0,0,0;0,-1,0,0];
2255mtk = [0,0,0,(y-1)/3;0,0,(y-1)/3,0;0,1,0,0;1,0,0,0];
2256mt = [matid(4),mti,mtj,mtk];
2257al = alginit(nf,mt*Mod(1,nf.pol));
2258
2259print("al_MAT over al_CSA");
2260setrand(1);
2261nf = nfinit(y^2-5);
2262mti = [0,-1,0,0;1,0,0,0;0,0,0,-1;0,0,1,0];
2263mtj = [0,0,(y-1)/2,0;0,0,0,-(y-1)/2;1,0,0,0;0,-1,0,0];
2264mtk = [0,0,0,(y-1)/2;0,0,(y-1)/2,0;0,1,0,0;1,0,0,0];
2265mt = [matid(4),mti,mtj,mtk];
2266al = alginit(nf,mt*Mod(1,nf.pol));
2267m1 = [algrandom(al,2)/2,algrandom(al,1)/3;algrandom(al,2)/5,algrandom(al,3)];
2268m2 = [[y,1/2,2-y,0]~,[3,y/2,-3*y,1]~;[7,0,0,y]~,[1-y,2,-y,0]~];
2269algnorm(al,algmul(al,m1,m2)) == algnorm(al,m1)*algnorm(al,m2)
2270algnorm(al,algmul(al,m1,m2),1) == algnorm(al,m1,1)*algnorm(al,m2,1)
2271algtrace(al,algadd(al,m1,m2)) == algtrace(al,m1)+algtrace(al,m2)
2272algtrace(al,algadd(al,m1,m2),1) == algtrace(al,m1,1)+algtrace(al,m2,1)
2273mid = [[1,0,0,0]~,vector(4)~;vector(4)~,[1,0,0,0]~];
2274t = -3/7;
2275testcp(m)=
2276{
2277  my(P,Q);
2278  P = algcharpoly(al,m);
2279  print(algnorm(al,m) == polcoeff(P,0));
2280  print(algtrace(al,m) == -polcoeff(P,3));
2281  print(algnorm(al,algsub(al,m,t*mid))==subst(P,x,t));
2282  Q = algcharpoly(al,m,,1);
2283  print(algnorm(al,m,1) == polcoeff(Q,0));
2284  print(algtrace(al,m,1) == -polcoeff(Q,31));
2285  print(algnorm(al,algsub(al,m,t*mid),1)==subst(Q,x,t));
2286}
2287testcp(m1);
2288testcp(m2);
2289print("algleftmultable");
2290M1 = algtomatrix(al,m1,1);
2291M2 = algtomatrix(al,m2,1);
2292#M1==32
2293#M1[,1]==32
2294#M2==32
2295#M2[,1]==32
2296M1+M2 == algtomatrix(al,algadd(al,m1,m2),1)
2297M1*M2 == algtomatrix(al,algmul(al,m1,m2),1)
2298
2299
2300
2301print("nfgrunwaldwang SEGV #1669");
2302nfgrunwaldwang(nfinit(y),[2,3],[1,2],Vecsmall(1))
2303nfgrunwaldwang(nfinit(x),[2,3],[1,2],Vecsmall(1))
2304
2305\\alg_model bug for 1-dim al_CSA
2306al = alginit(nfinit(y),[matid(1)]);
2307algtomatrix(al,[Mod(1,y)]~)
2308algtomatrix(al,[1]~)
2309algtomatrix(al,[1+y-y]~)
2310algtomatrix(al,[1/2]~)
2311algtomatrix(al,[Mod(1/2,y)]~)
2312
2313\\al_MATRIX sqr bug
2314al = algtableinit([matid(2),[0,0;1,0]],2);
2315m = [[1,1]~,[1,0]~;[0,1]~,[0,0]~];
2316type(algsqr(al,m))=="t_MAT"
2317
2318print("GW modified arguments");
2319L = [2,3];
2320Q = nfinit('z);
2321gw = nfgrunwaldwang(Q,L,[2 | p <- L],[1],'y);
2322print(L==[2,3]);
2323
2324\\algpoleval type check
2325setrand(1);
2326nf = nfinit(y^2-5);
2327al = alginit(nf,[y,-1]);
2328a = [1..8]~;
2329pol = algcharpoly(al,a);
2330algpoleval(al,pol,a)==0
2331algpoleval(al,pol,[;])
2332pol = algcharpoly(al,a,,1);
2333algpoleval(al,pol,a)==0
2334ma = algtomatrix(al,a,1);
2335algpoleval(al,pol,[a,ma])==0
2336
2337\\not implemented
2338al = algtableinit([Mat(1)]);
2339al2 = algtensor(al,al);
2340al = algtableinit([Mat(1)],2);
2341al2 = algtensor(al,al);
2342
2343
2344