From 9b8c3722b3e12846139cfd7195467dba20d39647 Mon Sep 17 00:00:00 2001
From: kalmar <kalmar@amu.edu.pl>
Date: Mon, 13 Mar 2017 11:11:09 +0100
Subject: [PATCH] Allow generation of SL(n,p) (matrices over Z/p)

---
 SL3Z.jl | 77 +++++++++++++++++++++++++++++++++++++++++++++++++--------
 1 file changed, 67 insertions(+), 10 deletions(-)

diff --git a/SL3Z.jl b/SL3Z.jl
index 67f552b..bfc37f6 100644
--- a/SL3Z.jl
+++ b/SL3Z.jl
@@ -1,25 +1,82 @@
 using JLD
 using JuMP
+import Primes: isprime
 import SCS: SCSSolver
 import Mosek: MosekSolver
 
+using Mods
+
 using Groups
 using ProgressMeter
 
 
-function SL₃ℤ_generatingset()
+function SL_generatingset(n::Int)
 
-    function E(i::Int, j::Int, N::Int=3)
-        @assert i≠j
-        k = eye(N)
-        k[i,j] = 1
-        return k
+    indexing = [(i,j) for i in 1:n for j in 1:n if i≠j]
+
+    S = [E(i,j,N=n) for (i,j) in indexing];
+    S = vcat(S, [convert(Array{Int,2},x') for x in S]);
+    S = vcat(S, [convert(Array{Int,2},inv(x)) for x in S]);
+    return unique(S)
+end
+
+function E(i::Int, j::Int; val=1, N::Int=3, mod=Inf)
+   @assert i≠j
+   m = eye(Int, N)
+   m[i,j] = val
+   if mod == Inf
+       return m
+   else
+       return [Mod(x,mod) for x in m]
+   end
+end
+
+function cofactor(i,j,M)
+    z1 = ones(Bool,size(M,1))
+    z1[i] = false
+
+    z2 = ones(Bool,size(M,2))
+    z2[j] = false
+
+    return M[z1,z2]
+end
+
+import Base.LinAlg.det
+
+function det(M::Array{Mod,2})
+    if size(M,1) ≠ size(M,2)
+        d = Mod(0,M[1,1].mod)
+    elseif size(M,1) == 2
+        d =  M[1,1]*M[2,2] - M[1,2]*M[2,1]
+    else
+        d = zero(eltype(M))
+        for i in 1:size(M,1)
+            d += (-1)^(i+1)*M[i,1]*det(cofactor(i,1,M))
+        end
     end
+#     @show (M, d)
+    return d
+end
 
-    S = [E(1,2), E(1,3), E(2,3)];
-    S = vcat(S, [x' for x in S]);
-    S = vcat(S, [inv(x) for x in S]);
-    return S
+function adjugate(M)
+    K = similar(M)
+    for i in 1:size(M,1), j in 1:size(M,2)
+        K[j,i] = (-1)^(i+j)*det(cofactor(i,j,M))
+    end
+    return K
+end
+
+import Base: inv, one, zero, *
+
+one(::Type{Mod}) = 1
+zero(::Type{Mod}) = 0
+zero(x::Mod) = Mod(x.mod)
+
+function inv(M::Array{Mod,2})
+    d = det(M)
+    d ≠ 0*d || thow(ArgumentError("Matrix is not invertible!"))
+    return inv(det(M))*adjugate(M)
+    return adjugate(M)
 end
 
 function prepare_Δ_sdp_constraints(identity, S)