{   prog5a.pas -- Solve Least Squares Problem by QR Factorization
    Written Mar 30, 2001 by Eric Olson for Mathematics 483 

    For more information see Step 27 in First Steps in Numerical 
    Analysis, 2nd Edition, by Hosking, Joe, Joyce and Turner. }

program main;

const Nc=200;
const m=4;

type matrix=array [1..m] of array [1..Nc] of real;
    nvector=array [1..Nc] of real;
    mvector=array [1..m] of real;

var n:integer;
    At:matrix;
    c:mvector;
    y:nvector;

function phi(x:real; i:integer):real;
begin
    case i of
1:
        phi:=1;
2:
        phi:=x;
3:
        phi:=x*x;
4:
        phi:=x*x*x;
else
        begin
            writeln('Unknown basis function number ',i,'!');
            halt(1)
        end
    end
end;

procedure writeC;
    var i:integer;
begin
    write('C=[ ');
    for i:=1 to m do begin
        write(C[i],' ')
    end;
    writeln(']')
end;

procedure QRsolve;
    var w:nvector;
        i,j,k,l:integer;
        s,t:real;
begin
    for l:=1 to m do begin
        s:=0.0;
        for j:=l to n do s:=s+At[l][j]*At[l][j];
        s:=sqrt(s);
        if At[l][l]>0 then s:=-s;
        s:=2*s;
        t:=sqrt(0.5-At[l][l]/s);
        w[l]:=t;
        t:=-t*s;
        for j:=l+1 to n do w[j]:=At[l][j]/t;
        for i:=l to m do begin
            s:=0.0;
            for j:=l to n do s:=s+w[j]*At[i][j];
            s:=2*s;
            for j:=l to n do At[i][j]:=At[i][j]-w[j]*s
        end;
        s:=0.0;
        for j:=l to n do s:=s+w[j]*Y[j];
        s:=2*s;
        for j:=l to n do Y[j]:=Y[j]-w[j]*s
    end;
    for l:=m downto 1 do begin
        s:=0.0;
        for k:=l+1 to m do s:=s+At[k][l]*C[k];
        C[l]:=(Y[l]-s)/At[l][l]
    end
end;

var i,j:integer;
    xn,yn:real;

begin
    read(n);
    if n>Nc then begin
        writeln('Capacity of ',N,' data points exceeded!');
        halt(2)
    end;
    writeln('Reading ',n,' data points:');
    for j:=1 to n do begin
        read(xn,yn);
        for i:=1 to m do At[i][j]:=phi(xn,i);
        Y[j]:=yn
    end;
    QRsolve;
    writeC
end.