IRREPS

<<C:\book\programs\BookPrograms\Class.m

Print["Group ",(Subscript[𝒞, 4 v])];
Print["GROUP ORDER = ",g=8];
(*Print["NUMBER OF GENERATORS = ",gen=3];
Print["GROUP GENERATORS: ",L={Range[4],{2,3,4,1},{2,1,4,3}}];*)
Print["GROUP ELEMENTS: ",L];
Print["MULTIPLICATION TABLE: ",m];
Print["INVERSE ELEMENTS: ",LI1];
Print["NUMBER OF CLASSES, nc= ",nc];
Print["CLASSES: ",LC];
Print["CLASS MULTIPLICATION MATRICES"];
Do[
     Print["H(",i,") =  ",MatrixForm[h[[i]]]],{i,2,nc}];



(* 
                                          DIXON'S METHOD 
*)

      (*   LC1 is the list of class exponents. nn is the LCM exponent.  *)

LC1={1};
Do[
     t=LC[[i,1]];ncc=1;While[t≠1,t=m[[1,LC[[i,1]],t]];ncc=ncc+1];AppendTo[LC1,ncc],
     {i,2,nc}
     ];
Print["CLASS EXPONENTS: ",LC1];
nn=LC1[[2]];
Do[
     nn=LCM[nn,LC1[[i]]],{i,3,nc}
     ];
Print["GROUP EXPONENT = ",nn];

    (*    pp is the chosen prime number. Z is the primitive root mod(p).    *)

i=2;
pp=Prime[2];
While[
           Mod[pp,nn]≠1,i++;pp=Prime[i]
           ];
Print["PRIME NUMBER = ",pp];
Z=2;
mm=1;
While[
           mm!=nn,Z++;mm=MultiplicativeOrder[Z,pp]
           ];
Print["PRIMITIVE ROOT = ",Z];
Do[
     zz[i,j]=PowerMod[Z,-1*(i-1)*(j-1),pp],{i,1,nn},{j,1,nn}
     ];
Print["POWERS OF PRIMITIVE ROOT IN THE FORM Z**[(I-1)*(J-1)] : ",MatrixForm[zst=Array[zz,{nn,nn}]]];
Do[
     ccl[i,1]=1;el=LC[[i,1]];ccl[i,2]=i;       
     Do[
         Switch[
                      i==1,True,ccl[i,j]=1,False,el=m[[1,el,LC[[i,1]]]];k=1;
                      While[
                                 FreeQ[LC[[k]],el],k=k+1
                                 ];ccl[i,j]=k
                      ],{j,3,nn}
          ],{i,1,nc}
     ];
Print["CLASS TO WHICH CLASS(I)**J BELONGS : ",MatrixForm[class=Array[ccl,{nc,nn}]]];

