# Hermite interpolation

We already have a challenge for polynomial interpolation: given a list of points, output the coefficients of the polynomial that passes through them.

Hermite interpolation is a generalization of polynomial interpolation. Instead of just specifying the points that the polynomial passes through, we also specify the higher-order derivatives at those points. If we have $$\n\$$ points, and we specify the value and the first $$\m-1\$$ derivatives at each point, then we can uniquely determine a polynomial of degree less than $$\mn\$$.

For example, the following data uniquely determines the polynomial $$\x^8+1\$$:

$$\x\$$ $$\y=x^8+1\$$ $$\y'=8x^7\$$ $$\y''=56x^6\$$
0 1 0 0
1 2 8 56
2 257 1024 3584

Given an $$\m \times n\$$ matrix $$\\mathbf{A}=(a_{ij})_{1\leq i\leq m, \; 1\leq j\leq n}\$$, find a polynomial $$\p\$$ of degree less than $$\mn\$$ such that $$\p^{(i)}(j) = a_{i+1,j+1}\$$ for all $$\0 \leq i < m\$$ and $$\0 \leq j < n\$$. Here $$\p^{(i)}(j)\$$ denotes the $$\i\$$th derivative of $$\p\$$ evaluated at $$\j\$$.

For example, if $$\\mathbf{A}=\begin{bmatrix}1&0&0\\2&8&56\\257&1024&3584\end{bmatrix}\$$, then $$\p(x)=x^8+1\$$.

You may take input and output in any convenient format.

You may transpose the input matrix if you wish. You may flatten it into a list. You may take $$\m\$$ or $$\n\$$ or both as additional inputs.

Here are some example formats for the output polynomial:

• a list of coefficients, in descending order, e.g. $$\24 x^4 + 96 x^3 + 72 x^2 + 16 x + 1\$$ is represented as [24,96,72,16,1];
• a list of coefficients, in ascending order, e.g. $$\24 x^4 + 96 x^3 + 72 x^2 + 16 x + 1\$$ is represented as [1,16,72,96,24];
• taking an additional input $$\k\$$ and outputting the coefficient of $$\x^k\$$;
• a built-in polynomial object.

When the degree of the polynomial is less than $$\mn-1\$$, you may pad the coefficient list with zeros.

When the coefficients are not integers, you may output them as rational numbers, floating-point numbers, or any other convenient format.

This is , so the shortest code in bytes in each language wins.

## Test Cases

Here I output lists of coefficients in descending order.

[[1,0,0],[1,0,0]] -> 
[[1,2],[0,0]] -> [4,-7,2,1]
[[1,2],[3,4]] -> [2,-2,2,1]
[[0,0,1],[3,11,33]] -> [1,1/2,1,1/2,0,0]
[[1,2],[3,4],[5,6]] -> [3,-14,21,-10,2,1]
[[1,2,3],[4,5,6]] -> [-3/2,9/2,-7/2,3/2,2,1]
[[1,16],[209,544],[1473,2224]] -> [24,96,72,16,1]
[[0,0,0],[0,1,0],[0,0,0]] -> [-1,7,-18,20,-8,0,0,0]
[[1,0,0],[2,8,56],[257,1024,3584]] -> [1,0,0,0,0,0,0,0,1]
[[1,0,0],[0,0,0],[0,0,0]] -> [3/2,-207/16,713/16,-1233/16,1083/16,-99/4,0,0,1]


# Mathematica, 76 bytes

Mathematica has built-in InterpolatingPolynomial to construct an interpolating polynomial that reproduces function values and derivatives.

Try it online!

