%% Copyright (c) 2025 Duong Ngoc Son <son@duongnotes.com>

%% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
%% "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
%% LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
%% A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
%% HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
%% SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
%% LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
%% DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
%% THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
%% (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
%% OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

%' @title Polynomial Interpolation
%'
%' @description
%' Compute the coefficients of the interpolating polynomial by solving the
%' Vandermonde matrix. 

% Hàm Vandermonde(x) tính định thức Vandermonde

function vdm = Vandermonde(x)
	n = length(x);
	vdm = ones(n, 1);
	c = x.';	% Chuyển vị của vectơ x
    t = c;
    
    % Xây dựng ma trận Vandermonde bằng cách thêm vào các cột của nó
    
	for i = 1:(n-1)
		vdm = [vdm t];
		t = t.*c;
	end
end

% Hàm InterpolatingPolynomial() tính các hệ số của đa thức nội suy

function p = InterpolatingPolynomial(x, y)
	n = length(x);
	m = length(y);
	if m ~= n
		p = NaN;	% Độ dài các dãy x và y khác nhau.
	elseif n == 1
		p = [y(1)];
	else
		A = Vandermonde(x);
		if det(A) == 0
			p = NaN;
		else
			p = (A \ (y.')).';	% Giải hệ phương trình tuyến tính
		end
	end
end

% Ví dụ áp dụng: x = (0, 2, 4, 5), y = (-2, 8, 115, 150)
% p(x) = -2 - (2911/60) * x + (1363/40) * x^2 - (439/120) * x^3

x = [0 2 4 5];
y = [-2 8 115 150];
xi = 1.2;

% Tính các hệ số của đa thức và hiển thị theo thứ tự bậc tăng dần

w = InterpolatingPolynomial(x, y);
disp(w);

