# Engel expansion

How do I do an Engel expansion in Mathematica?

And find the unique non-decreasing sequence of positive integers $\{ a_1 ,a_2,a_3,\dots \}$ , such $$\frac{e}{\pi}=\frac{1}{a_1}+\frac{1}{a_1 a_2}+\frac{1}{a_1 a_2 a_3}+\cdots$$

The MathWorld article you linked to has a link to a notebook containing Weisstein's code for generating Engel expansions.

In any case, here is a slight simplification of Weisstein's code:

engel[x_, n_Integer?Positive] :=
Nest[Append[#, Ceiling[1/(x - Total[1/FoldList[Times, #]])/Times @@ #]] &,
{}, n]


Test:

engel[E/π, 27]
{2, 2, 3, 3, 7, 23, 43, 58, 30503, 32703, 44621, 51291, 479922, 781181, 2136095,
2435396, 2600287, 2860451, 12557067, 35938396, 37552004, 75289454, 109334609,
321224695, 336383475, 451039452, 586390841}


Check:

Total[1/FoldList[Times, %]] - E/π // N
0.


The code at OEIS works for me:

 EngelExp[A_, n_] :=
Join[Array[1 &, Floor[A]],
First@Transpose@
NestList[{Ceiling[1/Expand[#[[1]] #[[2]] - 1]],
Expand[#[[1]] #[[2]] - 1]} &, {Ceiling[1/(A - Floor[A])],
A - Floor[A]}, n - 1]]

res = EngelExp[N[E/Pi, 500000], 27]


returns:

{2, 2, 3, 3, 7, 23, 43, 58, 30503, 32703, 44621, 51291, 479922, 781181, 2136095, 2435396, 2600287, 2860451, 12557067, 35938396, 37552004, 75289454, 109334609, 321224695, 336383475, 451039452, 586390841}

Test:

 Sum[1/Apply[Times, res[[1 ;; j]]], {j, 1, Length[res]}] - E/Pi // N


returns 0. (meaning it's numerically close to zero; for higher precision use N[(expression),(n_digit_precision)])