Algorithm for 1d spline interpolation suitable for 8 bit microcontroler What is a concise, fast, down to earth algorithm for doing (or closely approximating) spline interpolation on a 1d continuous stream of data?
(Edit1: The paragraph below equates to saying &quot;the data is uniform in interval.&quot; in an awkward way.)
The data is 1d in that sampling on x is at fixed, regular intervals (likely a power of 2 constant) known well ahead of time. This makes y the only variant, which should allow for quite a bit of simplification and precomputation (LUT?).
Here is a graphical example of approximately what we're trying to do. It's our algo applied to a coarsely sampled sine function.
(Edit2: Note that this is only an example of what our algo should do with coarsely sampled sine data, however, the data we would like to process will be random in the set [0,255], most likely 3 or 4 points)

(source: mathworks.com)
Assume high accuracy is not required, but that we must compute the number of results (red dots) between samples (blue circles) in less than 8ms (~120Hz). Also assume the signal processor available is limited in power and instruction set (PIC/AVR), so has only the following relevant instructions and limitations.

(signed + carry) Addition/subtraction instructions.
(unsigned 8x8-to-16) Multiplication instruction. (PIC18, megaAVR only)
Byte wide Boolean instructions (XOR, OR, NOR, AND, NAND, etc.)
Single bit left and right logical shifts. (no barrel shifter)
Can only execute at 2~4 MIPS


Additional Notes:

Responses would preferentially be in pseudocode, so they are more generally applicable.
Cheating is totally OK; it doesn't need to be perfect, just better than liner interpolation.
Bonus points for alternatives that don't need multiplication.
More bonus points for alternatives likely to complete in less than 1ms!

This is for an RGB mood lamp x-mas present for my sister and mommy :3, which I would do myself, but the maths for this are apparently beyond me.

Edit 12-21-2016: Better list formatting
Take a look at the cubic Hermite spline.  The interpolated function is continuous at the data points and the first derivative is also continuous.  Away from the data points all of the derivatives are continuous.

Let's say that the function $f(x)$ is defined by equally-spaced data points for all $x$ that is an integer.  This means you know the values of $f(0), f(1), f(2), ...$

Then separate $x$ into an integer and fractional parts:

$$ x \triangleq n+u $$

where 

$$ n = \lfloor x \rfloor = \operatorname{floor}(x) $$

and

$$ u = x - n \quad \text{  ,  } \quad 0 \le u &lt; 1  $$

$$ $$

$$\begin{align}
f(n+u) &amp; = 
\begin{bmatrix}
1 &amp; u &amp; u^2 &amp; u^3 \\
\end{bmatrix}
\cdot
\begin{bmatrix}
0 &amp; 1 &amp; 0 &amp; 0 \\
-\tfrac12 &amp; 0 &amp; \tfrac12 &amp; 0 \\
1 &amp; -\tfrac52 &amp; 2 &amp; -\tfrac12 \\
-\tfrac12 &amp; \tfrac32 &amp; -\tfrac32 &amp; \tfrac12 \\
\end{bmatrix}
\cdot
\begin{bmatrix} f(n-1) \\ f(n) \\ f(n+1) \\ f(n+2) \end{bmatrix} \\
\\
&amp; =  \frac12 \begin{bmatrix} -u^3 +2u^2 - u \\ 3u^3 - 5u^2 + 2 \\ -3u^3 + 4u^2 + u \\ u^3 - u^2 \end{bmatrix}^T \cdot \begin{bmatrix} f(n-1)\\f(n)\\f(n+1)\\f(n+2) \end{bmatrix} \\
\\
&amp; = \frac12 \begin{bmatrix} u ((2-u) u-1) \\ u^2 (3 u-5)+2 \\ u ((4-3 u) u+1) \\ u^2 (u-1) \end{bmatrix}^T \cdot \begin{bmatrix} f(n-1)\\f(n)\\f(n+1)\\f(n+2) \end{bmatrix} \\
\\
&amp; = \tfrac12 \bigg( (u^2(2-u)-u)f(n-1) \ + \ (u^2(3u-5)+2)f(n) \\
&amp; \quad \quad \quad \quad + \ (u^2(4-3u)+u)f(n+1) \ + \ u^2(u-1)f(n+2) \bigg) \\
\end{align}$$

Is this enough information for how to do this in your PIC?  You need to be able to separate into integer and fractional parts and you need to be able to multiply.  

in my opinion, Olli's method [now moved into its own answer] is not the best way of looking at it for the OP's case of implementing this simply in a PIC.  [his formulation] separates the four data points and computes four coefficients that are attached to powers of $u$.  that's the way to do it if your fractional ordinate is any arbitrary value that is $0 \le u &lt; 1$.  but the OP has only a few values like $u=0, \tfrac14, \tfrac12, \tfrac34$.  or maybe 8 multiples of $\tfrac18$.

so my recommendation is to compute the vaules of these four polynomials:

$$  c_{-1} = \tfrac12 (-u^3 +2u^2 - u) \\
    c_0 = \tfrac12 (3u^3 - 5u^2 + 2) \\
    c_1 = \tfrac12 (-3u^3 + 4u^2 + u) \\
    c_2 = \tfrac12 (u^3 - u^2)  $$

and do that for every fractional value of $u$ (such as $u=0, \tfrac14, \tfrac12, \tfrac34$) that you will use many, many times.

then the code in the PIC only needs to implement is a dot product between the 4 data points and the selected set of coefficients:

$$ f(x) = f(n+u) = c_{-1} f(n-1) + c_0 f(n) + c_1 f(n+1) + c_2 f(n+2) $$

since $c_{-1}$ and $c_2$ can be shown to always be negative for $0 &lt; u &lt; 1$, then put into the table their absolute values and subtract their terms:

$$ f(x) = f(n+u) = c_0 f(n) + c_1 f(n+1) - (-c_2) f(n+2) - (-c_{-1}) f(n-1) $$

the coefficients stored will be 256 times bigger than their actual value (and stored as 8-bit unsigned integers) then, after multiplying and accumulating your answer (that is 256 times too big), you add 128 (for rounding) and shift right 8 bits (which is the same as taking the answer out of the higher-order byte of the product).