In this task you'll be given a set of N points (x1,y1),…,(xN,yN) with distinct xi values and your task is to interpolate a polynomial through these points. If you know what Lagrange interpolation is you can skip this section.

The goal of a polynomial interpolation is to construct the (unique) polynomial p(x) with degree N-1 (for higher degrees there are infinite solutions) for which p(xi) = yi for all i = 1…N. One way to do so is to construct N Lagrange basis polynomials and form a linear combination. Such a basis polynomial is defined as follows:

                          l_i(x) := \cfrac{x-x_1}{x_i-x_1}\cdots\cfrac{x-x_{i-1}}{x_i-x_{i-1}}\cdot\cfrac{x-x_{i+1}}{x_i-x_{i+1}}\cdots\cfrac{x-x_N}{x_i-x_N}

As you can see if you evaluate li at the points x1,…,xi-1 ,xi+1,… ,xN it is 0 by construction and 1 for xi, multiplying by yi will only change the value at xi and set it to yi. Now having N such polynomials that are 0 in every point except one we can simply add them up and get the desired result. So the final solution would be:

                                                        p(x) = \overset{N}{\underset{i=1}{\sum}} y_i \cdot l_i(x)


  • the input will consist of N data points in any reasonable format (list of tuples, Points, a set etc.)
  • the coordinates will all be of integer value
  • the output will be a polynomial in any reasonable format: list of coefficients, Polynomial object etc.
  • the output has to be exact - meaning that some solutions will have rational coefficients
  • formatting doesn't matter (2/2 instead of 1 or 1 % 2 for 1/2 etc.) as long as it's reasonable
  • you won't have to handle invalid inputs (for example empty inputs or input where x coordinates are repeated)


These examples list the coefficients in decreasing order (for example [1,0,0] corresponds to the polynomial x2):

[(0,42)] -> [42]
[(0,3),(-18,9),(0,17)] -> undefined (you don't have to handle such cases)
[(-1,1),(0,0),(1,1)] -> [1,0,0]
[(101,-42),(0,12),(-84,3),(1,12)] -> [-4911/222351500, -10116/3269875, 692799/222351500, 12]

3 Answers 3


Pari/GP, 14 bytes


Takes two lists [x1, ..., xn] and [y1, ..., yn] as input. Outputs a polynomial.

Try it online!

Pari/GP, 35 bytes, without built-in


Takes two lists [x1, ..., xn] and [y1, ..., yn] as input. Outputs a list of coefficients.

Try it online!


Enlist, 25 bytes


Try it online!

Monadic link takes an array of length 2 (e.g., [[101,0,-84,1],[-42,12,3,12]]) as input. The first element of the array is list of x value, the second is list of y value. Output as list of coefficient in increasing degree.

Because @HyperNeutrino forgot to wrap numbers in rational wrapper if it is evaluated in list, (Try it online!) the argument must be hardcoded in the code.

The is necessary because HyperNeutrino make some mistake in implementing · I don't know (or is this a feature? Not sure)

Equivalent Mathematica code:

Expand[ Table[
 Times @@ ((z - #)/(xi - #) & /@ DeleteCases[x, xi])
 , {xi, x}] . y]

Try it online!

where x and y are list of coordinates. This prints the polynomial.

What Enlist have for this challenge:

  • Built in rational number support
  • Convolution æc (my feature request)
  • 1
    \$\begingroup\$ If you need a rational number version of Jelly, M might be worth a shot \$\endgroup\$ Dec 23, 2017 at 0:47

Python, with numpy, 301 bytes

See this notebook for how it works.

Golfed version. 301 bytes, it can be golfed much more. Attempt this online!

import numpy as np
from itertools import combinations as C
from numpy.polynomial import Polynomial as P
def l(k,x):return P([sum([np.prod(x) for x in C(np.delete(x,k),i)])*(-1)**(i+1)/np.prod(np.delete(x,k)-x[k]) for i in range(len(x))][::-1])
def f(x,y):return sum(l(k,x)*y[k] for k in range(len(y)))

Ungolfed version. Attempt this online!

import numpy as np
from itertools import combinations 
from numpy.polynomial import Polynomial

# get the Lagrange basis function
def l(k, x):
    n = len(x)
    assert (k < len(x))
    x_k = x[k]
    x_copy = np.delete(x, k)
    denominator = np.prod(x_copy - x_k)
    coeff = []
    for i in range(n):
        coeff.append(sum([np.prod(x) for x in combinations(x_copy, i)]) * (-1)**(i+1) / denominator)
    return Polynomial(coeff)

def f(x,y):return sum(l(k,x)*y[k] for k in range(len(y)))

print("y[0]*l(0,x)+y[1]*l(1,x)+y[2]*l(2,x)+y[3]*l(3,x)== f(x):")

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.