Wzór Rodriguesa 
http://mathworld.wolfram.com/RodriguesRotationFormula.html

omegam[w1_,w2_,w3_]:={{0,-w3,w2},{w3,0,-w1},{-w2,w1,0}}

MatrixForm[omegam[w1,w2,w3]]

MatrixForm[omegam[w1,w2,w3].omegam[w1,w2,w3]]

Wzór Rodriguesa podaje macierz dla obrotu o kąt phi przeciwnie do ruchu wskazówek 
zegara wokół osi zdefiniowanej przez wektor o długości 1 i składowych kartezjańskich (w1,w2,w3). 
BARDZO WAŻNE: OŚ OBROTU MUSI PRZECHODZIĆ PRZEZ POCZĄTEK UKŁADU WSPÓŁRZĘDNYCH !!!

R[w1_,w2_,w3_,phi_]:=IdentityMatrix[3]+Sin[phi]*omegam[w1,w2,w3]+(1-Cos[phi])*omegam[w1,w2,w3].omegam[w1,w2,w3]

FullSimplify[MatrixForm[R[w1,w2,w3,phi]],Element[w1,Reals]&&Element[w2,Reals]&&Element[w3,Reals]&&Element[phi,Reals]&&(w1*w1+w2*w2+w3*w3==1)]

(* Macierz R powinna być macierzą ortogonalną ! *)

FullSimplify[Inverse[R[w1,w2,w3,phi]]==Transpose[R[w1,w2,w3,phi]],Element[w1,Reals]&&Element[w2,Reals]&&Element[w3,Reals]&&Element[phi,Reals]&&(w1*w1+w2*w2+w3*w3==1)]

(* Wierzchołki kwadratu w płaszczyźnie xy *)

pktA={-1,1,0};

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

pktC={1,-1,0};

pktD={1,1,0};

e=IdentityMatrix[3];

vb={1,0,0};  (* wersor osi x *)

phib=Pi;

b=R[vb[[1]],vb[[2]],vb[[3]],phib] ;

vc={0,0,1};  (* wersor osi z *)

phic=Pi/2;

c=R[vc[[1]],vc[[2]],vc[[3]],phic] ;

(* Zachodzi b^2=c^4=(b.c)^2=e  *)

b.b==e

c.c.c.c==e

b.c.b.c==e

(* Teraz próbuję wypisać wszystkie elementy grupy ... *)

G={e,c,c.c,c.c.c,b,b.c,b.c.c,b.c.c.c};

(* ... i przy pomocy tej podwójnej pętli sprawdzam, czy któryś się nie powtarza.  *)

Do[Do[If[G[[i]]==G[[k]],Print["i= ",i," k=",k]],{i,1,k-1}],{k,2,8}]

(* Wygląda na to, że wszystkie elementy są istotnie różne ! *)

(* W takim razie, który element pokrywa się z c.b ? *)

Do[If[G[[i]]==c.b,Print["i= ",i]],{i,1,8}]

(* A który element jest odwrotny do G[[7]]= b.c.c ? *)

Do[If[(G[[7]].G[[i]]==e)&&(G[[i]].G[[7]]==e),Print["i= ",i]],{i,1,8}]

(* Szukamy centrum grupy D4 *)

Do[{l=0;Do[If[G[[i]].G[[k]]==G[[k]].G[[i]] ,l=l+1],{i,1,8}];If[l==8,Print["k= ",k]]},{k,1,8}]

(* Czyli w centrum są e i c.c ! *)

(* Zbudujmy klasy sprzężenia (a), szukając wszystkich g.a.(g)^(-1). Uwaga: Mathematica używa list, dla których ważna jest pozycja elementu. Dlatego, by móc porównywać różne klasy, sortuję wyniki, a po tej operacji mogę już usunąć elementy, które się powtarzają.   *)

Do[klasa[k]=Sort[DeleteDuplicates[Table[Simplify[G[[i]].G[[k]].Inverse[G[[i]]]],{i,1,8}]]],{k,1,8}]

KLASY=DeleteDuplicates[Table[klasa[i],{i,1,8}]]

(* KLASY jest listą list macierzy (listą klas) ! *)

Length[KLASY]

KLASY[[1]]

(* KLASY[[i]] jest listą macierzy (z klasy nr i) *)

Length[KLASY[[2]]]

KLASY[[2]][[2]]

(* KLASY[[k]][[j]] jest macierzą nr j w klasie nr k *)

Do[Do[Do[If[G[[i]]==KLASY[[k]][[j]],Print["Element i= ",i," jest w klasie nr=",k]],{i,1,8}],{j,1,Length[KLASY[[k]]]}],{k,1,Length[KLASY]}]