%% 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 Forward Newton Interpolation Formula
%'
%' @description
%' Compute value of the interpolating polynomial by divided difference
%' formula of Newton 

% Tính các hệ số trong nội suy Newton và giá trị của nó tại một điểm.

function [t, w] = NewtonInterpolation(x, y, xi)

	% Tham số: x và y là hai dãy dữ liệu cần nội suy, a là điểm cần
	% tính giá trị của đa thức nội suy.
	% Kết quả: t là dãy hệ số trong công thức nội suy Newton, w là giá
	% trị của đa thức nội suy tại x = xi.
	
    % Tính độ dài của các dãy x và y

	n = length(x);
	m = length(y);

    % Kiểm tra độ dài các dãy có phù hợp hay không

	if (m ~= n)
		t = 0;	% Hai dãy x và y khác độ dài
	elseif (n == 1)
	
		% Tính các hệ số A_i trong công thức nội suy Newton
		
		t = y(1);
	else
		t = zeros(1, n);	% Dãy t gồm n thành phần
		t(1) = y(1);
		for i = 2:1:n
			for j = 1:1:(n - i + 1)
				y(j) = (y(j + 1) - y(j)) / (x(j + i - 1) - x(j));
            end
			t(i) = y(1);
        end
    end

	% Tính giá trị của đa thức nội suy P(x) tại x = xi ở dạng lồng nhau
  
	w = t(n);
	while n > 1
		n = n - 1;
		w = t(n) + w * (xi - x(n));
	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 giá trị của đa thức nội suy đối với ví dụ trên

[t, w] = NewtonInterpolation(x, y, xi);

% In ra giá trị của đa thức nội suy tại a

fprintf('%d\n', w);
