Computing zillions of digits of the "derangement constant" This is a sort of inspired sequel to the following question:
Evaluating $\sum\limits_{x=2}^\infty \frac{1}{!x}$ in exact form.
where the question is the discussion of the "$e$-like constant"
$$\sum_{n=2}^{\infty} \frac{1}{!n}$$
which is at least visually similar to the series for $e$, i.e.
$$\sum_{n=0}^{\infty} \frac{1}{n!}$$
except there's also an indexing difference. I call the former the "derangement constant" ($ɘ$?), for lack of a better name. This number is about 1.638. I want more digits. Millions; billions maybe. Kind of like $\pi$, but $\pi$ is old stuff; need a new kid on the block. :)
Actually the reason is because that comments on that question were talking about how that the series for the derangement constant converges fast, superlinear actually in a digit-for-digit sense, which at first seems like we should be able to use it to compute millions of digits. But in reality, things aren't so simple - "fast convergence" alone by no means does a fast algorithm make for a computer.
On a computer, it's not just the amount of terms that matters (i.e. the convergence rate), but also how many digits you need to keep for each term, because the computer has to go through the digits and add and multiply them. Addition, of course, is cheap; multiplication is not - in this domain, it's typically done using tricked out and once written, often kept secret, code based on Fast Fourier Transform (FFT) methods, which for numbers of $n$ digits have complexities
$$O(n \log n\ a(n))$$
for some slow-grow function $a(n)$. The ideal FFT has $a(n) = 1$, but this is unachievable in practice because of precision limits of the computer processor. We can call the complexity above $M(n)$, or complexity of multiplication of $n$-digit numbers.
Fortunately though, for a naive implementation of $e$, we don't need this FFT multiplication. Naively, to sum a series of the type above, if we want millions of digits, we better evaluate each term to millions of digits precision. Given the factorial has a recurrence relation
$$(n+1)! = (n+1) \cdot n!$$
which gives
$$\frac{1}{n!} = \frac{1}{n} \frac{1}{(n-1)!}$$
we can get $n$ digits in time
$$O\left(n\ e^{W(\ln n)}\right)$$
which isn't too bad, but is still bad enough for large $n$ (millions or more) that we want something better. So we need to do something more clever. The usual way to overcome this is to use methods based on binary splitting - basically, consider the following. Note that we can write a factorial as
$$n! = \left(1 \cdot 2 \cdot 3 \cdot \cdots \cdot m\right) \cdot \left((m+1) \cdot (m+2) \cdot \cdots \cdot (n-1) \cdot n\right)$$
for some $m$ with $1 < m < n$. Defining
$$P(a, b) := (a+1)(a+2)\cdots(b-1)b$$
we have the identity
$$P(a, b) = P(a, m) P(m, b)$$
which allows us to recursively compute
$$P(0, n) = n!$$
using much smaller multiplications. Something similar can be done by splitting the sum for $e$. Take
$$e = \sum_{n=0}^{\infty} \frac{1}{n!}$$
and now write the equation
$$\frac{P(a, b)}{Q(a, b)} = \sum_{n=a+1}^{b} \frac{1}{(a+1)(a+2)\cdots (n-1)n}$$
Then note that
$$\frac{P(a, b)}{Q(a, b)} = \frac{P(a, m)}{Q(a, m)} + \frac{1}{Q(a, m)} \frac{P(m, b)}{Q(m, b)}$$
once we take
$$Q(a, b) := (a+1)(a+2) \cdots (b-1)b$$
like before. In particular, we have the recurrence equations
$$Q(a, b) = Q(a, m) Q(m, b)$$
and
$$P(a, b) = P(a, m) + Q(m, b) P(m, b)$$
of which the second can be found by multiplying out. Then,
$$e \approx 1 + \frac{P(0, N)}{Q(0, N)}$$
gives $N+1$ terms of the $e$-series, with, again, smaller multiplies. This saves a lot; and permits efficient computation of $e$ by leveraging fast multiplication - in particular, the complexity to generate $n$ digits is now on the order instead of $O(\log(n)\ M(n \log n))$; obviously much better at least if $M(n)$ is suitably good, which it is for FFT-based multiplication methods (giving $O(n \log^2(n)\ a(n))$).
However, this trick does not appear to generalize so easily to the derangements. In particular, while factorials expand nicely as a product, derangements $!n$ have a "branching" recurrence identity
$$!(n+1) = n(!n\ +\ !(n-1))$$
which means that, in particular, we cannot simply parcel out the computation of $!n$ into a "1 to $m$" and "$m+1$ to $n$" part as we did above. Worse, because the derangement is in the denominator, we cannot do something analogous to
$$\frac{1}{n!} = \frac{1}{n} \frac{1}{(n-1)!}$$
which was key to the naive method for $e$! Instead, we must do the full reciprocal, and now have horrible costs on the order of $O(e^{W(\ln n)} \lg(n)\ M(n))$!
So is there some trick that we can use to compute $ɘ$ efficiently like how we can use binary splitting + FFT multiplication to compute $e$ efficiently?

