1setrand(1)
2
3vec2mat(d,n,V) =
4{
5  my(M = matrix(d,d), x = 'x, c);
6  for(i=0,d^2*n-1,
7    c = V[i+1];
8    M[i\(n*d) + 1, (i\n)%d + 1] += c*x^(i%n)
9  );
10  c = polcoeff(M[1,1],0);
11  for(i=2,d,M[i,i]+=c);
12  M
13};
14
15mat2vec(d,n,M) =
16{
17  my(V = vector(d^2*n), c = polcoeff(M[1,1],0));
18  for(i=2,d,M[i,i]-=c);
19  for(i=1,d,
20    for(j=1,d,
21      for(k=0,n-1,
22        V[n*d*(i-1) + n*(j-1) + k + 1] = polcoeff(M[i,j],k);
23      )
24    )
25  );
26  V
27};
28
29matalg(d,n,p) =
30{
31  my(pol = ffinit(p,n), mtx, N=d^2*n, x, basis, mt=vector(N));
32  basis = vector(N,i,vec2mat(d,n,vector(N,k,k==i)));
33  basis = basis*Mod(1,p)*Mod(1,pol);
34  for(i=1,N,
35    x = basis[i];
36    mtx = matconcat(vector(N,j,mat2vec(d,n,liftall(x*basis[j]))~));
37    mt[i] = mtx;
38  );
39  algtableinit(mt,p)
40};
41
42smallchg(mt,p) =
43{
44  my(i,j,N=#mt,c);
45  if(N==1, return(mt));
46  if(p, c=random([1,p-1]), c=1);
47  i = random([2,N]);
48  j=i;
49  while(j==i,j=random([1,N]));
50  for(k=1,N,
51    mt[k][,i] += c*mt[k][,j];
52    mt[k][j,] -= c*mt[k][i,];
53  );
54  mt[i] += c*mt[j];
55  if(p,mt%p,mt)
56};
57
58chg(mt,nb=3*#mt,p=0) =
59{
60  for(c=1,nb,mt=smallchg(mt,p));
61  mt
62};
63
64test2(d,n,p,nb=1) =
65{
66  my(al,mt,mt2,res,map,x,y,Mx,My);
67  setrand(1);
68  al = matalg(d,n,p);
69  mt = algmultable(al);
70  mt2 = mt;
71  for(c=1,nb,
72    if(c>1, mt2 = chg(mt2,#mt\4 + 5,p));
73    al = algtableinit(mt2,p);
74    res = algsplit(al,'t);
75    map = res[1];
76    mapi = res[2];
77    x = algrandom(al,10);
78    y = algrandom(al,10);
79    Mx = map*x;
80    My = map*y;
81    if(Mx*My != map*algmul(al,x,y),  print("FAIL"); return(0));
82  );
83  [al,map,mapi]
84};
85
86print("M_2(F_4)");
87test2(2,2,2,200);
88
89print("M_2(F_9)");
90test2(2,2,3,100);
91
92q = 31;
93print("M_2(F_31^2)");
94test2(2,2,q,30);
95
96p = nextprime(10^20);
97print("M_2(F_p^2)");
98test2(2,2,p,30);
99
100print("M_d(F_31^n)");
101for(d=1,4,for(n=1,min(18\(d^2),10),print("d=",d," n=",n);if(!test2(d,n,q,10), break(2))));
102
103isblock(M,L) =
104{
105  my(d,mini,nxt);
106  d = #M[,1];
107  mini=L[1];
108  nxt=2;
109  for(j=1,d,
110    if(j>mini,
111      mini += L[nxt];
112      nxt += 1;
113    );
114    for(i=mini+1,d,
115      if(M[i,j]!=0,return(0));
116    );
117  );
118  1
119};
120matblock(L,n,p) =
121{
122  d = sum(i=1,#L,L[i]);
123  al0 = matalg(d,n,p);
124  basis = [b | b<-Vec(matid(d^2*n)), isblock(vec2mat(d,n,b),L)];
125  basis = Mat(basis);
126  algsubalg(al0,basis)[1];
127};
128
129print("examples from documentation");
130al0 = alginit(nfinit(y^2+7), [-1,-1]);
131al = algtableinit(algmultable(al0), 3); \\isomorphic to M_2(F_9)
132[map,mapi] = algsplit(al, 'a);
133x = [1,2,1,0,0,0,0,0]~; fx = map*x
134y = [0,0,0,0,1,0,0,1]~; fy = map*y
135map*algmul(al,x,y) == fx*fy
136map*mapi[,6]
137
138print("bad input");
139algsplit("toto");
140algsplit(alginit(nfinit('y),[-1,-1]));
141algsplit(algtableinit([matid(3),[0,0,0;1,1,0;0,0,1],[0,0,0;0,0,0;1,0,0]],2));
142algsplit(algtableinit([matid(2),[0,0;1,1]],2));
143algsplit(matblock([2,1,1,1],1,q));
144algsplit(matblock([1,2,3],1,q));
145{mt=[[1, 0, 0, 0, 0, 0, 0, 0, 0; 0, 1, 0, 0, 0, 0, 0, 0, 0; 0, 0, 1, 0, 0, 0, 0,
1460, 0; 0, 0, 0, 1, 0, 0, 0, 0, 0; 0, 0, 0, 0, 1, 0, 0, 0, 0; 0, 0, 0, 0, 0, 1, 0,
1470, 0; 0, 0, 0, 0, 0, 0, 1, 0, 0; 0, 0, 0, 0, 0, 0, 0, 1, 0; 0, 0, 0, 0, 0, 0, 0,
1480, 1], [0, 0, 0, 1, 0, 0, 0, 0, 0; 1, 0, 0, 0, 0, 0, 0, 0, 0; 0, 0, 0, 0, 0, 1,
1490, 0, 0; 0, 1, 0, 0, 0, 0, 0, 0, 0; 0, 0, 0, 1, 0, 0, 0, 0, 1; 0, 0, 0, 1, 0, 0,
1501, 0, 0; 0, 1, 1, 0, 0, 0, 0, 0, 0; 0, 0, 0, 0, 1, 0, 0, 0, 0; 0, 1, 0, 0, 0, 0,
1510, 1, 0], [0, 0, 0, 0, 0, 0, 1, 0, 1; 0, 0, 0, 0, 1, 1, 0, 0, 0; 1, 0, 0, 0, 1,
1520, 0, 0, 0; 0, 1, 0, 0, 0, 0, 0, 0, 0; 0, 0, 0, 0, 0, 0, 1, 0, 0; 0, 0, 0, 1, 0,
1530, 1, 0, 0; 0, 0, 0, 0, 0, 0, 0, 1, 0; 0, 0, 0, 0, 1, 0, 0, 0, 0; 0, 0, 1, 0, 0,
1540, 0, 0, 0], [0, 0, 1, 0, 0, 0, 0, 1, 0; 0, 0, 0, 0, 0, 0, 1, 0, 1; 0, 0, 0, 0,
1550, 0, 0, 0, 1; 1, 0, 0, 0, 0, 0, 0, 0, 0; 0, 1, 1, 0, 0, 0, 0, 0, 0; 0, 0, 1, 0,
1560, 0, 0, 0, 0; 0, 0, 0, 0, 1, 0, 0, 0, 0; 0, 0, 0, 1, 0, 0, 0, 0, 1; 0, 0, 0, 0,
1570, 1, 0, 0, 0], [0, 0, 0, 0, 1, 1, 0, 0, 0; 0, 0, 1, 0, 0, 0, 0, 1, 0; 0, 1, 0,
1580, 0, 0, 0, 1, 0; 0, 0, 0, 1, 0, 0, 0, 0, 0; 1, 0, 0, 0, 0, 1, 0, 0, 0; 0, 0, 0,
1590, 0, 1, 0, 0, 0; 0, 0, 0, 1, 0, 0, 0, 0, 1; 0, 0, 0, 0, 0, 0, 0, 1, 0; 0, 0, 0,
1601, 0, 0, 1, 0, 0], [0, 0, 0, 0, 1, 1, 0, 0, 0; 0, 1, 0, 0, 0, 0, 0, 0, 0; 0, 1,
1610, 0, 0, 0, 0, 1, 0; 0, 0, 0, 0, 0, 0, 1, 0, 1; 0, 0, 0, 0, 1, 0, 0, 0, 0; 1, 0,
1620, 0, 1, 0, 0, 0, 0; 0, 0, 0, 1, 0, 0, 0, 0, 1; 0, 1, 1, 0, 0, 0, 0, 0, 0; 0, 0,
1630, 0, 0, 0, 0, 0, 1], [0, 1, 0, 0, 0, 0, 0, 0, 0; 0, 0, 0, 0, 0, 0, 1, 0, 1; 0,
1640, 0, 1, 0, 0, 1, 0, 0; 0, 0, 0, 0, 1, 1, 0, 0, 0; 0, 1, 1, 0, 0, 0, 0, 0, 0; 0,
1651, 0, 0, 0, 0, 0, 1, 0; 1, 0, 0, 0, 0, 1, 0, 0, 0; 0, 0, 0, 0, 0, 0, 1, 0, 0; 0,
1660, 0, 0, 0, 1, 0, 0, 0], [0, 0, 0, 1, 0, 0, 0, 0, 0; 0, 0, 0, 0, 1, 1, 0, 0, 0;
1670, 0, 0, 0, 0, 1, 0, 0, 0; 0, 0, 1, 0, 0, 0, 0, 1, 0; 0, 0, 0, 0, 0, 0, 1, 0, 0;
1680, 0, 0, 0, 0, 0, 0, 0, 1; 0, 1, 1, 0, 0, 0, 0, 0, 0; 1, 0, 0, 0, 0, 1, 0, 0, 0;
1690, 0, 1, 0, 0, 0, 0, 0, 0], [0, 0, 1, 0, 0, 0, 0, 1, 0; 0, 0, 0, 1, 0, 0, 0, 0,
1700; 0, 0, 0, 0, 0, 0, 0, 0, 1; 0, 0, 0, 0, 1, 1, 0, 0, 0; 0, 0, 0, 0, 0, 0, 0, 1,
1710; 0, 1, 0, 0, 0, 0, 0, 1, 0; 0, 0, 0, 0, 1, 0, 0, 0, 0; 0, 0, 0, 0, 0, 0, 1, 0,
1720; 1, 0, 0, 0, 1, 0, 0, 0, 0]]};
173algsplit(algtableinit(mt,2));
174
175