(* ::Package:: *)

BeginPackage["MatrixManifolds`"];


G2::usage = "Lie Algebra g2 matrices of the Italians";
Gg2::usage = "Lie Algebra g2 matrices of Ruben Arenas PhD thesis";
sG2::usage = "s of the Italian paper";
\[Sigma]G2::usage = "\[Sigma] of the Italian paper";
uSO4::usage = "SO(4) parametrization of the Italian paper";
gG2::usage = "G2 parametrization of the Italian paper";
paramBooleG2::usage = "Parameter boundaries";
LOctMultMatrix::usage = "Octonian matrix multi LEFT";
ROctMultMatrix::usage = "Octonian matrix multi RIGHT";
OctMult::usage = "Octonion mult";
OctCommutator::usage = "Octonion commutator";
OctAssociator::usage = "Octonion Associator";
OctReal::usage = "Octonion Re as a single Real";
OctReal8::usage = "Octonion Re as 8 dim array with last 7 0";
realMult::usage = "Real multiplication stored in an array of Reals padded with 0s";
realConjugate::usage = "Real conjugate which is itself";
OctComplex8::usage = "Complex numbers stored as 2 Reals padded with 6 0s";
complexMult::usage ="complex numbers multiplied as array of 8 reals, 6 0s";
complexConjugate::usage = "complex conjugation with 6 padded 0s";
OctQuaterion8::usage = "Quaternion with 4 padded 0s";
quaternionMult::usage = "Quaternion multiplication with 4 padded 0s";
quaternionConjugate::usage = "conjugate quaternion with 4 padded 0s";
OctIm::usage = "Octonion Im as 8 dim array, with 0 for Re";
OctIm7::usage = "Octonion Im as 7 dim array, Re dropped";
OctExp::usage = "EXP for Octonions";
OctLog::usage = "Log for Octonions";
OctForget::usage ="forget representation in complex numbers";
MatrixCommutator::usage = "General matrix commutator";
OctConjugator::usage = "Conjugate of an Octonion";
OctConjugate::usage = "Conjugate of an Octonion";
JacobiIdentity::usage = "JacobiIdentity";
OctPower::usage = "positive integer powers of octonions";
OctExpSeries::usage = "Taylor series expansion of Exp for Octonions";
OctBCHsym::usage = "non-evaluated symbolic version of the OctBCH";
OctBCH::usage = "BCH formula best for unit vectors ";
BCH::usage = "General BCH formula";
actionG2::usage = "7x7 G2 matrix acting on the Im part of an octonion, leaving Re intact";
oneparamgroupmaker::usage = "1-param subgroup";
SO::usage = "SO(n) Lie Algebra basis";
SO4Basis::usage="Basis for SO(4)";

quaternionMatrices::usage = "init call for quaternion matrices";
Rq1:="Basis matrix for Quaternions, Right handside";
Rq2:="Basis matrix for Quaternions, Right handside";
Rq3:="Basis matrix for Quaternions, Right handside";
Rq4:="Basis matrix for Quaternions, Right handside";

Lq1:="Basis matrix for Quaternions, Left handside";
Lq2:="Basis matrix for Quaternions, Left handside";
Lq3:="Basis matrix for Quaternions, Left handside";
Lq4:="Basis matrix for Quaternions, Left handside";

C1::usage = "Lie Algebra g2 matrices of the Italians";
C2::usage = "Lie Algebra g2 matrices of the Italians";
C3::usage = "Lie Algebra g2 matrices of the Italians";
C4::usage = "Lie Algebra g2 matrices of the Italians";
C5::usage = "Lie Algebra g2 matrices of the Italians";
C6::usage = "Lie Algebra g2 matrices of the Italians";
C7::usage = "Lie Algebra g2 matrices of the Italians";
C8::usage = "Lie Algebra g2 matrices of the Italians";
C9::usage = "Lie Algebra g2 matrices of the Italians";
C10::usage = "Lie Algebra g2 matrices of the Italians";
C11::usage = "Lie Algebra g2 matrices of the Italians";
C12::usage = "Lie Algebra g2 matrices of the Italians";
C13::usage = "Lie Algebra g2 matrices of the Italians";
C14::usage = "Lie Algebra g2 matrices of the Italians";

G2Rarenas::usage = "initialize lie Algebra g2 matrices of Rarenas PhD";

