function A = polyfit(x,y,n)
    // Polynomial fitting function
    // Usage:
    //       A = polyfit(x,y,n)
    // where
    //       x: Vector containing values of the independent variable
    //       y: Vector containing values of the dependent variable
    //       n: Degree of the polynomial to fit the data
    //       A: Vector of n+1 components containing the coefficients
    //          of the polynomial of degree n that best fits the
    //          data in the sense that the sum of the squared errors
    //          are minimal (least square error). The coefficients 
    //          are given in increasing degree order.
    //
    // NOTES: 1) The polynomial can be computed as
    //           P = poly(A,"x","coeff")
    //           and its value in a given array of data x, as
    //           Y = horner(P,x)
    //        2) If n is non integer, floor(n) is used
    //        3) If n is negative a polynomial of 1/x is fitted
    //           In such case A is ordered in increasing absolute
    //           value of the degree
    //        4) x and y must have the same length
    //  
    // --------------------------------
    // Author: Federico Miyara
    // Date:   2020-03-13
    
    // Handle input arguments
    rhs = argn(2);
    
    // Test data --delete--
    if 1==2
        x = 1:20
        y = 1./x
        n=-1
    end

    if rhs<3
        disp("n is not specified, linear regression is used")
        n = 1;
    end
    if floor(n)~=n
        n = floor(n);
        msg = gettext("%s: n is non-integer, floor(n) is used.\n")
        mprintf(msg, "polyfit")        
    end
    if n<0 then
        if or(x==0)==%t
            msg = gettext("%s: n is negative, x cnnot be 0.\n")
            error(msprintf(msg, "polyfit"))
        else
            x = 1 ./ x;
            n = -n;
        end
    end    
         
    
    // Ensure x and y are column vectors
    x = x(:);
    y = y(:);
    
    // Initialize matrix X
    X = ones(length(x),n+1);
    
    tic
    // Compute Vandermonde's matrix
    for k =2:n+1
        X(:,k) = X(:,k-1).*x;
    end
    toc
    
    // Apply the Moore-Penrose pseudoinverse matrix and
    // multiply by the dependent data vector to get the
    // least squares approximation of the polynomial
    // coefficients
    A = inv(X'*X)*X'*y;
            
endfunction
