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!

• RGS: related
• Robin Ryder: May we take m as part of the input?
• RGS: @RobinRyder only if it is standard in your language to take the length of a vector as an argument alongside the vector (like it is in C, for example)
• agtoever: Some answers round the resulting coefficients to integers, although I don’t read this as a hard requirement. Could you please explicitly (not) require this?
• xnor: @RGS That's fair. I was misunderstanding the "lowest degree" rule to refer to how the polynomials are shown in the output (without extraneous higher-power zero terms) rather than as a mathematical property of the polynomials themselves, which is of course not affected by how they are shown.

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.

• Bubbler: -4 bytes by replacing m.reduce(...) with m.reduceRight((s,[v],i)=>v*g(0,m.map(([,...r])=>r).filter(_=>i--))-s,0).
• Arnauld: @Bubbler Very nice!
• Bubbler: Looks like reduce instead of reduceRight also works -- it would invert the sign of even-sized(?) matrices' determinants, but the same inversion happens on both numerator and denominator, cancelling each other. (It gives spurious -0s, but -0==0 anyway.)
• Arnauld: @Bubbler Yes, that's fine.

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.

• Giuseppe: 52 bytes
• Giuseppe: alternative 55 bytes
• Giuseppe: Note you need the a= in the seq statement for the first test case -2.
• Robin Ryder: @Giuseppe Thanks! I really need to remember about seq_along; this isn't the first time you've improved one of my answers with that!
• JDL: I have this niggling feeling that R has an inbuilt function to do this, but I cannot for the life of me think what it is (and it would probably be disqualified anyway)

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

• AviFS: Nice solution; wouldn't have thought of it! Pretty sure you can golf it down even further to 13 bytes with either of {⍵⌹m∘.*m←⍳⍴⍵} or, closer to your style, n⌹b∘.*b←⍳⍴n←⎕ :)
• Graham: @AviF.S. Thanks for your suggestion but the first option is not available in APL+WIN and the second the comma before quad is to force a vector if a scalar is entered and the 0⍕ is required to filter out the extremely small errors that can occur and display as nonexistent coefficients
• AviFS: You're totally right about the , being necessary to vectorize! I looked at the other test cases, but didn't notice the domain error with the scalar; you're totally right! About the 0⍕, I'm curious which test cases yielded errors. I haven't been able to reproduce any, but then again TIO is running Dyalog classic and not APL+WIN, could that be the difference?
• Graham: @AviF.S. Take [10, 20, 90, 280] -> [10, 0, 0, 10] for example without 0⍕ the result would be 10 ¯9.371726267E¯15 1.308691287E¯14 10 as opposed to 10 0 0 10
• Graham: @AviF.S. if you can golf my approach down using a modern rather than my ancient APL with your own entry by all means do so as I like to see APL being able to compete with golfing languages.

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.

• ZaMoC: 49 bytes
• Lukas Lang: According to the comments under the question you can also return the polynomial itself rather than its coefficients, which gets you to 38 bytes
• the default.: @LukasLang Nice. Reading the docs for Expand, I noticed Apart that achieves the same thing for one byte less in this case.

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.

• RGS: The output doesn't have to be reversed, you can give the coefficient list from highest degree to lowest, so 69 bytes.
• Command Master: Thanks, I edited the answer
• agtoever: I think you can skip the -1: this gives an extra leading 0, but this is allowed. Saves you 2 bytes.
• agtoever: If you use Python 3.8 with PEP572 (assignment expressions), you can save another 2 bytes with l:=len(x) for the first len and replace the second len(x) with just l.
• Command Master: I can't find numpy for pyhon 3.8

Answer 8

Charcoal, (削除) 68 (削除ここまで) 62 bytes

×ばつκιη»Iυ

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:

≔⟦1⟧η

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

FLθ«

Loop over the \$ m \$ terms.

⊞υ0

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

×ばつκXιλ∨ΠEι⊕κ1ζ

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 \$.

»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)

• Razetime: APL smashes the competition. Nice.

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)

• Grimmy: 24 without using matrices. Feels like just brute-forcing polynomials could be even shorter.
• Grimmy: 23
• Kevin Cruijssen: @Grimmy Ah, it's in reverse order of course. Still, the [4,5,62,733,4160,15869,47290,118997] test case lacks its 0 between the 1s. Same applies to some of the other test cases, like [-3,-1,5,15,29].

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 $$

• RGS: I don't understand how to use the "try it online" link you provide
• Adám: Shalom Uriel. Baruch mechayei hameitim!
• Uriel: @RGS wrong link, fixed :)
• Uriel: @Adám Thanks! Long time away
• Uriel: @RGS forget what I said - updated to use the Rationals ring that gives integer answers

Answer 14

Jelly, 14 bytes

×ばつæ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?

×ばつæ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)
 8 - chain's left argument, V
 ×ばつ - matrix-multiplication
 ær0 - round to zero decimal places (vectorises)