1) Set up variables for Koca 0906.2109v2 analysis

While this notebook will use octonions rather than quaternions, the quaternions will be shown to be seamlessly integrated into octonions. As long as we use only{(Subscript[ℯ, 0]),(Subscript[ℯ, 1]),(Subscript[ℯ, 2]),(Subscript[ℯ, 3])} in our computations and select a multiplication table with the upper left quadrant being quaternions with the first triad={1,2,3} where (Subscript[ℯ, 1])∘(Subscript[ℯ, 2])=(Subscript[ℯ, 3]), the rest of the octonion elements {(Subscript[ℯ, 4]),(Subscript[ℯ, 5]),(Subscript[ℯ, 6]),(Subscript[ℯ, 7])} are not referenced within the results.

We select one of 480 possible nontrivial octonion octonion multiplication tables . There are 48 possible with triad={1,2,3} where(Subscript[ℯ, 1])∘(Subscript[ℯ, 2])=(Subscript[ℯ, 3]). There are another 48 with the first triad={1,3,2}  where  (Subscript[ℯ, 2])∘(Subscript[ℯ, 1])=(Subscript[ℯ, 3]), this forms an alternate type of quaternion, which by convention is not generally used.

setFM[87,flip=True,split=0];
triads

fanoPlane
fmDispN

Define τ, σ, and p  as used in the paper in terms of my base MTM symbol φ and show p in symbolic and numeric forms

σ=-φ;
τ=Φ=1/φ;
(**)
Clear[φ,Φ,p]

p=(1/2)(τ+(Subscript[ℯ, 1])+σ (Subscript[ℯ, 3]))
%/.φ->((Sqrt[5])-1)/2.

The octonion power results can be symbolically messy, so it is useful to use numeric values and rounding for comparison .
We also can revert from numeric back to symbolic forms in some cases.

(* chop defines the rounding value *)
chop=(10^(-(*)6**)4(**)));

(* These are my customized round[] functions where rndMat performs rounding on lists of inputs *)
round@in_:=N@If[NumberQ[1.in],Chop[Round[N@in,chop],chop],in];
rndMat@in_List:=round@#&/@round@#&/@in;

round10@in_:=N@If[NumberQ[1.in],Chop[Round[N@in,10chop],chop],in];
rndMat10@in_List:=round10@#&/@round10@#&/@in;

(* This converts from symbolic to rounded numeric octonions *)rndOct@in_:=octonion@FullSimplify[rndMat[oct2List@in/.φ->((Sqrt[5])-1)/2.]];

(* This reverts rounded numeric octonion to Golden Ratio φ symbolic form  *)
octSym@in_:=rndOct@in/.{
N@round10[ (1+(Sqrt[5]))/2.]->τ,N@round10[ -(1+(Sqrt[5]))/2.]->-τ,
N@round[ (1-(Sqrt[5]))/2.]->σ,N@round[ ((Sqrt[5])-1)/2.]->-σ,
N@round10[ (1+(Sqrt[5]))/4.]->τ/2,N@round10[ -(1+(Sqrt[5]))/4.]->-τ/2,
N@round[ (1-(Sqrt[5]))/4.]->σ/2,N@round[ ((Sqrt[5])-1)/4.]->-σ/2,
.5->1/2,-.5->-1/2,
N@round[ 1/(Sqrt[2.])]->1/(Sqrt[2]),N@round[ -1/(Sqrt[2.])]->-1/(Sqrt[2]),
N@round10[ (Sqrt[2.])]->(Sqrt[2]),N@round10[ -(Sqrt[2.])]->-(Sqrt[2])};

Function to τ↔σ in the 48 IcosO where (p^5)=±1.

(* This works for rounding with chop=(10^(-4)) 
after converting the input Quaternion or Octonion to a numeric list *)
switchστ@in_:=rndOct[oct2List@in]/.
 {-0.809->x,0.809->-x,-0.309->0.809,0.309->-0.809}/.x->0.309;

Show the numeric powers of (p^(n=0-5)) & ((OverBar[p])^(n=0-5)) 
Note: I will occasionally use a factor of 2 to simplify output readability

Join[{{"i","(p^i)","N@(p^i)","((OverBar[p])^i)","N@((OverBar[p])^i)"}},
FullSimplify[{#,
2octPower[p,#],2octPower[p,#]/.φ->(1./2) ((Sqrt[5])-1),
2octPower[(SuperStar[p]),#],2octPower[(SuperStar[p]),#]/.φ->(1./2) ((Sqrt[5])-1)}&/@Range[0,5]]]//MatrixForm

2) Demonstrate Quaternion & Octonion JGM-Built-In functions

Octonion Power octPower[] and Conjugation (Subscript[ℯ, n]^( * )) 
Note: we use an asterisk not the MTM conjugate symbol.

octPower[(Subscript[ℯ, 3]),2]

{#,(Subscript[ℯ, 3]^(#)),octPower[(Subscript[ℯ, 3]^( * )),#]}&/@Range[-5,5]//MatrixForm

Need to define a few variables to play with....

Λ0=(*)octonion@**){1,1,1,1,0,0,0,0};
Λ1=(*)octonion@**){1/φ,1,1/φ,1/φ,0,0,0,0};
N@(SuperStar[Λ1])

Λ2=-(*)octonion@**){φ,-1,φ,φ,0,0,0,0};
N@(SuperStar[Λ2])

Λ0=(*)octonion@**){1,1,1,1,0,0,0,0}
Λ1=(*)octonion@**){1/φ,1,1/φ,1/φ,0,0,0,0}
N@(SuperStar[Λ1])
Λ2=-(*)octonion@**){φ,-1,φ,φ,0,0,0,0}
N@(SuperStar[Λ2])

