function nma_proj_1()
%
% function nma_proj_1
% By Nasser Abbasi.
%
% Project 1 for course MAE 185 by Dr Chaudhry
% UCI. April 30, 2003.
% 

TRUE=1;
FALSE=0;

[degree,a] = getInput;

% use these variable to store all the quadratic factors
% for table display
maxNumberOfFactors=degree/2 +1;
savedFactors=zeros(maxNumberOfFactors,3); 
savedInitialR=zeros(maxNumberOfFactors);
savedInitialS=zeros(maxNumberOfFactors);
savedErrorPer=zeros(maxNumberOfFactors);
savedNumberOfIter=zeros(maxNumberOfFactors);

%Save original poly before factorization.
originalPoly=a;

fprintf('\nInput polynomial is:\n');
displayPolynomial(a);
fprintf('\n');

moreReductionNeeded=TRUE;

nReduced=0;

%
% Main loop. Loop untill no more factorization needed
%
while(moreReductionNeeded)
    
    nReduced=nReduced+1;    
    fprintf('\nInput for round %d of reduction\n',nReduced);
    [r,s,tol]=getIterationInput;       
    
    [degree,a,q,nIter]=nma_reduce(degree,a,r,s,tol,nReduced);
    savedFactors(nReduced,:)=q(1,:);
    savedInitialR(nReduced)=r;
    savedInitialS(nReduced)=r;
    savedInitialPer(nReduced)=tol;    
    savedNumberOfIter(nReduced)=nIter;        
    
    if(degree<=2)
        moreReductionNeeded=FALSE;
        fprintf('\n*******************************************\n');
        fprintf('Reduction completed...... Displaying final results \n\n');
        fprintf('f(x)='); displayPolynomial(originalPoly); fprintf('\n\n');
        
        for(i=1:nReduced)        
            fprintf('Q(x)=('); displayPolynomial(savedFactors(i,:)); fprintf(')\n');
            fprintf('Initial R= %d\n',savedInitialR(i));
            fprintf('Initial S= %d\n',savedInitialS(i));
            fprintf('Error percentage= %6.6f\n',savedInitialPer(i));
            fprintf('Number of Iterations= %d\n\n',savedNumberOfIter(i));                      
        end
        fprintf('Remaining polynomial = ('); displayPolynomial(a); fprintf(')');
    end    
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function called by main loop to generate one
% quadratic factor using Bairstow algorithm.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [newDegree,new_a,q,nIter]=nma_reduce(degree,a,r,s,tol,nReduced)

b=zeros(1,degree+3);
c=zeros(1,degree+3);

nIter=0;

while(1);
    nIter=nIter+1;
    
    b(end)=0;
    b(end-1)=0;
    c(end)=0;
    c(end-1)=0;
    
    for(i=degree+1:-1:1)
        b(i) = a(i) + r*b(i+1) + s*b(i+2);    
    end
    
    for(i=degree+1:-1:2)
        c(i) = b(i) + r*c(i+1) + s*c(i+2);    
    end
    
    den=nma_det(c(3),c(4),c(2),c(3));
    
    new_r = r + nma_det(-b(2),c(4),-b(1),c(3)) / den;
    new_s = s + nma_det(c(3),-b(2),c(2),-b(1)) / den;
    
    fprintf('\nIteration= %d\t r= %6.15f\t s= %6.15f',nIter,new_r,new_s);
    
    if(  abs((new_r-r)/new_r) *100  <= tol  &  abs((new_s-s)/new_s)  *100  <= tol )
        fprintf('\n\nCompleted round %d of quadratic factorization. Result is:\n\n',nReduced);
        fprintf('Q(x)= x^2 %+6.6f x %+6.6f\n',-new_r,-new_s);        
        q=[-new_s -new_r 1];
        
        newDegree=degree-2;
        numOfTerms=newDegree+1;
        new_a=b(3:numOfTerms+2);        
        fprintf('f(x)= ');
        displayPolynomial(new_a);
        fprintf('\n');
        break;
    else
        r=new_r;
        s=new_s;
    end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function den=nma_det(a11,a12,a21,a22)
den=( a11*a22 - a21*a12 );

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% reads input. Degree and Coefficients.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [degree,a]=getInput()

INTEGER=0;
FLOAT=1;

upperLimit=12;
prompt=sprintf('degree of poynomial? [max=%d]  >',upperLimit);

prompt=sprintf('degree of polynomial? >');

degree = nma_inputNumeric(INTEGER,prompt,3,upperLimit);
a      = inputPolyCoefficients(degree);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Reads input for each iteration. The r,s, and error
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [r,s,tol]=getIterationInput()

INTEGER=0;
FLOAT=1;

prompt=sprintf('Initial value for r?  >');
r = nma_inputNumeric(FLOAT,prompt,-Inf,Inf);

prompt=sprintf('Initial value for s?  >');
s = nma_inputNumeric(FLOAT,prompt,-Inf,Inf);

prompt=sprintf('percentage error?  >');
tol = nma_inputNumeric(FLOAT,prompt,eps,100);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function that reads the polynomial coefficients from user
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function a = inputPolyCoefficients(degree)

FLOAT=1;

a=zeros(1,degree+1);

for(i=1:degree+1)
    prompt=sprintf('Coefficent a%d ?  >',i-1);
    a(i)=nma_inputNumeric(FLOAT,prompt,-Inf,Inf);    
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
% a pretty print of a polynomial
% (not the best code, but have no more time to improve it.
% it does the job though)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function displayPolynomial(a)

numOfTerms=length(a);
for(i=1:numOfTerms)
    
    rounded=round(a(i));            
    
    if(i==1)
        if(rounded ~= 0 )
            if( abs(rounded-a(i))<eps)
                fprintf('+%d  ',rounded);               
            else
                fprintf('%6.6f  ',a(i));                                   
            end        
        end
    else
        if(i==2)
            if(rounded ~= 0 )                
                if( abs(rounded-a(i))<eps)
                    if(abs(rounded)==1)
                        if(rounded<0)
                            fprintf('%- x  ',rounded);               
                        else
                            fprintf('%+ x  ',rounded);                                   
                        end
                    else
                        fprintf('%+d x  ',rounded);                                  
                    end
                else
                    fprintf('%+ 6.6f x  ',a(i));                               
                end
            end
        else
            if(rounded ~= 0 )
                
                if( abs(rounded-a(i))<eps)
                    if(abs(rounded)==1)
                        if(rounded<0)
                            fprintf('- x^%d  ',i-1);
                        else
                            fprintf('+ x^%d  ',i-1);
                        end
                    else
                        fprintf('%+d x^%d  ',rounded,i-1);
                    end                
                else            
                    fprintf('%+ 6.6f x^%d  ',a(i),i-1);                                 
                end        
            end
        end
    end
end