In which direction to accelerate to pass through the origin In which direction, $\hat{n}$ (unit vector), should an object accelerate at constant rate $a$, with initial position $\vec{p}$ and velocity $\vec{v}$, such that it passes through the origin? So solving for $t$ and $\hat{n}$ in the following equation,
$$
\vec{p}+\vec{v}\;t+\frac{a}{2}\,\hat{n}\;t^2=\vec{0}. \tag{1}
$$
Solving this for $\hat{n}$ yields,
$$
\hat{n} = \frac{-2}{a\;t^2} \left(\vec{p} + \vec{v}\; t\right). \tag{2}
$$
The next equation can be obtain by taking the dot product of equation $(2)$ with itself, and solve it for $t$ to insure that $\hat{n}$ is indeed a unit vector,
$$
\frac{a^2}{4} t^4 = \vec{v}\cdot\vec{v}\;t^2 + 2\;\vec{v}\cdot\vec{p}\;t+\vec{p}\cdot\vec{p}. \tag{3}
$$
After some trial and error I noticed that the roots of this polynomial have the following structure,
$$
t=r-\delta \enspace\vee\enspace t=-r-\delta \enspace\vee\enspace t=\delta+i\,\sigma \enspace\vee\enspace t=\delta-i\,\sigma, \quad\forall\; r,\delta,\sigma\in\mathbb{R}.
$$
Using substitution I was able to find expressions for $r$ and $\sigma$,
$$
r^2 = -\delta^2 + \frac{2\;\vec{v}\cdot\vec{v}}{a^2} - \frac{2\;\vec{v}\cdot\vec{p}}{a^2} \delta^{-1}, \tag{4}
$$
$$
\sigma^2 = \delta^2 - \frac{2\;\vec{v}\cdot\vec{v}}{a^2} - \frac{2\;\vec{v}\cdot\vec{p}}{a^2} \delta^{-1}. \tag{4}
$$
However I am now stuck solving for $\delta$, which after the substitution of $x=\delta^2$ needs to satisfy,
$$
x^3 - \frac{2\;\vec{v}\cdot\vec{v}}{a^2} x^2 + \left(\frac{\left(\vec{v}\cdot\vec{v}\right)^2}{a^4} + \frac{\vec{p}\cdot\vec{p}}{a^2}\right) x - \frac{\left(\vec{v}\cdot\vec{p}\right)^2}{a^4} = 0. \tag{5}
$$
The general analytical solution for a cubic equation is known, however I wondered if in my case, that the even power coefficient are always negative and the odd power coefficient are always negative. According to Descartes' rule of signs there should be no negative real roots, but potentially two complex roots. Would there in this case exist a simpler solution for real (and positive) $x$?

How much more difficult would the problem become if an additional acceleration/disturbing-force $\vec{f}$ is added? The disturbance alters equations $(1)$, $(2)$ and $(3)$ into,
$$
\vec{p}+\vec{v}\;t+\frac{1}{2}\left(a\hat{n} + \vec{f}\right)t^2=\vec{0}, \tag{6}
$$
$$
\hat{n} = \frac{-2}{a\;t^2} \left(\vec{p} + \vec{v}\; t\right) - \frac{\vec{f}}{a}, \tag{7}
$$
$$
\frac{a^2 - \vec{f}\cdot\vec{f}}{4} t^4 = \vec{v}\cdot\vec{f}\;t^3 + \left(\vec{v}\cdot\vec{v} + \vec{p}\cdot\vec{f}\right)t^2 + 2\;\vec{v}\cdot\vec{p}\;t+\vec{p}\cdot\vec{p}. \tag{8}
$$
Due to the introduction of $\vec{f}$ the coefficient of the cubic, quadratic and linear terms of equation $(8)$ can have any sign. Under the assumption that $a^2>\vec{f}\cdot\vec{f}$, then the coefficient of the quartic term is always positive and the coefficient of the constant term is always negative. This implies that if the polynomial would be split into two second order polynomials, then the constant coefficients of each polynomial would have opposite signs. This implies at least two, but potentially four real roots. I am however only interested in the smallest positive real root. Does that simplify the disturbed problem?
 A: After some more searching online I did came across multiple methods, but to avoid having to compute complex values I went for Descartes' solution.
