(* 2017 Oct 30
Minimal  3-tuple program
  Input Data:
	a : positive definite symmetric matrix.
	data1 : candidates of 1st columns of minimal bivectors of a.
	data2 : candidates of 2nd columns of minimal bivectors of a.
	data3 : candidates of 3rd columns of minimal bivectors of a.
	(data1, data2, data3 are computed by qfminim Pari/GP command.)

  Usage:
	min3v[a, data1, data2, data3] : 
		a complete system of representatives of equivalent classes of minimal 3-tuples of a.
	print3v[list] : input list <- min3v[a, data1, data2, data3]:
		print  ( 3-arithmetical mimimam, Rankin invariant, numerical Rankin invariant, cardinal number of 
		the list = min3vector[a, data1, data2, data3]).
	
	p3rank[a, list] : input list <- min3v[a, data1, data2, data3] :
		a pair of the 3-perfection rank of "a" and the dimension of the space of n * n symmetric matrices  
	        (Note: "a" should be converted to floating point to saving time.)                                                                                                *)
	
	
Clear["Global'*"]
min3v[a_, data1_, data2_, data3_] :=
	Module[{amatrix,index3, wedge3,listcolumns, data123,data123regular,position, pwedgesv0,
		     sv0,sv0rep, sv0wedgelist, wedgesv0},
		
		index3[n_] := Sort@Flatten[Table[{k,j+1, i+2}, {i, 1,n-2}, {j, 1, i}, {k,1,j}],2];
		wedge3[x_] := Table[Det[x[[All, i]]], {i, index3[Last[Dimensions[x]]]}];
		listcolumns[x_] :=  Table[x[[;;,i]],{i,1,Last[Dimensions[x]]}];
		(* note that elements of this list are row vectors *)
		amatrix[x_] := x.a.Transpose[x];
		 (* this gives 3 times 3 matrix *)

		data123= Tuples[{Union[listcolumns[data1],-listcolumns[data1]], 
				 Union[listcolumns[data2],-listcolumns[data2]], listcolumns[data3]}];
		(* create a list of 3 times n matrices whose row vectors are in +-data1, +-data2, data3 *)

		data123regular= Select[data123,  
					(MatrixRank[amatrix[#]] ==3 && amatrix[#][[1,2]]>= 0 && amatrix[#][[1,3]]>= 0 &)];
		(* select elements satisfying a part of Minkowski's reduction conditions *)

		sv0= MinimalBy[data123regular, Det[#.a.Transpose[#]] &];

		sv0wedgelist = wedge3 /@ sv0;
		(* old : wedgesv0= DeleteDuplicates@Table[wedge[sv0[[i]]],{i, 1, Length[sv0]}]; *)

		pwedgesv0= Union[DeleteDuplicates@sv0wedgelist, SameTest -> (#1 == #2 || #1 == -#2 &)];

		position = Flatten[Table[FirstPosition[sv0wedgelist, pwedgesv0[[i]], 
						FirstPosition[sv0wedgelist, -pwedgesv0[[i]]]], {i, 1, Length[pwedgesv0]}]];
		sv0rep = sv0[[#]]& @ position
		(* old : sv0rep = Flatten[Table[Select[sv0, (wedge[#] \[Equal] pwedgesv0[[i]] || wedge[#] \[Equal] -pwedgesv0[[i]] &),1], 
					{i, 1, Length[pwedgesv0]}],1] *)
		]

print3v[a_,list_] := 
	{Det[list[[1]].a.Transpose[list[[1]]]],Det[list[[1]].a.Transpose[list[[1]]]]/(Det[a]^(3/Length[a])),
	N[Det[list[[1]].a.Transpose[list[[1]]]]/(Det[a]^(3/Length[a]))],Length[list]}
 

p3rank[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}
		]