(* 2017 Oct 30
Minimal 2-tuple program
  Input Data:
	a : positive definite symmetric matrix.
	data1 : candidates of 1st columns of minimal 2-tuples of a, which is computing by qfminim Pari/GP command.
	data2 : candidates of 2nd columns of minimal 2-tuples of a, which is computing by qfminim Pari/GP command.

  Usage:
	minbivectors[a, data1, data2] : 
		give a complete system of representatives of Subscript[GL, 2](Z)-equivalent classes of minimal 2-tuples of a.
	printminbivectors[list] : input list �� minbivectors[a, data1, data2] :
		print a quintuplets ( 2-arithmetical mimimam, Rankin invariant, numerical Rankin invariant, 
		cardinal number and a list of matrix representations of elements in list = minbivectors[a, data1, data2]).
	p2rank[a, list] : input list �� minbivectors[a, data1, data2] :
		give a pair of the 2-perfection rank of "a" and the dimension of 
		the space of n * n symmetric matrices. (Note: "a" should be converted by floating point to saving time.)
	
	mpfd[a, data1, data2] :
		give all minimal 2-tuples v  modulo {�}I} of a such that v.a.transpose[v]  is 
		contained in Minowski's fundamental domain. 
	mpfdequiv[a, data1, data2] : 
		 give a list of equivalent classes of mpfd[a, data1, data2] .
	printmpfd[list1] : input list1 �� mpfd[a, data1, data2]
		print a list of  matrix representations of elements in list1 = mpfd[a, data1, data2].
	printmpfdequiv[list2] : input list2 �� mpfdequiv[a, data1, data2]:
		print a list of matrix representations of elements in list2 = mpfdequiv[a, data1, data2]. 
		*)
	
	
Clear["Global'*"]
minbivectors[a_, data1_, data2_] :=
	Module[{index2, wedge2,listcolumns, data12,data12regular,sv0, sv0equiv, sv0rep, sv0repreduced, wedgesv0,pwedgesv0},
		
		index2[n_] := Sort@Flatten[Table[{j, i+1}, {i, 1,n-1}, {j, 1, i}],1]; 
		(* index of basis of 2-exterior product *)

		wedge2[x_] := Table[Det[x[[All, i]]], {i, index2[Last[Dimensions[x]]]}]; 
		(* 2-exterior prodct of 2 row vectors of a 2 by n matrix x *)

		listcolumns[x_] :=  Table[x[[;;,i]],{i,1,Last[Dimensions[x]]}]; 
		(* list of column vectors of a given matrix *)

		data12= Tuples[{Union[listcolumns[data1],-listcolumns[data1]], listcolumns[data2]}];
		(* make 2 by n matrices whose first rows are in +- data1 and second rows are in data2 *)

		data12regular= Select[data12,  MatrixRank[#.a.Transpose[#]] ==2 &];
		(* select regular elements from data12 *)

		sv0= MinimalBy[data12regular, Det[#.a.Transpose[#]] &];
		(* select elements with minimal determiant from data12regular *)

		wedgesv0= DeleteDuplicates@Table[wedge2[sv0[[i]]],{i, 1, Length[sv0]}];
		(* list of 2-exterior products of elements in sv0 *)

		pwedgesv0= Union[wedgesv0, SameTest -> (#1 == #2 || #1 == -#2 &)];
		(* list of representatives of wedgesv0 mod +-1 *)

		sv0equiv = DeleteDuplicates@Table[
                           Select[sv0, wedge2[#] == pwedgesv0[[i]]||wedge2[#] ==-pwedgesv0[[i]]&], {i, 1, Length[pwedgesv0]}];
		(* GL(2,Z)-equivalent classes of sv0 *)

		sv0rep = Table[sv0equiv[[i]][[1]], {i, 1, Length[sv0equiv]}]
	        (* list of representatives of sv0equiv *) 
		]

printminbivectors[a_,list_] := 
	{Det[list[[1]].a.Transpose[list[[1]]]],Det[list[[1]].a.Transpose[list[[1]]]]/(Det[a]^(2/Length[a])),
	N[Det[list[[1]].a.Transpose[list[[1]]]]/(Det[a]^(2/Length[a]))],Length[list],
	Table[(Transpose@list[[i]])//MatrixForm, {i, 1, Length[list]}]}
 

p2rank[a_,list_] := 
	Module[{g,n, olist, r,s},
		n = First[Dimensions[a]];
		s = Length[list];
		g = CholeskyDecomposition[a];
		olist = Orthogonalize[#.Transpose[g]] & /@ list ;
		r = MatrixRank@Table[Flatten[Diagonal[Transpose[olist[[i]]].olist[[i]], #]& /@ Range[0,n-1]], {i, 1,s}];
		{r,n*(n+1)/2}
		]
 

mpfd[a_, data1_, data2_] :=
	Module[{listcolumns, data12,data12regular,sv0, sv1},		
		listcolumns[x_] :=  Table[x[[;;,i]],{i,1,Last[Dimensions[x]]}];
		data12= Tuples[{Union[listcolumns[data1],-listcolumns[data1]], listcolumns[data2]}];
		data12regular= Select[data12,  MatrixRank[#.a.Transpose[#]] ==2 &];
		sv0= MinimalBy[data12regular, Det[#.a.Transpose[#]] &];
		sv1 = Select[sv0, 
				(#.a.Transpose[#])[[1,1]] <= (#.a.Transpose[#])[[2,2]] && 0 <= 2*(#.a.Transpose[#])[[1,2]] <= (#.a.Transpose[#])[[1,1]]
				&]
		]


mpfdequiv[a_, data1_, data2_] :=
	Module[{index2, wedge2,listcolumns, data12,data12regular,sv0, sv1,sv1equiv, wedgesv1,pwedgesv1},		
		index2[n_] := Sort@Flatten[Table[{j, i+1}, {i, 1,n-1}, {j, 1, i}],1];
		wedge2[x_] := Table[Det[x[[All, i]]], {i, index2[Last[Dimensions[x]]]}];
		listcolumns[x_] :=  Table[x[[;;,i]],{i,1,Last[Dimensions[x]]}];
		data12= Tuples[{Union[listcolumns[data1],-listcolumns[data1]], listcolumns[data2]}];
		data12regular= Select[data12,  MatrixRank[#.a.Transpose[#]] ==2 &];
		sv0= MinimalBy[data12regular, Det[#.a.Transpose[#]] &];
		sv1 = Select[sv0, 
				(#.a.Transpose[#])[[1,1]] <= (#.a.Transpose[#])[[2,2]] && 0 <= 2*(#.a.Transpose[#])[[1,2]] <= (#.a.Transpose[#])[[1,1]]
				&];
		wedgesv1= DeleteDuplicates@Table[wedge2[sv1[[i]]],{i, 1, Length[sv1]}];
		pwedgesv1= Union[wedgesv1, SameTest -> (#1 == #2 || #1 == -#2 &)];
		sv1equiv = DeleteDuplicates@Table[
                           Select[sv1, wedge2[#] == pwedgesv1[[i]]||wedge2[#] ==-pwedgesv1[[i]]&], {i, 1, Length[pwedgesv1]}]
		]


printmpfd[list1_] := Table[Transpose@list1[[i]]//MatrixForm, {i, 1, Length[list1]}]


printmpfdequiv[list2_] := Column@Table[MatrixForm/@(Transpose/@list2[[i]]), {i, 1, Length[list2]}]