% Givens formula applied to the matrix from the textbook:

% More efficient

clc

% Set up the problem
t = (-1:.5:1)';
A = [1+0*t t t.^2];
b = [1 .5 0 .5 2]';
[m, n] = size(A);

list = [4 5 1 ; 3 4 1 ; 2 3 1 ; 1 2 1 ; 4 5 2 ;3 4 2; 2 3 2 ; 4 5 3 ; 3 4 3];

Qt = eye(m);

while ( 1 ) % Loop forever
    A
    i = list(1,1); j = list(1,2); k = list(1,3);
    list(1,:) = [];
 
    % Make g:
    a = A(:,k);
    ai = a(i);
    aj = a(j);
    tmp = sqrt(ai^2+aj^2);
    c = ai/tmp;
    s = aj/tmp;
    G  = [c s ; -s c];
    
    % Update A, rhs b, and Q':
    A([i j], :) = G*A([i j], :);
    b([i j]) = G*b([i j]);
    Qt([i j], :) = G*Qt([i j], :);
    
    if ( isempty(list) ), break, end
end

% At this stage, R should be upper-triangular
A
R = A(1:n,1:n)

% Solve the least squares problem:
x = R\b(1:n)

Q = Qt'

Q(:,1:n)*R