Both equation $(3)$ and $(8)$ can be written as a depressed quartic. Equation $(3)$ is already close to this form and only has to be normalized with respect to the forth order term. Equation $(8)$ also requires normalization, but also requires the following substitution, 
$$
t = x - \frac{\vec{v}\cdot\vec{f}}{\vec{f}\cdot\vec{f} - a^2}. \tag{1}
$$
So now the solutions to a equation of the following form need to be found,
$$
x^4 + A\, x^2 + B\, x + C = 0. \tag{2}
$$
For Descartes' solution the LHS of equation $(2)$ is split up into two quadratic terms,
$$
(x^2\! + \alpha x + \beta) (x^2\! + \gamma x + \delta) = x^4\! + (\alpha + \gamma)x^3\! + (\alpha\gamma + \beta + \delta)x^2\! + (\alpha\delta + \gamma\beta)x + \beta\delta. \tag{3}
$$
Since the coefficients should be the same as the LHS of equation $(2)$, so therefore the following constraints can be formulated,
$$
\begin{align}
\alpha + \gamma & = 0, \tag{4a} \\
\alpha\gamma + \beta + \delta & = A, \tag{4b} \\
\alpha\delta + \gamma\beta & = B, \tag{4c} \\
\beta\delta & = C. \tag{4d}
\end{align}
$$
Equation $(4a)$ can easily be solved for $\gamma$ as $\gamma = -\alpha$, which leaves,
$$
\begin{align}
\beta + \delta & = A + \alpha^2, \tag{5a} \\
\alpha(\delta - \beta) & = B, \tag{5b} \\
\beta\delta & = C. \tag{5c}
\end{align}
$$
By squaring both sides of equation $(5a)$ and $(5b)$, also multiplying the squared equation $(5a)$ by $\alpha^2$, and subtracting the resulting equation $(5b)$ from the resulting equation $(5b)$ yields,
$$
\begin{align}
\alpha^2 \left((\beta + \delta)^2 - (\delta - \beta)^2\right) & = \alpha^2 (A + \alpha^2)^2 - B^2, \tag{6a} \\
\alpha^2 \left(4\beta\delta\right) & = \alpha^6 + 2 A \alpha^4 + \alpha^2 A^2 - B^2. \tag{6b}
\end{align}
$$
After substituting equation $(5c)$ into equation $(6b)$ and also making the following substitution $\alpha^2 = \lambda$ yields the following cubic equation,
$$
\lambda^3 + 2 A \lambda^2 + (A^2 - 4C) \lambda - B^2 = 0. \tag{7}
$$
For solving this cubic equation I used a method from Numerical Recipes in C: The art of scientific computing, page 184-185. For this I needed to calculate two intermediate values,
$$
\begin{align}
Q & = \left(\frac{A}{3}\right)^2 + \frac{4C}{3}, \tag{8a} \\
R & = -\left(\frac{A}{3}\right)^3 + \frac{4 AC}{3} - \frac{B^2}{2}. \tag{8b}
\end{align}
$$
One positive real root can then be found with,
$$
\lambda = \left\{\begin{align}
-2 \sqrt{Q} \cos\left(\frac13 \left(2\pi + \cos^{-1}\frac{R}{\sqrt{Q^3}}\right)\right) - \frac{2A}{3}, & \quad\text{if}\quad Q^3 > R^2 \\
- \frac{2A}{3}, & \quad\text{if}\quad Q=R=0 \\
\frac{-\text{sign}(R)\, Q}{\sqrt[3]{|R| + \sqrt{R^2 - Q^3}}} - \text{sign}(R) \sqrt[3]{|R| + \sqrt{R^2 - Q^3}} - \frac{2A}{3}, & \quad\text{else}
\end{align}\right. \tag{9}
$$
The values for $\alpha$ and $\gamma$ can then be found by taking the square root of $\lambda$, and choosing $\alpha$ to be the positive solution and $\gamma$ the negative solution. The values for $\beta$ and $\delta$ can then be found using equations $(5a)$ and $(5b)$,
$$
\begin{align}
\beta & = \frac{A + \lambda - \frac{B}{\alpha}}{2}, \tag{10a} \\
\delta & = \frac{A + \lambda + \frac{B}{\alpha}}{2}. \tag{10b}
\end{align}
$$
The potential real roots can then be found by checking if the discriminant of each quadratic term, from equation $(3)$, is positive and then find the smallest positive solution from,
$$
\begin{align}
t_1 & = \frac{\gamma - \sqrt{\lambda - 4 \beta}}{2} - \frac{\vec{v}\cdot\vec{f}}{\vec{f}\cdot\vec{f} - a^2}, & \quad\text{if}\quad \lambda > 4 \beta \tag{11a} \\
t_2 & = \frac{2\beta}{\gamma - \sqrt{\lambda - 4 \beta}} - \frac{\vec{v}\cdot\vec{f}}{\vec{f}\cdot\vec{f} - a^2}, & \quad\text{if}\quad \lambda > 4 \beta \tag{11b} \\
t_3 & = \frac{\alpha + \sqrt{\lambda - 4 \delta}}{2} - \frac{\vec{v}\cdot\vec{f}}{\vec{f}\cdot\vec{f} - a^2}, & \quad\text{if}\quad \lambda > 4 \delta \tag{11c} \\
t_4 & = \frac{2\delta}{\alpha + \sqrt{\lambda - 4 \delta}} - \frac{\vec{v}\cdot\vec{f}}{\vec{f}\cdot\vec{f} - a^2}, & \quad\text{if}\quad \lambda > 4 \delta \tag{11d}
\end{align}
$$
I also implemented this into MATLAB and compared it to the general roots function and found that this method is about six times faster.