ADD: Not sure if this is useful. But I found on Wikipedia that
$$!n = n! \sum_{i=0}^{n} \frac{(-1)^i}{i!}$$
Thus
$$\frac{1}{!n} = \frac{1}{n!} \frac{1}{\sum_{i=0}^{n} \frac{(-1)^i}{i!}}$$
and then we can use the series reciprocal formula to get
$$\frac{1}{!n} = \frac{1}{n!} \sum_{k=0}^{\infty} (-1)^k d_{k,n}$$
where $d_{0,n} = 1$ and $d_{k,n} = -\sum_{l=\max(0, k-n)}^{k-1} \frac{d_{l,n}}{(k-l)!}$. See:
https://mathoverflow.net/questions/53384/power-series-of-the-reciprocal-does-a-recursive-formula-exist-for-the-coeffic
Thus
$$ɘ = \sum_{n=2}^{\infty} \sum_{k=0}^{\infty} \frac{(-1)^k d_{k,n}}{n!}$$
however; I'm not sure if this recursion for this $d_{k,n}$ is any better. At least though we're no longer recursing in a denominator!
ADD 2: The next thing I wonder about is the possibility that, because this is a nested sum, we should not consider a simple binary splitting, but a splitting into four "squares" or "rectangles", i.e. a quad splitting. In particular, if we want to sum the double sum
$$\sum_{n=0}^{\infty} \sum_{m=0}^{\infty} a_{n,m}$$
as an approximant
$$\sum_{n=0}^{N} \sum_{m=0}^{M} a_{n,m}$$
we should break into fours
$$a_{0,0} + a_{0,1} + a_{1,0} + \sum_{n=1}^{N} \sum_{m=1}^{M} a_{n,m} = \left(\sum_{n=1}^{n_m} \sum_{m=1}^{m_m} a_{n, m}\right) + \left(\sum_{n=n_m+1}^{N} \sum_{m=1}^{m_m} a_{n, m}\right) + \left(\sum_{n=1}^{n_m} \sum_{m=m_m+1}^{M} a_{n, m}\right) + \left(\sum_{n=n_m+1}^{N} \sum_{m=m_m+1}^{M} a_{n, m}\right)$$
viz.
$$a_{0,0} + a_{0,1} + a_{1,0} + \frac{P(0, 0, N, M)}{Q(0, 0, N, M)} = \frac{P(0, 0, n_m, m_m)}{Q(0, 0, n_m, m_m)} + \frac{P(n_m, 0, N, m_m)}{Q(n_m, 0, N, m_m)} + \frac{P(0, m_m, n_m, M)}{Q(0, m_m, n_m, M)} + \frac{P(n_m, m_m, N, M)}{Q(n_m, m_m, N, M)}$$
 A: Here I present forms of the constant. The first will use algebraic properties. Let the constant be denoted c for “constant”. Note I will use “!x” for the derangement as seen in the question. Also note the Pi-Product notation:
