I'm trying to find the inverse of a particular cubic equation over a certain range. I have a partial solution and I'm not sure what's missing to make it complete. I'd like to solve the following equation in terms of $u$:
$$ c=ku^3+(1-k)u $$
Where:
$$ -0.5 <= k <= 0.5 \\ 0 <= u <= 1 \\ \therefore 0 <= c <= 1 $$
I know that for these value ranges, the function is strictly increasing and so has a 1:1 mapping between $c$ and $u$.
The reason I'm doing this is the above equation is involved in fisheye lens distortion correction and I'm trying to do a de-correction. I thought inverting the equation would be easy algebra and it turns out it's way too much for my 41-year-old brain to handle.
Ultimately, I ended up throwing it at Wolfram Alpha, which produced (I've introduced $p$ to save typing):
$$ p = \sqrt[3]{\sqrt{729 k^4 c^2 - 108 (k-1)^3 k^3} - 27 k^2 c} \\ u = - \frac{\sqrt[3]{2} (k-1)}{p} - \frac{p}{3 \sqrt[3]{2} k} $$
Now, despite being a small monster, it mostly works (eq 1 is the original, eq 3 is the inverse over the original $c$ and so should be identity, eq 4+ should be ignored) but only for $k > 0$. In particular, when $k <= 0$ then the original function is no longer generally 1:1 over $u$, and the inverse seems to be undefined over the areas with multiple solutions for $u$.
| k=-0.5 | k=0.5 |
|---|---|
 |
 |
Blue: original, purple: inverse. Note the undefined gap below around $c=u=2$.
Only 0 to 1 are of interest; zoomed out to show gap.
Because the inverse is undefined over $c$'s with multiple $u$'s, I'm almost certain that I need a $\pm$ (or two) somewhere, but I can't figure it out. I did try randomly changing $+$'s to $-$'s but never arrived at a correct solution. I'm at least pretty sure that the undefinedness originates with a negative term in that square root (e.g. when $s$ is negative here).
What do I need to do to that partial solution to complete it over the given range?
Previous work (can be skipped):
The very first thing I did was type the function into some random online algebra calculator and ask it to solve for $u$. None of them could solve the equation (I forgot about Wolfram Alpha at this point).
So I was going to do it by hand. At first I forgot that I was looking for "inverse" and started digging around for solutions to cubics. I discovered depressed cubics and set about finding an equation for the solutions (a doomed approach since I was really looking for the inverse). I found Cardano's method here and tested it out on Desmos (it's actually eq 6 in the link above, kind of a neat curve if you zoom out), where I finally realized this wasn't what I was looking for.
So then I searched around for inverses of depressed cubics (I didn't really know how to handle the $k$ in front of the 3rd-degree term but figured I'd work it out eventually). I found a lot of algebra work but I wasn't able to understand how to apply it with non-constant terms ($k$), or how to take advantage of the constraints that I had.
Eventually I found this comment on a post here and tried just entering the equation into Wolfram Alpha and, to my pleasant surprise, it worked.
I threw the equation back into Desmos to test it, which is where I noticed it wasn't defined over $0<=c<=1$ for $k<=0$, so I messed around with $+$'s and $-$'s to no avail, got pretty burnt out, and that's the point I'm at now.
I finally figured it out (after spending the better part of the last 24 hours on it)! The equation that Wolfram Alpha spit out is fundamentally correct, I was just missing a few things.
So the easy part first is when $k=0$ the denominators become $0$. This is trivially manageable by treating $k=0$ as a special case, where the resulting equation is simply identity ($u=c$) since $\lim_{k \to 0} u=c$.
Moving past that, one huge issue was that none of the software I was using to test (Desmos in particular, and also my own, which this is for) could handle complex numbers. While the inverse has real solutions over that undefined interval, calculating them involved complex intermediates, and this was making the function appear to be undefined over that range when in actuality it has real values. This greatly contributed to my confusion and, unbeknownst to me at the time, was preventing me from finding any kind of solution.
Now, the part that was missing from my partial solution was using all three cube roots instead of just an arbitrary one. I essentially had the right idea with assuming a $\pm$ was needed (before one of the square roots), except I missed the mark a bit: I was so focused on the sign change itself† that I neglected to realize that it was about finding all $n$ roots, and didn't think to do that with the cube roots.
My hypothesis, which ultimately proved correct, was that each of the three cube roots would correspond to one of the three overlapping curves in the cubic inverse. This was supported by Wolfram Alpha's graph (see below) showing what appeared to be one segment of that "missing" section.
The next difficulty was figuring out how to find all three cube roots. I searched around for a while, finding a lot of math that I didn't understand, until I stumbled on this page which described de'Moivre's completely mind-blowing formula (seriously; math is beautiful sometimes), which basically states that all the $n$-th roots of any number are $n$ evenly spaced points on a circle around the origin of the complex plane (What! Also makes sense thinking about square roots now, where a 180° rotation is just a sign swap.)
Anyways, head-exploding aside, to find the three roots of
$$ p = \sqrt[3]{\sqrt{729 k^4 c^2 - 108 (k-1)^3 k^3} - 27 k^2 c} $$
All I needed to do was:
pow(z, 1.0/3.0)wherezisstd::complex<>])So I did that, ran it through my test software, and it was a hit with correct results over the entire $-0.5<=k<=0.5$ range.
As a performance optimization in software, I implemented the rotation as a transformation matrix multiplication (this is math I'm familiar with); the math is rather clean. 120 degrees in one direction is the following transformation:
$$ r_1 = \begin{bmatrix}-\frac{1}{2} & \frac{\sqrt{3}}{2} \\ -\frac{\sqrt{3}}{2} & -\frac{1}{2}\end{bmatrix} r_0 $$
And the other direction is
$$ r_2 = \begin{bmatrix}-\frac{1}{2} & -\frac{\sqrt{3}}{2} \\ -\frac{\sqrt{3}}{2} & \frac{1}{2}\end{bmatrix} r_0 $$
Which yields, for a given root $r_0$, the two other roots $r_1$ and $r_2$ as follows (let $\Im$ yield a real number in this case):
$$ r_1 = -\frac{1}{2} \Re(r_0) + \frac{\sqrt{3}}{2} \Im(r_0) + (-\frac{\sqrt{3}}{2} \Re(r_0) - \frac{1}{2} \Im(r_0)) i\\ r_2 = -\frac{1}{2} \Re(r_0) - \frac{\sqrt{3}}{2} \Im(r_0) + (-\frac{\sqrt{3}}{2} \Re(r_0) + \frac{1}{2} \Im(r_0)) i$$
Also, to confirm my hypothesis, I graphed the real solutions for all three roots, and it made a beautifully perfect piecewise inverse curve:
From which I was picking a solution on the yellow line, which was the part missing from my partial solution.
Whew.
† To the point where I literally wrote a program that generated a second program which contained tests for all 1024 non-redundant sign changes of every term, and found no solutions.
Regarding the realization that complex numbers were causing computational difficulties:
I noticed this because, while all the software I was using had "undefined" regions for $k<0$, Wolfram Alpha was able to plot points in the whole range:
Then I saw an incidental remark in this post here, where the author wrote:
Putting those together it seemed clear that it was a computational issue, with Wolfram Alpha not surprisingly being able to handle it.
Confirming this, after adding complex number support (and the $k=0$ case) to my own test software, I was able to get real results over the entire $-0.5<=k<=0.5$ range.
As an aside, after reading the above post and working through the cubic formula, I was also able to see how the generated inverse made sense.