Following are three possible ideas, the first two are not that satisfactory.

The third one is a modification of the second ideas which might work.

I hope they can inspire others to create something that is useful.

**As a Fourier series**

First, we can rewrite $\lfloor x \rfloor$ as a Fourier series.

$$\lfloor x \rfloor = x - \{ x \} = x - \frac12 + \sum_{m=1}^\infty \frac{\sin(2\pi m x)}{\pi m}\tag{*1}$$

Since the discontinuities of $\lfloor x \rfloor$ is contained inside $\mathbb{Z} \subset \mathbb{Q}$. the floor function is continuous at irrational $x$. As a result, RHS of $(*1)$ converges pointwisely to LHS for irrational $x$.

Substitute $x$ by $ek$ for $k = 1, \ldots, n$ and sum over $k$, we obtain.

$$\sum_{k=1}^n \lfloor ek \rfloor
= \frac{e}{2}n(n+1) - \frac{n}{2}
+ \underbrace{\frac{1}{2\pi}\sum_{m=1}^\infty \frac{\cos(\pi m e) - \cos(\pi m e(2n+1))}{m\sin(\pi m e)}}_{I}
$$
In principle, if we can approximate the series $I$ on RHS accurate enough,
we can round the RHS to nearest integer and it will give us the value of LHS.
The problem is when we approximate $I$ by its partial sums,
the $\sin(m \pi e)$ factor in denominator make it very hard to figure out the correct number of terms to keep!

**Recursive evaluation**

If we didn't insist for a closed formula, it is possible to evaluate the sum in a recursive manner.

For $\alpha \in (1,\infty)\setminus \mathbb{Q}$ and $n \in \mathbb{Z}$, define
$\displaystyle\;S(\alpha,n) \stackrel{def}{=} \sum_{k=1}^n\lfloor \alpha k \rfloor$. The sum we want is $S(e,n)$.

There are two branches in the recursion:

**Case I** - $\alpha > 2$.

Rewrite $\alpha$ as $\beta + m$ where $\beta \in (1,2)$ and $m = \lfloor \alpha - 1\rfloor$, we have
$$S(\alpha,n) = \sum_{k=1}^n \left( mk + \lfloor \beta k\rfloor\right)
= \frac{m}{2}n(n+1) + S(\beta,n)$$

**Case II** - $\alpha < 2$.

Let $\beta = \frac{\alpha}{\alpha-1} \in (2,\infty) \setminus \mathbb{Q}$,
we have
$$S(\alpha,n) = \sum_{k=1}^n \lfloor\alpha k\rfloor
= \sum_{0 < \alpha k \le \alpha n} \lfloor\alpha k\rfloor
= \sum_{0 < \alpha k < \lceil\alpha n\rceil} \lfloor\alpha k\rfloor\tag{*2}
$$
For any $r \in (0,\infty) \setminus \mathbb{Q}$, sequences of the form $\left( \lfloor r k \rfloor \right)_{k\in\mathbb{Z}_{+}}$ are known as
Beatty sequence.

Since $\frac{1}{\alpha} + \frac{1}{\beta} = 1$, the two Beatty sequences $\left( \lfloor \alpha k\rfloor \right)_k$ and $\left( \lfloor \beta k\rfloor \right)_k$ are complementary. Every positive integers belongs to exactly one of these two sequences. As a corollary, for any $N \in \mathbb{Z}_{+}$, we have
$$\sum_{0 < \alpha k < N} \lfloor \alpha k\rfloor + \sum_{0 < \beta k < N}\lfloor \beta k \rfloor = \frac12 N(N-1)$$

Apply this to $(*2)$, we obtain
$$S(\alpha,n) = \frac12\lfloor \alpha n\rfloor\lceil \alpha n\rceil -
S\left( \beta, \left\lfloor \frac{1}{\beta}\lceil \alpha n\rceil\right\rfloor\right)$$

Combine these two branches, we can evaluate $S(\alpha,n)$ recursively.

It turns out a similar question
about $\sum_{k=1}^n \lfloor \sqrt{2} k \rfloor$ has been asked before.
In an answer by @merico,
there is another derivation of the recurrence formula in a slightly different form. Comparing our answers, I notice the term
$\left\lfloor \frac{1}{\beta}\lceil \alpha n\rceil\right\rfloor$ here can be simplified to $\lfloor (\alpha-1)n\rfloor$.

