+  s1.elem(0,2)*s2.elem(2,1))*s3.elem(1,0);
+  s1.elem(0,2)*s2.elem(2,1))*s3.elem(1,0)  + (s1.elem(0,1)*s2.elem(1,2) -  s1.elem(1,1)*s2.elem(0,2) -  s1.elem(0,2)*s2.elem(1,1) +  s1.elem(1,2)*s2.elem(0,1))*s3.elem(2,0) + (s1.elem(0,2)*s2.elem(1,0) -  s1.elem(1,2)*s2.elem(0,0) -  s1.elem(0,0)*s2.elem(1,2) +  s1.elem(1,0)*s2.elem(0,2))*s3.elem(2,1);
colorContract(const PColorMatrix<T1,3>& s1, const PColorMatrix<T2,3>& s2, const PColorMatrix<T3,3>& s3){  typename TrinaryReturn<PColorMatrix<T1,3>, PColorMatrix<T2,3>, PColorMatrix<T3,3>, FnColorContract>::Type_t  d;  // Permutations: +(0,1,2)+(1,2,0)+(2,0,1)-(1,0,2)-(0,2,1)-(2,1,0)  // d = \epsilon^{i1,j1,k1}\epsilon^{i2,j2,k2} Q1^{i1,i2} Q2^{j1,j2} Q3^{k1,k2}  d.elem() = (s1.elem(0,0)*s2.elem(1,1)           -  s1.elem(1,0)*s2.elem(0,1)           -  s1.elem(0,1)*s2.elem(1,0)           +  s1.elem(1,1)*s2.elem(0,0))*s3.elem(2,2)           + (s1.elem(1,1)*s2.elem(2,2)           -  s1.elem(2,1)*s2.elem(1,2)           -  s1.elem(1,2)*s2.elem(2,1)           +  s1.elem(2,2)*s2.elem(1,1))*s3.elem(0,0)           + (s1.elem(2,2)*s2.elem(0,0)           -  s1.elem(0,2)*s2.elem(2,0)           -  s1.elem(2,0)*s2.elem(0,2)           +  s1.elem(0,0)*s2.elem(2,2))*s3.elem(1,1)           + (s1.elem(1,0)*s2.elem(2,1)           -  s1.elem(2,0)*s2.elem(1,1)           -  s1.elem(1,1)*s2.elem(2,0)           +  s1.elem(2,1)*s2.elem(1,0))*s3.elem(0,2)           + (s1.elem(1,2)*s2.elem(2,0)           -  s1.elem(2,2)*s2.elem(1,0)           -  s1.elem(1,0)*s2.elem(2,2)           +  s1.elem(2,0)*s2.elem(1,2))*s3.elem(0,1)           + (s1.elem(2,0)*s2.elem(0,1)           -  s1.elem(0,0)*s2.elem(2,1)           -  s1.elem(2,1)*s2.elem(0,0)           +  s1.elem(0,1)*s2.elem(2,0))*s3.elem(1,2)           + (s1.elem(2,1)*s2.elem(0,2)           -  s1.elem(0,1)*s2.elem(2,2)           -  s1.elem(2,2)*s2.elem(0,1)           +  s1.elem(0,2)*s2.elem(2,1))*s3.elem(1,0);  return d;}