# L-15 MCS 471 Mon 26 Sep 2022 : newdivdif.jl

# This scripts illustrates the computation of divided differences
# for Newton interpolation.

"""
    divdif(x::Array{Float64,1},f::Array{Float64,1})

Returns the vector of divided differences for the
Newton form of the interpolating polynomial.

On entry are x and f;
x contains the abscisses, given as a column vector;
f contains the ordinates, given as a column vector.

On return are the divided differences.
"""
function divdif(x::Array{Float64,1},f::Array{Float64,1})
    n = length(x)
    d = deepcopy(f)
    for i=2:n
        for j=1:i-1
            d[i] = (d[j] - d[i])/(x[j] - x[i])
        end
    end
    return d
end

"""
Evaluates the Newton form of the interpolating polynomial,
with abscisses in x and divided differences in d at xx.
"""
function newtonform(x::Array{Float64,1},d::Array{Float64,1},xx::Float64)
    n = length(d)
    result = d[n]
    for i=n-1:-1:1
        result = result*(xx - x[i]) + d[i]
    end
    return result
end

"""
    newton(x::Array{Float64,1},f::Array{Float64,1},xx::Float64)

Implements the interpolation algorithm of Newton

ON ENTRY :
  x       abscisses, given as a column vector;
  f       ordinates, given as a column vector;
  xx      point where to evaluate the interpolating
          polynomial through (x[i],f[i]).

ON RETURN :
  d       divided differences, computed from and f;
  p       value of the interpolating polynomial at xx.

EXAMPLE :
  x = [32.0, 22.2, 41.6, 10.1, 50.5]
  f = [0.52992, 0.37784, 0.66393, 0.17537, 0.63608]
  xx = 27.5
  d, p = newton(x,f,xx)
"""
function newton(x::Array{Float64,1},f::Array{Float64,1},xx::Float64)
    divided = divdif(x,f)
    result = newtonform(x,divided,xx) 
    return divided, result
end

"""
Runs the example in newton.
"""
function main()
    x = [32.0, 22.2, 41.6, 10.1, 50.5]
    f = [0.52992, 0.37784, 0.66393, 0.17537, 0.63608]
    xx = 27.5
    d, p = newton(x,f,xx)
    println(d)
    println(p)
end

main()