Since the recursion is a tail recursion, we can speed up implementation of $S(\alpha,n)$ by unwinding the recursion. Following is my implementation of $S(\alpha,n)$ in the CAS maxima.

```
S(a0,n0) := block(
[sum:0,sign:1,a:a0,n:n0,m],
while (n > 0) do
if( a > 2 ) then
(
m : floor(a-1),
sum : sum + sign*m*n*(n+1)/2,
a : a - m
) else
(
m : floor(a*n),
sum : sum + sign*m*(m+1)/2,
sign : -sign,
a : a/(a-1),
n : m-n
),
sum
);
```

Using the command `S(bfloat(%e),10^9)`

in maxima with $100$ digits accuracy,
above code evaluates the sum $S(e,10^9)$ in $44$ steps and returns
$$S(e,10^9) = 1359140915088663532$$

As a double check, we can compare this value with the approximation
$$S_{appx}(\alpha,n) = \frac{\alpha}{2}n(n+1) - \frac{n}{2}$$
Since $S_{appx}(e,10^9) \approx 1359140915088663531.9\ldots$, above value of $S(e,10^9)$ should be correct.

The basic problem of this approach is when $n$ is large, we need a very accurate value of $e$ as a seed. We also need to keep the precision all the way during computation. For example, if we compute the number using default precision
in maxima, the command `S(%e,10^9),numer`

returns a wrong number $1359140915088663452$. If we only use the `S(bfloat(%e),10^9)`

without bumping up the precision, we get another wrong number $1359140915088663538$.

**Something that should work?**

Inspired by Jack D'Aurizio's answer to yet another variant of this question, I investigated whether one can replace $e$ by one of its convergent as input to $S(\alpha,n)$. It does seem to work.

The basic ideas go like this.

For any $\alpha \in (1,\infty)\setminus\mathbb{Q}$, consider its representation as a CF (continued fraction):

$$\alpha = [a_0; a_1, a_2, \ldots ]$$

Let $\displaystyle\;\frac{p_k}{q_k} = [a_0;a_1,\ldots, a_k]\;$ be the $k^{th}$ convergent of $\alpha$. One property of the convergent is
$$\left| \alpha - \frac{p_k}{q_k} \right| < \frac{1}{q_k^2}$$
Using this, one can show that $\displaystyle\;\left\lfloor \frac{p_k}{q_k} n \right\rfloor = \lfloor \alpha n\rfloor$ for all $n < q_k$.

When we feed $\displaystyle\;\frac{p_k}{q_k} = [ a_0, a_1, a_2, \ldots, a_k ]\;$ as input to above implementation of $S(\alpha,n)$, the variables will be updated in following manner.

$$\overbrace{\begin{align}
\alpha &\leftarrow [1; a_1, a_2, \ldots, a_k ]\\
n &\leftarrow n
\end{align}}^{\alpha > 2}
\quad\text{ and }\quad
\overbrace{\begin{align}
\alpha &\leftarrow [ 1 + a_1; a_2, \ldots, a_k ]\\
n &\leftarrow \left\lfloor\frac{n}{ [ a_0 - 1; a_2, \ldots, a_k]} \right\rfloor
\end{align}}^{\alpha < 2}
$$

If one follow the steps in the while loop, the variables will be transformed in
essentially the same pattern.

All the finite CF appearing during this process are convergents of corresponding CF associated with $\alpha$. If the denominator of all these finite CF is larger than the $n$ they see in a step, they will produce the same result as if $\alpha$ is the input.

In short, if one feed a high enough order convergent of $\alpha$ to above implementation of $S(\alpha,n)$, one obtain the same result. The advantage
of this approach is we will be using exact rational number arithmetic
and we no longer need to worry about numerical error.

For the problem at hand, if one want to compute $S(e,n)$ for a large $n$,
we can estimate the required order of convergent of $e$ as follows.

Find the first $\ell$ such that $2^\ell \ell! > n$ and then set $k = 3\ell$.
For $n \approx 10^{4000}$, $k \approx 4011$ should be enough.

On my PC, I can compute $S(e,10^{4000})$ using maxima in less than a minute.
However, I've to admit I have no way to verify I got the right answer.