Lagrange Interpolation

Lagrange Interpolation is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.

The task is to implement the Lagrange Interpolation formula and use it to solve the example problem to find a polynomial P of degree<4 satisfying P(1)=1 P(2)=4 P(3)=1 P(4)=5 as described at [1]

Related task Curve that touches three points

F#

// Lagrange Interpolation. Nigel Galloway: September 5th., 2023
let symbol=MathNet.Symbolics.SymbolicExpression.Variable
let qi=MathNet.Symbolics.SymbolicExpression.FromInt32
let eval (g:MathNet.Symbolics.SymbolicExpression) x=let n=Map["x",MathNet.Symbolics.FloatingPoint.Real x] in MathNet.Symbolics.Evaluate.evaluate n g.Expression
let fN g=let x=symbol "x" in g|>List.fold(fun z c->(x-c)*z)(qi 1)
let fG(n,g)=let n,g=n|>List.map qi,g|>List.map qi in List.map2(fun i g->i,g,n|>List.except [i]) n g
let LIF n=fG n|>List.sumBy(fun(ci,bi,c)->bi*(fN c)/(c|>List.fold(fun z c->(ci-c)*z)(qi 1)))
printfn $"%s{LIF([1;2;3;4],[1;4;1;5]).Expand().ToString()}"

Output:

-21 + 215/6*x - 16*x^2 + 13/6*x^3

Julia

Note the Polynomials module prints polynomials in order from degree 0 to n rather than n to zero degree.

using Polynomials, SpecialPolynomials

const pts = [[1, 1], [2, 4], [3, 1], [4, 5]]
const xs = first.(pts)
const ys = last.(pts)
const p = Lagrange(xs, ys)

@show p Polynomial(p)

Output:

p = Lagrange(1⋅ℓ_0(x) + 4⋅ℓ_1(x) + 1⋅ℓ_2(x) + 5⋅ℓ_3(x))
Polynomial(p) = Polynomial(-21.0 + 35.83333333333333*x - 16.0*x^2 + 2.1666666666666665*x^3)

Rational base polynomial answer

using Polynomials

