Algorithm for 1d spline interpolation suitable for 8 bit microcontroler
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).