Rotation around line

format('v',15)

//rotation axis
u = [rand(), rand(), rand()];
u = u./ sqrt(u*u')

//angle
t = 2*%pi*rand();

c = cos(t);
s = sin(t);

T= [rand(); rand(); rand()]

R = [c+u(1)*u(1)*(1-c), u(1)*u(2)*(1-c)-u(3)*s, u(1)*u(3)*(1-c)+u(2)*s;
 u(1)*u(2)*(1-c)+u(3)*s, c+u(2)*u(2)*(1-c), u(2)*u(3)*(1-c)-u(1)*s;
 u(1)*u(3)*(1-c)-u(2)*s, u(2)*u(3)*(1-c)+u(1)*s, c+u(3)*u(3)*(1-c)]

M = [eye(3,3),T;0,0,0,1]*[R,[0;0;0];0,0,0,1]*[eye(3,3),-T;0,0,0,1]
//M = [eye(3,3),T;0,0,0,1]*[R,[0;0;0];0,0,0,1]

//identity
I = eye(4,4);

//system
A = I-M

//row reduction
S = rref(A)