Δmitry Baranovskiy

Maybe I have a tendency to reinvent the wheel every now and then, but this time, I promise, it was actually useful. While working on the easing implementation for a project that I can’t disclose, I noticed that the most popular library for easing calculations uses Newton’s method. This got me wondering: is there a way to compute it algebraically?

That became my challenge for the weekend. As an exercise, I decided to see if I could derive a solution myself. Here, I’d like to share what I discovered.

As a reminder, a Bézier cubic curve is a parametric curve. For easing calculations, it is positioned such that its start point is at $(0, 0)$ and its end point is at $(1, 1)$. While it is straightforward to calculate any point on the curve as a function of $t$, for easing to work correctly, we need to express $y$ as a function of $x$. This requires finding a function that converts $x$ into $t$. Essentially, we need to reverse $x(t)$.

$$x(t)=(1-t)^3x_1+3(1-t)^2tx_2+3(1-t)t^2x3+t^3x_4$$

Since $x_1=0$ and $x_4=1$, we could simplify it.

$$x(t)=3(1-t)^2tx_2+3(1-t)t^2x3+t^3$$

And after expanding the brackets we can bring it to standard polynomial form:

$$\frac{α}{3}t^3+\frac{β}{2}t^2+γt-x=0,$$

where

$$\begin{align*} α&=3(3x_1-3x_2+1)\\ β&=6(x_2-2x_1)\\ γ&=3x_1 \end{align*}$$

From here we could apply Cardano’s formula (named after Girolamo Cardano 1501–1576). To improve readability, let's introduce some intermediate constants and functions.

$$\begin{align*} q_1&=\frac{-β^3}{8α^3}+\frac{3βγ}{4α^2}\\ w&=\left(\frac{γ}{α}-\frac{β^2}{4α^2}\right)^3\\ δ&=\frac{β}{2α}\\ q(x)&=q_1+\frac{3x}{2α}\\ s(x)&=q(x)^2+w \end{align*}$$

Now, the final (for now) formula for $t(x)$ should look like this:

$$t(x)=\sqrt[3]{q(x)+\sqrt{s(x)}}+\sqrt[3]{q(x)-\sqrt{s(x)}}-δ$$

Great! We plug it in, and it works… well, most of the time. However, there are cases when $s(x) < 0$, causing the calculation to fail or “cough” and halt. At this point, we could stop here, concede defeat, and revert to Newton’s method… or we could venture into the realm of complex numbers. :).

Complex Realm

Complex numbers are usually presented in the form $x + iy$, where $x$ is the real part and $y$ is the imaginary part. If you look closely, you'll notice that if we define a complex number as $z = q(x) + i s(x)$, our function can be transformed into $\sqrt[3]{z} + \sqrt[3]{\bar{z}} - \delta$. Here, $\bar{z}$ denotes the conjugate of $z$—a number that is symmetrical with respect to the $x$-axis.

Complex numbers can also be treated as vectors, where the real part $x$ is along the horizontal axis and the imaginary part $y$ is along the vertical axis. When taking the cube root of a complex number, it's analogous to scaling down its length and dividing its angle by three.

If we write the complex number in polar form as $(r, \theta)$, then taking the cube root is expressed as:

$$\sqrt[3]{(r, \theta)} = \left( \sqrt[3]{r}, \frac{\theta + 2\pi k}{3} \right),$$

where $k = 0, 1, 2$, represents the three possible cube roots due to the $2\pi k$ term (the full circle). Although there are three possible answers, let’s set that aside for now.

After taking the cube root, we can add the number and its conjugate, $z_1 + z_2$, to obtain a real number. This is because their imaginary parts cancel out, leaving only the real part.

Lets do the calculations:

$$\begin{align*} a&=q(x)\\ b&=\sqrt{-s(x)}\\ r&=\sqrt[6]{a^2+b^2}\\ θ&=\tan^{-1}\left(\frac{b}{a}\right)\\ t(x)&=2r\cos\left(\frac{θ}{3}\right) \end{align*}$$

This is not too bad, however because of the multiple answers to cubic root and qualities of angles the final version will look like a massive piecewise function. To make it look simplier, lets find the result angle value first:

$$φ(x)=\left\{\begin{array}{rl} \frac{3}{2α}<0 & \left\{\begin{array}{rl} x>0 & 2π-θ\\ & π-θ \end{array}\right.\\ δ<0 & \left\{\begin{array}{rl} x>0 & 2π+θ\\ & θ-3π \end{array}\right.\\ & \left\{\begin{array}{rl} x>0 & θ\\ & π+θ \end{array}\right. \end{array}\right.$$

Our final function will look like this:

$$t(x)=\left\{\begin{array}{rl} α=0 & x\\ s(x)>0 & \sqrt[3]{q(x)+\sqrt{s(x)}}+\sqrt[3]{q(x)-\sqrt{s(x)}}-δ\\ s(x)\leq0 & 2r\cos\left(\frac{φ(x)}{3}\right)-δ\\ \end{array}\right.$$

This might look overcomplicated, but as benchmark shows, it’s almost three times faster than Newton’s method for this case and more precise. The only thing left is to submit pull request to the library.