newPackage(
     "Kronecker",
     Version => "0.2.1",
     Date => "March 10, 2008",
     Authors => {{Name => "Edward Carter",
               Email => "ecarter@math.berkeley.edu",
               HomePage => "http://math.berkeley.edu/~ecarter/"}},
     DebuggingMode => true
     )
    
export {kroneckerNormalForm, kroneckerIndices, rationalNormalForm}

-- columns(A)
--
-- Returns a list consisting of the columns of the matrix A.  The zero matrix
-- gives an empty list.
columns = method()
columns(Nothing) := N -> {}
columns(Matrix) := A -> (
     if A == 0 then {}
     else apply(numgens source A, i -> submatrix(A, {i}))
     )

-- joinColumns(M,L)
--
-- The opposite of the method columns(A).  The elements of the list L are
-- assumed to be columns which are joined to form a matrix.
--
-- The parameter M is only used in the case of L being an empty list.  In this
-- case, we return the zero map from M to itself.
joinColumns = method()
joinColumns(Module, List) := (M,L) -> (
     if #L == 0 then map(M, M, 0)
     else fold((v,w) -> v | w, L)
     )

-- gen(I)
--
-- Assert that the ideal I is principal and return its generator.
gen = method()
gen(Ideal) := I -> (
     if I == 0 then 0 else (
     	  assert(numgens I == 1);
     	  (gens I)_0_0
	  )
     )

-- monomialMatrix(M, n, x, y)
--
-- Returns the homogeneous, degree 0 map to the module M determined by the
-- matrix:
--
--     | y^n -xy^(n-1) . . . (-x)^n |
--
monomialMatrix = method()
monomialMatrix(Module, ZZ, RingElement, RingElement) := (M, deg, x, y) -> (
     map(M, deg+1, (i,j) -> y^(deg-j)*(-x)^j)
     )

-- getDegree(v)
--
-- Return the degree of the entries of a homogeneous matrix v.
getDegree = method()
getDegree(Matrix) := v -> (
     degs := degrees v;
     degs#1#0#0 - degs#0#0#0
     )

