$$\int \frac{5x^3+3x-1}{\left(x^3+3x+1\right)^3},dx$$
$x^3+3x+1$ Isn't going to factor cleanly, so we can use the cubic formula:
$$x = \sqrt[3]{\left(\frac{-b^3}{27a^3}+\frac{bc}{6a^2}-\frac{d}{2a}\right)+\sqrt{\left(\frac{-b^3}{27a^3}+\frac{bc}{6a^2}-\frac{d}{2a}\right)^2+\left(\frac{c}{3a}-\frac{b^2}{9a^2}\right)^3}}+\sqrt[3]{\left(\frac{-b^3}{27a^3}+\frac{bc}{6a^2}-\frac{d}{2a}\right)-\sqrt{\left(\frac{-b^3}{27a^3}+\frac{bc}{6a^2}-\frac{d}{2a}\right)^2+\left(\frac{c}{3a}-\frac{b^2}{9a^2}\right)^3}}-\frac{b}{3a}$$
$a=1,b=0,c=3,d=1$
$$\sqrt[3]{\left(-\frac{1}{2}\right)\pm\sqrt{\frac{1}{4}+1}}$$
$x = \pm\sqrt[3]{\frac{\sqrt{5}\pm1}{2}}$
This is the combination of third roots of the golden ratio $\varphi$:
$$x=-\sqrt[3]{\varphi}+\sqrt[3]{\varphi^{-1}}$$
$x_0 =-\sqrt[3]{\varphi}+\sqrt[3]{\varphi^{-1}}$
$$x^3+3x+1=\left(x-x_0\right)\left(x-x_1\right)\left(x-x_2\right)$$
$$x_{1,2} = \frac{-x_0\pm \sqrt{x_0^2+\frac{4}{x_0}}}{2}$$
$$x_1 = \frac{\sqrt[3]{\frac{\sqrt{5}+1}{2}}-\sqrt[3]{\frac{\sqrt{5}-1}{2}}}{2}+\frac{\sqrt{3}}{2}\left(\sqrt[3]{\frac{\sqrt{5}+1}{2}}+\sqrt[3]{\frac{\sqrt{5}-1}{2}}\right) i, \quad x_2 = \frac{\sqrt[3]{\frac{\sqrt{5}+1}{2}}-\sqrt[3]{\frac{\sqrt{5}-1}{2}}}{2}-\frac{\sqrt{3}}{2}\left(\sqrt[3]{\frac{\sqrt{5}+1}{2}}+\sqrt[3]{\frac{\sqrt{5}-1}{2}}\right) i$$
$$x_1 = \frac{\sqrt[3]{\varphi}-\sqrt[3]{\varphi^{-1}}}{2}+\frac{\sqrt{3}}{2}\left(\sqrt[3]{\varphi}+\sqrt[3]{\varphi^{-1}}\right) i, \quad x_2 = \frac{\sqrt[3]{\varphi}-\sqrt[3]{\varphi^{-1}}}{2}-\frac{\sqrt{3}}{2}\left(\sqrt[3]{\varphi}+\sqrt[3]{\varphi^{-1}}\right) i$$
$$x^3+3x+1=\prod_{k=0}^{2} \left(x-x_k\right)$$
$$\left(x-x_1\right)\left(x-x_2\right)= x^2 - \frac{2^{2/3}x\sqrt[3]{2\varphi}}{2}+\frac{2^{2/3}x \sqrt[3]{2\varphi^{-1}}}{2}+\frac{\sqrt[3]{2}\left(2\varphi^{-1}\right)^{2/3}}{2}+1+\frac{\sqrt[3]{2}\left(2\varphi\right)^{2/3}}{2}$$
So, we know that $\left(\left(x-x_0\right)\left(x-x_1\right)\left(x-x_2\right)\right)^3$ has three roots (one real, two complex) each with multiplicity three.
We'll apply the cubic formula to the numerator as well, which yields:
$$n_0 = \sqrt[3]{\frac{3\sqrt{5}+5}{50}}-\sqrt[3]{\frac{3\sqrt{5}-5}{50}}$$
$\Delta_{+} = \frac{3\sqrt{5}+5}{50}, \quad \Delta_{-} = \frac{3\sqrt{5}-5}{50}$
$$n_{1,2} = -\frac{n_0}{2} \pm \frac{\sqrt{3}}{2} \sqrt{\frac{3\left(1+n_0\right)}{5 n_0}} i$$
$$n_1 = -\frac{\sqrt[3]{\Delta_+}-\sqrt[3]{\Delta_-}}{2}+\frac{\sqrt{3}}{2} \sqrt{\frac{3\left(1+\sqrt[3]{\Delta_+}-\sqrt[3]{\Delta_-}\right)}{5\left(\sqrt[3]{\Delta_+}-\sqrt[3]{\Delta_-}\right)}} i, \quad n_2 = -\frac{\sqrt[3]{\Delta_+}-\sqrt[3]{\Delta_-}}{2}-\frac{\sqrt{3}}{2} \sqrt{\frac{3\left(1+\sqrt[3]{\Delta_+}-\sqrt[3]{\Delta_-}\right)}{5\left(\sqrt[3]{\Delta_+}-\sqrt[3]{\Delta_-}\right)}} i$$
$$\int\frac{5x^3+3x-1}{\left(x^3+3x+1\right)^3}\,dx = \int \frac{\left(x-n_0\right)\left(x-n_1\right)\left(x-n_2\right)}{\left[\left(x-x_0\right)\left(x-x_1\right)\left(x-x_2\right)\right]^3}\,dx\quad \left[\bigcup_{k\geqslant1} \lbrace x_k;n_k\rbrace \in \mathbb{C}\right] \wedge \left[\bigcap_{k\geqslant0} \lbrace x_k;n_k\rbrace =0\right]$$
Let $P\left(x\right)$ represent the numerator, and $Q\left(x\right)$ represent the denominator; then let $\mathbb{C}\left[x\right]$ be the polynomial ring over $\mathbb{C}$, with the ideal $I_k = \left<\left(x-x_k\right)^3\right>$; therefore $Q^3\left(x\right)$ is defined as the intersection of the ideals $I_k$ such that $\prod_{k\geqslant 0}\left(x-x_k\right)^3=\bigcap_{k\geqslant 0} I_k = I$. From here, we can utilize the Chinese remainder theorem, which states that every equivalence class (ie. two polynomials are congruent if their difference $\lambda \in I$) modulo $I$ decomposes (due to the inequivalence of roots) into unique components modulo $I_k$, such that:$$\mathbb{C}\left[x\right]/ I \cong \bigoplus_{k\geqslant0}\mathbb{C}\left[x\right]/I_k$$ Which allows for the decomposition of the polynomials via partial fractions.
From here, let us define a function $D_k\left(x\right)$ which outputs the product of all poles for which $x_k\neq x_j$ such that $D_k\left(x\right) = \prod_{\substack{j\geqslant 0 \\ j\neq k}}\left(x-x_j\right)^3$. For every pole $x_k$, partial fraction coefficients $A_k,B_k,C_k$ are derived from their congruence to the taylor expansion of $P\left(x\right)/D_k\left(x\right)$ centered at $x=x_k$, modulo $\left(x-x_k\right)^3$ up and equal to degree 2:
$$\frac{P\left(x\right)}{D_k\left(x\right)}\equiv C_k + B_k \left(x-x_k\right)+ A_k\left(x-x_k\right)^2 \pmod{\left[x-x_k\right]^3}$$
Now, we will develop indexing system, which we will define over a set of cyclic indices taken modulo three:
$$k+1 \to k+1 \pmod{3} \\\\k+2 \to k+2 \pmod{3}$$
Now, we can determine coefficients:
Residue $C_k$ is the constant term of the Laurent expansion $\left[x-x_k\right]^3 f\left(x\right) = C_K + B_k\left(x-x_k\right)+A_k\left(x-x_k\right)^2$, which is isolated via truncating modulo $\left[x-x_k\right]^3$.
$$\left[x-x_k\right]^3 f\left(x\right) \equiv C_k \pmod{\left[x-x_k\right]^3}$$
$$C_k = \lim_{x\to x_k} \left[x-x_k\right]f\left(x\right)$$
$$C_k = \frac{P\left(x_k\right)}{D_k\left(x_k\right)}$$
$$C_k = \frac{\left(x_k-n_0\right)\left(x_k-n_1\right)\left(x_k-n_2\right)}{\left(x_k-x_{k+1}\right)^3 \left(x_k - x_{k+2}\right)^3}$$
Coefficients residue $B_k$ and $A_k$ will be far easier to obtain:
$$\frac{d}{dx}\left[\left(x-x_k\right)^3 f\left(x\right)\right] \equiv B_k + 2 A_k\left(x-x_k\right) \pmod{\left[x-x_k\right]^3}$$
$$B_k = \lim_{x\to x_k} \frac{d}{dx}\left[\left(x-x_k\right)^3 f\left(x\right)\right]$$
$$B_k = \frac{P'\left(x_k\right)-3 P\left(x_k\right)\left(\frac{1}{x_k - x_{k+1}}+\frac{1}{x_k - x_{k+2}}\right)}{\left(x_k-x_{k+1}\right)^3\left(x_k-x_{k+2}\right)^3}$$
$$\frac{d^2}{dx^2}\left[\left(x-x_k\right)^3 f\left(x\right)\right] \equiv 2 A_k \pmod{\left[x-x_k\right]^3}$$
$$A_k = \lim_{x\to x_k} \frac{1}{2} \frac{d^2}{dx^2}\left[\left(x-x_k\right)^3 f\left(x\right)\right] $$
$$A_k = \frac{1}{2}\frac{P''\left(x_k\right)-6 P'\left(x_k\right) \left(\frac{1}{x_k - x_{k+1}}+\frac{1}{x_k - x_{k+2}}\right) + 9 P\left(x_k\right)\left(\frac{1}{\left(x_k - x_{k+1}\right)^2}+\frac{1}{\left(x_k - x_{k+2}\right)^2}\right)}{\left(x_k-x_{k+1}\right)^3\left(x_k-x_{k+2}\right)^3}$$
Now that we have the partial fraction coefficients, we can say:
$$\int \frac{P\left(x\right)}{Q^3\left(x\right)}\,dx= \sum_{k\geqslant 0}\int \frac{A_k}{x-x_k} + \frac{B_k}{\left(x-x_k\right)^2} + \frac{C_k}{\left(x-x_k\right)^3}\,dx$$
Which can easily be evaluated to:
$$\sum_{k\geqslant 0} A_k \log\left(x-x_k\right) - \frac{B_k}{x-x_k} - \frac{C_k}{2\left(x-x_k\right)^2}$$
$$\boxed{\int \frac{5x^3+3x-1}{\left(x^3+3x+1\right)^3}\,dx=\sum_{0 \leqslant k \leqslant 2} A_k \log\left(x-x_k\right) - \frac{B_k}{x-x_k} - \frac{C_k}{2\left(x-x_k\right)^2}}$$
