<s>[INST] Evaluating $\int_0^1 \frac{\mathrm dx}{(x^2+ax+1)^{n+1}}$ with real methods I would like to know if there are any (preferably easier) methods of evaluating
$$Q_n(a)=\int_0^1 \frac{\mathrm dx}{(x^2+ax+1)^{n+1}}$$
With real methods. Here's the way I did it.
Complete the square: 
$$Q_n(a)=4^{n+1}\int_0^1\frac{\mathrm dx}{((2x+a)^2+4-a^2)^{n+1}}$$
Then $u=2x+a$ gives
$$Q_n(a)=2^{2n+1}\int_a^{a+2}\frac{\mathrm du}{(u^2+4-a^2)^{n+1}}$$
Then consider the indefinite integral 
$$F_n^{w}(x)=\int\frac{\mathrm dx}{(x^2+w)^{n+1}}$$
Integration by parts yields the recurrence relation (for $n\in\Bbb Z\geq 1$)
$$F_n^{w}(x)=\frac{x}{2wn(x^2+w)^n}+\frac{2n-1}{2wn}F_{n-1}^{w}(x)$$
With the base case
$$F_0^{w}(x)=\frac1{\sqrt{w}}\arctan\frac{x}{\sqrt{w}}$$
And since this recurrence is in the form $$f_n=\alpha_n+\beta_nf_{n-1}$$
the solution to which is $$f_n=f_0\prod_{k=1}^{n}\beta_k+\sum_{k=0}^{n-1}\alpha_{n-k}\prod_{j=1}^{k}\beta_{n-j+1}$$
We have (I omit the simplification steps)
$$F_n^{w}(x)=\frac{{2n\choose n}}{2^{2n}w^{n+1/2}}\arctan\frac{x}{\sqrt{w}}+S_n^w(x)$$
With $$S_n^w(x)=\frac{x}{2w}\sum_{r=0}^{n-1}\frac{(x^2+w)^{r-n}}{(2w)^r}R_r^{(n)}$$
and $$R_r^{(n)}=\frac1{n-r}\prod_{j=1}^{r}\frac{2n-2j+1}{n-j+1}$$
So we have that 
$$Q_n(a)=2^{2n+1}\left[F_n^{4-a^2}(a+2)-F_n^{4-a^2}(a)\right]$$
Which is, after some simplification,
$$\begin{align}
Q_n(a)=&\frac{2{2n\choose n}}{(4-a^2)^{n+1/2}}\left[\arctan\sqrt{\frac{2+a}{2-a}}-\arctan\frac{a}{\sqrt{4-a^2}}\right]\\
&+\frac1{4-a^2}\sum_{r=0}^{n-1}\frac{2^rR_r^{(n)}}{(4-a^2)^r}\left[(2+a)^{r-n+1}-a\right]
\end{align}$$
 [/INST] $$Q_n(a)=\int_0^1\frac{dx}{(x^2+ax+1)^{n+1}}$$
and since $x^2+ax+1=\left(x+\frac a2\right)^2+\left(1-\frac{a^2}{4}\right)$. If we let $u=x+\frac a2$ and $\alpha^2=1-\frac{a^2}{4}$ we get similar to what you got:
$$Q_n(a)=\int_{a/2}^{1+a/2}\frac{du}{(u^2+\alpha^2)^{n+1}}$$
now by observing that: $u^2+\alpha^2=\alpha^2\left(\left[\frac u\alpha\right]^2+1\right)$ and letting $v=\frac u\alpha$ we obtain:
$$Q_n(a)=\int_{a/2}^{1+a/2}\frac{\alpha^{-(2n+1)}dv}{(v^2+1)^{n+1}}$$
now by letting $v=\tan\omega$ and $\beta=\arctan(1/2),\gamma=\arctan(1+a/2)$ we get:
$$Q_n(a)=\alpha^{-(2n+1)}\int_\beta^\gamma\frac{\sec^2\omega}{(\sec^2\omega)^{n+1}}d\omega=\alpha^{-(2n+1)}\int_\beta^\gamma\cos^{2n}(\omega)d\omega$$
Now by making the substitution $\sigma=\cos\omega$ we get:
$$Q_n(a)=-\alpha^{-(2n+1)}\int_{\cos\beta}^{\cos\gamma}\sigma^{2n}(1-\sigma^2)^{-1}d\sigma$$
Where we can now let $\epsilon=\sigma^2$, giving:
$$Q_n(a)=-\frac{\alpha^{-(2n+1)}}{2}\int_{\cos^2\beta}^{\cos^2\gamma}\epsilon^{\frac{2n-1}{2}}(1-\epsilon)^{-1}d\epsilon$$
Now we can use the standard function for the incomplete beta function which is:
$$B(x;a,b)=\int_0^xt^{a-1}(1-t)^{b-1}dt$$
And so we can now say that:
$$Q_n(a)=-\frac{\alpha^{-(2n+1)}}{2}\left[B\left(\cos^2\gamma;\frac{2n+1}{2},0\right)-B\left(\cos^2\beta;\frac{2n+1}{2},0\right)\right]$$
We can now go back and calculate that:
$$\cos^2\gamma=\frac{4}{(2+a)^2+4}$$
$$\cos^2\beta=\frac{4}{a^2+4}$$
$$-\frac{\alpha^{-(2n+1)}}{2}=-\frac{2^{4n+1}}{(4-a^2)^{-(2n+1)}}$$
Overall if we combine this together we get:
$$Q_n(a)=-\frac{2^{4n+1}}{(4-a^2)^{-(2n+1)}}\left[B\left(\frac{4}{(2+a)^2+4};\frac{2n+1}{2},0\right)-B\left(\frac{4}{a^2+4};\frac{2n+1}{2},0\right)\right]$$
I believe this is right but correct me if there are any mistakes
</s>