Generating generating expressions for sequences

Question

(yes, "generating generating" in the title is correct :) )

Context

In middle (?) school we are taught about sequences and, in particular, we are taught about linear sequences where the nth term is generated with an expression of the form an + b, where a and b are some coefficients. In this challenge, we will deal with sequences generated by polynomials of arbitrary degree.

Task

Given the first m terms of a sequence, find the coefficients of the polynomial of lowest degree that could have generated such a sequence.

A polynomial, and thus the generating expression you are looking for, is to be seen as a function \$p(n)\$ that takes n as an argument and returns

$$a_0 + a_1 n + a_2 n^2 + a_3 n^3 + \cdots + a_k n^k$$

where \$k \geq 0\$ and \$a_i, 0 \leq i \leq k\$ have to be found by you.

You will assume that the m terms you were given correspond to taking n = 0, n = 1, ..., n = m-1 in the generating polynomial above.

Examples

If I am given the sequence [2, 2, 2] then I realize this is a constant sequence and can be generated by a polynomial of degree 0: p(n) = 2.

If I am given the sequence [1, 2, 3] then I realize this cannot come from a constant polynomial but it could come from a linear polynomial p(n) = n + 1, so that is what my output should be. Notice how

p(0) = 1
p(1) = 2
p(2) = 3    # and NOT p(1) = 1, p(2) = 2, p(3) = 3

Input

Your input will be the first terms of a sequence, which you can take in any reasonable format/data type. A standard list is the most obvious choice.

You may assume the input sequence is composed of integers (positive, 0 and negative).

Output

The coefficients of the polynomial of lowest degree that could have generated the input sequence. The output format can be in any sensible way, as long as the coefficients can be retrieved unambiguously from the output. For this, both the value of each coefficient and the degree of each coefficient are important. (e.g. if using a list, [1, 0, 2] is different from [0, 1, 2]).

You can assume the polynomial you are looking for has integer coefficients.

Test cases

For these test cases, the input is a list with the first terms; the output is a list of coefficients where (0-based) indices represent the coefficients, so [1, 2, 3] represents 1 + 2x + 3x^2.

[-2] -> [-2]
[0, 0] -> [0]
[2, 2, 2] -> [2]
[4, 4] -> [4]
[-3, 0] -> [-3, 3]
[0, 2, 4, 6] -> [0, 2]
[2, 6] -> [2, 4]
[3, 7] -> [3, 4]
[4, 8, 12, 16] -> [4, 4]
[-3, -1, 5, 15, 29] -> [-3, 0, 2]
[0, 1, 4, 9] -> [0, 0, 1]
[3, 2, 3, 6, 11] -> [3, -2, 1]
[3, 4, 13, 30, 55] -> [3, -3, 4]
[4, 12, 28, 52, 84] -> [4, 4, 4]
[2, 4, 12, 32, 70] -> [2, 1, 0, 1]
[3, 6, 21, 54] -> [3, -1, 3, 1]
[4, 2, 12, 52, 140] -> [4, -2, -3, 3]
[10, 20, 90, 280] -> [10, 0, 0, 10]
[-2, 8, 82, 352, 1022, 2368, 4738] -> [-2, 4, -1, 4, 3]
[4, 5, 32, 133, 380] -> [4, -2, 0, 2, 1]
[1, 0, 71, 646, 2877, 8996, 22675] -> [1, -1, 0, -3, 0, 3]
[4, 2, 60, 556, 2540, 8094, 20692] -> [4, -2, -1, 0, -2, 3]
[1, 2, -17, 100, 1517, 7966, 28027, 78128, 186265] -> [1, 3, -2, 4, -3, -2, 1]
[4, 5, 62, 733, 4160, 15869, 47290, 118997] -> [4, 3, -1, -3, 1, 0, 1]

Test cases generated with this code

---

This is code-golf so shortest submission in bytes, wins! If you liked this challenge, consider upvoting it! If you dislike this challenge, please give me your feedback. Happy golfing!

Answer 1

JavaScript (ES7),  193 ... 154  145 bytes

Saved 9 bytes thanks to @Bubbler

Returns \$(a_0,a_1,...,a_k)\$ with some possible trailing zeros.