X1::usage = "Lie Algebra g2 matrices of the Rarenas PhD";
X2::usage = "Lie Algebra g2 matrices of the Rarenas PhD";
X3::usage = "Lie Algebra g2 matrices of the Rarenas PhD";
X4::usage = "Lie Algebra g2 matrices of the Rarenas PhD";
X5::usage = "Lie Algebra g2 matrices of the Rarenas PhD";
X6::usage = "Lie Algebra g2 matrices of the Rarenas PhD";
X7::usage = "Lie Algebra g2 matrices of the Rarenas PhD";
Y1::usage = "Lie Algebra g2 matrices of the Rarenas PhD";
Y2::usage = "Lie Algebra g2 matrices of the Rarenas PhD";
Y3::usage = "Lie Algebra g2 matrices of the Rarenas PhD";
Y4::usage = "Lie Algebra g2 matrices of the Rarenas PhD";
Y5::usage = "Lie Algebra g2 matrices of the Rarenas PhD";
Y6::usage = "Lie Algebra g2 matrices of the Rarenas PhD";
Y7::usage = "Lie Algebra g2 matrices of the Rarenas PhD";

hyperCube::usage = "generates a unit hypercube";
randomGrid2d::usage = "generates a 2D random grid of reals";
randomGrid3d::usage = "generates a 3D random grid of reals";
randomGrid4d::usage = "generates a 4D random grid of reals";
randomGrid::usage = "generates a general random grid of reals";
kolComplexity::usage = "Approximates Kolmogorov Complexity on delimited input";





Begin["`Private`"];