Quaternion to Octonion conversion

FromQuaternion@Quaternion[1,1,1,1]

Quaternion[1,2,0,0]**Quaternion[0,0,1,0]

MTM doesn't handle symbolic Quaternions, so we convert to reals

Λ1[[;;4]]
ToQuaternion[Sequence@Λ1[[;;4]]]
QuaternionQ@%

Λ1[[;;4]]/.φ->(1./2) ((Sqrt[5])-1)
Λ1Quat=Quaternion[Sequence@@%]
QuaternionQ@%

Using the Conjugate symbol for Quaternions and the asterisk for octonions.

FromQuaternion@(Conjugate[Λ1Quat])

...and show the results in the two styles (IJKL vs. (Subscript[ℯ, n]) )

IJKLstyle=True;
octStyle@(SuperStar[quat2oct[Λ1Quat]])

IJKLstyle=False;
octStyle@(SuperStar[quat2oct[Λ1Quat]])

Mixing and converting Quaternion to octonion math is included.
Note: I am using some of my predefined Split Real Even (SRE) E8 and D4 Lie groups here.

Quaternion[Sequence@@pE8[d4[[2]]][[;;4]]]

Λ1Quat**%

quat2oct@%

Show that octonion math with an appropriately selected multiplication table (see above), matches the native Quaternion math.

Λ1Quat∘octonion@pE8[d4[[2]]]

Λ1Quat**Quaternion[Sequence@@pE8[d4[[2]]][[;;4]]]

Octonion functions include integer exponential powers, which simply iterate simple octonion product (octProduct[] or ∘)

