Intriguing Indefinite Integral: $\int ( \frac{x^2-3x+1/3 }{x^3-x+1})^2 \mathrm{d}x$ 
Evaluate
$$\int \left( \frac{x^2-3x+\frac{1}{3}}{x^3-x+1}\right)^2 \mathrm{d}x$$

I tried using partial fractions but the denominator doesn't factor out nicely. I also substituted $x=\dfrac{1}{t}$ to get
$$\frac{-1}{9} \int \left(\frac{t^2-9t+3}{t^3-t^2+1}\right)^2 \, \mathrm{d}t $$
But I don't know how to solve this either.
Please Help.
Thanks in advance.
 A: Partial fractions also works, this may be a bit more work than trying to find a clever ad hoc solution, but the advantage is that it is a straightforward method. 
If $\alpha$ is a root of the denominator, then putting $x = \alpha + t$ and performing a so-called Laurent expansion around $t=0$ (which is analogous to a Taylor expansion, except that negative power terms may appear in the expansion) yields:
$$\frac{x^2-3x+\frac{1}{3}}{x^3 - x + 1} = \frac{\alpha^2 - 3\alpha + \frac{1}{3}}{3\alpha^2 -1}\frac{1}{t} +\left[\frac{2\alpha-3}{3\alpha^2-1} -\frac{3\alpha\left(\alpha^2-3\alpha+\frac{1}{3}\right)}{\left(3\alpha^2-1\right)^2}\right]+\mathcal{O}(t)$$
Squaring both sides yields:
$$\left(\frac{x^2-3x+\frac{1}{3}}{x^3 - x + 1}\right)^2 = \left( \frac{\alpha^2 - 3\alpha + \frac{1}{3}}{3\alpha^2 -1}\right)^2\frac{1}{t^2}+\mathcal{O}(1)$$
The $\frac{1}{t}$ term cancels out due to $\alpha$ satisfying the equation $\alpha^3 - \alpha +1 = 0$.
We thus have the partial fraction expansion:
$$\left(\frac{x^2-3x+\frac{1}{3}}{x^3 - x + 1}\right)^2 = \sum_{\alpha}\left( \frac{\alpha^2 - 3\alpha + \frac{1}{3}}{3\alpha^2 -1}\right)^2\frac{1}{(x-\alpha)^2}$$
(Proof: the difference between the l.h.s. and the r.h.s. is a rational function but which by construction has no singularities, so it must be a polynomial. But because the terms all tend to zero for $x\to\infty$, that polynomial must be identical to zero.)
This is obviously trivial to integrate, but how do we sum over the roots $\alpha$? Consider the function:
$$f(\alpha) = \frac{1}{\alpha^3 -\alpha +1}\frac{\left(\alpha^2 - 3\alpha + \frac{1}{3}\right)^2}{3\alpha^2 -1}\frac{1}{\alpha-x}$$
obtained by multiplying the summand of the answer by the logarithmic derivative of the denominator inside the square of the integrand. Then consider doing a partial fraction expansion using the Laurent expansion method as we did above. Since the function tends to zero for $\alpha \to \infty$, faster than $1/\alpha$, the sum of the coefficients of the terms with exponent -1 will sum to zero. By construction, the desired summation is part of this summation.
The answer is thus given by minus the sum of the coefficients the Laurent expansion terms with exponent -1 of $f(\alpha)$ around $\alpha = x$ and $\alpha = \pm \frac{1}{\sqrt{3}}$. Such coefficients are referred to as "residues" in the theory of complex functions.
The residue at $\alpha = x$ is obviously:
$$\frac{1}{x^3 -x +1}\frac{\left(x^2 - 3x + \frac{1}{3}\right)^2}{3x^2 -1}$$
The residue at $\alpha = \pm \frac{1}{\sqrt{3}}$ can be written as:
$$\frac{x+\frac{3}{2}}{1-3 x^2} \mp\frac{\frac{\sqrt{3}}{3}+\frac{3}{2}\sqrt{3}x}{1-3 x^2}$$
The integral, which is minus the sum of these residues, can thus be written as:
$$\frac{x^4+9x^3-\frac{35}{3}x^2+x+\frac{26}{9}}{\left(x^3-x+1\right)\left(3 x^2-1\right)}$$
A: I did not completely finished it but I think it is doable.
Consider differentiating $(x^3-x+1)^{-1}$, the derivative is $-3x^2(x^3-x^2+1)^{-2}$. The denominator is right, so quotient rule will be natural to use.
Next we try to guess what the antiderivative of the whole expression is. Consider quotient rule, $(\frac{u}{v})'=\frac{u'v-v'u}{v^2}$. Therefore our guess is that the numerator is $ax^2+bx+c$, since if we use quotient rule formula, $u'$ has highest order of 1, $v'$ has highest order of 2.
Use the formula to find the numerator, compare with $(x^2-3x+\frac{1}{3})^2$, you will get simultaneous equations w.r.t. $a,b,c$. Solve it, then you can get the final solution.
The answer can be simply gotten from Wolframalpha, check if you did correctly.
A: Ostrogradsky-Hermite method. As I have written in more detail here, we can use the Ostrogradski-Hermite method to integrate a rational fraction $P(x)/Q(x)$ without decomposing it into partial fractions and without finding the multiple roots of the denominator.  You can also find a simple description of this method here, subsection 4.5.2, as well as some exercises (1301-1304). 
Assume that $\deg P<\deg Q$. Then, there exist polynomials $$P_{1},\quad Q_{1}=\gcd \left\{ Q,Q^{\prime }\right\} ,\quad P_{2}\quad \text{and}\quad Q_{2}=Q/Q_{1},$$ with $\deg P_{1}<\deg Q_{1}$, $\deg P_{2}<\deg Q_{2}$, such that
\begin{equation*}
\int \frac{P}{Q}\, dx=\frac{P_{1}}{Q_{1}}+\int \frac{P_{2}}{
Q_{2}}\, dx.\tag{1}
\end{equation*}
Since the integrand is 
\begin{equation*}
\frac{P}{Q}=\frac{\left( x^{2}-3x+1/3\right) ^{2}}{\left(
x^{3}-x+1\right) ^{2}},\tag{2}
\end{equation*}
and
\begin{eqnarray*}
Q &=&\left( x^{3}-x+1\right) ^{2} \\
Q^{\prime } &=&2\left( x^{3}-x+1\right) \left( 3x^{2}-1\right)  \\
Q_{1} &=&\gcd \left\{ Q,Q^{\prime }\right\} =x^{3}-x+1,\\
Q_{2}&=&\frac{Q}{Q_{1}}=x^{3}-x+1,
\end{eqnarray*}
we can expand the given integral as
\begin{equation*}
\int \frac{\left( x^{2}-3x+1/3\right) ^{2}}{\left( x^{3}-x+1\right) ^{2}}\, dx=
\frac{Ax^{2}+Bx+C}{x^{3}-x+1}+\int \frac{Dx^{2}+Ex+F}{x^{3}-x+1}\, dx,\tag{3}
\end{equation*}
where $A$, $B$, $C$, $D$, $E$, $F$ are constants. They can be found by
differentiating $(3)$, reducing both sides to a common denominator and equating the coefficients of like powers of $x$ of the numerators. We thus have
\begin{eqnarray*}
\frac{\left( x^{2}-3x+1/3\right) ^{2}}{\left( x^{3}-x+1\right) ^{2}} &=&
\frac{\left( 2Ax+B\right) \left( x^{3}-x+1\right) -\left( 3x^{2}-1\right)
\left( Ax^{2}+Bx+C\right) }{\left( x^{3}-x+1\right) ^{2}}\\ &&+\frac{\left(Dx^{2}+Ex+F\right)\left(x^{3}-x+1\right)}{\left( x^{3}-x+1\right) ^{2}}.
\end{eqnarray*} 
Since  $\left( x^{2}-3x+1/3\right) ^{2}=x^{4}-6x^{3}+(29/3)x^{2}-2x+1/9$, after some algebra we obtain
\begin{align}
x^{4}-6x^{3}+(29/3)x^{2}-2x+1/9&= Dx^{5}+\left( -A+E\right) x^{4}+\left( -2B-D+F\right) x^{3}\\
&\qquad+\left(-A-3C+D-E\right) x^{2}+\left( 2A+E-F\right) x\\
&\qquad+\left( B+C+F\right),
\end{align}
which means that the constants satisfy the following system of equations
\begin{equation*}
\left\{
\begin{array}{l}
D =0,\\ -A+E=1,\\ -2B-D+F=-6, \\
-A-3C+D-E =\frac{29}{3},\\ 2A+E-F=-2, \\ B+C+F=\frac{1}{9},
\end{array}
\right.
\end{equation*}
whose solution is
\begin{equation*}
A=-1,\, B=3,\, C=-\frac{26}{9},\, D=0,\, E=0,\, F=0.
\end{equation*}
Hence
\begin{equation*}
\int \frac{\left( x^{2}-3x+1/3\right) ^{2}}{\left( x^{3}-x+1\right) ^{2}}\, dx=
\frac{-x^{2}+3x- 26/9 }{x^{3}-x+1}+C.\tag{4}
\end{equation*}
A: The goal is to represent
$$\dfrac{dt}{dx} = \dfrac19 \cdot \dfrac{(3x^2-9x+1)^2}{(x^3-x+1)^2} = \dfrac{d(u/v)}{dx} = \dfrac{vdu/dx - u dv/dx}{v^2}$$
Hence, it is tempting to choose $v=(x^3-x+1)$. This means we need $u$ such that
$$(x^3-x+1)u'(x) - (3x^2-1)u = (x^2-3x+1/3)^2$$
Hence, $u(x)$ must be a quadratic as $u(x) = ax^2+bx+c$. This gives us
$$(x^3-x+1)(2ax+b) - (3x^2-1)(ax^2+bx+c) = (x^2-3x+1/3)^2$$
Comparing coefficient of $x^4$, we obtain
$$2a - 3a = 1 \implies a = -1$$
Comparing coefficient of $x^2$, we obtain
$$-2a + a -3c = 9 + 2/3 \implies 3c = -a-29/3= 1-29/3 = -26/3 \implies c = -26/9$$
Comparing the constant term, we obtain
$$b +c = 1/9 \implies b = 3$$
Hence, we have
$$\dfrac{d}{dx}\left(\dfrac{-x^2+3x-26/9}{x^3-x+1}\right) = \dfrac19 \cdot \dfrac{(3x^2-9x+1)^2}{(x^3-x+1)^2}$$