function Lagrange(pts::Vector{Vector{Int}})
    xs = first.(pts)
    cs = last.(pts)
    n = length(xs)
    n == 1 && return Polynomial(cs[1]) 
    arr = ones(Int, n)
    for j in 2:n
        for k in 1:j
            arr[k] = (xs[k] - xs[j]) * arr[k]
        end
        arr[j] = prod(xs[j] - xs[i] for i in 1:(j - 1))
    end
    ws = 1 .// arr
    q = Polynomial(0 // 1)
    x = variable(q)
    for i in eachindex(ws)
        m = prod(x - xs[j] for j in eachindex(xs) if j != i)
        q += m * ws[i] * cs[i]
    end
    return q
end

const testpoints = [[1, 1], [2, 4], [3, 1], [4, 5]]
@show Lagrange(testpoints)

Output:

Lagrange(testpoints) = Polynomial(-21//1 + 215//6*x - 16//1*x^2 + 13//6*x^3)

Phix

Translation of: Wren

with javascript_semantics
function add_poly(sequence p1, p2)
    integer l1 = length(p1),
            l2 = length(p2)
    if l2>l1 then p1 &= repeat(0,l2-l1) end if
    for i=1 to l2 do
        p1[i] += p2[i]
    end for
    return p1
end function

function mul_poly(sequence p1,p2)
    integer m = length(p1),
            n = length(p2)
    sequence res = repeat(0,m+n-1)
    for i=1 to m do
        for j=1 to n do
            res[i+j-1] += p1[i]*p2[j]
        end for
    end for
    return res
end function

function scalar_multiply(sequence poly, atom x)
    return sq_mul(poly,x)
end function

function sclar_divide(sequence poly, atom x)
    return sq_div(poly,x)
end function

function eval_poly(sequence poly, atom x)
    integer c = length(poly)
    atom res = poly[$]
    for i=length(poly)-1 to 1 by -1 do
        res = res*x + poly[i]
    end for
    return res
end function

procedure show_poly(sequence poly)
    integer l = length(poly)
    for i=l to 1 by -1 do
        atom p = poly[i]
        if i<l then 
            puts(1,iff(p<0?" ":" +")) 
        end if
        printf(1,"%g",p)
        if i>1 then
            puts(1,"x")
            if i>2 then
                printf(1,"^%d",i-1)
            end if
        end if
    end for
    puts(1,"\n")
end procedure
constant show = show_poly

enum X, Y
function lagrange(sequence pts)
// Returns the Lagrange interpolating polynomial which passes through a list of points.
    integer c = length(pts)
    sequence polys = repeat(null,c)
    for i=1 to c do
        sequence poly = {1}
        for j=1 to c do
            if i!=j then
                poly = mul_poly(poly,{-pts[j][X],1})
            end if
        end for
        atom d = eval_poly(poly,pts[i][X])
        polys[i] = sclar_divide(poly,d)
    end for
    sequence res = {0}
    for i=1 to c do
        polys[i] = scalar_multiply(polys[i],pts[i][Y])
        res = add_poly(res,polys[i])
    end for
    return res
end function

constant pts = {{1,1},{2,4},{3,1},{4,5}}
sequence lip = lagrange(pts)
show(lip)

Output:

2.16667x^3 -16x^2 +35.8333x -21

RPL

Works with: HP version 49

[[1 2 3 4][1 4 1 5]] LAGRANGE

Output:

1: '(13*X^3-96*X^2+215*X-126)/6'

Wren

Library: Wren-dynamic

Library: Wren-math

Library: Wren-fmt

A polynomial is represented here by a list of coefficients from the lowest to the highest degree. However, the library methods which deal with polynomials expect the coefficients to be presented from highest to lowest degree so we therefore need to reverse the list before calling these methods.

import "./dynamic" for Tuple
import "./math" for Math
import "./fmt" for Fmt

var Point = Tuple.create("Point", ["x", "y"])

// Add two polynomials.
var add = Fn.new { |p1, p2|
    var m = p1.count
    var n = p2.count
    var sum = List.filled(m.max(n), 0)
    for (i in 0...m) sum[i] = p1[i]
    for (i in 0...n) sum[i] = sum[i] + p2[i]
    return sum
}

// Multiply two polynmials.
var multiply = Fn.new { |p1, p2|
    var m = p1.count
    var n = p2.count
    var prod = List.filled(m + n - 1, 0)
    for (i in 0...m) {
        for (j in 0...n) prod[i+j] = prod[i+j] + p1[i] * p2[j]
    }
    return prod
}

// Multiply a polynomial by a scalar.
var scalarMultiply = Fn.new { |poly, x| poly.map { |coef| coef * x }.toList }

// Divide a polynomial by a scalar.
var scalarDivide = Fn.new { |poly, x| scalarMultiply.call(poly, 1/x) }

// Returns the Lagrange interpolating polynomial which passes through a list of points.
var lagrange = Fn.new { |pts|
    var c = pts.count
    var polys = List.filled(c, null)
    for (i in 0...c) {
        var poly = [1]
        for (j in 0...c) {
            if (i == j) continue
            poly = multiply.call(poly, [-pts[j].x, 1])
        }
        var d = Math.evalPoly(poly[-1..0], pts[i].x)
        polys[i] = scalarDivide.call(poly, d)
    }
    var sum = [0]
    for (i in 0...c) {
        polys[i] = scalarMultiply.call(polys[i], pts[i].y)
        sum = add.call(sum, polys[i])
    }
    return sum
}

var pts = [
    Point.new(1, 1),
    Point.new(2, 4),
    Point.new(3, 1),
    Point.new(4, 5)
]
var lip = lagrange.call(pts)
Fmt.pprint("$f", lip[-1..0], "", "x")

Output:

2.166667x³ - 16.000000x² + 35.833333x - 21.000000