Yes, there is a polynomial time solution. It is not very pretty, so simplifications or alternative approaches are welcome.
TL;DR: the exact weights are rarely important, so we can round them up to a multiple of $K := \lceil \frac{3M}{n^3} \rceil$ or something similar (there are only $\approx n^3$ rounded item weights). Then, we can do dynamic programming that keeps track of rounded weights for two bins and exact weight for the third bin. If we somehow got a false negative result, then in all solutions two bins are almost full at the same time. Such cases turn out to be highly structured (all items have weights that are either close to $0$ or close to $3M$, the total weight of the former is small and there are exactly $3\ell + 2$ items of the latter type for some $\ell$) and can be solved by an algorithm that greedily moves items between bins.
Only the case when $M > n^4$ is interesting, otherwise simple dynamic programming will do the job. Similarly, assume that $n > 10^9$ to make things like $Kn^2$ negligible when compared to $M$. Do not worry, actual precise bounds are much better, I just chose these to not think about the details too much.
As mentioned above, the first step is to try finding a solution with "discretised" weights. Pick $K := \lceil \frac{3M}{n^3} \rceil$. For each item, its discrete weight is its weight rounded up to a multiple of $K$. For each bin, its discrete size is its size, rounded down to a nearest multiple of $K$. Items discretely fit into a bin, if the sum of their discrete weights does not exceed the discrete size of the bin. Clearly, if items discretely fit into a bin, then they do fit in the bin for real. The opposite implication is not true, but only because of rounding issues; if several items fit into a bin, but not discretely, then their total weight is at least $(S + M) - (n+1)K$ (at most $n$ items, each causes a rounding error of at most $K$; rounding of the bin size also causes a rounding error of at most $K$). Let us say that any bin with sum of weights at least $(S + M) - (n+1)K$ is almost full. Unless the bin is almost full, discrete fitting and real fitting is the same thing.
Consider any solution (by solution, I mean a correct distribution of items to bins). There can be at most two almost full bins; otherwise the total weight of items is at least $3(S + M) - 3(n+1)K = 3S + (3M - 3(n+1)K) > 3S$, as $K \approx \frac{3M}{n^3}$).
We can find all solutions with at most one full bin by simple dynamic programming, which answers the following question: suppose that we distributed some prefix of items to bins and we know the sums of discrete item weights in the second and the third bins (there are only $n \cdot \lceil \frac{3M}{K} \rceil \approx n^4$ possible discrete weights for each bin), what is the minimal possible total real weight of items in the first bin in this case? We will not miss out the solutions where there is at most one almost full bin (because we can track the weight in one bin precisely).
The only remaining problem is the case when all optimal solutions have exactly two almost full bins. Let us study the structure of the solutions in this case (this is not an algorithmic part, but rather only a mathematical setup; so, we are not interested in algorithmic aspects for the following few paragraphs).
Suppose that the first bin is filled up to $S + M - \varepsilon_1$ and the second is filled up to $S + M - \varepsilon_2$, where $\max(\varepsilon_1, \varepsilon_2) \leqslant (n+1)K$. Then, the third bin is filled up to $3S - (S + M - \varepsilon_1) + (S + M - \varepsilon_2) = S - 2M + \varepsilon_1 + \varepsilon_2$.
Let us try to move some items from the first bin to the third (in order to make them both not almost full). Why cannot we do that? Because all items in the first bin are either small (have weight at most $(n+1)K$) and the first bin will not stop being almost full, or because all items in the first bin are big (have weight at least $3M - 3K(n+1)$), so the third bin will become almost full (or, maybe, there is not even enough space for the item). About the latter part: if the weight of an item is less than $3M - 3K(n+1)$, then the third bin will be filled up to less than
\begin{equation*}
S - 2M + \varepsilon_1 + \varepsilon_2 + (3M - 3K(n+1)) = S + M - (3K(n+1) - \varepsilon_1 - \varepsilon_2) \leqslant S + M - K(n+1),
\end{equation*}
because $\max(\varepsilon_1, \varepsilon_2) \leqslant K(n+1)$.
Now, let us move all small items from the first bin to the third bin one-by-one. Then, either the first one stops being almost full at some point (while the third one is still not almost full, it still had almost $3M$ space left before all operations), or we move them all. Similarly, move all small items from the second bin to the third bin.
Now, all items in the first and the second bins are big, but they are still almost full. As before, denote the space left in the first and the second bins by $\varepsilon_1$ and $\varepsilon_2$. I claim that no items in the third bin have weight in the range $[10K(n+1), 3M - 10K(n+1)]$. Why? Because we can take any such item and swap its places with a big item from the first bin. Then, the total weight of items in the first bin will decrease by at least $(3M - 3K(n+1)) - (3M - 10K(n+1)) = 7K(n+1)$, hence it will not be almost full anymore. On the other hand, the total weight of items in the third bin was $S - 2M + \varepsilon_1 + \varepsilon_2 \leqslant S - 2M + 2K(n+1)$ and will increase by at most $3M - 10K(n+1)$. So, the third bin will be filled up to at most $S + M - 8K(n+1)$ and still will not be almost full.
In the end, there are only items of weight at most $10K(n+1)$ (now, forget about the previous terminology and call such items small) and items of weight at least $3M - 10K(n+1)$ (call them big). As per gnasher729 answer, small items can be easily distributed to the bins in the end, because $10K(n+1) < 3M/2$. Unfortunately, we cannot ignore
small items, because they affect $3S$ (the total weight of items) and, therefore, the sizes of the bins.
Fortunately, we know from the above manipulations that we can move all small items to the third (not almost full) bin. Hence, we are left with the following problem: there are $m$ big items with weights $3M - v_1$, $3M - v_2$, $\ldots$, $3M - v_m$, where $v_1 \leqslant v_2 \leqslant \ldots \leqslant v_m \leqslant 10K(n+1)$. There are also $n - m$ small items with total weight $R = 3S - (3Mm - V)$, where $V = v_1 + \ldots + v_m$. The exact weights of small items is not important.
There is barely any difference between weights of big items, because $10K(n+1)$ is negligible when compared to $M$. Moreover, the total weight of small items is also negligible. Because otherwise we could do something similar to the way we proved that all items are either big or small: we can pick small items (remember that they are all in the third bin) until their total weight exceeds $20K(n+1)$ for the first time (at such a moment, their total weight will exceed $20K(n+1)$ by at most a weight of a single small item; in other word, it will be in the range $[20K(n+1), 30K(n+1)]$) and then swap all these items with a single big item from the first or the second bin.
Hence, $R < 20K(n+1)$. Therefore, "from a bird's view", we have $m$ items of weight $3M$ and $n - m$ items of weight $0$. I claim that $m \bmod 3 = 2$. Indeed, consider the other cases.
When $m \bmod 3 = 0$, we can place arbitrary $m/3$ big items in each bin. Distribute small items arbitrarily as well. Then, the total weight
of items in each bin will not exceed $3M \cdot (m/3) + R < Mn + 20K(n+1)$. On the other hand,
\begin{equation*}
S + M - K(n+1) \geqslant (3M - 10K(n+1)) \cdot (m/3) + M - K(n+1) = Mm + (M - K(n+1) - 10Km(n+1) / 3),
\end{equation*} which is greater than $Mm + 20K(n+1)$, because $Km(n+1) \leqslant Kn(n+1)$ is negligible when compared to $M$. Hence, in this case, there is a solution with no almost full bins. Contradiction.
When, $m \bmod 3 = 1$, or $m = 3\ell + 1$, there is a bin with $\ell + 1$ items in any solution. Hence, their total weight is at least $(\ell + 1)(3M - 10K(n+1))$. On the other hand, the size of bin is $S + M \leqslant (R + (3\ell + 1) \cdot 3M) / 3 + M = R/3 + (3\ell + 2)/3 \cdot 3M$. Again, the multiplier before $3M$ is the deciding factor here, and $\ell + 1 > (3\ell + 2)/3$. Hence, there is no solution in this case, contradiction.
Finally, $m \bmod 3 = 2$, or $m = 3\ell + 2$ corresponds to the case when two bins have $\ell + 1$ big items in them and one has $\ell$. This case is actually nontrivial: not only there are inevitably exactly two almost full bins, but we also cannot distribute big items without considering their weights, like in the case $m \bmod 3 = 0$. The remaining text deals with this case.
The good thing is that we can distribute all small items and $\ell$ heaviest big items to the third bin and this is optimal. Intuitively, some bin has to contain at most $\ell$ big items, and this leaves a lot of free space left regardless of big items chosen. Formally, $\ell$ big items and all small items have total weight $3M\ell + R$, but $S + M \geqslant (3M - 10K(n+1)) \cdot (3\ell+2)/3 + M = 3(\ell + 1)m - \textrm{something small}$. Hence, there is almost $3M$ of free space left, even if we choose $\ell$ heaviest items to place in the bin (because the number of big items is more important than the actual items involved). On the other hand, other two bins have little free space left, because the total weights of items in them is $3M \cdot (\ell + 1) - \textrm{something small}$.
Hence, we need to actually compare small numbers and cannot get away with rough bounds, as we did in the first algorithmic part and in the arguments above.
Recall that we place $\ell$ heaviest big items (or $\ell$ items with the smallest values of $v_i$, which is the same) in the third bin. Hence, let $L$ be the multiset of $\ell$ smallest values of $v_i$ and $U$ be the multiset of $2(\ell + 1)$ largest values of $v_i$. Moreover, denote $u_i = v_{i+\ell}$ for $i := 1, 2, \ldots, 2(\ell + 1)$, so $U = \{u_1, u_2, \ldots, u_{2(\ell + 1)}\}$ (here, $u_i$ is non-strictly increasing).
Then, we want to distribute items of weights $3M - u_i$ for $i \in [1, 2(\ell + 1)]$ into two bins of size \begin{equation*}
S + M = (3M \cdot (3\ell + 2) - \mathrm{sum}(U) - \mathrm{sum}(L) + R))/3 + M = 3M \cdot (\ell + 1) - (\mathrm{sum}(U) + \mathrm{sum}(L) - R)/3, \end{equation*} where by $\mathrm{sum}(X)$ I mean the sum of elements in the multiset $X$.
Because we distribute exactly $\ell + 1$ items in both bins, we can focus on the space we saved when compared to $3M(\ell + 1)$. That is, we need to find a submultiset $X$ of $U$ with size $\ell + 1$, such that $\mathrm{sum}(X) \geqslant (\mathrm{sum}(U) + \mathrm{sum}(L) - R)/3$
and $\mathrm{sum}(U \setminus X) \geqslant (\mathrm{sum}(U) + \mathrm{sum}(L) - R)/3$. Here, $X$ represents the items that go to the first bin and $U \setminus X$ represents items that go the second bin.
In other words, we need to find $X$ with \begin{equation*}\mathrm{sum}(X) \in [(\mathrm{sum}(U) + \mathrm{sum}(L) - R)/3, \mathrm{sum}(U) - (\mathrm{sum}(U) + \mathrm{sum}(L) - R)/3] = [(\mathrm{sum}(U) + \mathrm{sum}(L) - R)/3, (2 \cdot \mathrm{sum}(U) - \mathrm{sum}(L) + R)/3].\end{equation*}
In the end, we need $\mathrm{sum}(X)$ to get in the gap of length
\begin{equation*}(2 \cdot \mathrm{sum}(U) - \mathrm{sum}(L) + R)/3 - (\mathrm{sum}(U) + \mathrm{sum}(L) - R)/3 + 1 = 1 + (2R + \mathrm{sum}(U) - 2\mathrm{sum}(L))/3.\end{equation*} This gap always has positive length, because $\mathrm{sum}(U) \geqslant 2\mathrm{sum}(L)$: $\min(U) \geqslant \max(L)$ by definition of $U$ and $L$ and there are $2(\ell + 1)$ elements in $U$, but only $\ell$ elements in $L$. Moreover, the gap is centered at $\mathrm{sum}(U)/2$ (both by formulas and its "physical" meaning, which must be symmetric with respect to swapping the first and the second bin).
Now, we can do the following: start from some submultiset $X$, such that $\mathrm{sum}(X) \leqslant \mathrm{sum}(U) / 2$ and gradually change it to its complement $U \setminus X$ by replacing its elements with elements of its complement, for which $\mathrm{sum}(U \setminus X) \geqslant \mathrm{sum}(U) / 2$. Then, either at some moment $\mathrm{sum}(\textrm{current subset})$ will get in the gap, leading to a valid way to place the remaining $2(\ell + 1)$ items in the first and the second bins), or we will somehow "jump over" the whole gap in a single move, meaning that the difference between two elements of $U$ is somehow greater than the length of the gap. Intuitively, the length of the gap is usually quite large. Hence, jumping over the gap is an extremely rare situation.
Not all ways to change $X$ into $U \setminus X$ work similarly well, so let us choose some fixed $X$ and some fixed way to transform $X$ into $U \setminus X$. Specifically, split all $2(\ell + 1)$ elements of $U$ into $\ell + 1$ pairs. One of the pairs is $u_{2\ell + 1} \leftrightarrow u_{2\ell + 2}$ (we pair up two largest elements), all other pairs are $u_{i} \leftrightarrow u_{i + \ell}$ (we pair up elements on the distance $\ell$, of course $i \in [1, \ell]$). Initially, $X$ contains the smaller element in each pair. If $\mathrm{sum}(X)$ is already in the gap, we won and do not have to do anything.
Otherwise, replace the elements of $X$ one-by-one by elements that they are paired up with, as long as we can do that without jumping over the gap. In the end, we are left with at least one pair that makes us jump over the whole gap (as mentioned before, we cannot reach $U \setminus X$ without either falling in the gap or jumping over it, because the gap is centered in $\mathrm{sum}(U)/2$).
I claim that only the pair $u_{2\ell + 1} \leftrightarrow u_{2\ell + 2}$ can lead to such a big jump. Assume the contrary, suppose that, for some $i \in [1, \ell]$, the jump $\Delta := u_{i + \ell} - u_i$ is more than $1 + (2R + \mathrm{sum}(U) - 2\mathrm{sum}(L))/3$ (the length of the gap, as shown before). Then, $u_{i + \ell} = u_i + \Delta \geqslant \max(L) + \Delta$. Because $i + \ell \leqslant 2\ell$, $u_{2\ell+2} \geqslant u_{2\ell+1} \geqslant u_{i + \ell} \geqslant \max(L) + \Delta \geqslant \Delta$. In the end,
\begin{equation*} \mathrm{sum}(U) \geqslant u_{2\ell + 1} + u_{2\ell + 2} + (2\ell) \cdot \max(L) + (u_{i + \ell} - \max(L)) \geqslant 2\ell \cdot \max(L) + 3\Delta \geqslant 2 \cdot \mathrm{sum}(L) + 3\Delta\end{equation*} (compared to "doubled $L$", $U$ has two extra elements, which are both at least $\Delta$, and an element that is at least $\Delta + \textrm{the corresponding element of "doubled $L$"}$). Hence, the length of the gap is \begin{equation*}1 + (2R + \mathrm{sum}(U) - 2\mathrm{sum}(L))/3 > 0 + (0 + 3\Delta)/3 = \Delta,\end{equation*} contradiction.
Therefore, the only way to "jump over" the gap is to replace $u_{2\ell + 1}$ with $u_{2\ell + 2}$. Let us do that immediately (as the first operation, while $X$ still contains $u_1$, $u_2$, $\ldots$, $u_{\ell}$). Did we jump over the gap? If we did not, we will never jump over it, so we will successfully find a solution. If we did, then the largest element of $U$ is so large, that even putting it in the same set with $\ell$ smallest elements of $U$ is already making the sum of the corresponding set too large. In this case, there is no solution, because the largest element of $U$ has to go somewhere (either to $X$, or to $U \setminus X$).
