an integer sum of products of tangents
This question arose from my initial attempts at answering this question. I later found a way to transform the desired sum into a sum of squares of tangents, but before I did, I found numerically that apparently
$$ \sum_{l=1}^n\tan\frac{jl\pi}{2n+1}\tan\frac{kl\pi}{2n+1}=m_{jkn}(2n+1) $$
with integer factors $m_{jkn}$, for which I haven't been able to find an explanation. If $j$ or $k$ is coprime to $2n+1$, we can sum over $jl$ or $kl$ instead, so most cases (in particular all for $2n+1$ prime) can be reduced to the case $j=1$. Here are the numerically determined factors $m_{1kn}$ for $n\le18$ (with $n$ increasing downward and $k$ increasing to the right):
$$ \begin{array}{r|rr} &1&2&3&4&5&6&7&8&9&10&11&12&13&14&15&16&17&18\\\hline1&1\\ 2&2&0\\ 3&3&-1&1\\ 4&4&0&1&0\\ 5&5&-1&1&1&1\\ 6&6&0&2&-2&0&0\\ 7&7&-1&2&-1&1&0&1\\ 8&8&0&2&0&2&2&2&0\\ 9&9&-1&3&1&1&-3&1&-1&1\\ 10&10&0&3&-2&2&-1&1&0&1&0\\ 11&11&-1&3&-1&1&-1&1&3&-1&1&1\\ 12&12&0&4&0&2&0&0&-4&0&0&0&0\\ 13&13&-1&4&1&3&0&1&-1&1&1&3&0&1\\ 14&14&0&4&-2&2&2&2&0&2&4&0&0&2&0\\ 15&15&-1&5&-1&3&-3&3&-1&3&-5&1&1&1&-1&1\\ 16&16&0&5&0&2&-1&2&0&1&-2&1&1&-2&-2&1&0\\ 17&17&-1&5&1&3&-1&2&-1&1&-1&1&5&1&0&1&1&1\\ 18&18&0&6&-2&4&0&2&0&2&-2&2&-6&0&0&4&2&0&0\\ \end{array} $$
(See also the table in this answer to the other question, which shows the case $j=k+1$; in that case the rows of the table sum to $0$ because of the identity that's the subject of the other question.)
The values $m_{11n}=n$ reflect the sum of squares of tangents that I determined in my answer to the other question. I have no explanation for the remaining values. I've tried using the product formula for the tangent; multiplying by a third tangent to use the triple tangent product formula; and finding a polynomial whose roots are the products being summed; but none of that worked out.
This vaguely reminds me of character theory; the values $\tan\frac{kl\pi}{2n+1}$ for fixed $k$ are like characters, and their dot products are integer multiples of the "group order" $2n+1$; though if they were characters the dot products couldn't be negative.
I'd appreciate any insight into this phenomenon, and of course ideally a way to calculate the $m_{jkn}$.
[Update:]
I've verified the periodicities that Brian observed in comments up to $n=250$:
$$m_{1,k,n+k} = m_{1kn}+[k \text{ odd}]\;,$$
$$m_{1,k+4d+2,k+4d+2+d}=m_{1,k,k+d}\;,$$
where the bracket is the Iverson bracket.
[Just as I was finishing this I saw user8268's answer. I suspect the explanations are related.]
I wasn't intending to answer my own question in this case, but I've now found an explanation. The character theory analogy turned out to be more useful than I expected. Thinking of the values $\tan\frac{kl\pi}{2n+1}$ for fixed $k$ as vectors composed of integer multiples of mutually orthogonal vectors made me wonder what these mutually orthogonal vectors might be. A natural choice was a Fourier-style set of sines or cosines, and indeed it turns out that
$$ \sum_{l=1}^{2n}\sin\frac{2jl\pi}{2n+1}\tan\frac{kl\pi}{2n+1}= \begin{cases} \pm(2n+1)&\gcd(k,2n+1)\mid j\;,\\ 0&\text{otherwise.} \end{cases} $$
I found this surprising at first, but it's actually not too difficult to explain. We have
$$ \def\ex#1{\mathrm e^{#1}}\def\exi#1{\ex{\mathrm i#1}}\def\exm#1{\ex{-\mathrm i#1}} \ex{2n\mathrm i\phi}-\ex{-2n\mathrm i\phi}=\left(\exi{\phi}+\exm{\phi}\right)\left(\ex{(2n-1)\mathrm i\phi}-\ex{(2n-3)\mathrm i\phi}+\dotso+\ex{-(2n-3)\mathrm i\phi}-\ex{-(2n-1)\mathrm i\phi}\right)\;, $$
so
$$ \sin(2j\phi)=2\cos\phi\left(\sin((2j-1)\phi)-\sin((2j-3)\phi)+\dotso+(-1)^{j+1}\sin\phi\right)\;. $$
Thus, for $k=1$ the cosine in the denominator of the tangent is cancelled, and the remaining sine picks out the last term in the sum of alternating sines with odd frequencies, which yields
$$ \sum_{l=1}^{2n}\sin\frac{2jl\pi}{2n+1}\tan\frac{l\pi}{2n+1}=(-1)^{j+1}(2n+1)\;. $$
But since the integers $jl$ and $kl$ in the arguments of both factors only matter $\bmod2n+1$, if $k$ is coprime to $2n+1$, we can sum over $kl$ instead of $l$ and will get the result for $k^{-1}j\bmod(2n+1)$, so for $k$ coprime to $2n+1$
$$ \sum_{l=1}^{2n}\sin\frac{2jl\pi}{2n+1}\tan\frac{kl\pi}{2n+1}=(-1)^{\sigma_k(j)+1}(2n+1)\;, $$
where $\sigma_k$ is the permutation effected by multiplication with $k^{-1}\bmod(2n+1)$. If $1\lt\gcd(k,2n+1)\mid j$, the sum reduces to $\gcd(k,2n+1)$ identical copies, whereas if $\gcd(k,2n+1)\nmid j$, cancellation lets the sum vanish.
Thus, the $m_{jkn}$ are integers because the vectors $\left\{\tan\frac{kl\pi}{2n+1}\right\}_l$ are integer linear combinations of vectors $\left\{\sin\frac{2jl\pi}{2n+1}\right\}_l$ whose dot products are all either $0$ or $\pm1$, and these values can be obtained using only elementary number theory, namely permutations induced by multiplicative inverses.
warning: I suppose that $p=2n+1$ is prime and only prove that $m_{jkn}$ is an integer, without actually computing it (the assumption can perhaps be removed, but the method is non-explicit).
We have $a_k:=i\tan\frac{k\pi}{p}=\frac{\alpha^k-1}{\alpha^k+1}$, where $\alpha=\exp \frac{2\pi i}{p}$. Your sum is in $K:=\mathbb{Q}(\alpha)$ and is invariant under the Galois group, so it is a rational number.
On the other hand, $a_k$'s are the roots of $(1+x)^{p}-(1-x)^{p}$. Now notice that $f(x):=((1+x)^{p}-(1-x)^{p})/2x$ is a monic polynomial with integer coefficients. In particular, your sum is an algebraic integer, and as it is rational, it is an integer.
As $f(0)=p$, we thus have $N_{K/\mathbb{Q}}a_k=p$. This means that $a_k$ is a unit times $1-\alpha$, so your sum is divisible by $1-\alpha$, which means (as it is rational) that your sum is a multiple of $p$.