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 https://brilliant.org/wiki/lagrange-interpolation or https://en.wikipedia.org/wiki/Lagrange_polynomial
- Related task Curve that touches three points
C++
#include <algorithm>
#include <cmath>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <string>
#include <vector>
// A double[] is used to represents a Polynomial
// with its coefficients reversed compared to the standard mathematical notation.
// For example, the polynomial 3x^2 + 2x + 1 is represented by the array [1, 2, 3].
std::vector<double> add(const std::vector<double>& one, const std::vector<double>& two) {
std::vector<double> sum(std::max(one.size(), two.size()), 0.0);
for ( uint32_t i = 0; i < one.size(); ++i ) {
sum[i] = one[i];
}
for ( uint32_t i = 0; i < two.size(); ++i ) {
sum[i] += two[i];
}
return sum;
}
std::vector<double> multiply(const std::vector<double>& one, const std::vector<double>& two) {
std::vector<double> product(one.size() + two.size() - 1, 0.0);
for ( uint32_t i = 0; i < one.size(); ++i ) {
for ( uint32_t j = 0; j < two.size(); ++j ) {
product[i + j] += one[i] * two[j];
}
}
return product;
}
std::vector<double> scalar_multiply(std::vector<double> vec, const double& value) {
std::vector<double> result(vec.size(), 0.0);
std::transform(vec.begin(), vec.end(), result.begin(), [value](double& d) { return d * value; });
return result;
}
std::vector<double> scalar_divide(std::vector<double> vec, const double& value) {
return scalar_multiply(vec, 1.0 / value);
}
double evaluate(const std::vector<double>& vec, const double& value) {
double result = 0.0;
for ( int32_t i = vec.size() - 1; i >= 0; --i ) {
result = result * value + vec[i];
}
return result;
}
void display(const std::vector<double> vec) {
const int32_t degree = vec.size() - 1;
if ( degree == 0 ) {
std::cout << std::fixed << std::setprecision(5) << vec[0] << std::endl;
return;
}
for ( int32_t i = degree; i >= 0; i-- ) {
if ( vec[i] == 0.0 ) {
continue;
}
const std::string sign = ( vec[i] < 0.0 && i == degree ) ?
"-" : ( vec[i] < 0.0 ) ? " - " : ( i < degree ) ? " + " : "";
std::cout << sign;
const double coeff = std::abs(vec[i]);
if ( coeff > 1.0 ) {
std::cout << std::fixed << std::setprecision(5) << coeff;
}
const std::string term = ( i > 1 ) ? "x^" + std::to_string(i) : ( i == 1 ) ?
"x" : ( coeff == 1.0 ) ? "1" : "";
std::cout << term;
}
std::cout << std::endl;
}
struct Point {
double x;
double y;
};
std::vector<double> lagrange_interpolation(const std::vector<Point>& points) {
std::vector<std::vector<double>> polys(points.size(), std::vector<double>(points.size(), 0.0));
for ( uint32_t i = 0; i < points.size(); ++i ) {
std::vector<double> poly(1, 1.0);
for ( uint32_t j = 0; j < points.size(); ++j ) {
if ( i != j ) {
poly = multiply(poly, { -points[j].x, 1.0 });
}
}
const double value = evaluate(poly, points[i].x);
polys[i] = scalar_divide(poly, value);
}
std::vector<double> sum(1, 0.0);
for ( uint32_t i = 0; i < points.size(); ++i ) {
polys[i] = scalar_multiply(polys[i], points[i].y);
sum = add(sum, polys[i]);
}
return sum;
}
int main() {
const std::vector<Point> points = { Point(1, 1), Point(2, 4), Point(3, 1), Point(4, 5) };
display(lagrange_interpolation(points));
}
- Output:
2.16667x^3 - 16.00000x^2 + 35.83333x - 21.00000
EasyLang
func[] add_poly p1[] p2[] .
l1 = len p1[]
l2 = len p2[]
if l2 > l1
len p1[] len p2[]
.
for i = 1 to l2
p1[i] += p2[i]
.
return p1[]
.
func[] mul_poly p1[] p2[] .
m = len p1[]
n = len p2[]
len res[] m + n - 1
for i = 1 to m
for j = 1 to n
res[i + j - 1] += p1[i] * p2[j]
.
.
return res[]
.
func[] scal_mul p[] x .
for e in p[]
r[] &= e * x
.
return r[]
.
func[] scal_div p[] x .
return scal_mul p[] (1 / x)
.
func eval_poly p[] x .
res = p[$]
for i = len p[] - 1 downto 1
res = res * x + p[i]
.
return res
.
proc show_poly p[] . .
l = len p[]
for i = l downto 1
p = p[i]
if i < l
if p < 0
write " - "
p = -p
else
write " + "
.
.
write p
if i > 1
write "x"
if i > 2
write "^" & i - 1
.
.
.
print ""
.
func[] lagrange pts[][] .
c = len pts[][]
for i = 1 to c
poly[] = [ 1 ]
for j = 1 to c
if i <> j
poly[] = mul_poly poly[] [ -pts[j][1] 1 ]
.
.
d = eval_poly poly[] pts[i][1]
polys[][] &= scal_div poly[] d
.
res[] = [ 0 ]
for i = 1 to c
polys[i][] = scal_mul polys[i][] pts[i][2]
res[] = add_poly res[] polys[i][]
.
return res[]
.
pts[][] = [ [ 1 1 ] [ 2 4 ] [ 3 1 ] [ 4 5 ] ]
lip[] = lagrange pts[][]
show_poly lip[]
- Output:
2.17x^3 - 16x^2 + 35.83x - 21
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
FreeBASIC
Sub AddPoly(p1() As Single, p2() As Single)
Dim As Integer i, l1 = Ubound(p1), l2 = Ubound(p2)
If l2 > l1 Then
Redim Preserve p1(l2)
For i = l1 + 1 To l2
p1(i) = 0
Next i
End If
For i = 0 To l2
p1(i) += p2(i)
Next i
End Sub
Sub MulPoly(res() As Single, p1() As Single, p2() As Single)
Dim As Integer i, j, m = Ubound(p1), n = Ubound(p2)
Redim res(m + n)
For i = 0 To m
For j = 0 To n
res(i + j) += p1(i) * p2(j)
Next j
Next i
End Sub
Sub ScalarMultiply(poly() As Single, x As Single)
For i As Integer = 0 To Ubound(poly)
poly(i) *= x
Next i
End Sub
Sub ScalarDivide(poly() As Single, x As Single)
For i As Integer = 0 To Ubound(poly)
poly(i) /= x
Next i
End Sub
Function EvalPoly(poly() As Single, x As Single) As Single
Dim As Integer i, c = Ubound(poly)
Dim As Single res = poly(c)
For i = c - 1 To 0 Step -1
res = res * x + poly(i)
Next i
Return res
End Function
Sub ShowPoly(poly() As Single)
Dim As Integer i, l = Ubound(poly)
Dim As Single p
For i = l To 0 Step -1
p = poly(i)
If i < l Then Print Iif(p < 0, " -", " +") & Abs(p); Else Print Abs(p);
If i > 0 Then
Print "x";
If i > 1 Then Print "^" & i;
End If
Next i
Print
End Sub
Type Punto
X As Single
Y As Single
End Type
Sub Lagrange(res() As Single, pts() As Punto)
Dim As Integer i, j, k, c = Ubound(pts)
Dim As Single basis(c, c), temp(1), mulRes(c)
For i = 0 To c
Redim As Single poly(0)
poly(0) = 1
For j = 0 To c
If i <> j Then
temp(0) = -pts(j).X
temp(1) = 1
MulPoly(mulRes(), poly(), temp())
Redim Preserve poly(Ubound(mulRes))
For k = 0 To Ubound(mulRes)
poly(k) = mulRes(k)
Next k
End If
Next j
Dim As Single d = 1
For j = 0 To c
If i <> j Then d *= (pts(i).X - pts(j).X)
Next j
ScalarDivide(poly(), d)
For k = 0 To Ubound(poly)
basis(i, k) = poly(k)
Next k
Next i
Redim res(c)
For i = 0 To c
res(i) = 0
Next i
For i = 0 To c
For k = 0 To c
res(k) += basis(i, k) * pts(i).Y
Next k
Next i
End Sub
Dim As Punto pts(3)
pts(0).X = 1 : pts(0).Y = 1
pts(1).X = 2 : pts(1).Y = 4
pts(2).X = 3 : pts(2).Y = 1
pts(3).X = 4 : pts(3).Y = 5
Dim As Single lip()
Lagrange(lip(), pts())
ShowPoly(lip())
Sleep
- Output:
2.166667x^3 -16x^2 +35.83334x -21
J
From the J wiki polynomial phrases page, we have a couple implementations:
lagrange=: ] +/ .* [ (] % p.~) 1 p.@;~&1\. [
lagrange1=: %.@(^/ i.@#)@[ +/ .* ]
Task example:
1 2 3 4 lagrange 1 4 1 5
_21 215r6 _16 13r6
1 2 3 4 lagrange1 1 4 1 5
_21 35.8333 _16 2.16667
(The J representation of a polynomial is the list of coefficients of the polynomial such that the index of the coefficient is the power applied to the polynomial's variable.)
Java
import java.awt.Point;
import java.util.Arrays;
public final class LagrangeInterpolation {
public static void main(String[] args) {
Point[] points = new Point[] { new Point(1, 1), new Point(2, 4), new Point(3, 1), new Point(4, 5) };
display(lagrangeInterpolation(points));
}
private static double[] lagrangeInterpolation(Point[] points) {
double[][] polys = new double[points.length][points.length];
for ( int i = 0; i < points.length; i++ ) {
double[] poly = new double[] { 1.0 };
for ( int j = 0; j < points.length; j++ ) {
if ( i != j ) {
poly = multiply(poly, new double[] { -points[j].x, 1.0 });
}
}
final double value = evaluate(poly, points[i].x);
polys[i] = scalarDivide(poly, value);
}
double[] sum = new double[] { 0.0 };
for ( int i = 0; i < points.length; i++ ) {
polys[i] = scalarMultiply(polys[i], points[i].y);
sum = add(sum, polys[i]);
}
return sum;
}
// A double[] is used to represents a Polynomial
// with its coefficients reversed compared to the standard mathematical notation.
// For example, the polynomial 3x^2 + 2x + 1 is represented by the array [1, 2, 3].
private static double[] add(double[] one, double[] two) {
double[] sum = Arrays.copyOf(one, Math.max(one.length, two.length));
for ( int i = 0; i < two.length; i++ ) {
sum[i] += two[i];
}
return sum;
}
private static double[] multiply(double[] one, double[] two) {
double[] product = new double[one.length + two.length - 1];
for ( int i = 0; i < one.length; i++ ) {
for ( int j = 0; j < two.length; j++ ) {
product[i + j] += one[i] * two[j];
}
}
return product;
}
private static double[] scalarMultiply(double[] array, double value) {
return Arrays.stream(array).map( d -> d * value ).toArray();
}
private static double[] scalarDivide(double[] array, double value) {
return scalarMultiply(array, 1.0 / value);
}
private static double evaluate(double[] array, double value) {
double result = 0.0;
for ( int i = array.length - 1; i >= 0; i-- ) {
result = result * value + array[i];
}
return result;
}
private static void display(double[] array) {
final int degree = array.length - 1;
if ( degree == 0 ) {
System.out.println(String.format("%2.5f", array[0]));
return;
}
for ( int i = degree; i >= 0; i-- ) {
if ( array[i] == 0.0 ) {
continue;
}
String sign = ( array[i] < 0.0 && i == degree ) ?
"-" : ( array[i] < 0.0 ) ? " - " : ( i < degree ) ? " + " : "";
System.out.print(sign);
final double coeff = Math.abs(array[i]);
if ( coeff > 1.0 ) {
System.out.print(String.format("%2.5f", coeff));
}
String term = ( i > 1 ) ? "x^" + String.valueOf(i) : ( i == 1 ) ?
"x" : ( coeff == 1.0 ) ? "1" : "";
System.out.print(term);
}
System.out.println();
}
}
- Output:
2.16667x^3 - 16.00000x^2 + 35.83333x - 21.00000
jq
Works with gojq, the Go implementation of jq
With minor modifications related to the module system, the program presented here also works with jaq, the Rust implementation of jq
The function `lagrange/1` is adapted from Wren.
This entry uses the "polynomials" module in which a polynomial is represented by the JSON array comprised of the polynomial's coefficients, with the entry at .[i] corresponding to the coefficient of x^i. The canonical form of a polynomial of degree n is represented by an array of length n+1.
include "polynomials" {search: "./"};
# Returns the Lagrange interpolating polynomial which passes through a list of points.
def lagrange($pts):
($pts|length) as $c
| reduce range(0;$c) as $i ({polys: []};
.poly = [1]
| reduce range(0;$c) as $j (.;
if ($i != $j)
then .poly = (multiply(.poly; [-$pts[$j][0], 1]))
else .
end )
| (.poly | eval($pts[$i][0])) as $d
| .polys[$i] = (.poly | scalarDivide($d)) )
| reduce range(0;$c) as $i (.sum = [0];
.polys[$i] = (.polys[$i] | scalarMultiply($pts[$i][1]))
| .sum = add(.sum; .polys[$i]) )
| .sum ;
def pts: [
[1, 1],
[2, 4],
[3, 1],
[4, 5]
];
lagrange(pts) | pp- Output:
2.1666666666666665x³-16.0x²+35.83333333333333x-21.0
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)
Perl
# 20241002 Perl programming solution
use strict;
use warnings;
use Math::Polynomial;
use Math::BigRat;
use PDL;
sub Lagrange {
my $pts = shift;
my ($xs, $cs) = map {
my $outer = $_;
[ map { Math::BigRat->new($_->[$outer]) } @$pts ]
} (0, 1);
my $n = scalar @$xs;
return Math::Polynomial->new($cs->[0]) if $n == 1;
my @arr = (Math::BigRat->bone) x $n;
for my $j (1 .. $n - 1) {
for my $k (0 .. $j - 1) { $arr[$k] *= $xs->[$k] - $xs->[$j] }
$arr[$j] = Math::BigRat->bone;
$arr[$j] *= $_ for map { $xs->[$j] - $xs->[$_] } (0 .. $j - 1);
}
my @ws = map { Math::BigRat->bone / $_ } @arr;
my ($q, $x) = (Math::Polynomial->new(0), Math::Polynomial->new(0, 1));
for my $i (0 .. $n - 1) {
my $m = Math::Polynomial->new(1);
for my $j (0 .. $n - 1) {
next if $j == $i;
$m *= $x - Math::Polynomial->new($xs->[$j]);
}
$q += $m * $ws[$i] * Math::Polynomial->new($cs->[$i]);
}
$q->string_config({ fold_sign => 1 });
return $q
}
my $testpoints = [[1, 1], [2, 4], [3, 1], [4, 5]];
print "Lagrange(testpoints) = ", Lagrange($testpoints), "\n";
- Output:
Lagrange(testpoints) = (13/6 x^3 - 16 x^2 + 215/6 x - 21)
Phix
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
Raku
# 20231020 Raku programming solution
class Point { has Rat ($.x, $.y) }
sub add(@p1, @p2) { return @p1 <<+>> @p2 } # Add two polynomials.
sub multiply(@p1, @p2) { # Multiply two polynomials.
my @prod;
for ^+@p1 X ^+@p2 -> ($i, $j) { @prod[$i + $j] += @p1[$i] * @p2[$j] }
return @prod;
}
# Multiply a polynomial by a scalar.
sub scalarMultiply(@poly, $x) { return @poly.map: * × $x }
# Divide a polynomial by a scalar.
sub scalarDivide(@poly, $x) { return scalarMultiply(@poly, 1/$x) }
# rosettacode.org/wiki/Horner%27s_rule_for_polynomial_evaluation#Raku
sub evalPoly(@coefs, $x) { return ([o] map { $_ + $x × * }, @coefs.reverse)(0) }
sub lagrange(@pts) {
my ($c, @polys) = @pts.elems;
for ^$c -> $i {
my @poly = 1;
for ^$c -> $j {
next if $i == $j;
@poly = multiply @poly, [ -@pts[$j].x, 1 ]
}
@polys[$i] = scalarDivide @poly, evalPoly(@poly.reverse, @pts[$i].x)
}
my @sum = 0;
for ^$c -> $i { @sum = add @sum, scalarMultiply @polys[$i], @pts[$i].y }
return @sum;
}
my @pts = map { Point.new: x=>.[0].Rat,y=>.[1].Rat }, [<1 1>,<2 4>,<3 1>,<4 5>];
say [~] lagrange(@pts).kv.rotor(2).reverse.map: -> ($expo, $coef) {
"{ '+' if $coef ≥ 0 }$coef.nude.join('/') " ~ do given $expo {
when 0 { " " }
when 1 { "𝑥 " }
default { "𝑥^$_ " }
}
}
You may Attempt This Online!
RPL
[[1 2 3 4][1 4 1 5]] LAGRANGE
- Output:
1: '(13*X^3-96*X^2+215*X-126)/6'
Wren
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