# Factoring orthonormal into Givens' : the good, the bad, ... 
#   text output U/S : try screen-grab! 
restart; # CMND-3 to zoom, Text input 
with(LinearAlgebra): 

RB := 10^10; # precision of generator ~ 1.0E-12 
randunit := proc () 1.0 - (rand()mod(2*RB))/RB end; # random real in (-1, 1] 

# Givens rotation matrix, given cos & sin as  t^2 + s^2 = 1 
given := proc(n, p, q, t, s) local i, j, R; 
  R := [seq([seq(`if`(i = j, 1, 0), j = 1..n)], i = 1..n)]; 
  R[p][p] := t; R[q][q] := t; R[p][q] := -s; R[q][p] := s; 
  Matrix(R) end; 

# Eliminate A_pq from A (ignoring overflow): Digits trampled (Maple BUG ??) 
elim := proc(A, n, p, q) local s, t, r, R; 
  s := A[p, q]; t := A[q, q]; r := sqrt(s^2 + t^2); 
  R := given(n, p, q, t/r, s/r); 
  Multiply(R, A) end; 

# Random special orthonormal matrix 
randorth := proc(n) local i, j, k, R, A, x, l; 
  A := Matrix([seq([seq(`if`(i = j, 1, 0), j = 1..n)], i = 1..n)]); 
  l := binomial(n, 2); # isometry freedom 
  for k from 1 to l do x := evalf(Pi*randunit()); 
    R := given(n, k mod(n-1)+1, k mod(n-1)+2, cos(x), sin(x)); 
    A := Multiply(R, A); od; 
  eval(A) end; 

n := 4; # matrix order 
Digits := 6; # ignored (Maple BUG ??) 
t,s := 't','s'; 
given(n, 2, 4, t, s); # test 

# Column-by-column downwards: success! 
A0 := randorth(n); 
A1 := elim(A0, n, 2, 1); 
A2 := elim(A1, n, 3, 1); 
A3 := elim(A2, n, 4, 1); 
A4 := elim(A3, n, 3, 2); 
A5 := elim(A4, n, 4, 2); 
A6 := elim(A5, n, 4, 3); 

# Triangular growing SE starting SW corner: col 1 corrupt? 
A0 := randorth(n); 
A1 := elim(A0, n, 4, 1); 
A2 := elim(A1, n, 3, 1); 
A3 := elim(A2, n, 4, 2); 
A4 := elim(A3, n, 2, 1); 
A5 := elim(A4, n, 3, 2); 
A6 := elim(A5, n, 4, 3);