function P3()
% Checks TwoVecRotate and Uniform

% Check TwoVecRotate
clc
format long
disp('Example 1:')
x = [3;4]; alfa = 2;
[c,s] = TwoVecRotate(x,alfa);
y = [c s;-s c]'*x
disp('Example 2:')
x = [4;3]; alfa = 2;
[c,s] = TwoVecRotate(x,alfa);
y = [c s;-s c]'*x

% Check uniform...
disp('Example 3:')
n = 9;
v = sort(abs(randn(n,1)));
v = v/norm(v);
[c,s,idx,jdx] = Uniform(v);
disp(' c s c^2+s^2 idx(i) jdx(i)')
disp('---------------------------------------------------')
for i=1:n-1
fprintf('%8.3f %8.3f %8.3f %2d %2d\n',c(i),s(i),c(i)^2+s(i)^2,idx(i),jdx(i))
end
for k=1:n-1
i = idx(k); j = jdx(k);
v([i j]) = [c(k) s(k);-s(k) c(k)]'*v([i j]);
end
final_v = v

disp('Example 4:')
disp(' ')
n = 9;
v = randn(n,1);
[~,idx] = sort(abs(v),'descend');
v = v(idx);
v = v/norm(v);
[c,s,idx,jdx] = Uniform(v);
disp(' c s c^2+s^2 idx(i) jdx(i)')
disp('---------------------------------------------------')
for i=1:n-1
fprintf('%8.3f %8.3f %8.3f %2d %2d\n',c(i),s(i),c(i)^2+s(i)^2,idx(i),jdx(i))
end
for k=1:n-1
i = idx(k); j = jdx(k);
v([i j]) = [c(k) s(k);-s(k) c(k)]'*v([i j]);
end
final_v = v

function [c,s,idx,jdx] = UniformCVL(v)
% v is a column n-vector with unit 2-norm.
% c and s are column (n-1)-vectors with c.^2 + s.^2 = ones(n-1,1) and
% idx is a column (n-1)-vector of integers with the property that
%
% G(c(n-1),s(n-1),n-1,idx(n-1))'*...*G(c(1),s(1),1,idx(1))*v = ones(n,1)/sqrt(n)

n = length(v);
[y,kdx] = sort(abs(v),'descend');
mu = 1/sqrt(n);
idx = zeros(n-1,1);
jdx = zeros(n-1,1);
c = zeros(n-1,1);
s = zeros(n-1,1);
for k=1:n-1
% Replace the kth largest value in v with mu.
if kdx(k)<kdx(k+1)
i = kdx(k); j = kdx(k+1);
% [ck sk;-sk ck]'*[v(i);v(j)] = [mu;beta]
% where mu^2+beta^2 = v(i)^2+v(j)^2
[ck,sk] = TwoVecRotate(v([i j]),mu);
v([i j]) = [ck sk;-sk ck]'*v([i j]);
c(k)=ck; s(k)=sk; idx(k)=i; jdx(k)=j;
else
i = kdx(k+1); j = kdx(k);
% [ck sk;-sk ck]'*[v(j);v(i)] = [beta;mu]
% where mu^2+beta^2 = v(i)^2+v(j)^2
[ck,sk] = TwoVecRotate(v([j i]),mu);
v([i j]) = [ck -sk; sk ck]'*v([i j]);
c(k)=ck; s(k)=-sk; idx(k)=i; jdx(k)=j;
end
end

function [c,s] = TwoVecRotateCVL(x,alfa)
% x is a nonzero column 2-vector and 0<=alfa<=norm(x,2)
% c and s satisfy c^2 + s^2 = 1 and have the property that if [c s,-s,c]'*x = y,
% then y(1) = alfa

tau = -x(1)/x(2); s1 = 1/sqrt(1+tau^2); c1 = s1*tau;
r = c1*x(1)-s1*x(2);
c2 = alfa/r;
s2 = sqrt(1-c2^2);
%[c s;-s c] = [c1 s1;-s1 c1]*[c2 s2; -s2 c2]
c = c1*c2 - s1*s2;
s = c1*s2 + s1*c2;