  function [c,s] = Rotate(x,y)
% x and y are unit 2-norm column n-vectors.
% |x(1)| <= ... <= |x(n)| and |y(1)|>= ... >= |y(n)|
% c and s are column (n-1)-vectors  with c(i)^2 + s(i)^2 = 1, i=1:n-1.
% y = G_{1}...G_{n-1}x  where G_{i} is the identity except
% G_{i}(i:i+1,i:i+1) =  [ c(i) s(i) ; -s(i) c(i)]
 
% Theorem: if u and v are m-vectors with the same norm satisfying
% |u(1)|<=...<= |u(m)| and |v(1)|>=...>=|v(m)|, then 
% |u(m)| >= |v(m)|.
% Proof: |v(m)|<= norm(v,2)/sqrt(m) = norm(u,2)/sqrt(m) <= |u(m)|

  n = length(x);
  c = zeros(n-1,1);
  s = zeros(n-1,1);
  for i = n-1:-1:1
%    Generate G_{i} 
     [c(i),s(i)] = twist(x(i),x(i+1),y(i),y(i+1));
%    Need the top component of the updated x(i:i+1)...
     x(i) = c(i)*x(i) + s(i)*x(i+1);
     
% Note that at this stage, x(1:i) and y(1:i) have the same 2-norm
% since x(i+1:n) = y(i+1:n). Moreover, since the new x(i) is larger
% than the original in absolute value, so |x(1))<= ... <= |x(i)|.
% Thus, because of the theorem the process will complete.
  end
  
  
 
 
 
  function [c,s] = twist(alpha,beta,delta,gamma)
% alpha, beta, delta, and gamma are scalars
% alpha^2 + beta^2 >= gamma^2
% c and s are a sine-cosine pair so that
% [c s;-s c]*[alpha;beta] = [mu;gamma] where sign(mu) = sign(delta)

  if alpha^2 + beta^2 == 0
%     Degenerate case
      c = 1; s = 0;
  else
      
% First determine (c1,s1) so [c1 s1;-s1 c1]*[alpha;beta] = [ r;0].
% This is the standard Givens rotation problem...

%    s1*alpha = c1*beta 
     if abs(alpha)> abs(beta)
        tau = beta/alpha; c1 = 1/sqrt(1+tau^2); s1 = c1*tau;
     else
        tau = alpha/beta; s1 = 1/sqrt(1+tau^2); c1 = tau*s1;
     end
     r = c1*alpha + s1*beta;
  
% Now find (c2,s2) so gamma is the second component of [c2 s2;-s2 c2]*[r;0]
% i.e., gamma = -s2*r:
     s2 = -gamma/r;
% This is <=1 since |r|^2 = alpha^2 + beta^2

% r*c2 should have same sign as delta..
     if delta*r>0
         c2 =  sqrt(1-s2^2);
     else
         c2 = -sqrt(1-s2^2);
     end

% [c s;-s c] = [c2 s2;-s2 c2]*[c1 s1; -s1 c1]
     c = c2*c1 - s2*s1;
     s = c2*s1 + s2*c1;
  end
  
  
% -3, -5 if divide by zero, related mistakes
% -5 if G_1 not correct