(* We deal with qubits as polynomials in x[i,0],x[i,1] with the
degree in x[i,0],x[i,1] exactly 1.  This is assumed in all of the
following code.*)

(* x[i,j],a[i],b[i],c[i],v[i,j],V,A,M,y[i,j]  are all global variables.  Handle
with care! *)

(* This calculates the element uf, u the sum of matrices
				a[i] b[i]
				c[i] -a[i]
each acting on the i-th factor. *)
Lie[f_] := Block[{i, g}, g = 0; 
    For[i = 0, i < 4, i++, 
      g = g + (f /. {x[i, 0] -> a[i]*x[i, 0] + c[i]*x[i, 1], 
                x[i, 1] -> b[i]*x[i, 0] - a[i]*x[i, 1]})];
    g]

(*This is the bracket op from pxp->k *)
Brack[f_,g_] := Block[{L,M,m,i,j,k,l,u,u1,v1,v,h,h1,h2},L={}; M=Array[m,3];
	For[i=0,i<4,i++,
	u=Coefficient[f,x[i,0]];
	v=Coefficient[g,x[i,0]];u1=u;v1=v;
	For[j=0,j<i,j++,u1=u1/.{x[j,0]->y[j,0],x[j,1]->y[j,1]};
	v1=v1/.{x[j,0]->y[j,0],x[j,1]->y[j,1]}];
	For[j=i+1,j<4,j++,u1=u1/.{x[j,0]->y[j-1,0],x[j,1]->y[j-1,1]};
	v1=v1/.{x[j,0]->y[j-1,0],x[j,1]->y[j-1,1]}];
    h=0;
	For[j=0,j<=1,j++,For[k=0,k<=1,k++,For[l=0,l<=1,l++,
		h = h+(-1)^(j+k+l)Coefficient[u1,y[0,j]*y[1,k]*y[2,l]]*
			Coefficient[v1,y[0,Mod[j+1,2]]*y[1,Mod[k+1,2]]*y[2,Mod[l+1,2]]]]]];
	m[1]=h;
	u=Coefficient[f,x[i,0]];
	v=Coefficient[g,x[i,1]];
    u1=u;v1=v;
	For[j=0,j<i,j++,u1=u1/.{x[j,0]->y[j,0],x[j,1]->y[j,1]};
	v1=v1/.{x[j,0]->y[j,0],x[j,1]->y[j,1]}];
	For[j=i+1,j<4,j++,u1=u1/.{x[j,0]->y[j-1,0],x[j,1]->y[j-1,1]};
	v1=v1/.{x[j,0]->y[j-1,0],x[j,1]->y[j-1,1]}];
    h1=0;
	For[j=0,j<=1,j++,For[k=0,k<=1,k++,For[l=0,l<=1,l++,
		h1 = h1+(-1)^(j+k+l)Coefficient[u1,y[0,j]*y[1,k]*y[2,l]]*
			Coefficient[v1,y[0,Mod[j+1,2]]*y[1,Mod[k+1,2]]*y[2,Mod[l+1,2]]]]]];
	u=Coefficient[f,x[i,1]];
	v=Coefficient[g,x[i,0]];
	 u1=u;v1=v;
    For[j=0,j<i,j++,u1=u1/.{x[j,0]->y[j,0],x[j,1]->y[j,1]};
	v1=v1/.{x[j,0]->y[j,0],x[j,1]->y[j,1]}];
	For[j=i+1,j<4,j++,u1=u1/.{x[j,0]->y[j-1,0],x[j,1]->y[j-1,1]};
	v1=v1/.{x[j,0]->y[j-1,0],x[j,1]->y[j-1,1]}];
    h2=0;
	For[j=0,j<=1,j++,For[k=0,k<=1,k++,For[l=0,l<=1,l++,
		h2 = h2+(-1)^(j+k+l)Coefficient[u1,y[0,j]*y[1,k]*y[2,l]]*
			Coefficient[v1,y[0,Mod[j+1,2]]*y[1,Mod[k+1,2]]*y[2,Mod[l+1,2]]]]]];
	m[2]=h1+h2;
	u=Coefficient[f,x[i,1]];
	v=Coefficient[g,x[i,1]];
	 u1=u;v1=v;
    For[j=0,j<i,j++,u1=u1/.{x[j,0]->y[j,0],x[j,1]->y[j,1]};
	v1=v1/.{x[j,0]->y[j,0],x[j,1]->y[j,1]}];
	For[j=i+1,j<4,j++,u1=u1/.{x[j,0]->y[j-1,0],x[j,1]->y[j-1,1]};
	v1=v1/.{x[j,0]->y[j-1,0],x[j,1]->y[j-1,1]}];
    h=0;
	For[j=0,j<=1,j++,For[k=0,k<=1,k++,For[l=0,l<=1,l++,
		h = h+(-1)^(j+k+l)Coefficient[u1,y[0,j]*y[1,k]*y[2,l]]*
			Coefficient[v1,y[0,Mod[j+1,2]]*y[1,Mod[k+1,2]]*y[2,Mod[l+1,2]]]]]];
	m[3]=h;
	L=Append[L,M]];
	L]

(* computes ad(u)f with u in sl(2)^4  given as in the output of
Brack[] *)

ad[u_,f_] := Block[{i,j,S},
				S = {};
				For[i=0,i<4,i++,
					  S = Append[S,b[i]->2u[[i+1,1]]];
					  S = Append[S,c[i]->-2u[[i+1,3]]];
					  S = Append[S,a[i]->-u[[i+1,2]]]];
			    Expand[Lie[f] /.S]]
					  
(*This makes a 3 qubit product vector *)
V = Product[v[i,0]*x[i,0] + v[i,1]*x[i,1],{i,1,3}];

(* This calculates the transformation in two qubits corresponding to
the 4 qubit vector corresponding to f.*)

Transform[f_] := Block[{L,l,i,j,p,q,g},
			L = Array[l,{4,4}];
			For[i=0,i<2,i++,For[j=0,j<2,j++,
				g = ((-1)^(i+j))*Coefficient[f,x[0,1-i]*x[1,1-j]];
				For[p=0,p<2,p++,For[q=0,q<2,q++,
					L[[i*2+j+1,p*2+q+1]]=Coefficient[g,x[2,p]*x[3,q]]]]]];
			Transpose[L]]

(* The transpose relative to the SL(2)xSL(2) invariant form is 
	A.Transpose[T].A *)
A = {{0,0,0,1},{0,0,-1,0},{0,-1,0,0},{1,0,0,0}};

M = {a[0],a[1],a[2],a[3],b[0],b[1],b[2],b[3],c[0],c[1],c[2],c[3]};

(* This leads to the example:
   x=|1000>+|0011>+|0101>+|0110>
   y=|1110>+|1101>+|1011>
Then [x,y]=-4f_1-4f_2-4f_3.  Take z=[x,y].  Then [z,x] = 
4|1100>+4|1010>+4|1001>+12|0111>=w.  [w,x]=6h_0 -4h_1 - 4h_2 - 4h_3.
So the pair (x,y) is not in the nullcone.  However, tx+sy is
in the nullcone as is [x,y]. *)

(*Completed the calculation in 4qubits-success.nb*)

B[i_,u_] := 2 u[[i,1]] u[[i,3]] - (u[[i,2]]^2)/2