I did it for the DCT-I and DCT-II.
At first I thought it is about circular convolution, but it is not. After desperate attempts to do it by myself I found the article: Convolution Using Discrete Sine and Cosine Transforms (Alternative Link).
In this paper there is a derivation of circular convolution using DTT's, but it is constructed from sum of $C_1^{-1}$ and $S_1^{-1}$, so, it's not what we are interested in.
What I suggest is to look at Fourier transform first, in order to get the idea:
Let we have
$\displaystyle F[t] = \sum_{k=0}^{N-1} f[k]e^{-\frac{i2π}{N}tk}$,
$\displaystyle G[t] = \sum_{k=0}^{N-1} g[k]e^{-\frac{i2π}{N}tk}$.
So then
$$\displaystyle ℱ\{f*g\}[t] = \sum_{x=0}^{N-1} f[x]e^{-\frac{i2π}{N}tx}\cdot\sum_{y=0}^{N-1} g[y]e^{-\frac{i2π}{N}ty} = \sum_{x=0}^{N-1}\sum_{y=0}^{N-1} f[x]g[y]e^{-\frac{i2π}{N}t(x+y)}$$
Now, let $x+y = m$, $m = 0\ldots2N-2$ and then group all $(x,y)$ which share the same $m$. So you will get:
$$\displaystyle \sum_{m=0}^{2N-2}\biggl(\sum_{0\leq x,y\leq N-1\\x+y=m}f[x]g[y] \biggl) e^{-\frac{i2π}{N}tm} \tag{1}$$
Because $x, y$ are both positive and less than $N$, then sum inside is just linear convolution at $m$. Let me denote it as $A[m]$, then everything is going to be rewritten in that form
$$\sum_{m=0}^{2N-2} A[m] e^{-\frac{i2π}{N}tm} \tag{2}$$
DFT kernel is N-periodic that is why we can group terms another time:
$$\displaystyle \sum_{m=0}^{N-1} (A[m] + A[m+N])e^{-\frac{i2π}{N}tm}$$
$A[m]$ captures $A[0]\ldots A[N-1]$, $A[m+N]$ captures $A[N]\ldots A[2N-1]$. But remember that $A[m]$ is linear convolution of two sequences of size $N$ and therefore size is $N + N - 1$ then $A[2N-1] = 0$. So these two sums are equivalent. Also it is a DFT of $A[m] + A[m+N]$. What is it? Right, circular convolution, as expected.
Now plug in DHT kernel instead of DFT into $(1)$ and you will get the same result ─ circular convolution. That is because DHT and DFT kernels are quite similar.
Well, if you have 2 signals $f$, $g$ with lengths $Q$ and $P$, and $A = f*g$, with $Q+P-1 = L$. Now, add zeros to $f$ up to $L$, do the same for $g$. Now $f$, $g$, and $A$ are $L$ periodic signals, so we can take the DFT of them. DFT of $A$ is going to be:
$$\sum_{m=0}^{L-1} A[m] e^{-\frac{i2π}{L}tm} = \sum_{m=0}^{L-1}\biggl(\sum_{0\leq x,y\leq L-1\\x+y=m}f[x]g[y] \biggl) e^{-\frac{i2π}{L}tm} = \sum_{x=0}^{L-1}\sum_{y=0}^{L-1}f[x]g[y]e^{-\frac{i2π}{L}t(x+y)} \tag{3}$$
So, via $(2)$ we can get the rule and expression for the convolution and $(3)$ shows that zero-padding up to convolution length is required to get linear convolution using the same rule. We're going to apply it to the DCT's and DST's.
Because there are many varieties of DCT's and DST's, I will mention what I used:
DCT-I (N+1) points:
$\displaystyle C_1\{f\}[n] = \sum_{k=0}^{N} c_k f[k]\\cos\Bigl[\frac{π}{N}nk\Bigl], \text{where}\
c_k =
\left\{
\begin{array}{l}
\tfrac12, \qquad &k (\mod N)=0, \\
1,\qquad &\text{otherwise} \\
\end{array}
\right. \tag{4}$
DST-I (N+1) points:
$$S_1\{f\}[n] = \sum_{k=0}^{N} f[k]\sin\Bigl[\frac{π}{N}nk\Bigl]=\sum_{k=0}^{N} c_k f[k] \sin\Bigl[\frac{π}{N}nk\Bigl] \tag{5}$$
(N-1) for data: $f[1] \ldots f[N-1]$ it doesn't matter what is placed in $f[0]$ and $f[N]$.
DCT-II (N+1) points:
$$C_2\{f\}[n] = \sum_{k=0}^{N} f[k]\cos\Bigl[\frac{π}{N}\Bigl(k+\tfrac12\Bigl)n\Bigl] \tag{6}$$
N for data: $f[0]\ldots f[N-1]$ and $f[N] = 0$
DST-II (N+1) points, N for data:
$f[0]\ldots f[N-1]$ and $f[N] = 0$:
$$S_2\{f\}[n] = \sum_{k=0}^{N} f[k]\sin\Bigl[\frac{π}{N}\Bigl(k+\tfrac12\Bigl)n\Bigl] \tag{7}$$
Inverse DCT-I:
$$C_1^{-1} = \frac{2}{N}C_1 \tag{8}$$
Inverse DST-I:
$$S_1^{-1} = \frac{2}{N}S_1 \tag{9}$$
Inverse DCT-II (N+1) points:
$$C_2^{-1}\{F_{C_2}\}[k] = \frac{2}{N}\sum_{n=0}^{N} c_n F_{C_2}[n]\cos\Bigl[\frac{π}{N}n\Bigl(k+\tfrac12\Bigl)\Bigl] \tag{10}$$
$C_2^{-1}\{C_2\{f\}\}[k] = f[k]$ except for $k=N$ there it is going to be $f[N-1]$.
Inverse DST-II (N+1) points:
$$C_2^{-1}\{F_{S_2}\}[k] = \frac{2}{N}\sum_{n=0}^{N} c_n F_{S_2}[n]\sin\Bigl[\frac{π}{N}n\Bigl(k+\tfrac12\Bigl)\Bigl] \tag{11}$$
$S_2^{-1}\{S_2\{f\}\}[k] = f[k]$ except for $k=N$ there it is going to be $-f[N-1]$. Also, they are not interchangeable.
Everything mentioned above is just Wikipedia ones, but modified to achieve more similar look for easier coding and derivation.
So, now again, applying the rule $(2)$ for the DCT-I we get:
$$\sum_{m=0}^{2N} A[m] c_m \cos\Bigl[\frac{π}{N}nm\Bigl] = \sum_{m=0}^{N} p_m(A[m] + A[2N-m]) c_m \cos\Bigl[\frac{π}{N}nm\Bigl]$$
$p_m = \tfrac12$ if $m = N$ else $p_m = 1$
from the other hand we have
$$\sum_{m=0}^{2N} \frac{A[m]}{c_m} c_m \cos\Bigl[\frac{π}{N}nm\Bigl] = \sum_{x=0}^{N}\sum_{y=0}^{N} f[x]g[y]\cos\Bigl[\frac{π}{N}n(x+y)\Bigl]$$
$$\sum_{x=0}^{N}\sum_{y=0}^{N} f[x]g[y]\cos\Bigl[\frac{π}{N}n(x+y)\Bigl] = \\\sum_{x=0}^{N}\sum_{y=0}^{N} f[x]g[y]\biggl(\cos\Bigl[\frac{π}{N}nx\Bigl]\cos\Bigl[\frac{π}{N}ny\Bigl] - \sin\Bigl[\frac{π}{N}nx\Bigl]\sin\Bigl[\frac{π}{N}ny\Bigl]\biggl) = \\ \biggl(\sum_{x=0}^{N} f[x]\cos\Bigl[\frac{π}{N}nx\Bigl]\biggl)\cdot\biggl(\sum_{y=0}^{N} g[y]\cos\Bigl[\frac{π}{N}ny\Bigl]\biggl) - \biggl(\sum_{x=0}^{N} f[x]\sin \Bigl[\frac{π}{N}nx\Bigl]\biggl)\cdot\biggl(\sum_{y=0}^{N} g[y]\sin\Bigl[\frac{π}{N}ny\Bigl]\biggl)$$
It is now easy to see that it can be rewriten in a such way:
$$C_1\{\frac{f[m]}{c_m}\}[n]\cdot C_1\{\frac{g[m]}{c_m}\}[n] - S_1\{f[m]\}[n]\cdot S_1\{g[m]\}[n]$$
Let $\frac{f[m]}{c_m} = f^{/}$ then we can write everything in this form:
$$\frac{p_m}{c_m}(A[m] + A[2N-m]) ⟷ F^{/}_{C_1}G^{/}_{C_1} - F_{S_1}G_{S_1}$$
alternative version
$$δ_{0m}(A[0] + A[2N]) + (A[m] + A[2N-m]) ⟷ F^{/}_{C_1}G^{/}_{C_1} - F_{S_1}G_{S_1}$$
or even
$$(A[m] + A[2N-m]) ⟷ F^{/}_{C_1}G^{/}_{C_1} - F_{S_1}G_{S_1} - \frac{f[0]g[0] + f[N]g[N]}{2}$$
Now, remember the $(3)$? If you add zeroes, so the size of f and g is equal to the convolution's length, then you have:
$$(f*g)^{/}⟷ F^{/}_{C_1}G^{/}_{C_1} - F_{S_1}G_{S_1}$$
It looks a little bit different to what is presented in the paper, but it is because I used different version of DCT.
For DCT-II we have the rule (remember that $f[N] = g[N] = 0$):
$$\sum_{m=0}^{2N} A[m] \cos\Bigl[\frac{π}{N}\Bigl(m+\tfrac12\Bigl)n\Bigl] = \sum_{m=0}^{2N-2} A[m] \cos\Bigl[\frac{π}{N}\Bigl(m+\tfrac12\Bigl)n\Bigl] = \sum_{m=0}^{N-1} (A[m] + A[2N-1-m]) \cos\Bigl[\frac{π}{N}\Bigl(m+\tfrac12\Bigl)n\Bigl] $$
And the expression is going to be:
$$\sum_{m=0}^{2N} A[m] \cos\Bigl[\frac{π}{N}\Bigl(m+\tfrac12\Bigl)n\Bigl] =
\sum_{m=0}^{2N} A[m] \Bigl(\cos\Bigl[\frac{πn}{2N}\Bigl]\cos\Bigl[\frac{π}{N}kn\Bigl] - \sin\Bigl[\frac{πn}{2N}\Bigl]\sin\Bigl[\frac{π}{N}kn\Bigl]\Bigl) = \\
\cos\Bigl[\frac{πn}{2N}\Bigl](F^{/}_{C_1}G^{/}_{C_1} - F_{S_1}G_{S_1}) - \sin\Bigl[\frac{πn}{2N}\Bigl](F^{/}_{C_1}G_{S_1} + F_{S_1}G^{/}_{C_1}) = \\
=F^{/}_{C_1}(\cos\Bigl[\frac{πn}{2N}\Bigl]G^{/}_{C_1} - \sin\Bigl[\frac{πn}{2N}\Bigl]G_{S_1}) - F_{S_1}(\cos\Bigl[\frac{πn}{2N}\Bigl]G_{S_1}+ \sin\Bigl[\frac{πn}{2N}\Bigl]G^{/}_{C_1}) = \\
=F^{/}_{C_1}G_{C_2} - F_{S_1}G_{S_2} = G^{/}_{C_1}F_{C_2} - G_{S_1}F_{S_2} $$
Now we finally can write it as:
$$A[m] + A[2N-1-m] ⟷ F^{/}_{C_1}G_{C_2} - F_{S_1}G_{S_2}$$
And again, if you add enough zeroes, then it will be:
$$f*g ⟷ F^{/}_{C_1}G_{C_2} - F_{S_1}G_{S_2}$$
You can play with it here: https://repl.it/@dogtyping/dct#main.py — to make sure that it does work.
Thanks for the question, it was a lot of fun.
