# Exponential transform of an integer sequence

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

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 , 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 =  * len(a)
b = 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]

• A worked example would be nice. Jul 24, 2022 at 8:17
• Is the solution allowed to fail due to floating-point errors? Jul 24, 2022 at 8:20
• @CommandMaster Yes, allowed. Jul 24, 2022 at 9:02
• @emanresuA Added a Python example. Jul 24, 2022 at 9:27
• Seems to be EXP on OEIS transforms (Maple) (search for "EXP:="), which can be accessed from this page. Not sure what the leading 1+ in the comment means though. There's also a wiki page for sequence transforms but many aren't documented yet. Jul 25, 2022 at 1:33

# 05AB1E, 212018 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.

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

• I think you can write the base case as _&l=l to save a byte
– xnor
Jul 25, 2022 at 10:20

# 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

• Shorter: BellY[#,k,Rest@#2]~Sum~{k,0,#}&. Or even shorter if you don't take the leading zero: BellY[#,k,#2]~Sum~{k,0,#}& Jul 25, 2022 at 1:39
• @alephalpha Thanks, you can save even two more bytes with BellY[#,,#2]~Sum~{,0,#}&. Jul 25, 2022 at 2:05

# VyxalR, 16 bytes

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


Try it Online!

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

# 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;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; \\ 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]));

• In PARI/GP, curly braces are only requires for multiline functions, and the last expression in automatically returned. So this can be f(a)=b=vector(#a,i,0);b=1;for(i=2,#a,b[i]=sum(j=1,i-1,binomial(i-2,j-1)*a[j+1]*b[i-j]));b. Apr 13 at 5:49
• @alephalpha thanks for your comment! :-) Apr 13 at 6:12

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

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