(*    GAUSS'S METHOD FOR DETERMINING THE EIGENVECTORS OF THE CLASS MATRICES IN MODULAR ARITHMATIC.    *)

nc2=2*nc
ii=1
q=False
evec={}
evec1={}
evec2={}

(*       nev is the number of non-degenerate eigenvectors found.     *)

nev=0
Do[last[i]=0,{i,1,nc}]


(*       The initial degeneracy is that of the identity class H(1), namely, nc.    *)

last[1]=nc

(*       augm is the augmenting matrix.    *)

Do[Do[augm[i,j]=0;c[i,j]=0,{j,1,nc}];augm[i,i]=1,{i,1,nc}]
Print["AUGMENTING MATRIX, AUGM = ",MatrixForm[cx=Array[augm,{nc,nc}]]]
While[
           And[nev<nc,ii<nc],
          Label[Begin]; ii++;
           b=h[[ii]]; Print["H(",ii,") = ",MatrixForm[b]];

(*      Use the transpose of the H-matrix so that the eigenvectors are row vectors.    *)

           b=Transpose[b];

(*      The ist'th subspace in the eigenspace of H(i) is designated by the label of the first row, ven, and the last row, last[ist]=ind.    *)

           ist=1;
            Label[Try];
           Switch[last[ist]==ist,
                        True,Goto[Fin]
                       ];
           ven=ist;ind=last[ist];
           Print["subspace starts at row = ",ist,", and ends at row = ",ind,", ven = ",ven];
           Do[
               Do[
                   c[i,j]=0;
                   Do[c[i,j]+=augm[i,k]b[[k,j]],{k,1,nc}
                        ],{j,1,nc}
                    ],{i,ist,ind}
               ];
            Print["c = ",MatrixForm[ccx=Array[c,{nc,nc}]]];ev=-1;

(*     Run over ev, (λ=-ev) the trial eigenvalue, from 0 to (prime-1).    *)

            While[
                       And[ev<pp-1,nev<nc,ven<=ind],
                       ev++;

(*     flag will give the degeneracy of the eigenvalue.    *)

                       flag=0;

(*     evd=-λ*augm + augm*H(transpose), λ=-ev, which is just the subspace orthogonal to the known eigenvectors.    *)

Do[evd[i,j]=ev augm[i,j]+c[i,j],{i,1,nc},{j,1,nc}];
                      Print["evd[",ev,"] = ",MatrixForm[xvd=Array[evd,{nc,nc}]]];

(*    Construct the augmented matrix d=[ augm | evd ].    *)

                       Do[
                            Do[
                                 k=j+nc;d[i,j]=augm[i,j];d[i,k]=evd[i,j],{j,1,nc}
                                 ],{i,1,nc}
                          ];

 (*    Set up the augmented matrix, modulo pp, with eigenvalue λ=-ev    *)
        
                       Do[
                           d[i,j]=Mod[d[i,j],pp],
                           {i,1,nc},{j,nc+1,nc2}
                           ];
(*     Print["ev = ",ev,", d = ",MatrixForm[d3=Array[d,{nc,nc2}]]];    *)
                       ix=ven;

(*     Begin the Gauss elimination over the degenerate subspace that starts at row ist, and ends at row ind.    *)

                       Do[
                            jj=nc+1;
                            While[
                                        And[d[i,jj]==0,jj<=nc2],jj=jj+1
                                        ];
                             Switch[jj<=nc2,
                                           True,inv=PowerMod[d[i,jj],-1,pp];
	                                  Do[
                                                d[i,j]=Mod[inv*d[i,j],pp],{j,1,nc2}
                                                ];
	                                  Do[
	                                       dd=pp-d[k,jj];
                                                Do[
                                                     d[k,l]=Mod[dd*d[i,l]+d[k,l],pp],
                                                     {l,1,nc2}
                                                     ],{k,i+1,ind}
                                                 ],
                                           False,flag++;ven++;

(*     If a row of zeros is found in evd, then the corresponding row in augm is an eigenvector.    *)

   Print["THE AUGMENTED MATRIX WITH ZERO ROW IS d3 = ",MatrixForm[d3=Array[d,{nc,nc2}]]];
                                           Do[
                    dv[j]=d[i,j];,{j,1,nc}
                                                 ];
                                            dvd=Array[dv,nc];
                                            Print["evec[ev=",ev,", ven=",ven,", flag=",flag,"] = ",MatrixForm[dvd]];

(*     The list evec1 is a temporary storage for the eigenvectors of the current H    *)

(*     The list evec2 is a temporary storage used to modify the matrix augm to include the non-degenerate eigenvectors.    *)

                                            AppendTo[evec1,dvd];AppendTo[evec2,dvd]
                                            ],{i,ist,ind}
                             ];Switch[flag==0,True,Print["ev = ",ev," is not an eigenvalue"]];

                               last[ix]=ven-1;Print["last[ix=",ix,"] = ",last[ix],
", last[ist=",ist,"] = ",last[ist]];

(*       If eigenvalue is non-degenrate, and if eigenvector, stored in temporary list evec1,
        is not already in permanent list evec, add it to evec. Then empty evec1.   *)

                               Switch[
                                             And[flag==1,FreeQ[evec,evec1]],
                                             True,AppendTo[evec,evec1];nev=nev+1;evec1={},
                                             False,evec1={}
                                             ];

(*       If all nc eigenvectors have been obtained, exit.    *)

                               Switch[
                                             nev==nc,True,q=True
                                             ]
                               ];
            Label[Fin]; Print["Fin, last[ist=",ist,"] = ",last[ist]];

(*        Update the augm matrix so that it contains the newly determined eigenvectors as rows from ist to ind.    *)

            Do[
                  augm[i,j]=evec2[[i,j]],{i,ist,ind},{j,1,nc}
                 ];
            Print["augm = ",MatrixForm[cx=Array[augm,{nc,nc}]]];evec2={}
       Switch[last[ist]==nc,
                     True,i=0;
                    While[i<nc,i++;
                               Switch[last[i]!=i,True,Goto[Begin]]
                               ];qq1=1
                    Print["ok"];Switch[qq1==1,True,Goto[Gd]],
                   False,ist=last[ist]+1;
                   Switch[last[ist]==0,
                                 True,Goto[Panic],
                                 False,Goto[Try]
                                 ]
                   ]
       ]


(* MatrixForm[evec]*)
Label[Gd]
Label[Panic]
dirac=Flatten[evec,1]
Print["first dirac",MatrixForm[dirac[[1]]]]

(* Do[
     Switch[Count[evec,evec[[i]]]==1,False,dirac=Delete[evec,i]]
      ,{i,1,Length[evec]}] *)

Do[Print["dirac[",i,",1] = ",dirac[[i,1]]];inv1=PowerMod[dirac[[i,1]],-1,pp];
    Do[dirac[[i,j]]=Mod[dirac[[i,j]]*inv1,pp] ,{j,1,nc}
         ],{i,1,nc}
    ]
Print["DIRAC CHARACTERS IN MODULO PRIME NUMBER:  ",MatrixForm[dirac]]

                                           (* COMPUTE THE DIMENSIONS OF IRREPS *)
(* d(i)**2=g/sum(ch(j)**2/r(j)), where ch(j) is the dirac character of the j'th class, and r(j) the corresponding class size. *)
dim2={}
Do[ d2=0;
     Do[
          ri=PowerMod[Length[LC[[j]]],-1,pp]; 
          d2=Mod[d2+dirac[[i,j]]^2*ri,pp],{j,1,nc}
          ];
     d2=PowerMod[d2,-1,pp];
     dm=Mod[g*d2,pp];  
     AppendTo[dim2,dm],{i,1,nc}
     ]
Print["SQUARES OF IRREP DIMENSIONS:  ",dim2]
Do[Switch[dim2[[i]]==1,False,j=1;dm=1;While[And[j<pp/2,dim2[[i]]≠dm],j=j+1;dm=Mod[j^2,pp]];dim2[[i]]=j],{i,1,nc}]
Print["IRREP DIMENSIONS, UNORDERED, dim =",dim2]

(* 
                                            COMPUTE CHARACTER MULTIPLICITIES 
*)

cip={}
chit={}
Do[Do[ri=PowerMod[Length[LC[[j]]],-1,pp]; ch=Mod[ri*dim2[[i]]*dirac[[i,j]],pp];AppendTo[chit,ch],{j,1,nc}];AppendTo[cip,chit];chit={},{i,1,nc}]
(*Print[cip]*)
nm1=PowerMod[nn,-1,pp]
Do[mult[i,j,k]=0,{i,1,nc},{j,1,nc},{k,1,nn}]
Do[mmx=0;Do[mmx+=cip[[i,ccl[j,l]]]*zz[k,l];mult[i,j,k]=Mod[mmx*nm1,pp],{l,1,nn}],{i,1,nc},{j,1,nc},{k,1,nn}]
Print["CHARACTER MULTIPLICITIES :",MatrixForm[m3=Array[mult,{nc,nc,nn}]]]

Do[char[i,j]=0,{i,1,nc},{j,1,nc}]
Do[Do[char[i,j]+=mult[i,j,k]*Exp[I(2*π*(k-1)/nn)],{k,1,nn}],{i,1,nc},{j,1,nc}]

(*      SET IRREPS IN ASCENDING ORDER OF THEIR DIMENSIONS    *)

tempo={}

(*      CREATE THE FILE GAMMA FOR THE IRREDUCIBLE REPRESENTATIONS OF THE GROUP ELEMENTS    *)

Γ={}
Do[charac[i,j]=0,{i,1,nc},{j,1,nc}]
Do[ban[i]=0,{i,1,nc}]
Do[dim[i]=0,{i,1,nc}]
dim[1]=1
nrep=1       
                   Do[
                      AppendTo[tempo,1],{j,1,g}
                      ];

(*      The array charac contains the characters of the irreps in ascending order of dimension. 
      The first irrep in charac is the identity.     *)                  

Do[
                    charac[1,j]=1,{j,1,nc}
                    ];

(*      Enter the identity irrep in the file GAMMA.    *)

                  AppendTo[Γ,tempo];tempo={};   

(*      Enter the remaining 1-dimensional irreps in the file GAMMA.    *)

  Do[Switch[
                    And[dim2[[i]]==1,ban[i]==0, Sum[char[i,j],{j,1,nc}]≠nc],True,nrep++;      
                    ban[i]=1;
                   dim[nrep]=dim2[[i]];
                   Do[
                        AppendTo[tempo,char[i,indc[j]]],{j,1,g}];
                  Do[
                     charac[nrep,j]=char[i,j],{j,1,nc}
                     ];
                  AppendTo[Γ,tempo];tempo={}
],{i,1,nc}
     ]
i1=1
i3=nrep+1
dmax=IntegerPart[(Sqrt[g])]
While[
           And[i1<=dmax,Sum[ban[ix],{ix,nc}]<nc,i3<=nc],i1=i1+1;i=1;
          While[
                      And[i≤nc,i3≤nc],
                      Switch[
                                    And[dim2[[i]]==i1,ban[i]==0],True,dim[i3]=dim2[[i]];ban[i]=1;
                                   Do[
                                        charac[i3,j]=char[i,j],{j,1,nc}
                                        ];
                                  i3=i3+1
                                 ];
                       i=i+1
                       ]
]
chi=Array[charac,{nc,nc}]
Print["IRREPS PRINTED IN ASCENDING ORDER OF DIMENSIONALITY:  "]
Print["IRREP DIMENSIONS:  dim = ",dimension=Array[dim,{nc}]]
Print["IRREDUCIBLE CHARACTERS:  "]
Do[Print[StyleForm["χ(",FontWeight->"Bold"],i,"] = ",Simplify[chi[[i]]]],{i,1,nc}]