(* Author Dara O Shayda, July 2012 *)
(* License http://www.apache.org/licenses/LICENSE-2.0.html *)




(* Basis for the Lie algebra g2 *)
(* C1 to C14 are golobals to avoid messy arguments to other functions *)
G2[]:=Module[{C = {}},
C1 = {{0,0,0,0,0,0,0},{0,0,0,0,0,0,0}, {0,0,0,0,0,0,0}, {0,0,0,0,0,0,-1},
{0,0,0,0,0,-1,0}, {0,0,0,0,1,0,0}, {0,0,0,1,0,0,0}};
C =Append[C, C1];

C2 = {{0,0,0,0,0,0,0},{0,0,0,0,0,0,0}, {0,0,0,0,0,0,0}, {0,0,0,0,0,1,0},
{0,0,0,0,0,0,-1}, {0,0,0,-1,0,0,0}, {0,0,0,0,1,0,0}};
C =Append[C, C2];


C3 = {{0,0,0,0,0,0,0},{0,0,0,0,0,0,0}, {0,0,0,0,0,0,0}, {0,0,0,0,-1,0,0},
{0,0,0,1,0,0,0}, {0,0,0,0,0,0,-1}, {0,0,0,0,0,1,0}};
C =Append[C, C3];


C4 = {{0,0,0,0,0,0,0},{0,0,0,0,0,0,1}, {0,0,0,0,0,1,0}, {0,0,0,0,0,0,0},
{0,0,0,0,0,0,0}, {0,0,-1,0,0,0,0}, {0,-1,0,0,0,0,0}};
C =Append[C, C4];


C5 = {{0,0,0,0,0,0,0},{0,0,0,0,0,-1,0}, {0,0,0,0,0,0,1}, {0,0,0,0,0,0,0},
{0,0,0,0,0,0,0}, {0,1,0,0,0,0,0}, {0,0,-1,0,0,0,0}};
C =Append[C, C5];


C6 = {{0,0,0,0,0,0,0},{0,0,0,0,1,0,0}, {0,0,0,-1,0,0,0}, {0,0,1,0,0,0,0},
{0,-1,0,0,0,0,0}, {0,0,0,0,0,0,0}, {0,0,0,0,0,0,0}};
C =Append[C, C6];


C7 = {{0,0,0,0,0,0,0},{0,0,0,-1,0,0,0}, {0,0,0,0,-1,0,0}, {0,1,0,0,0,0,0},
{0,0,1,0,0,0,0}, {0,0,0,0,0,0,0}, {0,0,0,0,0,0,0}};
C =Append[C, C7];


C8 = {{0,0,0,0,0,0,0},{0,0,-2,0,0,0,0}, {0,2,0,0,0,0,0}, {0,0,0,0,1,0,0},
{0,0,0,-1,0,0,0}, {0,0,0,0,0,0,-1}, {0,0,0,0,0,1,0}} / Sqrt[3];
C =Append[C, C8];


C9 = {{0,-2,0,0,0,0,0},{2,0,0,0,0,0,0}, {0,0,0,0,0,0,0}, {0,0,0,0,0,0,1},
{0,0,0,0,0,-1,0}, {0,0,0,0,1,0,0}, {0,0,0,-1,0,0,0}} / Sqrt[3];
C =Append[C, C9];


C10 = {{0,0,-2,0,0,0,0},{0,0,0,0,0,0,0}, {2,0,0,0,0,0,0}, {0,0,0,0,0,-1,0},
{0,0,0,0,0,0,-1}, {0,0,0,1,0,0,0}, {0,0,0,0,1,0,0}} / Sqrt[3];
C =Append[C, C10];


C11 = {{0,0,0,-2,0,0,0},{0,0,0,0,0,0,-1}, {0,0,0,0,0,1,0}, {2,0,0,0,0,0,0},
{0,0,0,0,0,0,0}, {0,0,-1,0,0,0,0}, {0,1,0,0,0,0,0}} / Sqrt[3];
C =Append[C, C11];

C12 = {{0,0,0,0,-2,0,0},{0,0,0,0,0,1,0}, {0,0,0,0,0,0,1}, {0,0,0,0,0,0,0},
{2,0,0,0,0,0,0}, {0,-1,0,0,0,0,0}, {0,0,-1,0,0,0,0}} / Sqrt[3];
C =Append[C, C12];

C13 = {{0,0,0,0,0,-2,0},{0,0,0,0,-1,0,0}, {0,0,0,-1,0,0,0}, {0,0,1,0,0,0,0},
{0,1,0,0,0,0,0}, {2,0,0,0,0,0,0}, {0,0,0,0,0,0,0}} / Sqrt[3];
C =Append[C, C13];


C14 = {{0,0,0,0,0,0,-2},{0,0,0,1,0,0,0}, {0,0,0,0,-1,0,0}, {0,-1,0,0,0,0,0},
{0,0,1,0,0,0,0}, {0,0,0,0,0,0,0}, {2,0,0,0,0,0,0}} / Sqrt[3];
C =Append[C, C14];

C 
]

paramBooleG2[r_] :=Boole[(0<=r[[1]]<=Pi) &&
(0<=r[[2]]<=Pi/2 )&& 
(0<=r[[3]]<=Pi/2) &&
(0<=r[[4]]<=2*Pi) &&
(0<=r[[5]]<=Pi/4) &&
(0<=r[[6]]<=Pi )&&
(0<=r[[7]]<=Pi/6) &&
(0<=r[[8]]<=Pi/2) &&
((3*r[[7]]) <= r[[8]]) &&
(0<=r[[9]]<=2*Pi) &&
(     0<=r[[10]]<=Pi/2) &&
(     0<=r[[11]]<=Pi) &&
(     0<=r[[12]]<=Pi) &&
(     0<=r[[13]]<=Pi/2) &&
(     0<=r[[14]]<=Pi) ]




sG2[x_, y_, z_]:= MatrixExp[x*C3].MatrixExp[y*C2].MatrixExp[z*C3]

\[Sigma]G2[x_,y_,z_] :=  MatrixExp[Sqrt[3]*x*C8].MatrixExp[Sqrt[3]*y*C9].MatrixExp[Sqrt[3]*z*C8]

uSO4[x1_, y1_, z1_, x2_, y2_,z2_]:= sG2[x1, y1, z1].\[Sigma]G2[x2,y2,z2] 

gG2 [r_] := N[\[Sigma]G2[r[[1]], r[[2]], r[[3]]].sG2[r[[4]], r[[5]], r[[6]]].N[MatrixExp[Sqrt[3]*r[[7]]*C11].MatrixExp[r[[8]]*C5]].uSO4[r[[9]], r[[10]], r[[11]],r[[12]], r[[13]], r[[14]]]]


G2Rarenas [] := Module[{X={}, Y={}},


X1 = {{0,0,0,0,0,0,0},{0,0,1,0,0,0,0}, {0,-1,0,0,0,0,0}, {0,0,0,0,-1,0,0},
{0,0,0,1,0,0,0}, {0,0,0,0,0,0,0}, {0,0,0,0,0,0,0}} ;

X =Append[X, X1];


Y1 = {{0,0,0,0,0,0,0},{0,0,0,0,0,0,0}, {0,0,0,0,0,0,0}, {0,0,0,0,1,0,0},
{0,0,0,-1,0,0,0}, {0,0,0,0,0,0,1}, {0,0,0,0,0,-1,0}} ;

Y =Append[Y, Y1];


X2 = {{0,0,-1,0,0,0,0},{0,0,0,0,0,0,0}, {1,0,0,0,0,0,0}, {0,0,0,0,0,-1,0},
{0,0,0,0,0,0,0}, {0,0,0,1,0,0,0}, {0,0,0,0,0,0,0}} ; 

X =Append[X, X2];


Y2 = {{0,0,0,0,0,0,0},{0,0,0,0,0,0,0}, {0,0,0,0,0,0,0}, {0,0,0,0,0,1,0},
{0,0,0,0,0,0,-1}, {0,0,0,-1,0,0,0}, {0,0,0,0,1,0,0}} ;

Y =Append[Y, Y2];


X3 = {{0,1,0,0,0,0,0},{-1,0,0,0,0,0,0}, {0,0,0,0,0,0,0}, {0,0,0,0,0,0,-1},
{0,0,0,0,0,0,0}, {0,0,0,0,0,0,0}, {0,0,0,1,0,0,0}} ;

X =Append[X, X3];


Y3 = {{0,0,0,0,0,0,0},{0,0,0,0,0,0,0}, {0,0,0,0,0,0,0}, {0,0,0,0,0,0,1},
{0,0,0,0,0,1,0}, {0,0,0,0,-1,0,0}, {0,0,0,-1,0,0,0}} ;

Y =Append[Y, Y3];

X4 = {{0,0,0,0,-1,0,0},{0,0,0,0,0,1,0}, {0,0,0,0,0,0,0}, {0,0,0,0,0,0,0},
{1,0,0,0,0,0,0}, {0,-1,0,0,0,0,0}, {0,0,0,0,0,0,0}} ;

X =Append[X, X4];


Y4 = {{0,0,0,0,0,0,0},{0,0,0,0,0,-1,0}, {0,0,0,0,0,0,1}, {0,0,0,0,0,0,0},
{0,0,0,0,0,0,0}, {0,1,0,0,0,0,0}, {0,0,-1,0,0,0,0}} ; 

Y =Append[Y, Y4];


X5 = {{0,0,0,1,0,0,0},{0,0,0,0,0,0,1}, {0,0,0,0,0,0,0}, {-1,0,0,0,0,0,0},
{0,0,0,0,0,0,0}, {0,0,0,0,0,0,0}, {0,-1,0,0,0,0,0}} ;

X =Append[X, X5];


Y5 = {{0,0,0,0,0,0,0},{0,0,0,0,0,0,1}, {0,0,0,0,0,1,0}, {0,0,0,0,0,0,0},
{0,0,0,0,0,0,0}, {0,0,-1,0,0,0,0}, {0,-1,0,0,0,0,0}} ;

Y =Append[Y, Y5];



X6 = {{0,0,0,0,0,0,-1},{0,0,0,1,0,0,0}, {0,0,0,0,0,0,0}, {0,-1,0,0,0,0,0},
{0,0,0,0,0,0,0}, {0,0,0,0,0,0,0}, {1,0,0,0,0,0,0}} ;

X =Append[X, X6];


Y6 = {{0,0,0,0,0,0,1},{0,0,0,0,0,0,0}, {0,0,0,0,1,0,0}, {0,0,0,0,0,0,0},
{0,0,-1,0,0,0,0}, {0,0,0,0,0,0,0}, {-1,0,0,0,0,0,0}} ; 

Y =Append[Y, Y6];


X7 = {{0,0,0,0,0,1,0},{0,0,0,0,1,0,0}, {0,0,0,0,0,0,0}, {0,0,0,0,0,0,0},
{0,-1,0,0,0,0,0}, {-1,0,0,0,0,0,0}, {0,0,0,0,0,0,0}} ;

X =Append[X, X7];


Y7 = {{0,0,0,0,0,0,0},{0,0,0,0,-1,0,0}, {0,0,0,1,0,0,0}, {0,0,-1,0,0,0,0},
{0,1,0,0,0,0,0}, {0,0,0,0,0,0,0}, {0,0,0,0,0,0,0}} ;

Y =Append[Y, Y7];

{X, Y, {X1, X2, X3, X4, X5, X6, X7, Y1, Y2, Y3, Y4, Y5, Y6, Y7}}

]



(* The Octonion arithmatic *)

LOctMultMatrix[a0_]:=Module[{a=a0, w},
(* Make the left-hand w matrix of Reals *)
w={{a[[1]], -a[[2]], -a[[3]],-a[[4]],-a[[5]], -a[[6]], -a[[7]], -a[[8]]},
{a[[2]], a[[1]], -a[[4]],a[[3]],-a[[6]], a[[5]], a[[8]], -a[[7]]},
{a[[3]], a[[4]], a[[1]],-a[[2]],-a[[7]], -a[[8]], a[[5]], a[[6]]},
{a[[4]], -a[[3]], a[[2]],a[[1]],-a[[8]], a[[7]], -a[[6]],a[[5]]},
{a[[5]], a[[6]], a[[7]],a[[8]],a[[1]], -a[[2]], -a[[3]], -a[[4]]},
{a[[6]], -a[[5]], a[[8]],-a[[7]],a[[2]], a[[1]], a[[4]], -a[[3]]},
{a[[7]], -a[[8]], -a[[5]],a[[6]],a[[3]], -a[[4]], a[[1]], a[[2]]},
{a[[8]], a[[7]], -a[[6]],-a[[5]],a[[4]], a[[3]], -a[[2]], a[[1]]}}
]

ROctMultMatrix[b0_]:=Module[{ b=b0, v},
(* Make the right-hand v matrix of Reals *)
v={{b[[1]], -b[[2]], -b[[3]],-b[[4]],-b[[5]], -b[[6]], -b[[7]], -b[[8]]},
{b[[2]], b[[1]],b[[4]],-b[[3]],b[[6]], -b[[5]], -b[[8]], b[[7]]},
{b[[3]], -b[[4]], b[[1]],b[[2]],b[[7]], b[[8]], -b[[5]],- b[[6]]},
{b[[4]], b[[3]], -b[[2]],b[[1]],b[[8]], -b[[7]], b[[6]],-b[[5]]},
{b[[5]], -b[[6]], -b[[7]],-b[[8]],b[[1]], b[[2]], b[[3]],b[[4]]},
{b[[6]], b[[5]], -b[[8]],b[[7]],-b[[2]], b[[1]], -b[[4]], b[[3]]},
{b[[7]], b[[8]], b[[5]],-b[[6]],-b[[3]], b[[4]], b[[1]],-b[[2]]},
{b[[8]], -b[[7]], b[[6]],b[[5]],-b[[4]], -b[[3]], b[[2]], b[[1]]}}
]



quaternionMatrices[]:= Module[{},

Rq1 = IdentityMatrix[4];
Rq2 ={{0,-1,0,0},{1,0,0,0},{0,0,0,1},{0,0,-1,0}};
Rq3 ={{0,0,-1,0},{0,0,0,-1},{1,0,0,0},{0,1,0,0}};
Rq4 = {{0,0,0,-1}, {0,0,1,0}, {0,-1,0,0},{1,0,0,0}};

Lq1 = IdentityMatrix[4];
Lq2 ={{0,-1,0,0},{1,0,0,0},{0,0,0,-1},{0,0,1,0}};
Lq3 ={{0,0,-1,0},{0,0,0,1},{1,0,0,0},{0,-1,0,0}};
Lq4 = {{0,0,0,-1}, {0,0,-1,0}, {0,1,0,0},{1,0,0,0}};


]


OctMult[a0_, b0_]:=Module[{a=a0,b=b0, v},
v = ROctMultMatrix[b];
v.a]

OctCommutator[a0_, b0_]:= Module[{a=a0,b=b0},

OctMult[a,b]-OctMult[b,a]]

OctAssociator[a0_, b0_,c0_]:=Module[{a=a0,b=b0,c=c0},
OctMult[OctMult[a,b],c]-OctMult[a,OctMult[b,c]]]

OctReal[o0_]:=Module[{o=o0},
o[[1]]]

OctReal8[o0_]:=Module[{o=o0},
{o[[1]], 0,0,0,0,0,0,0}]

OctComplex8[o0_]:=Module[{o=o0},
{o[[1]], o[[2]],0,0,0,0,0,0}]

OctQuaterion8[o0_]:=Module[{o=o0},
{o[[1]], o[[2]],o[[3]],o[[4]],0,0,0,0}]




realMult[a0_, b0_]:=Module[{a=a0, b=b0},

OctMult[OctReal8[a], OctReal8[b]]

]

realConjugate[a0_]:= Module[{a=a0},
a
]

complexMult[a0_, b0_]:=Module[{a=a0, b=b0},

OctMult[OctComplex8[a], OctComplex8[b]]

]

complexConjugate[a0_]:=Module[{a=a0},
OctConjugate[a]
]

quaternionConjugate[a0_]:=Module[{a=a0},
OctConjugate[a]
]

quaternionMult[a0_, b0_]:=Module[{a=a0, b=b0},

OctMult[OctQuaterion8[a], OctQuaterion8[b]]

]



OctIm[o0_]:=Module[{o=o0,w},
w=o;
w[[1]] = 0; 
w]

(* same as OctIm but returns a vector of 7 numbers *)
OctIm7 [o0_]:=Module[{o=o0},

Drop[o, 1]

]

OctExp[o0_]:=Module[{o=o0,imo},

imo = OctIm[o];

 Exp[OctReal[o]]*({Cos[Norm[imo]],0,0,0,0,0,0,0}+(Sin[Norm[imo]]/Norm[imo])*imo)
]

OctLog[o0_, k0_]:=Module[{o=o0,yu, ux, k=k0, y = o0},

yu = y/Norm[y];

ux = ArcCos[OctReal[yu]]+ 2*Pi*k; 

(* Try to reverse the exponentiation formular *)
N[Prepend[OctIm7[yu]*ux/Sin[ux],Log[Norm[y]]]]

]


MatrixCommutator[a0_, b0_]:= Module[{a=a0,b=b0},
Simplify[(a.b)-(b.a)]]

JacobiIdentity[a0_, b0_,c0_]:=Module[{a=a0,b=b0,c=c0, w1, w2, w3},
w1 = Simplify[MatrixCommutator[a,MatrixCommutator[b,c]]];
w2 = Simplify[MatrixCommutator[b,MatrixCommutator[c,a]]];
w3 = Simplify[MatrixCommutator[c,MatrixCommutator[a,b]]];

Simplify[w1+w2+w3]]


OctConjugator[a0_]:=Module[{a = a0},

Simplify[((2*a[[1,1]])*Inverse[a])-IdentityMatrix[8]]
]

OctConjugate[a0_]:=Module[{a=a0},
{a[[1]], -a[[2]], -a[[3]],-a[[4]],-a[[5]],-a[[6]],-a[[7]],-a[[8]]}
]

actionG2[g0_, x0_]:=Module[{x = x0, g = g0, b},

(* Simply mulitpl the Subscript[G, 2] element as if a 7x7 matrice, and do so from left *)
b=g.OctIm7[x];

(* take the 7 numbers from that multiplication and concatenate the original Real part
of the Octonion *)
PrependTo[b,OctReal[x]]]



OctPower[x0_, n0_]:=Module[{x=x0, n=n0, pwr={1, 0,0,0,0,0,0,0}},
Table[pwr=OctMult[x, pwr], {i, 1, n}];
pwr]

OctExpSeries[x0_, n0_]:=Module[{x=x0, n=n0}, Total[Table[OctPower[x, i]/i!, {i, 0, n}]]
]

OctBCHsym[x0_, y0_]:= Module[{xO=x0, yO=y0},
(OctCommutator[xO,yO]/2)+(OctCommutator[xO, OctCommutator[xO,yO]]/12)-(OctCommutator[yO, OctCommutator[xO,yO]]/12)-(OctCommutator[yO,OctCommutator[xO, OctCommutator[xO,yO]]]/24)-
((OctCommutator[OctCommutator[OctCommutator[OctCommutator[xO,yO],yO ],yO] ,yO]+ OctCommutator[OctCommutator[OctCommutator[OctCommutator[yO,xO],xO ],xO] ,xO])/720)+  
((OctCommutator[OctCommutator[OctCommutator[OctCommutator[xO,yO],yO ],yO] ,xO]+ OctCommutator[OctCommutator[OctCommutator[OctCommutator[yO,xO],xO ],xO] ,yO])/360)+
((OctCommutator[OctCommutator[OctCommutator[OctCommutator[xO,yO],xO ],yO] ,xO]+ OctCommutator[OctCommutator[OctCommutator[OctCommutator[yO,xO],yO ],xO] ,yO])/120)    
]

OctBCH[x0_, y0_]:= Module[{xO=x0, yO=y0},
N[(OctCommutator[xO,yO]/2)+(OctCommutator[xO, OctCommutator[xO,yO]]/12)-(OctCommutator[yO, OctCommutator[xO,yO]]/12)-(OctCommutator[yO,OctCommutator[xO, OctCommutator[xO,yO]]]/24)-
((OctCommutator[OctCommutator[OctCommutator[OctCommutator[xO,yO],yO ],yO] ,yO]+ OctCommutator[OctCommutator[OctCommutator[OctCommutator[yO,xO],xO ],xO] ,xO])/720)+  
((OctCommutator[OctCommutator[OctCommutator[OctCommutator[xO,yO],yO ],yO] ,xO]+ OctCommutator[OctCommutator[OctCommutator[OctCommutator[yO,xO],xO ],xO] ,yO])/360)+
((OctCommutator[OctCommutator[OctCommutator[OctCommutator[xO,yO],xO ],yO] ,xO]+ OctCommutator[OctCommutator[OctCommutator[OctCommutator[yO,xO],yO ],xO] ,yO])/120)    ]
]



OctForget[x0_]:=Module[{x=x0},

Exp[x[[1]]]*(Cos[Norm[OctIm7[x]]]+ I*Sin[Norm[OctIm7[x]]])
]


p[n_]:=p[n]=(F=Table[1/(j-i)!,{i,n+1},{j,n+1}];G=Table[1/(j-i)! Product[s[k],{k,i,j-1}],{i,n+1},{j,n+1}];qthpower=IdentityMatrix[n+1];FGm1=F.G-qthpower;Expand[-Sum[qthpower=qthpower.FGm1;(-1)^q/q qthpower,{q,n}][[1,n+1]]])

translated[n_]:=(temp=Expand[Product[s[k]^2,{k,n}] p[n]];Sum[term=Apply[List,temp[[i]]];term[[1]] Apply[Dot,Take[term,-n]/.{s[i_]^2->Symbol["x"],s[i_]^3->Symbol["y"]}],{i,Length[temp]}])


(* This algorithm would not work for the Octonions since there is no associativity *)
(* It works for the matrices since their products are associative *)
(* Alogorithm by Matthias W. Reinsch *)
BCH[A0_, B0_, n0_]:= Module[{A=A0, B= B0, n = n0, bch, i},

bch=translated[2]/.{Symbol["x"]->A, Symbol["y"]->B};

(* It has to be Symbol["x"] or it will get mixed up with local symbols as vars *)

Table[
bch = (translated[i]/.{Symbol["x"]->A, Symbol["y"]->B})+bch

, {i, 3, n}];

bch]






SO4Basis[x0_]:= Module[{x=x0},

{{{0,Symbol[x],0,0},{-Symbol[x],0,0,0},{0,0,0,Symbol[x]},{0,0,-Symbol[x],0}}, {{0,0,0,Symbol[x]},{0,0,-Symbol[x],0},{0,Symbol[x],0,0},{-Symbol[x],0,0,0}},{{0,0,Symbol[x],0},{0,0,0,-Symbol[x]},{-Symbol[x],0,0,0},{0,Symbol[x],0,0}},{{0,0,-Symbol[x],0},{0,0,0,-Symbol[x]},{Symbol[x],0,0,0},{0,Symbol[x],0,0}},{{0,0,0,Symbol[x]},{0,0,-Symbol[x],0},{0,Symbol[x],0,0},{-Symbol[x],0,0,0}},{{0,-Symbol[x],0,0},{Symbol[x],0,0,0},{0,0,0,Symbol[x]},{0,0,-Symbol[x],0}}}
]



SO[n0_]:=Module[{n=n0, i,j,k, matlist},

matlist=Table[ ConstantArray[0,{n,n}], {i, 1, n*(n-1)/2}];

indeces=Flatten[Table[Table[{i, j}, {j, 1,i-1}], {i, 2, n}],1];

Table[matlist[[k]][[indeces[[k]][[1]]]][[indeces[[k]][[2]]]] = 1;
matlist[[k]][[indeces[[k]][[2]]]][[indeces[[k]][[1]]]] = -1,{k, 1, n*(n-1)/2}];

matlist
]


oneparamgroupmaker[x0_,\[Theta]0_,exptext0_,dim0_, a0_, cylrad0_,n0_, range0_,arrowconfig0_]:=Module[
{x = x0,\[Theta]=\[Theta]0, exptext=exptext0,dim=dim0, a = a0, cylrad=cylrad0,n=n0,range=range0,arrowconfig=arrowconfig0,b, pixs,oneparamgroup,distances1,distances2,distances3,poly1,poly2,path,flow,
flow2,m,t,arrowheadsize ,arrowheadoffset ,arrowheadtextoffset ,separation,vector1,vector2,vectortext,vectortext2 ,exptextvector,exptextvector2,cylheight,cyl,expedgefactor,liealgebra,liegroup  },

pixs = 4;
(* Some factor needs to scale m in order for the strips not to be too long *)
m =Min[Round[Abs[Norm[MatrixExp[1*x] - MatrixExp[IdentityMatrix[dim]]]]]+5,20];  

(*m = If[m < 20,20,m]; *)
(* how much of the tail is used *)
b=2*m;
(* these are really tied to the a size *)
arrowheadsize = arrowconfig[[1]];
arrowheadoffset = arrowconfig[[2]];
arrowheadtextoffset =arrowconfig[[3]]; 
separation = arrowconfig[[4]]; (* Scale the Norm between two consecutive matrices *)
expedgefactor =  arrowconfig[[5]];
liealgebra = arrowconfig[[6]]; (* Font sizes must be changed at the source-code leve using the Mathematica GUI *)
liegroup = arrowconfig[[7]];

flow =Table[ArrayPlot[MatrixExp[t*x],PixelConstrained->pixs,ColorRules->{y_/;y<(-0.005*range)->Darker[Red,Abs[y/range]],y_/;y>(0.005*range)->Darker[Green,Abs[y/range]]},Frame->False],{t,1, 0,-1/(1*n)}] ;

(*oneparamgroup=GraphicsRow[flow,Spacings->0 ];*)
flow2=Drop[flow,n-m];
distances1 = Table[separation*Norm[MatrixExp[t*x] -MatrixExp[(t-(1/n))*x] ],{t,1, 0,-1/(1*n)}] ;
distances2 = Drop[distances1,n-m-1];
distances3  = Drop[distances2, -1];  (* Last one is bogus since we already hit t = 0 *)

oneparamgroup=GraphicsRow[flow2,{distances3} ];

poly1 = {{-b-cylrad ,-a,0},{-cylrad ,-a,0},{-cylrad ,a,0},{-b-cylrad ,a,0}};
vector1 = {{-b-cylrad,0,0}, {-b-cylrad-arrowheadoffset,0,0}};
vectortext={ {-b-cylrad-arrowheadoffset-arrowheadtextoffset,0,0}};
exptextvector = {{(-b-cylrad)/2, expedgefactor*a, 0}};

(*FIXME: rotation matrix to be multiplied from left*)
poly2=poly1.RotationMatrix[\[Theta],{0,1,0}];
vector2=vector1.RotationMatrix[\[Theta],{0,1,0}];
vectortext2=vectortext.RotationMatrix[\[Theta],{0,1,0}];
exptextvector2=exptextvector.RotationMatrix[\[Theta],{0,1,0}];

path={Graphics3D[{Texture[oneparamgroup],Opacity[.4],Polygon[poly2,VertexTextureCoordinates->{{0,0},{1,0},{1,1},{0,1}}]},Lighting->"Neutral",Boxed->False], Graphics3D[{Red,Arrowheads[arrowheadsize],Arrow[Tube[vector2,0.05]]}],Graphics3D[{Text[Style[Subscript[exptext[[2]],liealgebra], Black, Bold],vectortext2]}], Graphics3D[{Text[Style[Subscript[exptext[[1]],liegroup], Black, Bold],exptextvector2]}]        };

cylheight = 2*a;
cyl=Graphics3D[{Yellow,Text[Style[Subscript["1",liegroup], Black, Bold],{0,a,0}],Cylinder[{{0,-a,0},{0,a,0}},cylrad]},Boxed->False];

{path, cyl}
]



hyperCube[n_]:=Table[{-1,1}, {i, 1, n}]

randomGrid2d[a1_, a2_, n_]:=Table[{RandomReal[a1],RandomReal[a2]} , {i, 1, n}]

randomGrid3d[a1_, a2_, a3_,n_]:=Table[{RandomReal[a1],RandomReal[a2],RandomReal[a3]} , {i, 1, n}]

randomGrid4d[a1_, a2_, a3_,a4_,n_]:=Table[{RandomReal[a1],RandomReal[a2],RandomReal[a3],RandomReal[a4]} , {i, 1, n}]

randomGrid[intervals0_, n0_]:=Module[{intervals=intervals0, n = n0, m},

m = Length[intervals];
Table[
Table[RandomReal[intervals[[j]]], {j, 1, m}]

, {i, 1, n}]

]


(* Approximates the Kolmogorov Complexity function on arbitrary input *)
kolComplexity[x0_]:= Module[{x = x0, xdelim,compxdelim},

(* delimit the input x *)
xdelim = {x, ByteCount[x]};

(* Compress xdelim and count the bytes of the reuslting compression *)
compxdelim=ByteCount[Compress[xdelim]];

(* Return the tuple {K(delimx),|delimx|}*) 
{compxdelim,ByteCount[xdelim]}
]





End[];
EndPackage[];