v=>v.map((_,i)=>(g=(i,m=v.map((n,y)=>v.map((_,x)=>x==i?n:y**x)))=>+m||m.reduce((s,[v],i)=>v*g(0,m.map(([,...r])=>r).filter(_=>i--))-s,0))(i)/g())

Try it online!

(removed the penultimate test case, which requires more precision than IEEE-754 provides)

How?

We use Cramer's rule to solve a system of linear equations based on a square Vandermonde matrix:

1. Given an input vector of length \$n\$, we build a Vandermonde matrix \$V_n\$ of size \$n\times n\$ with coefficients \$\alpha_i=i,0\le i <n\$:

   $$Vn=\begin{pmatrix} 1&0&0&...&0\\ 1&1&1&...&1\\ 1&2&4&...&2^{n-1}\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ 1&n-1&(n-1)^2&...&(n-1)^{n-1} \end{pmatrix}$$

2. Using Cramer's rule, the coefficient \$a_i\$ of the polynomial is computed by taking the determinant of the matrix obtained by replacing the \$i\$-th column of \$V_n\$ with the input vector, and dividing by the determinant of \$V_n\$.

Example for \$(4,2,12,52,140)\$

The constant coefficient \$a_0\$ is given by:

$$a_0=\begin{vmatrix} \color{blue}4&0&0&0&0\\ \color{blue}2&1&1&1&1\\ \color{blue}{12}&2&4&8&16\\ \color{blue}{52}&3&9&27&81\\ \color{blue}{140}&4&16&64&256 \end{vmatrix}/|V_5|=\frac{1152}{288}=4$$

The coefficient \$a_1\$ is given by:

$$a_1=\begin{vmatrix} 1&\color{blue}4&0&0&0\\ 1&\color{blue}2&1&1&1\\ 1&\color{blue}{12}&4&8&16\\ 1&\color{blue}{52}&9&27&81\\ 1&\color{blue}{140}&16&64&256 \end{vmatrix}/|V_5|=\frac{-576}{288}=-2$$

And so on.

Answer 2

R, 55 52 bytes

-3 bytes thanks to Giuseppe.

round(solve(outer(n<-seq(a=u<-scan())-1,n,"^"))%*%u)

Try it online!

Outputs \$(a_0, a_1,\ldots,)\$ with possible trailing zeros.

Let \$u\$ be the output sequence, and \$X\$ be the \$m\times m\$ matrix such that \$X_{i,j}=i^j\$ (0-indexed), i.e.

\$ X=\begin{pmatrix} 1&0&0&\ldots&0\\ 1&1&1&\ldots&1\\ 1&2&4&\ldots&2^{m-1}\\ 1&3&9&\ldots&3^{m-1}\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ 1&m-1&(m-1)^2&\ldots&(m-1)^{m-1} \end{pmatrix}. \$

Then in matrix notation, \$u=Xa\$, hence \$a=X^{-1}u\$.

The code implements this: n is the vector (0, 1, ..., m-1) where m is the length of u; this is used to construct X = outer(n, n, "^"). The function solve performs matrix inversion, and the round is there to avoid numerical errors.

Answer 3

APL+WIN, 16 bytes

Index origin = 0

Prompts for input as a vector and outputs coefficients from a0 to an-1 where n is the length of the vector. The order of the polynomial can be obtained by summing the number of coefficients up to the last none zero coefficient:

0⍕n⌹m∘.*m←⍳⍴n←,⎕

Try it online! Courtesy of Dyalog Classic

Answer 4

Wolfram Language (Mathematica), 50 49 37 bytes

Returns a polynomial.

Mathematica is so awesome x+1 can be used as a variable in this context. Apart is a weird built-in that, quoting from the docs, seems to attempt to rewrite an expression as a sum of terms with minimal denominators, and also happens to expand polynomials (that are returned in a weird collapsed form by default) into something more sane.

