Limit of an arctangent summation 
Evaluate:
  $$S=\tan\left(\sum_{n=1}^{\infty}\arctan\left(\frac{16\sqrt{2557}n}{n^4+40916}\right)\right)$$

For $x$ small enough, $\arctan(x)<x$, and by the comparison test we know that the inner series is convergent (but we may end up with something like $\tan(\pi/2)$ ).
My best idea so far has been to use the formula for $\tan(x+y)$, to define a recursive function $f(n)$:
$$f(1)=\frac{16\sqrt{2557}}{40917}$$
$$f(n)=\frac{f(n-1)+\frac{16\sqrt{2557}n}{n^4+40916}}{1-f(n-1)\frac{16\sqrt{2557}n}{n^4+40916}}$$
And then calculate $S=\lim_{n\to\infty}f(n)$, but $f(n)$ seems to behave wildly. This is due to the fact that
$$g(x)=\sum_{n=1}^{x}\arctan\left(\frac{16\sqrt{2557}n}{n^4+40916}\right)$$
Passes near $\pi/2$.
Update:
Numerical calculation seems to hint that $\lim_{x\to\infty}g(x)=\pi$,  Some ideas?
 A: Method 1 - The quick but ugly way.
Consider following sequence $$a_n = \frac{8n(n-1)+12}{n^4-2n^3+3n^2-2n-40910}$$
By brute force, one can verify
$$\frac{a_{n} - a_{n+1}}{1 + 2557 a_n a_{n+1}} = \frac{16n}{n^4+40916}$$
From this, we can conclude up to some integral multiples of $\pi$, we have
$$\tan^{-1}\left(\frac{16\sqrt{2557}n}{n^4+40916}\right)
= \tan^{-1}(a_n\sqrt{2557}) - \tan^{-1}(a_{n+1}\sqrt{2557}) + N_n \pi \quad\text{ with } N_n \in \mathbb{Z}$$
The whole mess is a telescoping series and
$$\sum_{n=1}^\infty\tan^{-1}\left(\frac{16\sqrt{2557}n}{n^4+40916}\right) = \tan^{-1}(a_1\sqrt{2557}) + N\pi \quad\text{ with } N \in \mathbb{Z}\\
\implies S = a_1\sqrt{2557} = -\frac{6\sqrt{2557}}{20455}
$$
Note
In case people wonder how can one discover such a horrible looking sequence $a_n$. I first compute the series by expressing it in terms of a bunch of gamma function (see method 2 below). I then compute a few values of the series start at different points and then look for pattern which allow me to rewrite the whole mess as a telescoping series.

Method 2 - The hard way using Gamma functions.
To compute $S$ directly, we will use following representation of $\tan^{-1} x$ over the real axis:
$$\tan^{-1} x = \Im \log (1 + i x),\quad\forall x \in \mathbb{R}$$
The $\log(\cdots)$ here stands for the principal branch of the logarithm function over $\mathbb{C}$ with a branch cut along the negative real axis. For any $\alpha, \beta > 0$, we have
$$\begin{align}
  \tan\left[\sum_{n=1}^\infty \tan^{-1}\left( \frac{\alpha n}{n^4 + \beta}\right) \right]
= & \tan\Im\left[\sum_{n=1}^\infty\log\left(1 + i\frac{\alpha n}{n^4+\beta}\right)\right]\\
= &\tan\Im\left[\sum_{n=1}^\infty\left(\log(n^4 + i\alpha n + \beta) - \log(n^4+\beta)\right) \right]\\
= & \tan\Im\left[\sum_{n=1}^\infty\left(\log(n^4 + i\alpha n + \beta) - \log(n^4)\right) \right]\\
= &\tan\Im\left[\sum_{n=1}^\infty\sum_{k=1}^4\left(\log(n-\lambda_i)-\log n\right)\right]\\
= &\tan\Im\left[\sum_{k=1}^4\sum_{n=1}^\infty\log\left(1-\frac{\lambda_i}{n}\right)\right]
\end{align}\tag{*1}
$$
where $\lambda_1,\ldots\lambda_4$ are the $4$ roots of the polynomial
$x^4 + i\alpha x + \beta$. Since the coefficient for the cubic term of this polynomial 
is zero, we have
$$\sum_{i=1}^4 \lambda_i = 0$$
Recall the infinite product expansion of Gamma function:
$$\frac{1}{\Gamma(z)} = z e^{\gamma z}\prod_{n=1}^\infty\left(1 + \frac{z}{n}\right)e^{-\frac{z}{n}}$$
We can rewrite the last expression of $(*1)$ as
$$\tan\Im\left\{\sum_{k=1}^4\left[
-\gamma\lambda_i +\sum_{n=1}^\infty
\left(\log\left(1-\frac{\lambda_i}{n}\right) + \frac{\lambda_i}{n}\right)
\right]\right\}
=-\tan\Im\left[\sum_{i=1}^4\log\Gamma(1-\lambda_i)\right]$$
This leads to 
$$(*1) = -\tan\Im\log\Delta = -\frac{\Im\Delta}{\Re\Delta}
\quad\text{ where }\quad\Delta = \prod_{i=1}^4 \Gamma(1-\lambda_i)
$$
For the question at hand where $\alpha = 16\sqrt{2557}$ and $\beta = 40916$, the roots of
$x^4 + i\alpha x + \beta$ has the form
$$-1 + u, -1 -u, 1 - \bar{u}, 1 + \bar{u}\quad\text{ with }\quad u^2 = -1 + 4\sqrt{2557}i$$
Substitute this into above expression of $\Delta$, we find
$$\begin{align}
\Delta &= \Gamma(2-u)\Gamma(2+u)\Gamma(\bar{u})\Gamma(-\bar{u})
= |\Gamma(u)\Gamma(-u)|^2 u^2(u^2-1)\\
&= -|\Gamma(u)\Gamma(-u)|^2 (40910 + 12\sqrt{2557}i)
\end{align}
$$
As a result, the $S$ we seek is equal to
$$S = -\frac{\Im \Delta}{\Re\Delta} = - \frac{12\sqrt{2557}}{40910} = -\frac{6\sqrt{2557}}{20455}$$
A: If we consider that:
$$\arctan(x)-\arctan(y)=\arctan\frac{x-y}{1+xy}, $$
by taking $x_n=\frac{1}{n^2-n+1}$ and $y_n=x_{n+1}=\frac{1}{n^2+n+1}$ we have that:
$$ \sum_{n\geq 1}\arctan\frac{2n}{2+n^2+n^4}=\sum_{n\geq 1}\left(\arctan\frac{1}{n^2-n+1}-\arctan\frac{1}{n^2+n+1}\right)=\arctan 1=\frac{\pi}{4}$$
and by taking $x_n=\frac{1}{n^2-n}$ we also have:
$$ \sum_{n\geq 1}\arctan\frac{2n}{1-n^2+n^4}=\sum_{n\geq 1}\left(\arctan\frac{1}{n^2-n}-\arctan\frac{1}{n^2+n}\right)=\frac{\pi}{2}$$
with the same telescoping trick. I was not able to adapt this trick to the series we are dealing with, but maybe it is the good way, and for sure it is worth mentioning.
