Let $p(x)/q(x)$ be a real-valued rational function in $\mathbb{R}(x)$. I am wondering what the complete story is for integration of these functions. I would expect some algorithm to work like this:
1) Factor $q(x)$ into a product $\prod_{i = 1}^n q_i(x)^{n_i}$. This is actually impossible for polynomials $q(x)$ of high enough degree. So suppose we started with this factorization for $q(x)$ so we can expect to have a closed formula at the end.
2) Use the euclidian algorithm to express $p(x) / q(x) = \sum_{i =1}^n p_i(x) / q_i(x)^{n_i}$. Using linearity of the integral, it is enough to calculate the integral of $p_i(x)/q_i(x)^{n_i}$.
3) I would expect an algorithm to reduce to rational functions of the form $1/r_i(x)^{n_i}$ ($r_i$ of degree $2$ or $1$), and then to rational functions of the form $1/r_i(x)$ ($r_i$ of degree $2$ or $1$). I don't actually know how this would be done though.
4) For the integration of $1/r(x)$, where $r(x)$ is irreducible of degree $2$, we may use substitution, shifting, scaling, etc. to put $r(x)$ in the form $x^2 + 1$. This has the known integral $\text{arctan}(x) + C$. For the integration of $1/x$, that is simply $\text{ln}(x) + C$.
This is indeed the algorithm for integrating rational functions if you don't want to work over $\mathbb{C}$.
Step 1 can always be done, the catch is that you don't necessarily have an explicit formula for the roots (or the coefficients of $q_i$).
Step 2: You forgot to first compute the quotient polynomial to make $\deg p<\deg q$.
Step 3: If $q_i$ has degree $2$, you will have $$ \frac{p_i(x)}{q_i(x)^{n_i}}=\frac{a_1x+b_1}{q_i(x)}+\frac{a_2x+b_2}{q_i(x)^2}+\dots+\frac{a_{n_i}x+b_{n_i}}{q_i(x)^{n_i}} $$ (assuming $\deg p_i<2n_i$).
Step 4: You know an antiderivative of $1/x^n$. For $\dfrac{ax+b}{(x^2+1)^n}$, use the reduction formulae $$ \int\frac{\mathrm{d}x}{(x^2+1)^n}=\frac{2x+1}{4(n-1)(x^2+1)^{n-1}}+\frac{2n-3}{2(n-1)}\int\frac{\mathrm{d}x}{(x^2+1)^{n-1}} $$ and the obvious (dropping the $+C$) $$ \int\frac{x\,\mathrm{d}x}{(x^2+1)^n}= \begin{cases} \frac12\log(x^2+1) & n=1\\ -\dfrac{1}{2(n-1)(x^2+1)^{n-1}} & n>1. \end{cases} $$