function [p, S, mu] = polyfit(x, y, n)
%POLYFIT Convert roots to polynomial.
%   P = POLYFIT(X,Y,N)
%   [P,S] = POLYFIT(X,Y,N)
%   [P,S,MU] = POLYFIT(X,Y,N)
%   
%   Example:
%   spmd
%       x = codistributed.linspace(0, 3*pi, 1000);
%       y = sin(x);
%       p = polyfit(x, y, 5);
%       plot(x, y, 'r', x, polyval(p,x), 'b');
%   end
%   
%   See also POLYFIT, SIN, POLYVAL, CODISTRIBUTED, CODISTRIBUTED.LINSPACE.


%   Copyright 2014-2023 The MathWorks, Inc.

% Only first two inputs can be distributed
if nargin>2
    n = gather(n);
end

if ~iscodistributed(x) && ~iscodistributed(y)
    n = gather(n);
    switch nargout
        case {0, 1}
            p = polyfit(x, y, n);
        case 2
            [p, S] = polyfit(x, y, n);
        otherwise
            [p, S, mu] = polyfit(x, y, n);
    end
    return
end

if numel(x) ~= numel(y)
    error(message('MATLAB:polyfit:XYSizeMismatch'))
end

x = subsref(x, substruct('()', {':'}));
y = subsref(y, substruct('()', {':'}));
outputClass = superiorfloat(x, y);

if nargout > 2
    [sx, mx] = std(x);
    x = (x - mx) / sx;
    mu = [mx; sx];
end

% Construct the Vandermonde matrix V = [x.^n ... x.^2 x ones(size(x))]
V = codistributed.zeros(length(x), n+1, outputClass);
% V(:, n+1) = ones(length(x), 1, 'like', x);
V = subsasgn(V, substruct('()', {':', n+1}), ones(length(x), 1, 'like', x));
for j = n:-1:1
    % V(:, j) = x.*V(:, j+1);
    V = subsasgn(V, substruct('()', {':', j}), x.*subsref(V, substruct('()', {':', j+1})));
end
[nrowV, ncolV] = size(V);

% Convert y to the same class as V (leaving "remoteness" alone)
y = cast(full(y), outputClass);

% Solve least squares problem p = V\y to get polynomial coefficients p.
f = parallel.internal.decomposition.codistributed.DenseQR(V, []);
p = f.solve(y);

rankV = f.Rank;

% Override rank 0 solution to be NaN
if rankV == 0
    p = codistributed.NaN(ncolV, 1, outputClass);
end

if size(V, 1) < size(V, 2)
    warning(message('MATLAB:polyfit:PolyNotUnique'))
elseif rankV < size(V, 2)
    if nargout > 2
        warning(message('MATLAB:polyfit:RepeatedPoints'));
    else
        warning(message('MATLAB:polyfit:RepeatedPointsOrRescale'));
    end
end

if nargout > 1
    if ~anynan(V)
        QRfactor = f.getLocalQr();
        perm = f.getPiv();
        fCodistributor = f.getCodistributorA();
        QRfactor = codistributed.build(QRfactor, fCodistributor);
        QRfactor = redistribute(QRfactor, codistributor1d());
        perm = redistribute(perm, codistributor1d());
    else
        % Shortcut to return NaN
        QRfactor = codistributed.NaN(nrowV, ncolV, codistributor1d());
        perm = codistributed.colon(1, ncolV, codistributor1d());
    end

    r = y - V*p;
    minmn = min(size(V));

    perminv = codistributed.zeros(length(perm), 1);
    % R = triu(QRfactor(minmn, :))
    R = triu(subsref(QRfactor, substruct('()', {1:minmn, ':'})));
    % perminv(perm) = 1:length(perm);
    perminv = subsasgn(perminv, substruct('()', {perm}), 1:length(perm));

    % R = R(:,perminv);
    R = subsref(R, substruct('()', {':', perminv}));
    df = max(0, length(y) - (n+1));
    normr = norm(r);
    rsquared = 1 - (normr/norm(y - mean(y)))^2;
end

p = p.';

% The second output is a struct and all fields must be replicated, not distributed.
if nargout>=2
    [S.R, S.df, S.normr, S.rsquared] = gather(R, df, normr, rsquared);
end