$$c\mathop=^\text{def}\sum_2^\infty \frac 1{!x}\implies e^c=e^{\sum_2^\infty \frac 1{!x}}=e^{\frac1{!0}+ \frac1{!1} +\frac1{!2} +…}=e^{\frac1{!0} }e^{\frac1{!1}}e^ \frac1{!2}…= \prod_2^\infty e^\frac{1}{!x}$$
Now you see

various possibilities for a series expansion for $e^y$:

$$e^c= \prod_2^\infty e^\frac{1}{!x}= \prod_{x=2}^\infty \sum_{n=0}^\infty\frac{1}{(!x)^n n!} $$
There is also an option with multinomial theorem as it will converge because a finite sum would be taken with the theorem.
Another nice, but less interesting expansion would be using the Regularized Gamma functions and Poisson random variables section in the linked article:
$$c=\sum_2^\infty \frac1{!x}=e\sum_3^\infty\frac{1}{Γ(x,-1)}= e\sum_3^\infty\frac{1}{Q(x,-1)Γ(x)} $$
Let there be the usage of one of the reciprocal gamma function expansions and Euler-Mascheroni Constant:
$$c= e\sum_3^\infty\frac{1}{Q(x,-1)Γ(x)}= e\sum_3^\infty\frac{x \prod\limits_{n=1}^\infty\frac{\frac xn+1 }{\left(\frac1n+1\right)^n}}{Q(x,-1)}= \sum_3^\infty\frac{xe^{γx+1}\prod_\limits{n=1}^\infty\frac{\frac xn+1}{e^\frac xn}}{1-P(x,-1)} $$
Using the links, Q and P are both the Upper and Lower Regularized Incomplete Gamma functions and:
$Q(x,-1)$, because $x\in\Bbb Z$, is also the Cumulative Distribution Function
for Poissan random variables for Probability$(x>X\in\text{Poissan(-1) random variable})$.
Similarly, $P(x,-1)$ is the Cumulative Distribution Function for Gamma random variables with shape parameter $x$ and scale parameter of $1$ evaluated at $-1$.
I have little experience with either of these statistical interpretations, but they help put context to alternative forms of the constant. There are many more possible manipulations we can do of the constant, but not all of them are practical. Please correct me and give me feedback!
It can be shown that the following is true using @Jack Barber’s method in
$$\sum_\Bbb N \text{erfc}(x)$$
Here is an integral representation using the linearity of the Floor function and the Meijer G function:
$$\sum_2^\infty \frac{1}{!x}=-\int_2^\infty \lfloor x-1\rfloor \frac{d}{dx} \frac{1}{!x}dx=\int_2^\infty \frac{d}{dx} \frac{1}{!x} dx-\int_2^\infty\lfloor x\rfloor  \frac{d}{dx} \frac{1}{!x} dx=\frac1{!\infty}-\frac1{!2}-\int_2^\infty \lfloor x \rfloor \frac{d}{dx} \frac{1}{!x} =-1-\int_2^\infty \lfloor x \rfloor \frac{d}{dx} \frac{1}{!x}dx= -1-\int_2^\infty \lfloor x \rfloor\left(-\frac{\text G_{2,3}^{3,0}\big(-1\left|_{0,0,x+1}^{\ \ \ \ 1,1}\right)}{e(!x)^2}-\frac{i\pi}{!x}\right)dx =\frac1e\int_2^\infty \frac{\lfloor x\rfloor\text G_{2,3}^{3,0}\left(-1\big|_{0,0,x+1}^{\ \ \ \ 1,1}\right)}{(!x)^2} dx+i\pi\int_2^\infty\frac{\lfloor x\rfloor}{!x}dx-1$$
Please correct me and give me feedback!
There is also the classic Riemann Sum which you will use for approximation.