f@m_:=InterpolatingPolynomial[MapIndexed[Join[{#2[]-1},#1]&,m],x]//Expand


# Python3, 604 bytes:

import math
E=enumerate
def M(p,*P):
for i in P:
D={}
for a,b in i.items():
for A,B in p.items():D[a+A]=D.get(a+A,0)+b*B
p=D
return p
def O(p,*P):
for i in P:
for j,k in i.items():p[j]=p.get(j,0)+k
return p
def f(A):
z,l,d=zip(*[(i,a,a[1:])for i,a in E(A)for _ in a])
L=[[(a,[i])for i,a in E(l)]]
c=1
while len(U:=L[-1])>1:L+=[[(d[min(U[i])][c-1]/math.factorial(c)if 0==(D:=z[max(U[i+1]+U[i])]-z[min(U[i+1]+U[i])])else(U[i+1]-U[i])/D,U[i]+U[i+1])for i in range(len(U)-1)]];c+=1
return O(*[M({0:a},*[{1:1,0:-z[j]}for j in range(i)])for i,a in E(L)])


Try it online!

• Wow, Python is so verbose that import math and math.factorial is the golfiest way to take the factorial?
– Neil
Oct 27 at 10:14
• @Neil it was late and I could not be bothered to implement it myself :) Oct 27 at 14:26
• I=dict.items lets you do I(i) instead of i.items() and saves 2 bytes overall. Oct 27 at 18:55

# Charcoal, 99 94 bytes

ＦＬθＦ⌊θ⊞υι≔Ｅυ¬κη≔×⁰ηζＦＬυ«≔Ｅυ⎇⁼κ§υ⁺λι∕§§θκι∨Π…·¹ι¹∕⁻§ε⊕λ§ελ⁻§υ⁺λικεＵＭζ⁺κ×§ε⁰§ηλＵＭη⁺§η⊖λ×±§υικ»Ｉζ


Try it online! Link is to verbose version of code. Outputs the polynomial coefficients in ascending order of degree. Explanation:

ＦＬθＦ⌊θ⊞υι


Create a list of x-coordinates xᵢ where each coordinate is repeated by the number of times equal to the number of differentials provided.

≔Ｅυ¬κη


Create a unit polynomial.

≔×⁰ηζ


Create a zero polynomial, which will be the running total.

ＦＬυ«


Loop over each row of the divided difference table.

≔Ｅυ⎇⁼κ§υ⁺λι∕§§θκι∨Π…·¹ι¹∕⁻§ε⊕λ§ελ⁻§υ⁺λικε


Calculate the row. For each cell, if the x-coordinates of the triangle that generates this cell are all equal then divide the appropriate derivative of that coordinate by the factorial of the row number, otherwise divide the difference of the two cells from the previous row by the difference of the x-coordinates. (Some additional bogus entries are also calculated but none of them affect the final result.)

ＵＭζ⁺κ×§ε⁰§ηλ


Multiply the (originally unit) polynomial by the first entry in the row and add it to the running total.

ＵＭη⁺§η⊖λ×±§υικ


Multiply the (originally unit) polynomial by x-xᵢ. (This generates a bogus result on the last pass but since it won't be used this is irrelevant.)

»Ｉζ


Output the final total polynomial.

97 bytes for a version that accepts a list of x-coordinates as an additional argument, rather than assuming that they're 0-indexed:

ＦηＦ⌊θ⊞υι≔Ｅυ¬κη≔×⁰ηζＦＬυ«≔Ｅυ⎇⁼κ§υ⁺λι∕§§θ÷λＬ⌊θι∨Π…·¹ι¹∕⁻§ε⊕λ§ελ⁻§υ⁺λικεＵＭζ⁺κ×§ε⁰§ηλＵＭη⁺§η⊖λ×±§υικ»Ｉζ


Try it online! Link is to verbose version of code.

# PARI/GP, 52 bytes

f(a,m)=a/matrix(#a,,k,i,derivn(x^k--,i--%m)%(x-i\m))


Attempt This Online!

Takes input as a flatten list and $$\m\$$. Outputs a list of coefficients in ascending order.

Using the fact that Hermite interpolation is a linear map. We can construct a $$\mn\times mn\$$ matrix, whose entries are the $$\i\$$th derivative of $$\x^k\$$ evaluated at $$\x=j\$$, where $$\0\leq i< m\$$, $$\0\leq j< n\$$, and $$\0\leq k< mn\$$. Then we can multiply the flattened input by the inverse of this matrix.

# R, 149 bytes

\(y,n,m,o=m*n-1,z=y,v={},^=c){for(i in 0:o){v=v^0+(z=ifelse(x<-i:o%/%n-(j=0:(o-i)%/%n),diff(z)/x,y[i+1,j+1]/factorial(i)))*T;T=0^T-T^0*i%/%n};v}


Attempt This Online!

A function taking the matrix, n and m as arguments and returning a vector of coefficients, starting from the zeroth.

# Jelly, 71 63 bytes

I÷ḢṛṪḢ?@"Z}
FJḶṚ1ị+",Ʋ:L_/żṪ©ị"¥ɗ÷"J’!Ɗ$Żç\Ḣ€Ḋ×"®Ḣ1;$Ż_×ɗ\Ṗ¤F€S


Try it online!

A pair of monadic links that takes a list of lists and returns a list of polynomial coefficients.

The argument holding the input is transposed and then rotated both vertically by 1 such that the 0th column appears at the end. This is done to best fit with Jelly’s 1-indexing and the fact that zero refers to the last column. The output starts with the zeroth coefficient and works upwards.

Updated to a more efficient algorithm, but still based on the one on the Wikipedia page.

## Explanation

I÷ḢṛṪḢ?@"Z}                                          # ‎⁡Helper link: takes the current list of values for the divided difference algorithm on the left, and a list of the x differences and potential relevant derivatives on the right and returns the next column for the divided difference algorithm
I                                                    # ‎⁢Increments
÷                                                   # ‎⁣Divided by the right argument
Ḣ                                                  # ‎⁤Head (effectively gives the differences divided by the x differences)
@"Z}                                          # ‎⁢⁡Take each pair of x differences and possible relevant derivatives as the left argument and the relevant divided difference as the right and do the following for each:
Ḣ?                                              # ‎⁢⁢- If the x-difference was non-zero:
ṛ                                                 # ‎⁢⁣- Then the divided difference
Ṫ                                                # ‎⁢⁤- Else the relevant derivative

FJḶṚ1ị+",Ʋ:L_/żṪ©ị"¥ɗ÷"J’!Ɗ$Żç\Ḣ€Ḋ×"®Ḣ1;$Ż_×ɗ\Ṗ¤F€S  # ‎⁣⁡Main link: takes a matrix as described above and returns a list of coefficients
F                                                    # ‎⁣⁢Flatten
J                                                   # ‎⁣⁣Sequence along flattened matrix from 1 to length
Ḷ                                                  # ‎⁣⁤For each of these, sequence from 0 to value minus 1
Ṛ                                                 # ‎⁤⁡Reverse
Ʋ                                           # ‎⁤⁢Following as a monad (argument z):
1ị                                               # ‎⁤⁣- First list within z
+"                                             # ‎⁤⁤- Add zipped to z (so each member of first list in z increases by zero, each member of second list of z increases by 1, etc.)
,                                            # ‎⁢⁡⁡- Pair with z
:L                                         # ‎⁢⁡⁢Integer divide by length of matrix (m)
ɗ                                # ‎⁢⁡⁣Following as a dyad (see below for right argument details, referred to as y)
_/                                       # ‎⁢⁡⁤- Reduce by subtraction (gives x differences)
ż                                      # ‎⁢⁢⁡- Zip with
¥                                 # ‎⁢⁢⁢- Following as a dyad:
Ṫ©                                    # ‎⁢⁢⁣  - Tail (effectively z from above but after the integer divide bit), also copying this for later
ị"                                  # ‎⁢⁢⁤  - Index zipped into right argument
$# ‎⁢⁣⁡y, which will be the result of the following applied as a monad to the original argument (the matrix) ÷" Ɗ # ‎⁢⁣⁢- Divide zipped by the result of the following applied as a monad to the original matrix: J # ‎⁢⁣⁣ - Sequence along original matrix rows ’ # ‎⁢⁣⁤ - Decrease by 1 ! # ‎⁢⁤⁡ - Factorial Ż # ‎⁢⁤⁢Prepend zero ç\ # ‎⁢⁤⁣Reduce using the helper link, collecting all intermediate values Ḣ€ # ‎⁢⁤⁤Head of each Ḋ # ‎⁣⁡⁡Remove the original zero ×"® ¤ # ‎⁣⁡⁢Multiply zipped by the result of applying the following to the register value stored earlier Ḣ # ‎⁣⁡⁣- Head 1;$            # ‎⁣⁡⁤- Prepend 1 (\$ used to prevent the ¤ only tracking back to the 1)
ɗ\       # ‎⁣⁢⁡- Reduce using the following, collecting up intermediate values
Ż           # ‎⁣⁢⁢  - Prepend zero
_×         # ‎⁣⁢⁣  - Subtract the result of multiplying the previous list by the incoming x value
Ṗ      # ‎⁣⁢⁤- Remove the last list member
F€   # ‎⁣⁣⁡Flatten each (deals with the 1 introduced above that is not a list)
S  # ‎⁣⁣⁢Sum
💎


Created with the help of Luminespire.