-- expandPolynomialRelation(v, x, y)
--
-- Assuming v is a homogeneous column vector of degree n, returns a matrix
-- with columns w0...wn where
-- 
--     v = y^n*w0 - xy^(n-1)*w1 + . . . + (-x)^n*wn
--
expandPolynomialRelation = method()
expandPolynomialRelation(Matrix, RingElement, RingElement) := (v, x, y) -> (
     deg := getDegree v;
     dual(dual v // monomialMatrix(target dual v, deg, x, y))
     )
	
-- complementaryForm(f,x,y)
--
-- Here, f is some linear form, and we return a complementary linear form.  The parameters x
-- and y should be set to the two variables of the ring.
complementaryForm = method()
complementaryForm(RingElement, RingElement, RingElement) := (f,x,y) -> (
     if f % ideal(y) == 0 then x
     else y
     )

-- findGen(M,N,P)
--
-- Given modules N and P which are both submodules of M, return
-- some element of N which is not in P, in the form of a column
-- vector.  If no such element exists, we return 0.
findGen = method()
findGen(Module,Module,Module) := (M,N,P) -> (
     f := inducedMap(M,P);
     G := select(columns mingens N, w -> (
	       w2 := map(M, , w);
	       v := w2 // f;
	       w2-f*v != 0
	       )
	  );
     if #G > 0 then G#0 else map(M, 1, (i,j) -> 0)
     )

-- kroneckerBases(P,x,y)
--
-- Returns automorphisms {A1,A2} on the target and source of P, respectively,
-- such that, relative to the ordered bases which are the columns of A1 and A2,
-- P is in Kronecker normal form.
kroneckerBases = method()
kroneckerBases(Matrix,RingElement,RingElement) := (P,x,y) -> (
     -- The following variables will be referred to frequently in the comments.
     -- Our pencil is equal to Ax+By.
     R := ring P;
     M := source P;
     N := target P;
     A := map(N, , P ** R^1/(x-1,y));
     B := map(N, , P ** R^1/(x,y-1));
     
     -- Each homogeneous generator of the kernel corresponds to a block of the form
     --
     --     | x y 0 0 . . . 0 0 |
     --     | 0 x y 0 . . . 0 0 |
     --     | 0 0 x y . . . 0 0 |
     --     | . . . . . . . . . |
     --     | 0 0 0 0 . . . x y |
     --
     -- in the final Kronecker normal form.  The degree, d, of the element of the 
     -- kernel is equal to the number of rows.  We use expandPolynomialRelation 
     -- to produce d suitable basis elements in M.  Multiplying all of these by A 
     -- gives d basis elements in N, and multiplying B times the last one gives
     -- a (d+1)-st basis element in N.  Then P relative to these vectors is the
     -- above block.  Linear independence follows from the fact that "ker P"
     -- gives us *minimal* homogeneous generators of the kernel.
     --
     -- If d=0, using that kernel vector as a basis element in M produces a
     -- zero column, as desired.
     --
     -- The variables *SourceList are lists of basis elements in M, and the
     -- variables *TargetList are lists of basis elements in N.
     kerSourceList := apply(columns sort(gens ker P, DegreeOrder => Ascending), v -> (
	       map(M, , expandPolynomialRelation(v, x, y))
	       ));
     kerTargetList := apply(select(kerSourceList, v -> numgens source v > 1), v -> (
	       map(N, , A*submatrix(v, {0..(numgens source v - 2)}))
	       ));
     
     -- The second main type of block in the Kronecker normal form corresponds to 
     -- elements in the kernel of dual P.  (Or, after sheafification, to the locally
     -- free part of the cokernel of P)
     L := columns sort(gens ker dual P, DegreeOrder => Ascending);
     L2 := select(L, v -> getDegree(v) == 0);
     cokerTargetList := if #L2 == 0 then {} else (
	  X := dual joinColumns(N, L2);
	  {map(N, , id_(target X) // X)}
	  );
     cokerSourceList := {};
     -- Each homogeneous generator of the kernel of dual P corresponds to a
     -- block of the form
     --
     --     | x 0 0 . . . 0 |
     --     | y x 0 . . . 0 |
     --     | 0 y x . . . 0 |
     --     | 0 0 y . . . 0 |
     --     | . . . . . . . |
     --     | 0 0 0 . . . x |
     --     | 0 0 0 . . . y |
     --
     -- where the degree, d, is equal to the number of columns.  An element of
     -- this kernel is a row matrix w' such that w'P=0.  By using 
     -- expandPolynomialRelation on dual w', we get a matrix w with d+1 rows
     -- such that 
     --
     --     | (-y)^d  (-y)^(d-1)*x  (-y)^(d-2)*x^2  . . .  x^d | * w * P = 0
     --
     -- Therefore the columns of w*P are all linear combinations of columns
     -- of the above block.  Since "ker dual P" gives minimal generators of
     -- the kernel, we also know that all these columns are in the span of
     -- of the columns of w*P.  Therefore we can pull the above block back
     -- through w*P to get suitable basis vectors in M, then multiply by A
     -- and B to get suitable basis vectors in N.
     L2 = select(L, v -> getDegree(v) > 0);
     if #L2 > 0 then (
	  L3 := apply(L2, v -> (getDegree(v), 
		    map(R^(getDegree(v)+1), N, dual expandPolynomialRelation(v, x, y))));
	  wList := apply(L2, v -> map(R^(getDegree(v)+1), N, dual expandPolynomialRelation(v, x, y)));
	  tvList := apply(#wList, i -> (
		    curr := if i == 0 then first columns gens target wList#0 
		    else map(target wList#0, 1, (i,j) -> 0);
		    for j from 1 to #wList-1 do (
			 curr = curr || if i == j then first columns gens target wList#i 
			 else map(target wList#j, 1, (i,j) -> 0);
			 );
		    (-1 + numgens target wList#i, curr)
		    )
	       );
	  w := fold(wList, (v1,v2) -> v1 || v2);
	  scan(tvList, (deg, tmpv) -> (
		    tv := A*(map(target w, , tmpv) // (w*A));
		    cokerTargetList = cokerTargetList | {tv};
	       	    for i from 1 to deg do (
		    	 sv := tv // A;
		    	 tv = B*sv;
		    	 cokerTargetList = cokerTargetList | {tv};
		    	 cokerSourceList = cokerSourceList | {sv};
		    	 )
		    )
	       );
	  );
     
     -- What's left at this point is the regular part of the pencil (square blocks
     -- in the normal form with nonzero determinant, analogous to Jordan or rational
     -- normal form).  For this, we will have to generate basis elements in N
     -- which will span Q.
     Q := coker joinColumns(N, {P} | kerTargetList | cokerTargetList);
     if #(kerTargetList | cokerTargetList) > 0 then (
	  -- We're supposed to be building bases for M and N, so the list
	  -- of vectors in N at each point should be linearly independent.
	  assert(ker joinColumns(N, kerTargetList | cokerTargetList) == 0);
	  );
     regularTargetList := {};
     regularSourceList := {};
     if Q != 0 then (
	  -- The cokernel of the regular part of the pencil is annihilated
	  -- by the determinant.  Each irreducible factor of the determinant
	  -- will give rise to blocks analogous to the blocks of the rational
	  -- normal form.
	  Q2 := coker P;
	  scan(select(factor gen ann Q, t -> t#0 % ideal(x,y) == 0), t -> (
		    -- f is an irreducible polynomial, n = deg f, and e is
		    -- the degree to which f occurs in the minimal polynomial
		    -- of the matrix (or the size of the biggest Jordan block
		    -- with that eigenvalue)
		    f := t#0;
		    e := t#1;
		    n := (degree f)#0;

     	       	    -- We will produce basis vectors corresponding to
		    -- a block of the form diagForm*I+compForm*J, where
		    -- J is a matrix in rational normal form.
		    diagForm := if n == 1 then f else x;
		    compForm := complementaryForm(diagForm, x, y);
		    A = map(N, , P ** R^1/(diagForm-1, compForm));
		    B = map(N, , P ** R^1/(diagForm, compForm-1));
		    
		    rels := joinColumns(Q2, kerTargetList | cokerTargetList);
		    for i from 0 to e-1 do (
			 -- Find a block in the rational normal form of
			 -- size n*(e-i).
			 N2 := Q2/ker(f^(e-i-1)*id_Q2);
			 u := map(N2, , findGen(N2, ker(f*id_N2), image map(N2, , rels)));
			 while u != 0 do (
			      -- First we need a pair (v,w) of vectors in
			      -- M and N, respectively, such that:
			      --
			      --   1. w = A*v
			      --   2. w in Q is killed by f^(e-i) but
			      --      not by f^(e-i-1)
			      --
			      -- If both n and e-i are 1, we want to
			      -- produce a 1x1 block.  That means v should
			      -- be in the kernel of B, since compForm should
			      -- not appear in that column in the normal form.
			      v := if n*(e-i) == 1 then (
				   tmpv := u // (inducedMap(N2,N)*A*inducedMap(source A, ker B));
				   inducedMap(source A, ker B)*tmpv
				   )
			      else (
				   u // (inducedMap(N2, N)*A)
				   );
			      w := A*v;
			      
			      if n == 1 then (
				   -- If e-i>1, this block has the form
				   -- 
				   --     | x y . . 0 |
				   --     | 0 x . . 0 |
				   --     | . . . . . |
				   --     | 0 0 . . y |
				   --     | 0 0 . . x |
				   --
				   -- where x=diagForm, y=compForm for brevity.
				   -- Our original (v,w) pair gives us the x in the
				   -- lower right corner.  Then each new basis
				   -- vector w in N can be obtained by B*v from the
				   -- previous basis vector v in M.  Then the new
				   -- v can be obtained by pulling back w through A.
				   -- We just need to make sure that the last vector
				   -- v in M is in the kernel of B.
				   v2 := v;
				   for j from 2 to e-i do (
					w2 := map(N, , B*v2);
					v2 = if j == e-i then (
					     map(M, , (w2 || map(N, 1, (i,j) -> 0)) // (A || B))
					     )
					else (
					     map(M, , w2 // A)
					     );
					w = w2 | w;
					v = v2 | v;
					);
				   )
			      else (
				   -- This block has a form like
 				   --
				   --     | x-a1*y y 0 . . 0 |
				   --     |  -a2*y x y . . 0 |
				   --     |  -a3*y 0 x . . 0 |
				   --     |    .   . . . . . |
				   --     |  -an*y 0 0 . . y |
				   --     | -an'*y 0 0 . . x |
				   --
				   -- where n'=n+1.  We proceed as above, except for
				   -- the leftmost column, where we don't want the
				   -- last v to be in the kernel of B.
				   atmp := f // map(R^1, n+1, (i,j) -> (-compForm)^j*diagForm^(n-j));
				   a := map(R^n, 1, (i,j) -> atmp_j_(i+1) // atmp_j_0);
				   
				   w2 := {w};
				   v3 := {v};
				   for k from 1 to e-i do (
					for j from 1 to n do (
					     if (1 < j and j < n) or (j == n and k < e-i) then (
						  w2 = {map(N, , B*first v3)} | w2;
						  v3 = {map(M, , (first w2) // A)} | v3;
						  )	
					     else (
						  if 1 == j and k > 1 then (
						       a2 := joinColumns(N, select(n, w2, b -> true))*a;
						       w2 = {map(N, , a2+B*first v3)} | w2;
						       v3 = {map(M, , (first w2) // A)} | v3;
						       )
						  else if j == n and k == e-i then (
						       w2 = {map(N, , B*first v3)} | w2;
						       a3 := joinColumns(N, select(n, w2, b -> true))*a;
						       v3 = {map(M, , ((first w2) || (-a3)) // (A || B))} | v3;
						       );
						  );
					     );
					);
				   w = joinColumns(N, w2);
				   v = joinColumns(M, v3);
				   );
			      regularTargetList = {w} | regularTargetList;
			      regularSourceList = {v} | regularSourceList;
			      assert(numgens N   
			      	   >= #columns joinColumns(N, kerTargetList | cokerTargetList | regularTargetList));
			      
			      rels = rels | w;
			      u = map(N2, , findGen(N2, ker(f*id_N2), image map(N2, , rels)));
     			      );
			 );
		    ));
	  );
     
     -- Return bases as automorphisms on N and M.
     (joinColumns(N, kerTargetList | cokerTargetList | regularTargetList), 
	  joinColumns(M, kerSourceList | cokerSourceList | regularSourceList))
     )

returnNormalForm = method()
returnNormalForm(Sequence, List) := (m,cm) -> (
     (A,P,Q) := (m#0, m#1, m#2);
     (tchg,schg) := (cm#0, cm#1);
     if tchg then (
	  if schg then (A, P, Q) else (A, P)
	  )
     else (
	  if schg then (A, Q) else A
	  )
     )

checkPencil = method()
checkPencil(Module) := M -> (
     if isFreeModule M and isHomogeneous M then true else (
	  print "expected a map between homogeneous free modules";
	  false
	  )
     )
checkPencil(Matrix,RingElement,RingElement) := (P,x,y) -> (
     if P != 0 then (
     	  R := ring P;
     	  M := source P;
     	  N := target P;
     	  
     	  if not checkPencil(M) then return false;
     	  if not checkPencil(N) then return false;
	  if not isHomogeneous P then (
	       print "expected a homogeneous map";
	       return false;
	       );
     	  if 1 + (degrees N)#0#0 != (degree P)#0 + (degrees M)#0#0 then (
	       print "expected a matrix of linear forms";
	       return false;
	       );
	  );
     true
     )

shortBlock = method()
shortBlock(ZZ,RingElement,RingElement) := (n,x,y) -> (
     R := ring x;
     if n == 0 then map(R^0, R^1, 0)
     else map(R^n, n+1, (i,j) -> if i == j then x else if i+1 == j then y else 0)
     )

tallBlock = method()
tallBlock(ZZ,RingElement,RingElement) := (n,x,y) -> (
     R := ring x;
     if n == 0 then map(R^1, R^0, 0)
     else map(R^(n+1), n, (i,j) -> if i == j then x else if i == j+1 then y else 0)
     )

regularBlock = method()
regularBlock(Power,RingElement,RingElement) := (t,x,y) -> (
     f := t#0;
     e := t#1;
     n := (degree f)#0;
     
     diagForm := if n == 1 then f else x;
     compForm := complementaryForm(f,x,y);
     
     R := ring x;
     if n == 1 then map(R^e, e, (i,j) -> (
	       if i == j then diagForm
	       else if i+1 == j then compForm
	       else 0))
     else (
	  m := e*n;
	  atmp := f // map(R^1, n+1, (i,j) -> (-compForm)^j*diagForm^(n-j));
	  a := map(R^n, 1, (i,j) -> atmp_j_(i+1) // atmp_j_0);
	  map(R^m, m, (i,j) -> (
		    if i+1 == j then compForm
		    else if j % n == 0 then (
			 if j <= i and i < j+n then (
			      if j == i then diagForm - compForm*a_0_0
			      else compForm*(-a_0_(i-j))
			      )
			 else 0
			 )
		    else if i == j then diagForm else 0
		    ))
	  )
     )

normalFormFromIndices = method()
normalFormFromIndices(RingElement,RingElement,List) := (x,y,L) -> (
     if ring x != ring y then (
	  print "expected elements of the same ring";
	  )
     else (
	  ind := L#0;
     	  (kerList, dualKerList, divisors) := ind;
	  R := ring x;
	  P := fold(apply(kerList, n -> shortBlock(n,x,y))
	       | apply(dualKerList, n -> tallBlock(n,x,y))
	       | apply(divisors, t -> regularBlock(t,x,y)), (v,w) -> v ++ w);
	  map(R^(numgens target P), , P)
     	  )
     )
	
kroneckerNormalForm = method(
     Options => {
	  ChangeMatrix => {true, true}		    -- target, source
	  }
     )
kroneckerNormalForm(Matrix) := o -> f -> (
     v := vars ring f;
     kroneckerNormalForm(f, v_0_0, v_1_0, ChangeMatrix => o.ChangeMatrix)
     )
kroneckerNormalForm(Matrix,RingElement,RingElement) := o -> (ff,x,y) -> (
     R := ring ff;
     d := (degrees target ff)#0#0;
     f := map(R^(numgens target ff) ** R^{-d}, R^(numgens source ff) ** R^{-d-1}, 
	  (i,j) -> ff_j_i);
     if checkPencil(f,x,y) then (
	  cm := o.ChangeMatrix;
	  if cm#0 or cm#1 then (
     	       (P,Q) := kroneckerBases(f,x,y);
     	       assert(det P != 0);
     	       assert(det Q != 0);
     	       tmat := map(target f, target f, P^-1);
     	       returnNormalForm((map(target f, , tmat*f*Q), tmat, Q), o.ChangeMatrix)
	       )
	  else map(target f, , normalFormFromIndices(x, y, {kroneckerIndices(f,x,y)}))
	  )
     )

factorList = method()
factorList(RingElement) := r -> (
     p := factor r;
     apply(#p, i -> p#i)
     )

elementaryDivisors = method()
elementaryDivisors(Matrix, Ideal, ZZ, RingElement) := (P, J, n, prevGen) -> (
     r := gen saturate(minors(n, P), J);
     if r == 0 then {} else (
     	  if r == 1 then elementaryDivisors(P, J, n+1, prevGen)
     	  else factorList(r // prevGen) | elementaryDivisors(P, J, n+1, r)
	  )
     )

kroneckerIndices = method()
kroneckerIndices(Matrix, RingElement, RingElement) := (P,x,y) -> (
     checkPencil(P,x,y);
     (apply(columns sort(gens ker P, DegreeOrder => Ascending), getDegree),
     	  apply(columns sort(gens ker dual P, DegreeOrder => Ascending), getDegree),
	  elementaryDivisors(P, ideal(x,y), 1, promote(1, ring P)))
     )

rationalNormalForm = method(
     Options => {
	  ChangeMatrix => {true, true}
	  }
     )
rationalNormalForm(Matrix) := o -> f -> (
     R := ring f;
     if target f == source f then (
     	  n := numgens target f;
     	  S := R[rv0,rv1];
     	  use S;
     	  N := S^n;
     	  M := N ** S^{-1};
     	  P := map(N, M, rv0*id_M + rv1*(f**S));
     	  proj := map(R, S, matrix{{0_R,1_R}});
     	  returnNormalForm(apply(kroneckerNormalForm(P, rv0, rv1), g ->
	       	    map(target f, source f, (i,j) -> proj(g_j_i))), o.ChangeMatrix)
	  )
     else print "expected a square matrix"
     )

beginDocumentation()

document {
     Key => { (kroneckerNormalForm,Matrix,RingElement,RingElement),
	  kroneckerNormalForm,
	  [kroneckerNormalForm,ChangeMatrix]
	  },
     Headline => "Kronecker normal form of a matrix of linear forms",
     Usage => concatenate("(B,P,Q) = kroneckerNormalForm A\n",
	  "(B,P) = kroneckerNormalForm(A,ChangeMatrix=>{true,false})\n",
	  "(B,Q) = kroneckerNormalForm(A,ChangeMatrix=>{false,true})\n",
	  "B = kroneckerNormalForm(A,ChangeMatrix=>{false,false})\n",
	  "(B,P,Q) = kroneckerNormalForm(A,x,y)\n",
	  "(B,P) = kroneckerNormalForm(A,x,y,ChangeMatrix=>{true,false})\n",
	  "(B,Q) = kroneckerNormalForm(A,x,y,ChangeMatrix=>{false,true})\n",
	  "B = kroneckerNormalForm(A,x,y,ChangeMatrix=>{false,false})"),
     Inputs => {
	  "A" => Matrix => "a matrix of linear forms in two variables",
	  "x" => RingElement,
	  "y" => RingElement
	  },
     Outputs => {
	  "B" => Matrix => {"the Kronecker normal form of ", TT "A"},
	  "P" => Matrix => "invertible (left) change of basis matrix",
	  "Q" => Matrix => "invertible (right) change of basis matrix"
	  },
     "This function produces a matrix of linear forms ", TT "B", ", and invertible matrices ",
     TT "P", " and ", TT "Q", " such that ", TT "B = P*A*Q", " and ", TT "B", 
     " is in Kronecker normal form.  The optional parameters ", TT "x", " and ",
     TT "y", ", if given, should be complementary linear forms. ",
     "They are used to fill in nonzero entries of all blocks in the ",
     "Kronecker normal form for which there is a choice of which linear forms to use. ",
     "If they are not given, they are assumed to be the first two variables of the ring.",
     EXAMPLE lines ///
         R = ZZ/101[x,y]
         A = matrix{{x,y,x,y},{y,x,y,x},{x,y,x,y},{y,y,y,y},{x,x,y,y}}
	 (B,P,Q) = kroneckerNormalForm A
	 B == P*A*Q
	 kroneckerNormalForm(A, ChangeMatrix => {true,false})
	 kroneckerNormalForm(A, ChangeMatrix => {false,true})
	 kroneckerNormalForm(A, ChangeMatrix => {false,false})
	 (B,P,Q) = kroneckerNormalForm(A,x+y,x-y)
	 B == P*A*Q
     ///
     }

document {
     Key => { (rationalNormalForm,Matrix),
	  rationalNormalForm,
	  [rationalNormalForm,ChangeMatrix]
	  },
     Headline => "rational normal form of a matrix",
     Usage => concatenate("(B,P,Q) = rationalNormalForm A\n",
	  "(B,P) = rationalNormalForm(A,ChangeMatrix=>{true,false})\n",
	  "(B,Q) = rationalNormalForm(A,ChangeMatrix=>{false,true})\n",
	  "B = rationalNormalForm(A,ChangeMatrix=>{false,false})"),
     Inputs => {
	  "A" => Matrix => "a square matrix with constant entries"
	  },
     Outputs => {
	  "B" => Matrix => {"the rational normal form of ", TT "A"},
	  "P" => Matrix => "invertible (left) change of basis matrix",
	  "Q" => Matrix => "invertible (right) change of basis matrix"
	  },
     "This function produces a matrix ", TT "B", " in rational normal form, and ",
     "invertible matrices ", TT "P", " and ", TT "Q", " such that ",
     TT "P*Q = I", " and ", TT "B = P*A*Q", ".",
     EXAMPLE lines ///
         R = ZZ/101[x]
	 M = R^4
	 A = random(M,M)
	 factor det(x*id_M - A)
	 (B,P,Q) = rationalNormalForm A
	 B == P*A*Q
	 P*Q == id_M
     ///
     }

document {
     Key => { (kroneckerIndices,Matrix,RingElement,RingElement),
	  kroneckerIndices
	  },
     Headline => "data which classify a matrix pencil up to strict equivalence",
     Usage => "(L, dualL, d) = kroneckerIndices(P,x,y)",
     Inputs => {
	  "P" => Matrix => "a matrix of linear forms",
	  "x" => RingElement,
	  "y" => RingElement
	  },
     Outputs => {
	  "L" => List => "degrees of minimal generators of the kernel",
	  "dualL" => List => "degrees of minimal generators of the kernel of the dual",
	  "d" => List => "elementary divisors"
	  },
     "This function gives data which classify a matrix pencil up to strict ",
     "equivalence.  That is, these data determine the Kronecker normal form ",
     "of the pencil, up to rearranging the blocks.",
     EXAMPLE lines ///
         R = ZZ/101[x,y]
	 A = matrix{{x,y,x,y},{y,x,y,x},{x,y,x,y},{y,y,y,y},{x,x,y,y}}
	 kroneckerNormalForm(A, ChangeMatrix => {false,false})
	 kroneckerIndices(A,x,y)
     ///
     }

-- We compare the Kronecker normal form results for the same matrices,
-- asking for change of basis matrices on the one hand, and asking
-- for them to be omitted on the other.  This results in entirely
-- separate code paths being taken.  The code path where the change
-- of basis matrices are omitted can be easily verified to produce
-- a valid Kronecker normal form, whereas the code path where we ask
-- for change of basis matrices can be easily verified to give
-- a sequence whose product is equal to the input.  Then we compare
-- the two results to test both paths.
--
-- For test matrices, we have chosen examples with regular blocks
-- corresponding to only one eigenvalue or conjugate collection of
-- eigenvalues, or no regular blocks at all, for easier comparison
-- of results.  Otherwise we'd have to worry about rearranging the 
-- blocks.
TEST ///
     obfuscate = method()
     obfuscate(Matrix) := A -> (
	  P := random(target A, target A);
	  while det P == 0 do (
	       P = random(target A, target A);
	       );
	  Q := random(source A, source A);
	  while det Q == 0 do (
	       Q = random(source A, source A);
	       );
	  P*A*Q
	  )
     Q = ZZ/101[x,y]
     use Q
     maxDiff = 8
     maxN = 4
     LL = apply(maxDiff*(maxN+1), i -> (i // maxDiff, (i%maxDiff) + 1))
     L = apply(LL, (n,diff) -> random(Q^(n+diff), Q^n ** Q^{-1})) |
	       apply(LL, (n,diff) -> random(Q^n, Q^(n+diff) ** Q^{-1})) |
	       apply({matrix{{x,0,0,0},{y,x,0,0},{0,y,0,0},{0,0,x-y,0},{0,0,0,x-y}},
			 matrix{{x,y,0,0,0,0},{0,0,x,2*y,0,0},
			      {0,0,-y,x,y,0},{0,0,0,0,x,-2*y},{0,0,0,0,y,x}},
			 matrix{{y,x,0},{0,y,0},{0,0,y}},
			 matrix{{x,y,0,0,0,0},{-2*y,x,y,0,0,0},{0,0,x,y,0,0},
			      {0,0,-2*y,x,0,0},{0,0,0,0,x,y},{0,0,0,0,-2*y,x}}},
		    P -> obfuscate(map(Q^(numgens target P), , P))) |
	       {matrix{{0,x},{0,x-y}}}
     scan(L, P -> assert(kroneckerNormalForm(P, ChangeMatrix => {false, false}) == (kroneckerNormalForm P)#0))
///