Apart@InterpolatingPolynomial[#,x+1]&

Try it online!

Sledgehammer, 8 bytes

^{(it will try to deceive you into thinking it's actually 7.5, but it's actually not)}

⣕⢤⣏⠛⡪⣊⠵⢼

Explanation: It's Apart@InterpolatingPolynomial[Input[], x+1], but compressed via an awesome Mathematica compressor (it is so awesome that, as far as I understand, it translates Mathematica to an intermediate stack-based language).

Unfortunately, running this is fairly painful.

Answer 5

J, 10 bytes

%.^/~@i.@#

Try it online!

Obligatory J answer on a matrix-related challenge. Takes input as a vector of extended integers (otherwise the answer may have small floating-point errors), and gives the polynomial's coefficients in lowest-first order, possibly with some extra zeroes at the end.

How it works

%.^/~@i.@#  NB. Input: a vector V of extended integers.
         #  NB. Length of V
      i.@   NB. Generate 0..(len(V)-1)
  ^/~@      NB. Self outer product by ^(exponentiation)
%.          NB. Matrix-divide V by the matrix above,
            NB.   i.e. solve a linear system of equations

Answer 6

Pari/GP, 38 bytes

a->Vecrev(polinterpolate([0..#a-1],a))

Try it online!

Answer 7

Python 3 + Numpy, 69 bytes

lambda x:polyfit(range(len(x)),x,len(x)-1).round()
from numpy import*

Try it online!

May have leading zeros.

Answer 8

Charcoal, 68 62 bytes

≔⟦¹⟧ηＦＬθ«⊞υ⁰≔÷⁻§θιΣＥυ×κＸιλ∨ΠＥι⊕κ¹ζＵＭυ⁺κ×ζ§ηλ⊞η⁰≔Ｅη⁻§η⊖λ×κιη»Ｉυ

Try it online! Link is to verbose version of previous version of code that excludes trailing zeros, but apparently it isn't necessary to do that, thus saving 6 bytes. Outputs the terms in power order i.e. the constant term is printed first. Explanation:

≔⟦¹⟧η

Start by creating a helper polynomial \$ h(x) = 1 \$.

ＦＬθ«

Loop over the \$ m \$ terms.

⊞υ⁰

Add a \$ 0x^i \$ term to the result polynomial \$ u(x) \$.

≔÷⁻§θιΣＥυ×κＸιλ∨ΠＥι⊕κ¹ζ

Subtract the value of \$ u(i) \$ from the input term and divide that by \$ i! \$.

ＵＭυ⁺κ×ζ§ηλ

Multiply \$ h \$ by that value and add the result to \$ u \$. This doesn't change the values of \$ u(0) ... u(i-1) \$ but the value of \$ u(i) \$ is now the input term.

⊞η⁰≔Ｅη⁻§η⊖λ×κιη

Multiply \$ h \$ by \$ x - i \$.

»Ｉυ

Print the coefficients of \$ u \$, which may include trailing zeros.

Answer 9

Haskell, 77 bytes

h%(a:t)=h-a:a%t
h%_=[h]
f(h:t)=h:foldr(%)[](f$zipWith((/).(-h+))t[1..])
f e=e

Try it online!

Answer 10

APL (Dyalog Unicode), 10 bytes^SBCS

⊢⌹∘.*⍨∘⍳∘≢

Try it online!

A port of Graham's APL+WIN solution into a modern APL, which happens to work exactly the same (and have the same byte count) as my own J solution.

How it works

⊢⌹∘.*⍨∘⍳∘≢  ⍝ Input: V, result of a polynomial evaluated at 0..m-1
       ⍳∘≢  ⍝ Generate 0..m-1
  ∘.*⍨∘     ⍝ Self outer product by * (exponentiation)
⊢⌹          ⍝ Matrix divide V by above (solve linear system of equations)

Answer 11

05AB1E, 48 47 bytes

g≠iā<DδmUεXøINǝ}Xšεā<sUœε©2.ÆíÆ.±Xε®Nèè}«P}O}ć÷

Sometimes 05AB1E's lack of almost all matrix builtins is pretty annoying.. ;)
Inspired by @Arnauld's JavaScript answer.

Try it online or verify almost all test cases (removed the last two largest ones, since they time out on TIO).

Explanation:

First handle the edge case of a single-element input-list (would cause issues with the « later on in the code):

g                # Get the length of the (implicit) input-list
 ≠i              # And if it is NOT 1, continue with:
                 #  ... (see below)
                 # (implicit else:)
                 #  (output the implicit input-list as implicit output)

