[a b]   [E F]   [aE+bG aF+bH]
[c d] * [G H] = [cE+dG cF+dH]


---------------------------------------
X:
	y¹ = y⁰*cosR⁰ - z⁰*sinR⁰
	z¹ = y⁰*sinR⁰ + z⁰*cosR⁰
	x¹ = x⁰

	x⁰	y⁰	z⁰
x¹	[1	0	0	]
y¹	[0	cosR⁰	-sinR⁰	]
z¹	[0	sinR⁰	cosR⁰	]
----------------------------------------
Y:
	z² = z¹*cosR¹ - x¹*sinR¹
	x² = z¹*sinR¹ + x¹*cosR¹
	y² = y¹

	x¹	y¹	z¹
x²	[cosR¹	0	sinR¹	]
y²	[0	1	0	]
z²	[-sinR¹	0	cosR¹	]
---------------------------------------
Z:
	x³ = x²*cosR² - y²*sinR²
	y³ = x²*sinR² + y²*cosR²
	z³ = z²

	x²	y²	z²
x³	[cosR²	-sinR²	0	]
y³	[sinR²	cosR²	0	]
z³	[0	0	1	]
---------------------------------------
R = X * Y * Z

x	[cosR¹*cosR²	cosR¹*-sinR²	sinR¹*1]
y	[-sinR⁰*-sinR¹*cosR²+cosR⁰*sinR²	-sinR⁰*-sinR¹*-sinR²+cosR⁰*cosR²	-sinR⁰*cosR¹*1]
z	[cosR⁰*-sinR¹*cosR²+sinR⁰*sinR²	cosR⁰*-sinR¹*-sinR²+sinR⁰*cosR²	cosR⁰*cosR¹*1]

rotation around three axis , if you want to rotate around arbitrary axis , youll need to rotate back , rotate , then re rotate , about three big rotation

---------------------------------------
euler
XYZ	unit vector of arbitrary axis
		x⁰	y⁰	z⁰
x¹	[	(1-cosR)*X*X+cosR	(1-cosR)*X*Y-sinR*Z	(1-cosR)*X*Z+sinR*Y	0	]
y¹	[	(1-cosR)*X*Y-sinR*Z	(1-cosR)*Y*Y+cosR	(1-cosR)*Y*Z-sinR*X	0	]
z¹	[	(1-cosR)*X*Z-sinR*Y	(1-cosR)*Y*Z+sinR*X	(1-cosR)*Z*Z+cosR	0	]
	[	0			0			0			1	]

c = cosR
s = sinR
t = 1-c
tx = t*X
ty = t*Y
tz = t*Z
txy = tx*Y
txz = tx*Z
tyz = ty*Z
sx = s*X
sy = s*Y
sz = s*Z
		x⁰	y⁰	z⁰
x¹	[	tx*X+c	txy-sz	txz+sy	0	]
y¹	[	txy+sz	ty*Y-c	tyz-sx	0	]
z¹	[	txz-sy	tyz+sx	tz*z+c	0	]
	[	0	0	0	1	]

--------------------------------------------------
quaternion		[	w 	x 	y 	z	]
magnitude		sqrt(w^2 + x^2 + y^2 + z^2)
normalized quaternion	[	w/magnitude	
				x/magnitude	
				y/magnitude	
				z/magnitude	]
multiply
q1*q2	=	[	w1*w2 - x1*x2 - y1*y2 - z1*z2
			w1*x2 + x1*w2 + y1*z2 - z1*y2
			w1*y2 - x1*z2 + y1*w2 + z1*x2
			w1*z2 + x1*y2 - y1*x2 + z1*w2	]
with angle and axis flag make a temporary quaternion 
RR = R/2
tmpW = cosRR
tmpX = Xaxis * sinRR
tmpY = Yaxis * sinRR
tmpZ = Zaxis * sinRR
tmpQ = [tmpW tmpX tmpY tmpZ]
we have a real quaternion who have already received previous rotation, in case not [1 0 0 0]
newQ = tmpQ * oldQ	//order matter with matrix
quaternion rotation
[	(w*w)+(x*x)-(y*y)-(z*z)	(2*x*y)-(2*w*z)	(2*x*z)+(2*w*y)	0	]
[	(2*x*y)+(2*w*z)	(w*w)-(x*x)+(y*y)-(z*z)	(2*y*z)-(2*w*x)	0	]
[	(2*x*z)-(2*w*y)	(2*y*z)-(2*w*x)	(w*w)-(x*x)-(y*y)+(z*z)	0	]
[	0		0		0			1	]
optimized
[	1-2*y*y-2*z*z	2*x*y-2*w*z	2*x*z+2*w*y	0	]
[	2*x*y+2*w*z	1-2*x*x-2*z*z	2*y*z-2*w*x	0	]
[	2*x*z-2*w*y	2*y*z-2*w*x	1-2*x*x-2*y*y	0	]
[	0		0		0		1	]
more optimized
xx = 2*x*x
yy = 2*y*y
zz = 2*z*z
xy = 2*x*y
xz = 2*x*z
yz = 2*y*z
wx = 2*w*x
wy = 2*w*y
wz = 2*w*z
[	1-yy-zz	xy-wz	xz+wy	0	]
[	xy+wz	1-xx-zz	yz-wx	0	]
[	xz-wy	yz-wx	1-xx-yy	0	]
[	0	0	0	1	]