octPower[(SuperStar[quat2oct[Λ1Quat]]),#]&/@Range@5//MatrixForm

...with the same results as using native Quaternions.

((Conjugate[Λ1Quat])^(#))&/@Range@5//MatrixForm
FromQuaternion/@%//MatrixForm

3) Define some variables used in Koca papers.

Setup for an earlier 2007 paper:
(HyperlinkTemplate[Group theoretical analysis of 600-cell and 120-cell 4D polytopes with quaternions, https://www.academia.edu/14977836, https://www.academia.edu/14977836, HyperlinkActionRecycled, HyperlinkActive, , ]),

Simple Roots from 2007 paper after eq. 9

α1=octonion@{-1,0,0,0,0,0,0,0};
α2=octonion@{1,1,1,1,0,0,0,0}/2;
α3=octonion@{0,-1,0,0,0,0,0,0};
α4=octonion@{0,1,-σ,-τ,0,0,0,0}/2;
rα={α1,α2,α3,α4};
%//MatrixForm
%/.φ->((Sqrt[5])-1)/2.//MatrixForm

α0=-Total@rα
FullSimplify[%/.φ->((Sqrt[5])-1)/2.]

Setup for 2007 paper eq. 13

SU5=FullSimplify@Total[-rα{4,3,2,1}/5]
rndOct[SU5]
octNorm@%
rndOct[SU5/%]

octonion@{τ-σ,0,-τ,σ,0,0,0,0}/(Sqrt[8])
rndOct@%
octNorm@rndOct@%

rndOct[SU5/octNorm@SU5]==rndOct[octonion@{τ-σ,0,-τ,σ,0,0,0,0}/(Sqrt[8])]

Setup for 2007 paper eq. 15

(* c=octonion@{0,0,1,-1,0,0,0,0}/(Sqrt[2]) *)
cq=Quaternion[Sequence@@(#/(Sqrt[2.])&/@{0,0,1,-1})]
QuaternionQ@%
c=octonion@Chop@oct2List@quat2oct@%%
octSym@c

(* α=octonion@{-σ,1,τ,0,0,0,0,0}/2 *)
αQ=Quaternion[Sequence@@(#/2.&/@{-σ,1,τ,0}/.φ->((Sqrt[5])-1)/2.)]
QuaternionQ@%
α=quat2oct@%%
octSym@α

(* β=octonion@{-τ,1,σ,0,0,0,0,0}/2 *)
βq=Quaternion[Sequence@@(#/2.&/@{-τ,1,σ,0}/.φ->((Sqrt[5])-1)/2.)]
QuaternionQ@%
β=quat2oct@%%
octSym@β

(* This shows an error in an earlier 2007 paper (HyperlinkTemplate[Group theoretical analysis of 600-cell and 120-cell 4D polytopes with quaternions, https://www.academia.edu/14977836, https://www.academia.edu/14977836, HyperlinkActionRecycled, HyperlinkActive, , ]) ,
this occurs between eq. 9 and eq. 10 and suggests β=-c∘((SuperDagger[(OverBar[α])])∘c) , but  (Subscript[ℯ, 2])-> (Subscript[ℯ, 3]),
this may be a typo in the defintion for β *)
rndOct@β==rndOct[-c∘(switchστ@(SuperStar[α])∘c)]

(* Same result using non-commutative multiply in Quaternions *)
IJKLstyle=True;
-cq**ToQuaternion[octStyle@switchστ@(SuperStar[α])/.{𝕛 ->J,𝕜 ->K}]**cq
quat2oct@%
IJKLstyle=False;

(* Assuming a typo, redefine β=octonion@{-τ,1,0,σ,0,0,0,0}/2 *)
βq=Quaternion[Sequence@@(#/2.&/@{-τ,1,0,σ}/.φ->((Sqrt[5])-1)/2.)]
QuaternionQ@%
β=quat2oct@%%
octSym@β

rndOct@β==rndOct[-c∘(switchστ@(SuperStar[α])∘c)]

Setup for the 2009 paper analysis.

Define t1∘t2 with p being defined above in Section 1 to be used below in Section 7.

t1=(1+(Subscript[ℯ, 1])-(Subscript[ℯ, 2])-(Subscript[ℯ, 3]))/2;
t2=(1+(Subscript[ℯ, 1])+(Subscript[ℯ, 2])-(Subscript[ℯ, 3]))/2;
t1∘t2

FullSimplify[octPower[p,2]∘((SuperStar[t2])∘(SuperStar[t1]))]
FullSimplify[%/.φ-> ((Sqrt[5])-1)/2.]

(* Simple Roots (*)
α1=octonion@{0,0,1,-1,0,0,0,0};
α2=octonion@{0,1,0,1,0,0,0,0};
α3=octonion@{0,0,-1,-1,0,0,0,0};
α4=octonion@{1,-1,0,0,0,0,0,0};
rα={α1,α2,α3,α4};
%//MatrixForm

Weights

ω1=octonion@{1,1,1,-1,0,0,0,0}/2
ω2=(SuperStar[α4])
ω3=octonion@{1,1,-1,-1,0,0,0,0}/2
ω4=octonion@{1,0,0,0,0,0,0,0}
rω={ω1,ω2,ω3,ω4};

4) Define 120 vertex H4Φ Icosians (I in Table 1) using Quaternion(IcosQ) & Octonion(IcosO) forms.

This uses h4Φ from my E8->H4 folding matrix model.

IcosQ=Quaternion[Sequence@@N[pC600[#][[;;4]]]]&/@h4Φ;
IcosO=quat2oct@#&/@IcosQ;
{#,FullSimplify@FromQuaternion@IcosQ[[#]],
IcosO[[#]]}&/@Range@Length@IcosO//MatrixForm

rndIcosO=rndOct[FullSimplify[#/2]]&/@IcosO;

5) Define the Unit Quaternions using my Octonion functions (eq. 4)

F4(48)=T,which is made up of 1 component points (Subscript[V, 0])=±(Subscript[ℯ, j]) 
and 4 component elements (SubPlus[V])=Vp & (SubMinus[V])=Vm
   ⊕ T'=Tp=(Sqrt[2])V123 from {V1,V2,V3},
             which are made up of 2 component vectors ±(Subscript[ℯ, 0])±(Subscript[ℯ, i]) or 2 component spinor/anti-spinors ±(Subscript[ℯ, i])±(Subscript[ℯ, j]))

These are shown in symbolic & rounded numeric forms w/Norms.

Vp=Sort[octonion@Join[#,{0,0,0,0}]&/@Eperms4[1,1,1,1]]/2;
Length@%
%%//ColumnForm

Vm=Sort[octonion@Join[#,{0,0,0,0}]&/@Operms4[1,1,1,1]]/2;
Length@%
%%//ColumnForm

V0=Sort[octonion@pE8@#&/@excluded[[Join[Range[5,8],Range[13,16]]]]];
Length@%
%%

V123=Sort[octonion@pE8@#&/@d4]/2;
Length@%
%%
Total[{%%,%%%%%,%%%%%%%%,%%%%%%%%%%%}]

T=Join[V0,Vp,Vm]
Length@%

...and show the Octonion Norms octNorm[]=1.

octNorm@#&/@%

Tp=(Sqrt[2])V123
octNorm@#&/@%

v1={1,2,3,4,17,18,23,24};
v2={5,6,9,10,15,16,21,22};
v3={7,8,11,12,13,14,19,20};
Length@Union@Flatten@Join[v1,v2,v3]
{V1,V2,V3}={V123[[v1]],V123[[v2]],V123[[v3]]};
Flatten[2{V1,V2,V3}]//MatrixForm

(* Define a String version of T & T' for future formatting of outputs *)
TStr=ToString[2#,TraditionalForm]&/@T
TpStr=ToString[2#,TraditionalForm]&/@Tp

6)Apply σ scaling factor to T&V123 as (T,0)⊕(0,T)⊕((Subscript[V, 1]),(Subscript[V, 3]))⊕((Subscript[V, 2]),(Subscript[V, 1]))⊕((Subscript[V, 3]),(Subscript[V, 2]))(eq.5)

Where (A, B) = A + φB (vs . from E8 -> H4 folding) .
      This creates 72 of 120 H4 φ scaled elements of 72 = 2×T (24) + T' (24) or 96 if you add the 24 from (-T') . 
      This is not the 240 E8 -> H4 + H4φ as suggested but ...  
       using the predefined 120 Icosians (I, I) = I + φI it does trivially create a folded E8 -> H4 + H4φ set.

These are shown in symbolic & rounded numeric forms w/Norms.

Tgroup=2Join[T,σ T];
Length@%
Tgroup//MatrixForm
{out=rndOct@#,round[.5 octNorm@out]}&/@%%%;
%//MatrixForm;
TgroupInt=FullSimplify[.5%%[[All,1]]];
%//MatrixForm

TPgroup=2Join[Tp,σ Tp];
Length@%
TPgroup//MatrixForm
{out=rndOct@#,round[.5 octNorm@out]}&/@%%%;
%//MatrixForm;
TPgroupInt=FullSimplify[.5%%[[All,1]]];
%//MatrixForm

VgroupP=Flatten@Join[FullSimplify[{
(* On the (V1,V3) pairing for V3 Norms to 1 not φ unless negated *)
2V1[[#]]+If[#==2||#==3,-2,2]σ V3[[#]],
2V2[[#]]+2σ V1[[#]],
2V3[[#]]+2σ V2[[#]]}]&/@Range@8];
Length@%
%%//MatrixForm
{out=rndOct@#,round[.5 octNorm@out]}&/@%%%;
%//MatrixForm
VgroupPInt=FullSimplify[.5%%[[All,1]]];
%//MatrixForm;

VgroupM=Flatten@Join[FullSimplify[{
(* On the (V1,V3) pairing for V3 Norms to φ not 1 unless negated *)
2V1[[#]]-If[#==2||#==3,-2,2]σ V3[[#]],
2V2[[#]]-2σ V1[[#]],
2V3[[#]]-2σ V2[[#]]}]&/@Range@8];
Length@%
%%//MatrixForm
{out=rndOct@#,round[.5 octNorm@out]}&/@%%%;
%//MatrixForm
VgroupMInt=FullSimplify[.5%%[[All,1]]];
%//MatrixForm;

Vgroup=Flatten@Join[VgroupP,VgroupM]
Length@Union@%
%%//MatrixForm
{out=rndOct@#,round[.5 octNorm@out]}&/@%%%;
%//MatrixForm
VgroupInt=FullSimplify[.5%%[[All,1]]];
%//MatrixForm;

Flatten@Join[VgroupM,Tgroup,VgroupP,TPgroup]
Length@Union@%
{out=rndOct@#,round[.5 octNorm@out]}&/@%%;
%//MatrixForm
VTgroupInt=FullSimplify[.5%%[[All,1]]];
%//MatrixForm;

rndS=SgroupInt=Sort@Complement[rndIcosO,TgroupInt];
Length@%

7) Define the 12 vertex Icosahedron set (IcosSet) using a specific combinations of {p,t1,t2}

Define the base Icosahedron set functions in string form (for formatting output).

IcosSet=2FullSimplify[{
p,(SuperStar[p]),
p∘t1,p∘((SuperStar[t2])∘t1),
(SuperStar[p])∘t1,(SuperStar[p])∘octPower[t1,2],(SuperStar[p])∘(t1∘t2),
octPower[p,2],octPower[(SuperStar[p]),2]∘octPower[t1,2],octPower[(SuperStar[p]),2]∘(t2∘t1),octPower[(SuperStar[p]),2]∘((SuperStar[t2])∘t1),octPower[(SuperStar[p]),2]∘(t1∘(t2∘t1))
}];
IcosSet//MatrixForm

IcosSetStr={"p","\!\(\*FormBox[OverscriptBox[\"p\", \"_\"], TraditionalForm]\)","\!\(\*FormBox[RowBox[{\"p\", \"∘\", \"t1\"}], TraditionalForm]\)","\!\(\*FormBox[RowBox[{\"p\", \"∘\", RowBox[{\"(\", RowBox[{OverscriptBox[\"t2\", \"_\"], \"∘\", \"t1\"}], \")\"}]}], TraditionalForm]\)","\!\(\*FormBox[RowBox[{OverscriptBox[\"p\", \"_\"], \"∘\", \"t1\"}], TraditionalForm]\)","\!\(\*FormBox[RowBox[{OverscriptBox[\"p\", \"_\"], \"∘\", SuperscriptBox[\"t1\", \"2\"]}], TraditionalForm]\)","\!\(\*FormBox[RowBox[{OverscriptBox[\"p\", \"_\"], \"∘\", RowBox[{\"(\", RowBox[{\"t1\", \"∘\", \"t2\"}], \")\"}]}], TraditionalForm]\)","\!\(\*FormBox[SuperscriptBox[\"p\", \"2\"], TraditionalForm]\)","\!\(\*FormBox[RowBox[{SuperscriptBox[OverscriptBox[\"p\", \"_\"], \"2\"], \"∘\", SuperscriptBox[\"t1\", \"2\"]}], TraditionalForm]\)","\!\(\*FormBox[RowBox[{SuperscriptBox[OverscriptBox[\"p\", \"_\"], \"2\"], \"∘\", RowBox[{\"(\", RowBox[{\"t2\", \"∘\", \"t1\"}], \")\"}]}], TraditionalForm]\)","\!\(\*FormBox[RowBox[{SuperscriptBox[OverscriptBox[\"p\", \"_\"], \"2\"], \"∘\", RowBox[{\"(\", RowBox[{OverscriptBox[\"t2\", \"_\"], \"∘\", \"t1\"}], \")\"}]}], TraditionalForm]\)","\!\(\*FormBox[RowBox[{SuperscriptBox[OverscriptBox[\"p\", \"_\"], \"2\"], \"∘\", RowBox[{\"(\", RowBox[{\"t1\", \"∘\", RowBox[{\"(\", RowBox[{\"t2\", \"∘\", \"t1\"}], \")\"}]}], \")\"}]}], TraditionalForm]\)"}

Show the numeric and a final reduced symbolic form in terms of φ only as well as verify the octonion Norm=1.

rndOct/@IcosSet//MatrixForm
octSym/@%//MatrixForm

8)Find the 48 p∈I where (p^5)=±1 and show ((OverBar[p])^1)==±(p^4) && ((OverBar[p])^4)==±p , ((OverBar[p])^2)==(p^3) && ((OverBar[p])^3)==(p^2) (eq.8)

Show ((OverBar[p])^1)==±(p^4) && ((OverBar[p])^4)==±p , ((OverBar[p])^2)==(p^3) && ((OverBar[p])^3)==(p^2) when ((OverBar[p])^5)==±(p^5)==±((OverBar[p])^0)==±(p^0)==±1.

out={#,
(* 120 H4Φ Icosians (I[[#]]=p) *)
IcosO[[#]],
(* (P^5) *)
rndOct[2octPower[ IcosO[[#]]/2,5]],
(*  ((OverBar[p])^1)==±(p^4) && ((OverBar[p])^4)==±p *)
(rndOct@octPower[ (SuperStar[IcosO[[#]]])/2,1]===rndOct@octPower[ IcosO[[#]]/2,4]||
rndOct@octPower[ (SuperStar[IcosO[[#]]])/2,1]===-rndOct@octPower[ IcosO[[#]]/2,4])&&
(rndOct@octPower[ (SuperStar[IcosO[[#]]])/2,4]===rndOct@octPower[ IcosO[[#]]/2,1]||
rndOct@octPower[ (SuperStar[IcosO[[#]]])/2,4]===-rndOct@octPower[ IcosO[[#]]/2,1]),
(* ((OverBar[p])^2)==(p^3) && ((OverBar[p])^3)==(p^2) *)
(rndOct@octPower[ (SuperStar[IcosO[[#]]])/2,2]===rndOct@octPower[ IcosO[[#]]/2,3]||
rndOct@octPower[ (SuperStar[IcosO[[#]]])/2,2]===-rndOct@octPower[ IcosO[[#]]/2,3])&&
(rndOct@octPower[ (SuperStar[IcosO[[#]]])/2,3]===rndOct@octPower[ IcosO[[#]]/2,2]||
rndOct@octPower[ (SuperStar[IcosO[[#]]])/2,3]===-rndOct@octPower[ IcosO[[#]]/2,2])
}&/@Range@Length@IcosO;
Join[{{"#","120 H4Φ Icosians (I[[#]]=p)","(P^5)",Column[{"((OverBar[p])^1)=±(p^4)","((OverBar[p])^4)=±p"}],Column[{"((OverBar[p])^2)=(p^3)","((OverBar[p])^3)=(p^2)"}]}},
out]//MatrixForm
p50=Flatten@Position[out[[All,4]],True]
Length@%

Trim the ±1 from p∈I where (p^5)=±1 (which fail the power conjugate test).

(* Position in I *)
drop2I=Select[out,#[[2]]==2.||#[[2]]==-2.&][[All,1]]

(* Position in p50 *)
drop2p=Flatten[Position[p50,#]&/@drop2I,1]

p48=Delete[p50,drop2p]
Length@%

9)Create  I(120) using  T (24) and (p^n), where p∈I && (p^5)=±1  with S(96)=I-T.

We are going to perform the I(120) process for each of 50 p∈I && (p^5)=±1  (including the two trivial cases) 
  in order to demonstrate the issue between p50 vs. p48.
T(24)×1((p^5)=±1)×5((p^(n=0-4)))=120 permutations with no dups.
By definition (p^(j=0))T=1.T, so the (p^(n=1-4))T will be I-T=S(96).

in1=rndMat10@oct2List[#][[;;4]]&/@Sort@SgroupInt;
Join[{{"indx","Generated","Unique",Column[{"Naive","Match","Count","__________"},Center]}},
Table[
outTbl=TableForm@
Table[Row@{inp,"  ",inT,"  ",
MatrixForm[{#,out=rndOct[octPower[inp,#]∘inT],
(* Powers (p^(n=1-4)) or (p^(n=0-4)) *)
MemberQ[rndIcosO,out]}]&/@Range[1,4]},
(* These are individual elements of the 50 p∈I where (p^5)=±1 *)
{inp,rndIcosO[[{p50[[indx]]}]]},
(* Only the first 24 are the original T Norm=1, the rest are scaled by σ *)
{inT,TgroupInt[[;;24]]}];

outTblFlat=Flatten[outTbl[[1,All,All,1,5,All,1,2]],2];
Iout=DeleteDuplicates@outTblFlat;
{indx,
(* These Length's will be 96 when using Powers (p^(n=1-4)) or 120 for (p^(n=0-4)) *)
Length@outTblFlat,
Length@Iout,
(* This will be True when using Powers (p^(n=1-4)) on all but the 2 trivial (p^5)=±1 cases *)
in2=rndMat10@oct2List[#][[;;4]]&/@Sort@Iout;
Count[in1[[#]]==in2[[#]]&/@Range@Length@Iout,True]},
{indx,Length@p50}]]//MatrixForm

Notice the positions 7 & 44 have 72 duplicates leaving only 24 vertices matching. 
This coincides with drop2p = {{7}, {44}} obtained above .

10) Define and test the [p,q]:r function a prq[] (eq.6)

[p,q]:r is defined for combinations of inputs as elements or lists of equal length

prq[p_,r_,q_]:=FullSimplify[p∘(r∘q)];
prq[p_List,r_,q_]:=prq[p[[#]],r,q]&/@Range@Length@p;prq[p_,r_List,q_]:=prq[p,r[[#]],q]&/@Range@Length@r;prq[p_,r_,q_List]:=prq[p,r,q[[#]]]&/@Range@Length@q;prq[p_List,r_,q_List]:=prq[p[[#]],r,q[[#]]]&/@Range@Length@p;prq[p_List,r_List,q_List]:=prq[p[[#]],r[[#]],q[[#]]]&/@Range@Length@p;

Show Power Conjugations (Subscript[p, c]^j)=q∘(p^(j=0-4))∘(OverBar[q]) where p&q∈I using [p,q]:r->[q,(OverBar[q])]:(p^j)

(* Select a random element q=I[[12]] *)
inq=rndIcosO[[p48[[12]]]];
TableForm@Table[Row@{in,"  ",
MatrixForm[{#,out=rndOct[prq[inq,octExp[in,#],(SuperStar[inq])]],
MemberQ[rndIcosO,out]}]&/@Range[0,4]},
(* These are the p∈I where (p^5)=±1 *)
{in,rndIcosO[[p48]]}]

11) Create  600 vertices (J) of the 120 Cell ( the dual of I(120)=600 Cell) with  (∑_(i=0)^4)(p^i)∘((SuperDagger[((OverBar[p])^i)])∘I)  (eq.15).

For reference, these are the “canonical” 600 Cell values based on permutations of Φ &  φ based on the Cartesian coordinates of shown in Wikipedia.

cell120base={
{2 dn[4][[All,;;4]]},
perms[(Sqrt[5]),1,1,1],
perms[Φ,Φ,Φ,(φ^2)],
perms[(Φ^2),φ,φ,φ],
Union@Eperms[(Φ^2),(φ^2),1,0],
Union@Eperms[(Sqrt[5]),φ,Φ,0],
Eperms[2,1,Φ,φ]};
Cell120=DeleteDuplicates@Flatten[cell120base,2];
(* This adds 4 zero dimensions for use as Octonions *)
cell120=Join[#,{0,0,0,0}]&/@Cell120;

Round the 120 Cell vertices for future comparison of values to be generated . These have a 4D Norm=(Sqrt[8]).

rndCell120=Sort@rndMat10[Cell120/.{φ->((Sqrt[5])-1)/2.,Φ->((Sqrt[5])+1)/2.}];
rndScaledCell120=rndMat10[rndCell120/(Sqrt[8.])];
Length@%

(* Show the Norms are 2(Sqrt[2])=2.828...` *)
Union[rndMat10[Norm@#&/@rndCell120]]=={2.828`}

Visualize the canonical 120 Cell vertices.

ConvexHullMesh@rndCell120[[All,2;;]]

From the 2007 paper, validate eq. 15 (∑_(a=0)^4)(α^a)∘(c∘(β^a))==(∑_(a=0)^4)(α^a)∘(((OverBar[α])^(†a))∘c)

Join[{{"a","(α^a)∘(c∘(β^a))","(α^a)∘((SuperDagger[((OverBar[α])^a)])∘c)","Match?"}},
Table[{a,
o1=FullSimplify@rndMat@oct2List[prq[octPower[α,a],c,octPower[β,a]]];
octonion@o1,
o2=FullSimplify@rndMat@oct2List[prq[octPower[α,a],octPower[switchστ@(SuperStar[α]),a],c]];octonion@o2,
rndMat10@Abs[o1-o2]==Array[0.&,8]},
{a,0,4}]]//MatrixForm

From the 2007 paper, attempt to validate 120 cell {5,3,3} (∑_(a=0)^4)⊕I∘((α^a)∘(((OverBar[α])^(†a))∘c))

out1=Flatten[Join[{{{"a","I∘((α^a)∘(((OverBar[α])^(†a))∘c))","octNorm",
Column[{"Naive","Canonical","Cell120","Match","__________"},Center]}}},
Table[{a,
o1=FullSimplify@rndMat@oct2List[inI∘prq[octPower[α,a],octPower[switchστ@(SuperStar[α]),a],c]];
octonion@o1,
round10@octNorm@o1,
(* Determine if the vertex matches one of the rounded and scaled canonical values *)
MemberQ[rndScaledCell120,rndMat10@o1[[;;4]]]},
{a,0,4},
{inI,rndIcosO}]],1];
out1//MatrixForm
Length@out1-1

cell120Gen1=Sort[rndMat10[oct2List@#]&/@out1[[2;;,2]]];
cell120Gen1//MatrixForm;

From the 2007 paper, attempt to validate 120 cell {5,3,3}(∑_(a,b=0)^4)⊕(α^a)∘(T'∘(β^b))

(**)
Flatten[Join[{{"(α^a)∘(T'∘(β^b))"}},
out2=Table[Row@{"a=",a,",b=",b,
rndOct/@prq[octPower[α,a],TPgroupInt[[;;24]],octPower[β,b]]//MatrixForm},
{a,0,4},
{b,0,4}]],1]//MatrixForm;

cell120Gen2=Sort[rndMat10[oct2List@#]&/@Flatten@out2[[All,All,1,5,1]]];
Length@%
cell120Gen2//MatrixForm;

Show these are equivalent to each other.

A naive Sort of the 4D quaternion (or actually 8D octonion) vertices show mismatches.

out3={#,cell120Gen1[[#]],cell120Gen2[[#]],
diff=Chop[Abs[
rndMat@oct2List[cell120Gen1[[#]]]-
rndMat@oct2List[cell120Gen2[[#]]]],10chop],
0==Norm@diff}&/@Range@600;
%//MatrixForm

Evaluating the Norm vs. the Abs difference being  0 still miss 36 vertices.

outSel=Select[out3,#[[5]]==True&];
Length@outSel
(* Evaluating the Total Abs difference being within ±.01 of 0 *)
outSel=Select[out3,Total@#[[4]]<.01&][[All,1]];
Length@outSel
(* Matches vs. mismatches *)
600-%

The 36 mismatches are due to Sort errors caused by rounding errors in (Subscript[ℯ, 0])

outfail=out3[[Complement[Range@600,outSel]]];
%//MatrixForm

Match up the mismatching vertices by resorting them to find the nearest 4D Quaternion (or 8D octonion) distanced vertex, verifying they are all the same vertex within ±.01 of 0

failmatch={#,outfail[[#,1]],in=cell120Gen1[[outfail[[#,1]]]],
out=Sort[Table[{Abs[in-cell120Gen2[[outfail[[i,1]]]]],i},{i,Length@outfail}]][[1,2]],
cell120Gen2[[outfail[[out,1]]]],
Abs[in-cell120Gen2[[outfail[[out,1]]]]],
Total@Abs[in-cell120Gen2[[outfail[[out,1]]]]]}&/@Range@Length@outfail;
%//MatrixForm

Show that these also match the rounded and scaled canonical 120 Cell vertices.

The naive MemberQ matching, 464 of the 600 quaternion generated vertices match the canonical values.

Length@out1-1
Count[out1[[All,4]],True]

Match up the 136 unmatched generated 120 Cell  with canonical 120 Cell  vertices by resorting to find the nearest 4D Quaternion (or 8D octonion) distanced vertex, verifying they are all the same vertex within ±.01 of 0

not120Cell=Sort@Select[out1,#[[4]]==False&];
outfail1={#,rndMat10[oct2List@not120Cell[[#,2]]][[;;4]]}&/@Range@Length@not120Cell;

failmatch1={#,in=outfail1[[#,2]],
out=Sort[Table[{Abs[in-rndScaledCell120[[i]]],i},{i,600}]][[1,2]],
rndScaledCell120[[out]],
Abs[in-rndScaledCell120[[out]]],
Total@Abs[in-rndScaledCell120[[out]]]}&/@Range@Length@outfail1;
%//MatrixForm

Validate the 120 Cell with 2009 paper (∑_(i,j=0)^4)(p^i)∘(T'∘((OverBar[q])^j)) (eq.15)

(* Define the input as I(120) *)
inList=TPgroupInt[[;;24]];
(* Define the 48 p∈I where (p^5)=±1 *)
inqList=rndIcosO[[p48]];
(* Set ii as the max limit of i *)
ii=4;
(* Generate the annotated outputs *)
Join[{{"indx","Generated","Unique","Matching"}},
Table[
outTbl=TableForm@
Table[
(* c∈T'=TPgroupInt *)
inc=inList[[inP]];
(* q∈I where (q^5)=±1 or IcosO[[outp48]] *)
inq=inqList[[indx]];
(* p=(OverBar[c])∘((SuperDagger[(OverBar[q])])∘c) *)
inp=FullSimplify[(SuperStar[inc])∘(switchστ[(SuperStar[inq])]∘inc)];
Row@{MatrixForm@{inc,inq,inp},
Table[
(* (∑_(i,j=0)^4)(p^i)∘(T'∘((OverBar[q])^j)) *)
MatrixForm[{inTp,
i,octPower[inp,i],
j,inTp∘octPower[(SuperStar[inq]),j],
out=rndOct@prq[octPower[inp,i],inTp,octPower[(SuperStar[inq]),j]],
rndMat@oct2List[out][[;;4]],
(* Determine if the vertex matches one of the rounded and scaled canonical values *)
MemberQ[rndScaledCell120,rndMat10@oct2List[out][[;;4]]]}],
{i,0,ii},
{j,0,4}]},
(* 
inc=inList[[inP]]
{18}=Flatten@Position[rndMat10@oct2List@#&/@inList,rndMat10@oct2List@c],
 *)
{inP,Flatten@Position[rndMat10@oct2List@#&/@inList,rndMat10@oct2List@c]},
(* T' *)
{inTp,inList}];
outTblFlat=Flatten[outTbl[[1,All,All,1,2,All,All,1,6]],3];
Iout=DeleteDuplicates@outTblFlat;
{indx,
Length@outTblFlat,
Length@Iout,
Count[MemberQ[rndScaledCell120,rndMat10@oct2List[Iout[[#]]][[;;4]]]&/@Range@Length@Iout,True]},
(* Evaluate for indx up to inqList=rndIcosO[[p48]] or
{39}=Flatten@Position[rndMat10@oct2List@#&/@inqList,rndMat10@oct2List@β]
*)
{indx,Flatten@Position[rndMat10@oct2List@#&/@inqList,rndMat10@oct2List@β]}]]//MatrixForm

ConvexHullMesh[oct2List[#][[;;3]]&/@outTblFlat]

12) Create and visualize the 96 vertices of  S'=(∑_(i=1)^4)(p^i)∘((SuperDagger[((OverBar[p])^i)])∘T') , 
along with M=M1⊕M2 where M1=(∑_(i=1)^4)(p^i)∘T' and M2=(∑_(i=1)^4)(SuperDagger[((OverBar[p])^i)])∘T' .

Create S'=(∑_(i=1)^4)(p^i)∘((SuperDagger[((OverBar[p])^i)])∘T')

SPgroupInt=Flatten[Table[
rndOct/@prq[octPower[α,a],octPower[switchστ@(SuperStar[α]),a],TPgroupInt[[;;24]]],
{a,1,4}],1];
Length@SPgroupInt
SPgroupInt;

Visualize S'

ConvexHullMesh[oct2List[#][[;;3]]&/@SPgroupInt]

Tally and visualize S’ 3D hull Norms using varying Opacity

hulls3D[oct2List[#][[;;3]]&/@SPgroupInt]

Create M1=(∑_(i=1)^4)(p^i)∘T'

M1groupInt=Flatten[Table[
rndOct[octPower[α,a]∘inTp],
{a,1,4},
{inTp,inList}],1];
Length@M1groupInt
M1groupInt;

Visualize M1

ConvexHullMesh[oct2List[#][[;;3]]&/@M1groupInt]

Tally and visualize M1 3D hull Norms using varying Opacity

hulls3D[oct2List[#][[;;3]]&/@M1groupInt]

Create M2=(∑_(i=1)^4)(SuperDagger[((OverBar[p])^i)])∘T'

M2groupInt=Flatten[Table[
rndOct[octPower[switchστ@(SuperStar[α]),a]∘inTp],
{a,1,4},
{inTp,inList}],1];
Length@M2groupInt
M2groupInt;
MgroupInt=Join[M1groupInt,M2groupInt];

Visualize M2

ConvexHullMesh[oct2List[#][[;;3]]&/@M2groupInt]

Visualize M(192)

Length@Union@MgroupInt
ConvexHullMesh[oct2List[#][[;;3]]&/@MgroupInt]

Tally and visualize M(192) 3D hull Norms using varying Opacity

hulls3D[oct2List[#][[;;3]]&/@MgroupInt]

13) Visualize the  144 vertices of the dual snub 24 Cell =T⊕T'⊕ S'=(∑_(i=1)^4)(p^i)∘((SuperDagger[((OverBar[p])^i)])∘T') .

Visualize T⊕T'⊕ S'

dualSnub24Cell=oct2List[#][[;;4]]&/@Join[
TgroupInt[[;;24]],
TPgroupInt[[;;24]],
SPgroupInt];
Length@%

(* Visualize the dualSnub24Cell *)
ConvexHullMesh@dualSnub24Cell[[All,;;3]]

Tally and visualize dualSnub24Cell 3D hull Norms using varying Opacity

hulls3D[dualSnub24Cell[[All,;;3]]]

An alternate version of  S'=(∑_(i=1)^4)(p^i)∘(T'∘(SuperDagger[((OverBar[p])^i)]))

SPgroupInt2=Flatten[Table[
rndOct/@prq[octPower[α,a],TPgroupInt[[;;24]],octPower[switchστ@(SuperStar[α]),a]],
{a,1,4}],1];

ConvexHullMesh[oct2List[#][[;;3]]&/@SPgroupInt2]

Tally and visualize dualSnub24Cell 3D hull Norms using varying Opacity

hulls3D[oct2List[#][[;;3]]&/@SPgroupInt2]