12
\$\begingroup\$

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 , 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]
\$\endgroup\$
8
  • 1
    \$\begingroup\$ A worked example would be nice. \$\endgroup\$
    – emanresu A
    Jul 24 at 8:17
  • 1
    \$\begingroup\$ Is the solution allowed to fail due to floating-point errors? \$\endgroup\$ Jul 24 at 8:20
  • \$\begingroup\$ @CommandMaster Yes, allowed. \$\endgroup\$
    – alephalpha
    Jul 24 at 9:02
  • \$\begingroup\$ @emanresuA Added a Python example. \$\endgroup\$
    – alephalpha
    Jul 24 at 9:27
  • 1
    \$\begingroup\$ 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. \$\endgroup\$
    – Bubbler
    Jul 25 at 1:33

6 Answers 6

12
+100
\$\begingroup\$

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.

\$\endgroup\$
3
\$\begingroup\$

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

\$\endgroup\$
1
  • \$\begingroup\$ I think you can write the base case as _&l=l to save a byte \$\endgroup\$
    – xnor
    Jul 25 at 10:20
2
\$\begingroup\$

Vyxal R, 16 bytes

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

Try it Online!

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

\$\endgroup\$
2
\$\begingroup\$

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

\$\endgroup\$
2
  • \$\begingroup\$ 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,#}& \$\endgroup\$
    – alephalpha
    Jul 25 at 1:39
  • 1
    \$\begingroup\$ @alephalpha Thanks, you can save even two more bytes with BellY[#,,#2]~Sum~{,0,#}&. \$\endgroup\$
    – Steffan
    Jul 25 at 2:05
1
\$\begingroup\$

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.

\$\endgroup\$
1
\$\begingroup\$

Charcoal, 36 bytes

FLθ⊞υ∨¬ιΣEυ××κ§θ⁻ιλ∨¬λ÷Π…⁻ιλιΠ…·¹λIυ

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.

\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.