Simplest All Digital GMSK Modulator
I don't actually have an answer, just ideas, but since this is a riddle, I might of contribute these as (hidden as spoiler text) hints for others, which might have ideas that I miss:

So, MSK is often introduced as


    represented by a half-symbol time-offset BPSK on the I and Q components. That comes in very handy here (avoiding multiplicators) – BPSK can be implemened as $\pm 1$ on the respective axis. The idea is that at the boundary between symbols periods, the phase is continuous, which w.l.o.g cannot be realized in any other way than by having zeros exactly at the point where the other component carries the BPSK constellation point. Thus, the minimal viable MSK modulation has


&nbsp;


   $$\begin{align} I&amp;= [- {s_1} &amp;&amp; 0 &amp;&amp; -{s_2} &amp;&amp; 0 &amp;&amp; -{s_3}&amp;\dots ]\\Q&amp;= [ 0 &amp;&amp; s_1 &amp;&amp; 0 &amp;&amp; s_2 &amp;&amp; 0 &amp; \dots ] \end{align}$$


&nbsp;


   as samples, with $s_n$ being the bits to be transmitted represented as $\pm 1$ differentially. That doesn't look very much like CPFSK at all, but one has to realize that the advance between two consecutive samples takes two different values:




   For samples that have an even index $n$, the following sample always is -90° further, no matter which value $s_{n/2}$ has


   For samples that have an odd index $n$, the phase advance is -90° if $s_{\lfloor n/2\rfloor}=s_{(n+1)/2}$, and +90° in case the bits differ.



&nbsp;


   Considering the phase difference of consecutive samples as the instantaneous frequency, we thus see that we get an alternation of a constant $-\frac \pi2$ between samples and either $0$ or $-\frac \pi2$, carrying the information. I'm not even convinced any more this fulfills the criteria for being MSK.


&nbsp;


   Since we'd probably want spectral symmetry, we'd employ the usual tricks to shift the signal in frequency by half the sampling rate (i.e. "multiplication" with $[-1\,+1\,-1\,+1\,\dots]$).


Oooh and I'd have an approach for the Gaussian pulse shaper, but that simply applies


   the fact that repeated convolution of something with itself converges against a Gaussian – same reason the CLT works.


 That means I can just cascade a set of moving average filters – all of which work without multipliers to get an approximated Gauss filter.


Sadly, that doesn't inherently solve the issue


   that my MSK approach above basically only allows 1 sample per symbol – but we might at least recreate the inter-symbol dependence characteristic of GMSK trellises.