1#############################################################################
2##
3##  This file is part of GAP, a system for computational discrete algebra.
4##  This file's authors include A. Storjohann, R. Wainwright, F. Gähler, D. Holt, A. Hulpke.
5##
6##  Copyright of GAP belongs to its developers, whose names are too numerous
7##  to list here. Please refer to the COPYRIGHT file for details.
8##
9##  SPDX-License-Identifier: GPL-2.0-or-later
10##
11##  This file contains functions that compute Hermite and Smith normal forms
12##  of integer matrices, with or without the HNF/SNF  expressed as the linear
13##  combination of the input.
14##
15
16########################################################
17##
18##      auxiliary + main code for all in one function
19##
20##  MATINTsplit
21##  MATINTrgcd
22##  MATINTmgcdex
23##  MATINTbezout
24##  SNFofREF
25##  NormalFormIntMat
26##
27
28########################################################
29#
30# MATINTsplit(<N>,<a>) - returns product of prime factors of N which are not factors of a.
31#
32BindGlobal("MATINTsplit",function(N,a)
33local x,t;
34
35  x:=a;
36  t:=N;
37  while x<>1 do
38    x:=GcdInt(x,t);
39    t:=QuoInt(t,x);
40  od;
41  return t;
42end);
43
44################################################
45#
46# MATINTrgcd(<N>,<a>) - Returns smallest nonnegative c such that gcd(N,a+c) = 1
47#
48BindGlobal("MATINTrgcd",function(N,a)
49local k,r,d,i,c,g,q;
50
51  if N=1 then return 0; fi;
52  k := 1;
53  r:=[(a-1) mod N];
54  d:=[N];
55  c:=0;
56  while true do
57    for i in [1..k] do r[i]:=(r[i]+1) mod d[i]; od;
58    i:=PositionProperty(r,x->x<=0);
59    if i=fail then
60      g:=1;i:=0;
61      while g=1 and i<k do
62        i:=i+1;
63        g:=GcdInt(r[i],d[i]);
64      od;
65      if g=1 then return c; fi;
66      q:=MATINTsplit(QuoInt(d[i],g),g);
67      if q>1 then
68        k:=k+1;
69        r[k]:=r[i] mod q;
70        d[k]:=q;
71      fi;
72      r[i]:=0;
73      d[i]:=g;
74    fi;
75    c:=c+1;
76  od;
77
78end);
79
80#######################################################
81#
82#  MATINTmgcdex(<N>,<a>,<v>) - Returns c[1],c[2],...c[k] such that
83#
84#   gcd(N,a+c[1]*b[1]+...+c[n]*b[k]) = gcd(N,a,b[1],b[2],...,b[k])
85#
86BindGlobal("MATINTmgcdex", function(N,a,v)
87local h,g,M,c,i,d,b,l;
88  l:=Length(v);
89  c:=[]; M:=[];
90  h := N;
91  for i in [1..l] do
92    g := h;
93    h:=GcdInt(g,v[i]);
94    M[i]:=QuoInt(g,h);
95  od;
96  h:=GcdInt(a,h);
97  g:=QuoInt(a,h);
98
99  for i in [l,l-1..1] do
100    b:=QuoInt(v[i],h);
101    d:=MATINTsplit(M[i],b);
102    if d=1 then
103      c[i]:=0;
104    else
105      c[i]:=MATINTrgcd(d,g/b mod d);
106      g:=g+c[i]*b;
107    fi;
108  od;
109
110  return c;
111
112end);
113
114
115
116
117#####################################################
118#
119#  MATINTbezout(a,b,c,d) - returns row transform , P, to transform, A, to hnf :
120#
121#  PA=H;
122#
123#  [ s  t ] [ a  b ]   [ e  f ]
124#  [      ] [      ] = [       ]
125#  [ u  v ] [ c  d ]   [    g ]
126#
127BindGlobal("MATINTbezout", function(a,b,c,d)
128local e,f,g,q;
129
130  e := Gcdex(a,c);
131  f := e.coeff1*b+e.coeff2*d;
132  g := e.coeff3*b+e.coeff4*d;
133  if g<0 then
134    e.coeff3 := -e.coeff3;
135    e.coeff4 := -e.coeff4;
136    g := -g;
137  fi;
138  if g>0 then
139    q := QuoInt(f-(f mod g),g);
140    e.coeff1 := e.coeff1-q*e.coeff3;
141    e.coeff2 := e.coeff2-q*e.coeff4;
142  fi;
143  return e;
144
145end);
146
147#####################################################
148##
149## SNFofREF - fast SNF of REF matrix
150##
151##
152InstallGlobalFunction(SNFofREF , function(R,destroy)
153local k,g,b,ii,m1,m2,t,tt,si,n,m,i,j,r,jj,piv,d,gt,tmp,A,T,TT,kk;
154
155  Info(InfoMatInt,1,"SNFofREF - initializing work matrix");
156  n := NrRows(R);
157  m := NrCols(R);
158
159  piv := List(R,x->PositionProperty(x,y->y<>0));
160  r := PositionProperty(piv,x->x=fail);
161  if r=fail then
162    r := Length(piv);
163  else
164    r := r-1;
165    piv := piv{[1..r]};
166  fi;
167  Append(piv,Difference([1..m],piv));
168
169  if destroy then
170    T:=R;
171    ##  Need to be careful: we are trying to permute the cols in place
172    for i in [1..r] do
173      T[i]{[1..m]} := T[i]{piv};
174    od;
175  else
176    T := NullMat(n,m);
177    for j in [1..m] do
178      for i in [1..Minimum(r,j)] do
179        T[i,j]:=R[i,piv[j]];
180      od;
181    od;
182  fi;
183
184  si := 1;
185  A := [];
186  d := 2;
187  for k in [1..m] do
188    Info(InfoMatInt,2,"SNFofREF - working on column ",k);
189    if k<=r then
190      d := d*AbsInt(T[k,k]);
191      Apply(T[k],x->x mod (2*d));
192    fi;
193
194    t := Minimum(k,r);
195    for i in [t-1,t-2..si] do
196      t := MATINTmgcdex(A[i],T[i,k],[T[i+1,k]])[1];
197      if t<>0 then
198        AddRowVector(T[i],T[i+1],t);
199        Apply(T[i],x->x mod A[i]);
200      fi;
201    od;
202
203    for i in [si..Minimum(k-1,r)] do
204      g := Gcdex(A[i],T[i,k]);
205      T[i,k] := 0;
206      if g.gcd<>A[i] then
207        b := QuoInt(A[i],g.gcd);
208        A[i] := g.gcd;
209        for ii in [i+1..Minimum(k-1,r)] do
210          AddRowVector(T[ii],T[i],-g.coeff2*QuoInt(T[ii,k],A[i]) mod A[ii]);
211          T[ii,k] := b*T[ii,k];
212
213          Apply(T[ii],x->x mod A[ii]);
214        od;
215        if k<=r then
216          t := g.coeff2*QuoInt(T[k,k],g.gcd);
217          AddRowVector(T[k],T[i],-t);
218          T[k,k]:=b*T[k,k];
219        fi;
220        Apply(T[i],x->x mod A[i]);
221        if A[i]=1 then si := i+1; fi;
222      fi;
223    od;
224
225    if k<=r then
226      A[k] := AbsInt(T[k,k]);
227      Apply(T[k],x->x mod A[k]);
228    fi;
229
230  od;
231
232  for i in [1..r] do T[i,i] := A[i]; od;
233
234  return T;
235
236end);
237
238BindGlobal("BITLISTS_NFIM", MakeImmutable(
239  [ [ false, false, false, false, false ], [ true, false, false, false, false ],
240    [ false, true, false, false, false ], [ true, true, false, false, false ],
241    [ false, false, true, false, false ], [ true, false, true, false, false ],
242    [ false, true, true, false, false ], [ true, true, true, false, false ],
243    [ false, false, false, true, false ], [ true, false, false, true, false ],
244    [ false, true, false, true, false ], [ true, true, false, true, false ],
245    [ false, false, true, true, false ], [ true, false, true, true, false ],
246    [ false, true, true, true, false ], [ true, true, true, true, false ],
247    [ false, false, false, false, true ], [ true, false, false, false, true ],
248    [ false, true, false, false, true ], [ true, true, false, false, true ],
249    [ false, false, true, false, true ], [ true, false, true, false, true ],
250    [ false, true, true, false, true ], [ true, true, true, false, true ],
251    [ false, false, false, true, true ], [ true, false, false, true, true ],
252    [ false, true, false, true, true ], [ true, true, false, true, true ],
253    [ false, false, true, true, true ], [ true, false, true, true, true ],
254    [ false, true, true, true, true ], [ true, true, true, true, true ] ] ));
255
256
257###########################################################
258#
259# DoNFIM(<mat>,<options>)
260#
261# Options bit values:
262#
263# 1  - Triangular / Smith
264# 2  - No / Yes  Reduce off diag entries
265# 4  - No / Yes  Row Transforms
266# 8  - No / Yes  Col Transforms
267# 16 - change original matrix in place (The rows still change) -- save memory
268#
269# Compute a Triangular, Hermite or Smith form of the n x m
270# integer input matrix A.  Optionally, compute n x n / m x m
271# unimodular transforming matrices which satisfy Q C A = H
272# or  Q C A B P = S.
273#
274# Triangular / Hermite :
275#
276# Let I be the min(r+1,n) x min(r+1,n) identity matrix with r = rank(A).
277# Then Q and C can be written using a block decomposition as
278#
279#             [ Q1 |   ]  [ C1 | C2 ]
280#             [----+---]  [----+----]  A  =  H.
281#             [ Q2 | I ]  [    | I  ]
282#
283# Smith :
284#
285#  [ Q1 |   ]  [ C1 | C2 ]     [ B1 |   ]  [ P1 | P2 ]
286#  [----+---]  [----+----]  A  [----+---]  [----+----] = S.
287#  [ Q2 | I ]  [    | I  ]     [ B2 | I ]  [ *  | I  ]
288#
289# * - possible non-zero entry in upper right corner...
290#
291#
292BindGlobal("DoNFIM", function(mat, options)
293local opt, sig, n, m, A, C, Q, B, P, r, c2, rp, c1, j, k, N, L, b, a, g, c,
294      t, tmp, i, q, R, rank, signdet;
295
296  if not (IsMatrix(mat)
297         or (IsList(mat) and Length(mat)=1
298             and IsList(mat[1]) and Length(mat[1])=0))
299     or not IsInt(options) then
300    Error("syntax is DoNFIM(<mat>,<options>)");
301  fi;
302
303  #Parse options
304  opt := BITLISTS_NFIM[options+1];
305  #List(CoefficientsQadic(options,2),x->x=1);
306  #if Length(opt)<4 then
307  #  opt{[Length(opt)+1..4]} := List([Length(opt)+1..4],x->false);
308  #fi;
309
310  sig:=1;
311
312  #Embed matrix in 2 larger "id" matrix
313  n := NrRows(mat)+2;
314  m := NrCols(mat)+2;
315  k:=ListWithIdenticalEntries(m,0);
316  if opt[5] then
317    # change the matrix
318    A:=mat;
319    for i in [n-1,n-2..2] do
320      A[i]:=ShallowCopy(k);
321      A[i]{[2..m-1]}:=A[i-1];
322    od;
323  else
324    A := [];
325    for i in [2..n-1] do
326      #A[i] := [0];
327      #Append(A[i],mat[i-1]);
328      #A[i,m] := 0;
329      A[i]:=ShallowCopy(k);
330      A[i]{[2..m-1]}:=mat[i-1];
331    od;
332  fi;
333  A[1]:=ShallowCopy(k);
334  A[n]:=k;
335  A[1,1] := 1;
336  A[n,m] := 1;
337
338  if opt[3] then
339    C := IdentityMat(n);
340    Q := NullMat(n,n);
341    Q[1,1] := 1;
342  fi;
343
344  if opt[1] and opt[4] then
345    B := IdentityMat(m);
346    P := IdentityMat(m);
347  fi;
348
349  r := 0;
350  c2 := 1;
351  rp := [];
352  while m>c2 do
353    Info(InfoMatInt,2,"DoNFIM - reached column ",c2," of ",m);
354    r := r+1;
355    c1 := c2;
356    rp[r] := c1;
357    if opt[3] then Q[r+1,r+1] := 1; fi;
358
359    j := c1+1;
360    while j<=m do
361      k := r+1;
362      while k<=n and A[r,c1]*A[k,j]=A[k,c1]*A[r,j] do k := k+1; od;
363      if k<=n then c2 := j; j := m; fi;
364      j := j+1;
365    od;
366    #Smith with some transforms..
367    if opt[1] and (opt[4] or opt[3]) and c2<m then
368      N := Gcd(Flat(A{[r..n]}[c2]));
369      L := [c1+1..c2-1];
370      Append(L,[c2+1..m-1]);
371      Add(L,c2);
372      for j in L do
373        if j=c2 then
374          b:=A[r,c2];a:=A[r,c1];
375          for i in [r+1..n] do
376            if b<>1 then
377              g:=Gcdex(b,A[i,c2]);
378              b:=g.gcd;
379              a:=g.coeff1*a+g.coeff2*A[i,c1];
380            fi;
381          od;
382          N:=0;
383          for i in [r..n] do
384            if N<>1 then N:=GcdInt(N,A[i,c1]-QuoInt(A[i,c2],b)*a);fi;
385          od;
386        else
387          c := MATINTmgcdex(N,A[r,j],A{[r+1..n]}[j]);
388          b := A[r,j]+c*A{[r+1..n]}[j];
389          a := A[r,c1]+c*A{[r+1..n]}[c1];
390        fi;
391        t := MATINTmgcdex(N,a,[b])[1];
392        tmp := A[r,c1]+t*A[r,j];
393        while tmp=0 or tmp*A[k,c2]=(A[k,c1]+t*A[k,j])*A[r,c2] do
394          t := t+1+MATINTmgcdex(N,a+t*b+b,[b])[1];
395          tmp := A[r,c1]+t*A[r,j];
396        od;
397        if t>0 then
398          for i in [1..n] do A[i,c1] := A[i,c1]+t*A[i,j]; od;
399          if opt[4] then B[j,c1] := B[j,c1]+t; fi;
400        fi;
401      od;
402      if A[r,c1]*A[k,c1+1]=A[k,c1]*A[r,c1+1] then
403        for i in [1..n] do A[i,c1+1] := A[i,c1+1]+A[i,c2]; od;
404        if opt[4] then B[c2,c1+1] := 1; fi;
405      fi;
406      c2 := c1+1;
407    fi;
408
409    c := MATINTmgcdex(AbsInt(A[r,c1]),A[r+1,c1],A{[r+2..n]}[c1]);
410    for i in [r+2..n] do
411      if c[i-r-1]<>0 then
412        AddRowVector(A[r+1],A[i],c[i-r-1]);
413        if opt[3] then
414          C[r+1,i] := c[i-r-1];
415          AddRowVector(Q[r+1],Q[i],c[i-r-1]);
416        fi;
417      fi;
418    od;
419
420    i := r+1;
421    while A[r,c1]*A[i,c2]=A[i,c1]*A[r,c2] do i := i+1; od;
422    if i>r+1 then
423      c := MATINTmgcdex(AbsInt(A[r,c1]),A[r+1,c1]+A[i,c1],[A[i,c1]])[1]+1;;
424      AddRowVector(A[r+1],A[i],c);
425      if opt[3] then
426        C[r+1,i] := C[r+1,i]+c;
427        AddRowVector(Q[r+1],Q[i],c);
428      fi;
429    fi;
430
431    g := MATINTbezout(A[r,c1],A[r,c2],A[r+1,c1],A[r+1,c2]);
432    sig:=sig*SignInt(A[r,c1]*A[r+1,c2]-A[r,c2]*A[r+1,c1]);
433    A{[r,r+1]} := [[g.coeff1,g.coeff2],[g.coeff3,g.coeff4]]*A{[r,r+1]};
434    if opt[3] then
435      Q{[r,r+1]} := [[g.coeff1,g.coeff2],[g.coeff3,g.coeff4]]*Q{[r,r+1]};
436    fi;
437
438    for i in [r+2..n] do
439      q := QuoInt(A[i,c1],A[r,c1]);
440      AddRowVector(A[i],A[r],-q);
441      if opt[3] then AddRowVector(Q[i],Q[r],-q); fi;
442      q := QuoInt(A[i,c2],A[r+1,c2]);
443      AddRowVector(A[i],A[r+1],-q);
444      if opt[3] then AddRowVector(Q[i],Q[r+1],-q); fi;
445    od;
446
447  od;
448  rp[r+1] := m;
449  Info(InfoMatInt,2,"DoNFIM - r,m,n=",r,m,n);
450  if n=m and r+1<n then sig:=0;fi;
451
452  #smith w/ NO transforms - farm the work out...
453  if opt[1] and not (opt[3] or opt[4]) then
454    #R:=rec(normal:=SNFofREF(A{[2..n-1]}{[2..m-1]}),rank:=r-1);
455    for i in [2..n-1] do
456      A[i-1]:=A[i]{[2..m-1]};
457    od;
458    Unbind(A[n-1]);
459    Unbind(A[n]);
460    R:=rec(normal:=SNFofREF(A,opt[5]),rank:=r-1);
461    if n=m then R. signdet:=sig;fi;
462    return R;
463  fi;
464
465  # hermite or (smith w/ column transforms)
466  if (not opt[1] and opt[2]) or (opt[1] and opt[4]) then
467    for i in [r, r-1 .. 1] do
468      Info(InfoMatInt,2,"DoNFIM - reducing row ",i);
469      for j in [i+1 .. r+1] do
470        q := QuoInt(A[i,rp[j]]-(A[i,rp[j]] mod A[j,rp[j]]),A[j,rp[j]]);
471        AddRowVector(A[i],A[j],-q);
472        if opt[3] then AddRowVector(Q[i],Q[j],-q); fi;
473      od;
474      if opt[1] and i<r then
475        for j in [i+1..m] do
476          q := QuoInt(A[i,j],A[i,i]);
477          for k in [1..i] do A[k,j] := A[k,j]-q*A[k,i]; od;
478          if opt[4] then P[i,j] := -q; fi;
479        od;
480      fi;
481    od;
482  fi;
483
484  #Smith w/ row but not col transforms
485  if opt[1] and opt[3] and not opt[4] then
486    for i in [1..r-1] do
487      t := A[i,i];
488      A[i] := List([1..m],x->0);
489      A[i,i] := t;
490    od;
491    for j in [r+1..m-1] do
492      A[r,r] := GcdInt(A[r,r],A[r,j]);
493      A[r,j] := 0;
494    od;
495  fi;
496
497  #smith w/ col transforms
498  if opt[1] and opt[4] and r<m-1 then
499    c := MATINTmgcdex(A[r,r],A[r,r+1],A[r]{[r+2..m-1]});
500    for j in [r+2..m-1] do
501      A[r,r+1] := A[r,r+1]+c[j-r-1]*A[r,j];
502      B[j,r+1] := c[j-r-1];
503      for i in [1..r] do P[i,r+1] := P[i,r+1]+c[j-r-1]*P[i,j]; od;
504    od;
505    P[r+1] := List([1..m],x->0);
506    P[r+1,r+1] := 1;
507    g := Gcdex(A[r,r],A[r,r+1]);
508    A[r,r] := g.gcd;
509    A[r,r+1] := 0;
510    for i in [1..r+1] do
511      t := P[i,r];
512      P[i,r] := P[i,r]*g.coeff1+P[i,r+1]*g.coeff2;
513      P[i,r+1] := t*g.coeff3+P[i,r+1]*g.coeff4;
514    od;
515    for j in [r+2..m-1] do
516      q := QuoInt(A[r,j],A[r,r]);
517      for i in [1..r+1] do P[i,j] := P[i,j]-q*P[i,r]; od;
518      A[r,j] := 0;
519    od;
520    for i in [r+2..m-1] do
521      P[i] := List([1..m],x->0);
522      P[i,i] := 1;
523    od;
524  fi;
525
526  #row transforms finisher
527  if opt[3] then for i in [r+2..n] do Q[i,i]:= 1; od; fi;
528
529  for i in [2..n-1] do
530    A[i-1]:=A[i]{[2..m-1]};
531  od;
532  Unbind(A[n-1]);
533  Unbind(A[n]);
534  R:=rec(normal:=A);
535
536  if opt[3] then
537    R.rowC:=C{[2..n-1]}{[2..n-1]};
538    R.rowQ:=Q{[2..n-1]}{[2..n-1]};
539  fi;
540
541  if opt[1] and opt[4] then
542    R.colC:=B{[2..m-1]}{[2..m-1]};
543    R.colQ:=P{[2..m-1]}{[2..m-1]};
544  fi;
545
546  R.rank:=r-1;
547  if n=m then R.signdet:=sig;fi;
548  return R;
549
550end);
551
552#############################################################################
553##
554#F  NormalFormIntMat(<mat>,<options>)
555##
556InstallGlobalFunction(NormalFormIntMat,
557function(mat,options)
558  local r,opt;
559  r:=DoNFIM(mat,options);
560  opt := BITLISTS_NFIM[options+1];
561  #opt := List(CoefficientsQadic(options,2),x->x=1);
562  #if Length(opt)<4 then
563  #  opt{[Length(opt)+1..4]} := List([Length(opt)+1..4],x->false);
564  #fi;
565
566  if opt[3] then
567    r.rowtrans:=r.rowQ*r.rowC;
568    #Unbind(r.rowQ);
569    #Unbind(r.rowC);
570  fi;
571
572  if opt[1] and opt[4] then
573    r.coltrans:=r.colC*r.colQ;
574    #Unbind(r.colQ);
575    #Unbind(r.colC);
576  fi;
577  return r;
578end);
579
580#############################################################################
581##
582#O  TriangulizedIntegerMat(<mat>);
583##
584InstallMethod(TriangulizedIntegerMat,"dispatch",true,[IsMatrix],0,
585function(mat)
586  return DoNFIM(mat,0).normal;
587end);
588
589InstallOtherMethod(TriangulizedIntegerMat,"empty matrix",true,[IsList],0,
590function(mat)
591  if NrRows(mat)<>1 or (not IsList(mat[1])) or NrCols(mat)<>0 then
592    TryNextMethod();
593  fi;
594  return DoNFIM(mat,0).normal;
595end);
596InstallOtherMethod(TriangulizedIntegerMat,"empty",true,[IsEmpty],0,Immutable);
597
598#############################################################################
599##
600#O  TriangulizedIntegerMatTransform(<mat>);
601##
602InstallMethod(TriangulizedIntegerMatTransform,"dispatch",true,[IsMatrix],0,
603function(mat)
604  return NormalFormIntMat(mat,4);
605end);
606
607InstallOtherMethod(TriangulizedIntegerMatTransform,"empty matrix",true,[IsList],0,
608function(mat)
609  if NrRows(mat)<>1 or (not IsList(mat[1])) or NrCols(mat)<>0 then
610    TryNextMethod();
611  fi;
612  return NormalFormIntMat(mat,4);
613end);
614InstallOtherMethod(TriangulizedIntegerMatTransform,"empty",true,[IsEmpty],0,
615function(mat)
616  return rec(normal:=Immutable(mat),rowtrans:=Immutable([[1]]));
617end);
618
619#############################################################################
620##
621#O  TriangulizeIntegerMat(<mat>);
622##
623InstallMethod(TriangulizeIntegerMat,"dispatch",true,[IsMatrix and IsMutable],0,
624function(mat)
625  DoNFIM(mat,16);
626end);
627
628InstallOtherMethod(TriangulizeIntegerMat,"empty",true,[IsEmpty],0,Immutable);
629
630#############################################################################
631##
632#O  HermiteNormalFormIntegerMat(<mat>);
633##
634InstallMethod(HermiteNormalFormIntegerMat,"dispatch",true,[IsMatrix],0,
635function(mat)
636  return DoNFIM(mat,2).normal;
637end);
638
639InstallOtherMethod(HermiteNormalFormIntegerMat,"empty matrix",true,[IsList],0,
640function(mat)
641  if NrRows(mat)<>1 or (not IsList(mat[1])) or NrCols(mat)<>0 then
642    TryNextMethod();
643  fi;
644  return DoNFIM(mat,2).normal;
645end);
646InstallOtherMethod(HermiteNormalFormIntegerMat,"empty",true,[IsEmpty],0,
647  Immutable);
648
649#############################################################################
650##
651#O  HermiteNormalFormIntegerMatTransform(<mat>);
652##
653InstallMethod(HermiteNormalFormIntegerMatTransform,"dispatch",true,[IsMatrix],0,
654function(mat)
655  return NormalFormIntMat(mat,6);
656end);
657
658InstallOtherMethod(HermiteNormalFormIntegerMatTransform,"empty matrix",
659  true,[IsList],0,
660function(mat)
661  if NrRows(mat)<>1 or (not IsList(mat[1])) or NrCols(mat)<>0 then
662    TryNextMethod();
663  fi;
664  return NormalFormIntMat(mat,6);
665end);
666InstallOtherMethod(HermiteNormalFormIntegerMatTransform,"empty",true,
667  [IsEmpty],0,
668function(mat)
669  return rec(normal:=Immutable(mat),rowtrans:=Immutable([[1]]));
670end);
671
672#############################################################################
673##
674#O  SmithNormalFormIntegerMat(<mat>);
675##
676InstallMethod(SmithNormalFormIntegerMat,"dispatch",true,[IsMatrix],0,
677function(mat)
678  return DoNFIM(mat,1).normal;
679end);
680
681InstallOtherMethod(SmithNormalFormIntegerMat,"empty matrix",true,[IsList],0,
682function(mat)
683  if NrRows(mat)<>1 or (not IsList(mat[1])) or NrCols(mat)<>0 then
684    TryNextMethod();
685  fi;
686  return DoNFIM(mat,1).normal;
687end);
688InstallOtherMethod(SmithNormalFormIntegerMat,"empty",true,[IsEmpty],0,
689  Immutable);
690
691#############################################################################
692##
693#O  SmithNormalFormIntegerMatTransforms(<mat>);
694##
695InstallMethod(SmithNormalFormIntegerMatTransforms,"dispatch",true,[IsMatrix],0,
696function(mat)
697  return NormalFormIntMat(mat,13);
698end);
699
700InstallOtherMethod(SmithNormalFormIntegerMatTransforms,"empty matrix",
701  true,[IsList],0,
702function(mat)
703  if NrRows(mat)<>1 or (not IsList(mat[1])) or NrCols(mat)<>0 then
704    TryNextMethod();
705  fi;
706  return NormalFormIntMat(mat,13);
707end);
708InstallOtherMethod(SmithNormalFormIntegerMatTransforms,"empty",true,
709  [IsEmpty],0,
710function(mat)
711  return
712  rec(normal:=Immutable(mat),rowtrans:=Immutable([[1]]),
713      coltrans:=Immutable([[1]]));
714end);
715
716InstallGlobalFunction( DiagonalizeIntMat, function ( mat )
717  DoNFIM(mat,17);
718end);
719
720#############################################################################
721##
722#M  DiagonalizeMat(<integers>,<mat>)
723##
724InstallMethod( DiagonalizeMat, "over the integers",
725  [ IsIntegers, IsMatrix and IsMutable ],
726function(I,mat)
727  DiagonalizeIntMat(mat);
728end );
729
730#############################################################################
731##
732#M  ElementaryDivisorsTransformationsMat(<integers>,<mat>)
733##
734InstallMethod( ElementaryDivisorsTransformationsMat, "over the integers",
735  [ IsIntegers, IsMatrix and IsMutable ],
736function(I,mat)
737  return SmithNormalFormIntegerMatTransforms(mat);
738end );
739
740#############################################################################
741##
742#M  BaseIntMat(<mat>)
743##
744InstallMethod(BaseIntMat,"use HNF",true,
745  [IsMatrix and IsCyclotomicCollColl],0,
746function( mat )
747local norm;
748  norm := NormalFormIntMat( mat, 2 );
749  return norm.normal{[1..norm.rank]};
750end);
751
752InstallOtherMethod(BaseIntMat,"empty",true,
753  [IsEmpty],0,Immutable);
754
755#############################################################################
756##
757#M  BaseIntersectionIntMats(<m1>,<m2>)
758##
759InstallMethod(BaseIntersectionIntMats,"use HNF",true,
760  [IsMatrix and IsCyclotomicCollColl,IsMatrix and IsCyclotomicCollColl],0,
761function( M1, M2 )
762local M, Q, r, T;
763  M := Concatenation( M1, M2 );
764  r := NormalFormIntMat( M, 4 );
765  T := r.rowtrans{[r.rank+1..Length(M)]}{[1..Length(M1)]};
766  if not IsEmpty( T ) then T := T * M1; fi;
767  return BaseIntMat( T );
768end);
769
770InstallOtherMethod(BaseIntersectionIntMats,"emptyL",true,
771  [IsEmpty,IsObject],0,
772function(L,R)
773  return Immutable(L);
774end);
775
776InstallOtherMethod(BaseIntersectionIntMats,"emptyR",true,
777  [IsObject,IsEmpty],0,
778function(L,R)
779  return Immutable(R);
780end);
781
782#############################################################################
783##
784#M  ComplementIntMat(<m1>,<m2>)
785##
786InstallMethod(ComplementIntMat,"use HNF and SNF",true,
787  [IsMatrix and IsCyclotomicCollColl,IsMatrix and IsCyclotomicCollColl],0,
788function( full,sub )
789local F, S, M, r, T, R;
790  F := BaseIntMat( full );
791  if IsEmpty( sub ) or IsZero( sub ) then
792    return rec( complement := F, sub := [], moduli := [] );
793  fi;
794  S := BaseIntersectionIntMats( F, sub );
795  if S <> BaseIntMat( sub ) then
796    Error( "sub must be submodule of full" );
797  fi;
798  # find T such that BaseIntMat(T*F) = S
799  M := Concatenation( F, S );
800  T := NormalFormIntMat( M, 4 ).rowtrans^-1;
801  T := T{[Length(F)+1..Length(T)]}{[1..Length(F)]};
802
803  # r.rowtrans * T * F = r.normal * r.coltrans^-1 * F
804  r := NormalFormIntMat( T, 13 );
805  M := r.coltrans^-1 * F;
806  R := rec( complement := BaseIntMat( M{[1+r.rank..Length(M)]} ),
807            sub := r.rowtrans * T * F,
808            moduli := List( [1..r.rank], i -> r.normal[i,i] ) );
809  return R;
810end);
811
812InstallOtherMethod(ComplementIntMat,"empty submodule",true,
813  [IsMatrix and IsCyclotomicCollColl,IsList and IsEmpty],0,
814function( full, sub )
815  return rec( complement := BaseIntMat( full ), sub := [], moduli := [] );
816end );
817
818#############################################################################
819##
820#M  NullspaceIntMat(<mat>)
821##
822InstallMethod(NullspaceIntMat,"use HNF",true,
823  [IsMatrix and IsCyclotomicCollColl],0,
824function( mat )
825local norm, kern;
826  norm := NormalFormIntMat( mat, 4 );
827  kern := norm.rowtrans{[norm.rank+1..Length(mat)]};
828  return BaseIntMat( kern );
829end);
830
831#############################################################################
832##
833#M  SolutionIntMat(<mat>,<vec>)
834##
835InstallMethod(SolutionIntMat,"use HNF",true,
836  [IsMatrix and IsCyclotomicCollColl,
837   IsList and IsCyclotomicCollection],0,
838function( mat,v )
839local norm, rs, t, M, r;
840  if IsZero(mat) then
841    if IsZero(v) then
842      return ListWithIdenticalEntries( NrRows(mat), 0 );
843    else
844      return fail;
845    fi;
846  fi;
847  norm := NormalFormIntMat( mat, 4 );
848  t := norm.rowtrans;
849  rs :=  norm.normal{[1..norm.rank]};
850  M := Concatenation( rs, [v] );
851  r := NormalFormIntMat( M, 4 );
852  if r.rank = Length(r.normal) or
853     r.rowtrans[Length(M),Length(M)] <> 1 then
854    return fail;
855  fi;
856  return -r.rowtrans[Length(M)]{[1..r.rank]} * t{[1..r.rank]};
857end);
858
859InstallOtherMethod(SolutionIntMat,"empty",true,[IsEmpty,IsObject],0,
860function(mat,v)
861  return fail;
862end);
863
864#############################################################################
865##
866#M  SolutionNullspaceIntMat(<mat>,<vec>)
867##
868InstallMethod(SolutionNullspaceIntMat,"use HNF",true,
869  [IsMatrix and IsCyclotomicCollColl,
870   IsList and IsCyclotomicCollection],0,
871function( mat,v )
872local norm, rs, t, M, r, kern, len;
873  if IsZero(mat) then
874    len := Length(mat);
875    if IsZero(v) then
876      return [ListWithIdenticalEntries(len,0), IdentityMat(len)];
877    else
878      return [fail, IdentityMat(len)];
879    fi;
880  fi;
881  norm := NormalFormIntMat( mat, 4 );
882  kern := norm.rowtrans{[norm.rank+1..Length(mat)]};
883  kern := BaseIntMat( kern );
884  t := norm.rowtrans;
885  rs :=  norm.normal{[1..norm.rank]};
886  M := Concatenation( rs, [v] );
887  r := NormalFormIntMat( M, 4 );
888  if r.rank = Length(r.normal) or
889     r.rowtrans[Length(M),Length(M)] <> 1 then
890    return [fail,kern];
891  fi;
892  return [-r.rowtrans[Length(M)]{[1..r.rank]} * t{[1..r.rank]}, kern];
893end);
894
895
896#############################################################################
897##
898#F  DeterminantIntMat(<mat>)
899##
900InstallGlobalFunction(DeterminantIntMat,function(mat)
901local sig, n, m, A, r, c2, c1, j, k, c, i, g, q;
902
903  sig:=1;
904
905  #Embed mat in 2 larger "id" matrix
906  n := Length(mat)+2;
907  # Crossover point roughly 20x20 matrices, so farm the work if smaller..
908  if n<22 then return DeterminantMat(mat);fi;
909  m := NrCols(mat)+2;
910
911  if not n=m then
912    Error( "DeterminantIntMat: <mat> must be a square matrix" );
913  fi;
914
915  A := [List([1..m],x->0)];
916  for i in [2..n-1] do
917    A[i] := [0];
918    Append(A[i],mat[i-1]);
919    A[i,m] := 0;
920  od;
921  A[n] := List([1..m],x->0);
922  A[1,1] := 1;      A[n,m] := 1;
923
924  r := 0;    c2 := 1;
925  while m>c2 do
926    Info(InfoMatInt,2,"DeterminantIntMat - reached column ",c2," of ",m);
927    r := r+1;
928    c1 := c2;
929
930    j := c1+1;
931    while j<=m do
932      k := r+1;
933      while k<=n and A[r,c1]*A[k,j]=A[k,c1]*A[r,j] do k := k+1; od;
934      if k<=n then c2 := j; j := m; fi;
935      j := j+1;
936    od;
937
938    c := MATINTmgcdex(AbsInt(A[r,c1]),A[r+1,c1],A{[r+2..n]}[c1]);
939    for i in [r+2..n] do
940      if c[i-r-1]<>0 then
941        AddRowVector(A[r+1],A[i],c[i-r-1]);
942      fi;
943    od;
944
945    i := r+1;
946    while A[r,c1]*A[i,c2]=A[i,c1]*A[r,c2] do
947      i := i+1;
948    od;
949
950    if i>r+1 then
951      c := MATINTmgcdex(AbsInt(A[r,c1]),A[r+1,c1]+A[i,c1],[A[i,c1]])[1]+1;;
952      AddRowVector(A[r+1],A[i],c);
953    fi;
954
955    g := MATINTbezout(A[r,c1],A[r,c2],A[r+1,c1],A[r+1,c2]);
956    sig:=sig*SignInt(A[r,c1]*A[r+1,c2]-A[r,c2]*A[r+1,c1]);
957    if sig=0 then return 0;fi;
958    A{[r,r+1]} := [[g.coeff1,g.coeff2],[g.coeff3,g.coeff4]]*A{[r,r+1]};
959
960    for i in [r+2..n] do
961      q := QuoInt(A[i,c1],A[r,c1]);
962      AddRowVector(A[i],A[r],-q);
963      q := QuoInt(A[i,c2],A[r+1,c2]);
964      AddRowVector(A[i],A[r+1],-q);
965    od;
966  od;
967
968  for i in [2..r+1] do
969    sig:=sig*A[i,i];
970  od;
971
972  return sig;
973
974end);
975
976#############################################################################
977##
978#M  AbelianInvariantsOfList( <list> ) . . . . .  abelian invariants of a list
979##
980InstallMethod( AbelianInvariantsOfList,
981    [ IsCyclotomicCollection ],
982function ( list )
983    local   invs, elm;
984
985    invs := [];
986    for elm  in list  do
987        if elm = 0  then
988            Add( invs, 0 );
989        elif 1 < elm  then
990            Append( invs, List( Collected(Factors(elm)), x->x[1]^x[2] ) );
991        elif elm < -1 then
992            Append( invs, List( Collected(Factors(-elm)), x->x[1]^x[2] ) );
993        fi;
994    od;
995    Sort(invs);
996    return invs;
997end );
998
999InstallOtherMethod( AbelianInvariantsOfList,
1000    [ IsList and IsEmpty ],
1001    list -> [] );
1002
1003
1004# Reduce a list of abelianized relations: Heuristic reduction without
1005# making big vectors, iterate three times. Does not aim to do full HNF
1006InstallGlobalFunction(ReducedRelationMat,function(mat)
1007local n,zero,nv,new,pip,piv,i,v,p,w,g,nov,pin,now,rat,extra,clean,assign,try;
1008
1009  nv:=v->v*SignInt(v[PositionNonZero(v)]);
1010  assign:=function(p,v)
1011  local a,i,w,wn;
1012    a:=v[p];
1013    for i in [1..Length(pip)] do
1014      if i<>p and IsInt(pip[i]) and mat[pip[i]][p]<>0 then
1015        w:=mat[pip[i]]-QuoInt(mat[pip[i]][p],a)*v;
1016        wn:=w*w;
1017        if wn<=rat*pin[i] then
1018          mat[pip[i]]:=nv(w);
1019          pin[i]:=wn;
1020        fi;
1021      fi;
1022    od;
1023    mat[pip[p]]:=v;
1024    # also try to reduce extra vectors
1025    for i in [1..Length(extra)] do
1026      w:=extra[i];
1027      if not IsZero(extra[i]) then
1028        wn:=w*w;
1029        w:=w-QuoInt(w[p],a)*v;
1030        if w*w<=rat*wn then
1031          extra[i]:=w;
1032        fi;
1033      fi;
1034    od;
1035  end;
1036
1037  n:=NrCols(mat);
1038  rat:=2; # growth ratio
1039  zero:=ListWithIdenticalEntries(n,0);
1040  mat:=Filtered(mat,x->not IsZero(x));
1041  new:=Set(mat,nv); # kill duplicates
1042  Info(InfoMatInt,1,"Reduce ",Length(mat)," to ",Length(new));
1043  pip:=ListWithIdenticalEntries(n,fail);
1044  piv:=[];
1045  pin:=[];
1046  mat:=[];
1047  extra:=[];
1048
1049  # we once reduce and then go over the remainders again in case they were
1050  # nice and short
1051  for try in [1..3] do
1052    SortBy(new, x -> - x*x); # reversed norm sort
1053    i:=Length(new);
1054    while i>0 do
1055      v:=ShallowCopy(new[i]);
1056      Info(InfoMatInt,3,"Process ",i);#" Norm:",v*v,"\n");
1057      Unbind(new[i]); # take off stack
1058      i:=i-1;
1059      clean:=true;
1060      p:=PositionNonZero(v);
1061      while p<=n and pip[p]<>fail do
1062        if v[p] mod piv[p]=0 then
1063          # divides, reduce
1064          #v:=v-QuoInt(v[p],piv[p])*mat[pip[p]];
1065          AddRowVector(v,mat[pip[p]],-QuoInt(v[p],piv[p]));
1066          p:=PositionNonZero(v,p);
1067        elif clean and piv[p] mod v[p]=0 then
1068          # swap and clean out
1069          v:=nv(v);
1070          Info(InfoMatInt,2,"Replace pivot ",piv[p],"@",p," to ",v[p]);
1071          w:=mat[pip[p]];
1072          #mat[pip[p]]:=v;
1073          assign(p,v);
1074          pin[p]:=v*v;
1075          piv[p]:=v[p];
1076          v:=w;
1077          #v:=v-QuoInt(v[p],piv[p])*mat[pip[p]];
1078          AddRowVector(v,mat[pip[p]],-QuoInt(v[p],piv[p]));
1079          p:=PositionNonZero(v,p);
1080        else
1081          g:=Gcdex(v[p],piv[p]);
1082          # form new pivot with gcd
1083          #w:=g.coeff2*mat[pip[p]]+g.coeff1*v; # automatically normed by Gcdex
1084          w:=g.coeff2*mat[pip[p]];
1085          AddRowVector(w,v,g.coeff1); # automatically normed by Gcdex
1086          now:=w*w;
1087          if (not clean) or now>rat*pin[p] then
1088            # only reduce a bit, not full gcd
1089            #v:=v-QuoInt(v[p],piv[p])*mat[pip[p]];
1090            AddRowVector(v,mat[pip[p]],-QuoInt(v[p],piv[p]));
1091            p:=PositionNonZero(v,p);
1092            clean:=false;
1093          else
1094            # replace with cgd pivot
1095            Info(InfoMatInt,2,"Reduce pivot ",piv[p],"@",p," to ",g.gcd);
1096            new[i+1]:=v; # keep old vectors to process
1097            new[i+2]:=mat[pip[p]];
1098            i:=i+2;
1099            #mat[pip[p]]:=w;
1100            assign(p,w);
1101            piv[p]:=w[p];
1102            pin[p]:=now;
1103            p:=fail; # to bail out while loop
1104          fi;
1105
1106        fi;
1107      od;
1108      if not clean then
1109        # only reduced, did not do gcd
1110        Add(extra,v);
1111      elif p<=n then
1112        # new pivot position
1113        v:=nv(v); # norm (so we can compare with <)
1114        pip[p]:=Length(mat)+1;
1115        #Add(mat,v);
1116        assign(p,v);
1117        Info(InfoMatInt,1,"Added @",Length(mat));
1118        piv[p]:=v[p];
1119        pin[p]:=v*v;
1120      fi;
1121    od;
1122    # now we've processed all. Clean out extra
1123    new:=List(Filtered(Set(extra),x->not IsZero(x)),nv);
1124    Info(InfoMatInt,1,"After ",try,": ",Length(extra)," to ",Length(new),
1125      " new ones");
1126    extra:=[];
1127
1128  od;
1129
1130  mat:=Filtered(Concatenation(mat,new),x->not IsZero(x));
1131
1132  # need to keep one line.
1133  if Length(mat)=0 then mat:=[zero];fi;
1134
1135  return mat;
1136
1137end);
1138