Exponential transform of an integer sequence

Question

The exponential generating function (e.g.f.) of a sequence \$a_n\$ is defined as the formal power series \$f(x) = \sum_{n=0}^{\infty} \frac{a_n}{n!} x^n\$.

When \$a_0 = 0\$, we can apply the exponential function \$\exp\$ on this formal power series:

\$\begin{align} \exp(f(x)) &= \sum_{n=0}^{\infty} \frac{1}{n!} f(x)^n \\ &= \sum_{n=0}^{\infty} \frac{1}{n!} \left(\sum_{m=1}^{\infty} \frac{a_m}{m!} x^m\right)^n \\ &= \sum_{n=0}^{\infty} \frac{b_n}{n!} x^n \\ \text{where} \\ b_0 &= 1 \\ b_n &= \sum_{k=1}^n \binom{n-1}{k-1} a_k b_{n-k} \text{ when } n>0 \end{align}\$

This is the exponential generating function of the sequence \$b_n\$. If all \$a_n\$ are integers, then all \$b_n\$ are also integers. So this is a tranformation of integer sequences.

It seems that OEIS call this the Exponential transform. But I can't find a reference to its definition.

Here are some examples on OEIS:

• A001477 (the nonnegative integers, \$0,1,2,3,4,\dots\$) -> A000248 (\$1,1,3,10,41,\dots\$)
• A001489 (the nonpositive integers, \$0,-1,-2,-3,-4,\dots\$) -> A292952 (\$1, -1,-1,2,9,\dots\$)
• A000045 (Fibonacci numbers, \$0,1,1,2,3,\dots\$) -> A256180 (\$1,1,2,6,21,\dots\$)
• A160656 (\$0\$ and the odd primes, \$0,3,5,7,11,\dots\$) -> A353079 (\$1,3,14,79,521,\dots\$)
• A057427 (\$0,1,1,1,1,\dots\$) -> A000110 (Bell numbers, \$1,1,2,5,15,\dots\$)

Task

Given a finite integer sequence, compute its Exponential transform.

The length of the input sequence is always greater than \$0\$. Its \$0\$th term is always \$0\$. You can omit the leading \$0\$ in the input.

If the input sequence has length \$n\$, you only need to output the first \$n\$ terms of the output sequence.

Input and output can be in any reasonable format, e.g., a list, an array, a polynomial, a function that takes \$i\$ and returns the \$i\$th term (0-indexed or 1-indexed), etc.

You may also take the input sequence and an integer \$i\$, and output the \$i\$th term (0-indexed or 1-indexed) of the output sequence.

This is code-golf, so the shortest code in bytes wins.

Example Python code

This example code uses the above recurrence formula. There are other formulas that might give shorter answers.

import math

def exponential_transform(a):
    b = [0] * len(a)
    b[0] = 1
    for i in range(1, len(a)):
        b[i] = sum(math.comb(i-1, j-1) * a[j] * b[i - j] for j in range(1, i + 1))
    return b

Testcases

[0, 0, 0, 0, 0] -> [1, 0, 0, 0, 0]
[0, 1, 0, -1, 0, 1, 0, -1] -> [1, 1, 1, 0, -3, -8, -3, 56]
[0, 1, 2, 3, 4, 5, 6, 7] -> [1, 1, 3, 10, 41, 196, 1057, 6322]
[0, -1, -2, -3, -4, -5, -6, -7] -> [1, -1, -1, 2, 9, 4, -95, -414]
[0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55] -> [1, 1, 2, 6, 21, 86, 404, 2121, 12264, 77272, 525941]
[0, 3, 5, 7, 11, 13, 17, 19] -> [1, 3, 14, 79, 521, 3876, 31935, 287225]
[0, 1, 1, 1, 1, 1] -> [1, 1, 2, 5, 15, 52]
[0, 1, 2, 5, 15, 52] -> [1, 1, 3, 12, 60, 358]
[0, 1, 3, 12, 60, 358] -> [1, 1, 4, 22, 154, 1304]
[0, 1, 4, 22, 154, 1304] -> [1, 1, 5, 35, 315, 3455]
[0, 1, 5, 35, 315, 3455] -> [1, 1, 6, 51, 561, 7556]

Answer 1

05AB1E, 21 20 18 16 bytes

L©ææʒ˜{®Q}€€gèPO

Try it online!

Returns the \$ i \$-th term, zero indexed. It runs in complexity \$ \Omega(2^{2^i}) \$, so it can only work in a reasonable time for \$ i \leq 4 \$.

Explanation

We will define \$ c_n = \frac{a_n x^n}{n!} \$.

In the expression $$ \left(\sum_{m=1}^{\infty} \frac{a_m}{m!} x^m\right)^n = \left(\sum_{m=1}^{\infty} c_m\right)^n $$

we know from the multinomial theorem that the coefficient of \$ \prod {c_m ^ {k_m}} \$ is $$ \frac{n!}{\prod{k_m !}} $$ The \$ n! \$ cancels with the \$ \frac{1}{n!} \$ in \$ \sum_{n=0}^{\infty} \frac{1}{n!} \left(\sum_{m=1}^{\infty} \frac{a_m}{m!} x^m\right)^n \$, so we get $$ \exp(f(x)) = \sum \frac{\prod{c_m ^ {k_m}}}{\prod{k_m!}} $$ If the exponent of \$ x \$ in the equation is \$ j \$, we get that \$ \sum k_m m = j \$. If we create a sequence \$ d_m \$ by repeating \$ m \$ \$ k_m \$ times, we get that \$ d_m \$ is a partition of \$ j \$.

