# Generating generating expressions for sequences

(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.

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

• related
– RGS
Mar 29, 2020 at 10:12
• May we take m as part of the input? Mar 29, 2020 at 14:42
• @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)
– RGS
Mar 29, 2020 at 15:25
• 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? Mar 29, 2020 at 21:12
• @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.
– xnor
Mar 29, 2020 at 23:07

# 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 :

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

• -4 bytes by replacing m.reduce(...) with m.reduceRight((s,[v],i)=>v*g(0,m.map(([,...r])=>r).filter(_=>i--))-s,0). Mar 29, 2020 at 23:55
• @Bubbler Very nice! Mar 30, 2020 at 1:48
• 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.) Mar 30, 2020 at 2:06
• @Bubbler Yes, that's fine. Mar 30, 2020 at 2:11

# 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.

• 52 bytes Mar 30, 2020 at 16:45
• alternative 55 bytes Mar 30, 2020 at 16:53
• Note you need the a= in the seq statement for the first test case -2. Mar 30, 2020 at 18:18
• @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! Mar 30, 2020 at 19:43
• 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)
– JDL
Mar 31, 2020 at 15:48

# 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

• 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←⎕ :) Mar 29, 2020 at 15:37
• @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 Mar 29, 2020 at 15:53
• 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? Mar 29, 2020 at 16:49
• @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 Mar 29, 2020 at 17:02
• @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. Mar 29, 2020 at 17:06

# Wolfram Language (Mathematica), 5049 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.

• 49 bytes Mar 29, 2020 at 13:06
• According to the comments under the question you can also return the polynomial itself rather than its coefficients, which gets you to 38 bytes Mar 31, 2020 at 13:26
• @LukasLang Nice. Reading the docs for Expand, I noticed Apart that achieves the same thing for one byte less in this case. Mar 31, 2020 at 14:40

# 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


# Pari/GP, 38 bytes

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


Try it online!

# Python 3 + Numpy, 69 bytes

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


Try it online!

• The output doesn't have to be reversed, you can give the coefficient list from highest degree to lowest, so 69 bytes.
– RGS
Mar 29, 2020 at 16:38
• Thanks, I edited the answer Mar 29, 2020 at 16:41
• I think you can skip the -1: this gives an extra leading 0, but this is allowed. Saves you 2 bytes. Mar 29, 2020 at 21:13
• 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. Mar 29, 2020 at 21:19
• I can't find numpy for pyhon 3.8 Mar 30, 2020 at 5:20

# 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.

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! # APL (Dyalog Unicode), 10 bytesSBCS ⊢⌹∘.*⍨∘⍳∘≢  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)  • APL smashes the competition. Nice. Oct 11, 2020 at 10:42 # 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)  • 24 without using matrices. Feels like just brute-forcing polynomials could be even shorter. Apr 7, 2020 at 9:52 • 23 Apr 7, 2020 at 10:02 • @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]. Apr 8, 2020 at 7:19 # MATL, 12 bytes n:qGyz3$ZQYo


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

### 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]


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

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

# 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)
`