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