Now that we have a mathematical expression for how much time each partition of \$ i \$ occurs, we can think of a combinatorical one.

By playing with equations, we can find out that the coefficient of \$ \prod {a_{d_i}} \$ is the number of unique ways to divide a set with \$ j \$ elements to sets of sizes \$ d_i \$. Therefore, we can generate a set of size \$ j \$, list all of its separations to other sets, take the size of each set in those separations, index that size to \$ a \$, take the product of each, and take the total sum.

L       push the list 1...input index
©       save it in register ® (without popping)
æ       take its powerset
æ       take the powerset's powerset
ʒ       and only keep sets such that
 ˜       if you flatten them
 {       and sort the result
 ®       and then push the original list
 Q       they compare equal
}
€        for each of those sets
 €       map each of its sets
  g      to its size
è        index that size into the input list
P        take the product of each list
O        and take the sum

If the index is 0, L returns [1, 0], and no array sorted can be equal to it, so the filter returns an empty list. The map and index then do nothing, then P multiplies it to get 1, which O preserves.

Answer 2

Haskell, 59 bytes

f@(a:b)&g@(c:d)=a*c:zipWith(+)(b&g)(f&d)
_&_=[]
e b=1:e b&b

Try it online!

e accepts the sequence as a list of integers without the leading 0.

If \$f(x)\$ is the e.g.f. of \$a_n\$, then \$f(0) = a_0\$ and \$f'(x)\$ is the e.g.f. of \$a_{n + 1}\$. So this program is a direct encoding of the product rule \$[f(x)g(x)]' = f'(x)g(x) + f(x)g'(x)\$ and the chain rule \$[e^{f(x)}]' = e^{f(x)}f'(x)\$.

Answer 3

Wolfram Language (Mathematica), 37 24 bytes

BellY[#,,#2]~Sum~{,0,#}&

Try it online!

Based on the Mathematica code given by Vladimir Reshetnikov on OEIS A256180. Takes the 0-indexed item to output and the list.

-13 bytes thanks to alephalpha

Answer 4

Vyxal R, 16 bytes

ṗṗ'fs¹ɾ⁼;vvLİvΠ∑

Try it Online!

Port of Command Master's 05AB1E answer, upvote that!

Answer 5

PARI/GP, 103 92 bytes

saved 11 bytes thanks to the comment of @alephalpha

Golfed vesrion, try it online!

f(a)=b=vector(#a,i,0);b[1]=1;for(i=2,#a,b[i]=sum(j=1,i-1,binomial(i-2,j-1)*a[j+1]*b[i-j]));b

Ungolfde version, modified from the provided Python example code

exponential_transform(a) = {
    b = vector(#a, i, 0); \\ Initialize b as a zero vector of the same length as a
    b[1] = 1; \\ Set the first element of b to 1

    for(i=2, #a,
        b[i] = sum(j=1, i-1, binomial(i-2, j-1) * a[j+1] * b[i-j])
    );
    return(b);
}
print("Test case 1: ", exponential_transform([0, 0, 0, 0, 0]));
print("Test case 2: ", exponential_transform([0, 1, 0, -1, 0, 1, 0, -1]));
print("Test case 3: ", exponential_transform([0, 1, 2, 3, 4, 5, 6, 7]));
print("Test case 4: ", exponential_transform([0, -1, -2, -3, -4, -5, -6, -7]));
print("Test case 5: ", exponential_transform([0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55]));
print("Test case 6: ", exponential_transform([0, 3, 5, 7, 11, 13, 17, 19]));
print("Test case 7: ", exponential_transform([0, 1, 1, 1, 1, 1]));
print("Test case 8: ", exponential_transform([0, 1, 2, 5, 15, 52]));
print("Test case 9: ", exponential_transform([0, 1, 3, 12, 60, 358]));
print("Test case 10: ", exponential_transform([0, 1, 4, 22, 154, 1304]));
print("Test case 11: ", exponential_transform([0, 1, 5, 35, 315, 3455]));

Answer 6

Python, 85 bytes

import math
f=lambda n,l:sum(math.comb(n-1,i)*l[i]*f(n-i-1,l)for i in range(n))+(n<1)

Attempt This Online!

Doesn't take the leading zero in the input. If that's not allowed, I will change. Based on the Maple code given by Alois P. Heinz in A007446 and A353079.

Answer 7

Charcoal, 36 bytes

ＦＬθ⊞υ∨¬ιΣＥυ××κ§θ⁻ιλ∨¬λ÷Π…⁻ιλιΠ…·¹λＩυ

Try it online! Link is to verbose version of code. Explanation: A translation of the recurrence formula given in the question, except the subscript of b is used as the loop index, and the combination has to be written out as a ratio of products of ranges.

33 bytes using the newer version of Charcoal on ATO:

ＦＬθ⊞υ∨¬ιΣＥυ××κ§θ⁻ιλ÷Π…⁻ιλιΠ…·¹λＩυ

Attempt This Online! Link is to verbose version of code.