So it is very easy. To perform an arithmetic operation among two Simple Continued Fractions (SCF) we have to express the result or the current state which will eventually become the result, as a function $z(x,y)$ where $x$ and $y$ are the operand SCFs and let $z$ be the resulting SCF.
$z(x,y)=\frac{axy+bx+cy+d}{exy+fx+gy+h}$ where;
$x=[p_0,p_1,p_2,...,p_i]$ at hand
$y=[q_0,q_1,q_2,...,q_j]$ at hand
$z=[r_0,r_1,r_2,...,r_k]$ to be constructed
Right now I suggest you to read the remaining part first and then go over it backwards once more.
Before we start ingesting coefficients from $x$ or $y$ i have to explain how to choose whether to ingest from $x$ or $y$. Apparently it just doesn't matter as long as you deplete all $x$ and $y$ coefficients. I will actually do some tests on different strategies but for the time being we may use Ross Dempsey's suggestion.
If $|\frac{b}{f}-\frac{d}{h}|>|\frac{c}{g}-\frac{d}{h}|$ then choose $x$;
Edit: I have done some tests and it really doesn't matter. My choice is to avoid a comparison and interlace from $x$ and $y$ up until one of them finishes and then continue with the operand which has remaining coefficients.
Let's assume that we start ingesting from $x$. Remember that $x=p_0+\frac{1}{x_1}$. Which means we assign $(p_0+\frac{1}{x_1})$ for $x$ in the $z(x,y)=\frac{axy+bx+cy+d}{exy+fx+gy+h}$ formula to get the next state like;
$z(x,y)=\frac{a(p_0+\frac{1}{x_1})y+b(p_0+\frac{1}{x_1})+cy+d}{e(p_0+\frac{1}{x_1})y+f(p_0+\frac{1}{x_1})+gy+h}\implies\frac{(ap+c)xy+(bp+d)x+ay+b}{(ep+g)xy+(fp+h)x+ey+f}$
Which means an ingession from $x$ always yields a transformation like;
$\left\langle\begin{matrix}a&b&c&d\\e&f&g&h\end{matrix}\right\rangle\implies\left\langle\begin{matrix}ap+c&bp+d&a&b\\ep+g&fp+h&e&f\end{matrix}\right\rangle$
Now let's assume that we ingest from $y$. Remember that $y=q_0+\frac{1}{y_1}$. Which means we assign $(q_0+\frac{1}{y_1})$ for $y$ in the $z(x,y)=\frac{axy+bx+cy+d}{exy+fx+gy+h}$ formula to get the next state like;
$z(x,y)=\frac{ax(q_0+\frac{1}{y_1})+bx+c(q_0+\frac{1}{y_1})+d}{ex(q_0+\frac{1}{y_1})+fx+g(q_0+\frac{1}{y_1})+h}\implies\frac{(aq+b)xy+ax+(cq+d)y+c}{(eq+f)xy+ex+(gq+h)y+g}$
Which means an ingession from $y$ always yields a transformation like;
$\left\langle\begin{matrix}a&b&c&d\\e&f&g&h\end{matrix}\right\rangle\implies\left\langle\begin{matrix}aq+b&a&cq+d&c\\eq+f&e&gq+h&g\end{matrix}\right\rangle$
We keep ingesting from $x$ and $y$ up until $\lfloor\frac{a}{e}\rfloor=\lfloor\frac{b}{f}\rfloor=\lfloor\frac{c}{g}\rfloor=\lfloor\frac{d}{h}\rfloor$. When this condition is satisfied it means we have discovered a coefficient of $z$ such as $r_0=\lfloor\frac{a}{e}\rfloor$. Accordingly we subtract $r_0$ from $z$ and reciprocate the result.
Remember that $z=r_0+\frac{1}{z_1}$ and we have just discovered $r_0$ which means we have to subtract $r_0$ from $z$ and calculate the egession transformation for $z_1$.
$z=r_0+\frac{1}{z_1}=\frac{axy+bx+cy+d}{exy+fx+gy+h}\implies z-r_0=\frac{1}{z_1}\implies z_1=(z-r_0)^{-1}\implies\frac{exy+fx+gy+h}{(a\bmod{e})xy+(b\bmod{f})x+(c\bmod{g})y+(d\bmod{h})}$
Remember that $r_0=\lfloor\frac{a}{e}\rfloor$ is an integer and we are interested in the remainder of the division to reciprocate. This actually means $a\bmod{e}$. The transformation of egession is;
$\left\langle\begin{matrix}a&b&c&d\\e&f&g&h\end{matrix}\right\rangle\implies\left\langle\begin{matrix}e&f&g&h\\a\bmod{e}&b\bmod{f}&c\bmod{g}&d\bmod{h}\end{matrix}\right\rangle$
After egession we may still need further egessions if the above egession requirements satisfy. Once we reach to a state at which no more egession is possible then continue with ingession from $x$ and $y$ exactly as we have done before, up until we need another egession. After some time both of the $x$ and $y$ coefficients will deplete ($x_n=p_n+\frac{1}{\infty}$) and we will eventually end up with no $x$ and $y$ terms to carry on but a simple fraction. At this moment we already have obtained some $z$ coefficients like $[r_0,r_1,..]$ and a fraction. We simply convert the fraction into SCF and concatenate it's coefficients to the previously obtained $z$ coefficients to finalize the result.
So how do we add two SCFs like $x+y$? This addition is in fact turns our $z(x,y)$ function into $z(x,y) = \frac{x+y}{1}$.
Going back to the original $z(x,y)$ function; $z(x,y) = \frac{axy+bx+cy+d}{exy+fx+gy+h}$ we can now insert $a,b,c,d,e,f,g,h$ values to start with the addition like
$\frac{0xy+1x+1y+0}{0xy+0x+0y+1}\implies\left\langle\begin{matrix}0&1&1&0\\0&0&0&1\end{matrix}\right\rangle$
This is in fact a powerful tool which can combine multiple simple arithmetic operations into one operation, all at once. Consider this $\frac{(x + 3)(y + 4)}{(x - y)} = \frac{1xy+4x+3y+12}{0xy+1x-1y+0}$. All you need to do is to start with
$\left\langle\begin{matrix}1&4&3&12\\0&1&-1&0\end{matrix}\right\rangle$
The above matrix like thingy is in fact nothing but a binary operator like $+$ or $\times$ for SFCs. However it can be morphed to take the shape of the required composition of individual arithmetic operators. I would rather call it a "Composable Binary Operator" for SCFs. Once you embed multiple operations, it just does it's things at the same number of iterations. This is very powerful.
Note: Depending on the language of your choice, the egession formula can yield wrong results especially for subtraction or other operations if negative numbers are involved. As i am told while it is still OK for Mathematica language, for the wast majority of languages it's best to replace all modulus operators with remainder logic. In short, Replace $a\bmod{e}$ with $a-r_n.e$ where $r_n$ is $\lfloor\frac{a}{e}\rfloor$.