Next we'll get the exponentiation matrix of the list [0, input-length):

ā                #  Push a list in the range [1, (implicit) input-length] (without popping)
 <               #  Decrease each value by 1 to make the range [0, input-length)
  Dδ             #  Apply double-vectorized on itself by first duplicating:
    m            #   Take the power of the two values
     U           #  Pop and store this exponentiation matrix in variable `X`

Next we'll create a list of this matrix, with every column one by one replaced with the input-list:

ε     }          #  Map over the input-list that was still on the stack
 X               #   Push the exponentiation matrix from variable `X`
  ø              #   Zip/transpose it; swapping rows/columns
     ǝ           #   Replace the transposed row of the exponentiation matrix
    N            #   at the current map-index
   I             #   with the input-list

We'll prepend the original exponentiation matrix to this list:

Xš               #  Prepend the matrix `X` in front of this list

And we'll calculate the determinant of each inner matrix in this list:

ε              } #  Map over the list of matrices:
 ā               #   Push a list in the range [1, matrix-length] (without popping)
  <              #   Decrease it by 1 to make the range [0, matrix-length)
   sU            #   Swap to get the matrix again, and pop and store it in variable `X`
     œ           #   Get all permutations of the [0, matrix-length) list
      ε          #   Inner map over each permutation:
       ©         #    Store the current permutation in variable `®` (without popping)
        2.Æ      #    Get all 2-element combinations of this permutation
           í     #    Reverse each inner pair
            Æ    #    Reduce it by subtracting
             .±  #    And get it's signum (-1 if a<0; 0 if a==0; 1 if a>0)
       X         #    Push the matrix from variable `X`
        ε        #    Map over each of its rows:
         ®       #     Push the current permutation of variable `®`
          Nè     #     Get the value in the permutation at the current map-index
            è    #     And use that to index into the current matrix-row
        }«       #    After the map of rows: merge it together with the signum list
          P      #    And take the product of this entire list
      }O         #   After the map of permutations: sum all values

Now that we have all determinants of the matrices, we get the default one again to divide all others by it:

ć                #  Extract head: pop and push remainder-list and first item separated
 ÷               #  Integer-divide each value in the remainder-list by this head
                 #  (after which the result is output implicitly)

Answer 12

MATL, 12 bytes

n:qGyz3$ZQYo

The result is given with higher-order coefficients first, and may contain leading zeros.

Try it online! Or verify all test cases

Explanation

Consider input [-3, -1, 5, 15, 29] as an example.

n:q    % Implicit input. Number of elements. Range. Subtract 1, element-wise
       % STACK: [0, 1, 2, 3, 4]
G      % Push input again
       % STACK: [0, 1, 2, 3, 4], [-3, -1, 5, 15, 29]
yz     % Duplicate from below. Number of non-zero elements
       % STACK: [0, 1, 2, 3, 4], [-3, -1, 5, 15, 29], 4
3$ZQ   % Fit polynomial with inputs x, y, degree
       % STACK: [3.7536e-16, -3.1637e-15, 2.0000, -8.8363e-15, -3]
Yo     % Round, element-wise. Implicit display
       % STACK: [0, 0, 2, 0, -3]

Answer 13

SageMath, 63 48 bytes

lambda v:QQ[x].lagrange_polynomial(enumerate(v))

Try it online!

Outputs the polynomial as

$$a_k n^k + \cdots + a_3 n^3 + a_2 n^2 + a_1 n + a_0 $$

Answer 14

Jelly, 14 bytes

J’*þ`æ*-⁸æ×ær0

A monadic Link accepting a list of integers which yields a list of the exponents (floats and/or integers) with the lowest degree on the left of the same length as the input (with trailing zeros if need be).

Try it online!

How?

J’*þ`æ*-⁸æ×ær0 - Link: list of integers, V
J              - range of length (V)
 ’             - decrement (vectorises)
    `          - use as both arguments of:
   þ           -   outer-product using:
  *            -     exponentiation
       -       - minus one
     æ*        - matrix-exponentiation (i.e. inverse)
        ⁸      - chain's left argument, V
         æ×    - matrix-multiplication
           ær0 - round to zero decimal places (vectorises)