I think I get it
Jack's spent over three decades working on his root finder. At long last, he finally understands why it ever worked in the first place.
A long time ago, in a land far, far away (Philadelphia), our company had just bought a program to compute interplanetary trajectories. I was told that it was based upon Lambert's theorem.
That bothered me a lot. I thought I was something of a trajectory weenie, but I had never even heard of Lambert's theorem. So, in my inimitable fashion, I set out to learn everything there was to know about it.
Weeks later, my desktop was littered by reference books, notes, and scribblings. I had filled a threering binder with analyses, derivations, and proofs of Lambert's theorem, complete with carefully drawn figures in four colors and several shades of graynot easy when you're doing them in ink.
My boss began to worry about all the money I was burning. One day he walked into my office and asked, “Um . . . don't you think you've already learned enough about Lambert's theorem?”
“Are you kidding me?” I asked. “I've only begun to scratch the surface.”
As it turned out, the boss was quite right to worry. Today, 42 years later, I still have the notes from that study, but I'm still not finished analyzing Lambert's theorem. There's still lots more to learn.
I've been feeling the same way about the math behind The World's Best Root Finder (TWBRF). As I mentioned in my June column (June 2003, p.11), this algorithm bubbled back to the surface again, when I set out to apply it to a current problem. I suspected it wasn't working properly. Also, reader John Lull sent me a function that TWBRF failed to solve. That made me worry that I'd missed something. Perhaps I'd made a mistake translating the code from Fortran to Pascal to C to C++. Perhaps I'd missed something in the analysis. Perhaps IBM had missed something in their original implementation of the algorithm, a Fortran function called RTMI.
In June, I reported that I'd completed the analysis, satisfied myself that it was okay, and studied John's problem to see why the algorithm failed to catch it. I was going to suggest some changes to at least make it catch John's test case.
Still I worried. I do that, sometimes. So I went back and studied my original analysis34 years old, and almost as detailed as the one on Lambert's theoremmulticolor graphs and all. Two things impressed me immediately. First, I was a pretty clever young squirt. The analysis was well done and extremely thorough.
Second, it was wrong.
Oh, don't worry. The algorithm's still okay, and still TWBRF, in my studiedand now quite informedopinion. But after all these years, I still didn't understand it to the deep, profound level that I like to understand things. I'd made a mistake trying to justify a decision IBM's developers made those many years ago, and it had propagated through everything done since. So I'm here today to set the record straight.
There's getting it, and getting it
At this point, I feel as though I've done enough research on the RTMI algorithm for at least a master's degree (honorary ones will be gracefully accepted). This whole exercise brings me right back to the notion of Getting It. I think it's important for people to Get It. I try to help them do so. I'm distressed when people (including me) don't get it, and I'm even more distressed when they don't seem to want to.
I'm sure you know the type. They don't want to know the details; they just want a canned solution. These are the folks like my classmates in engineering, who said, “I don't need to memorize the equation; I just need to know where to find it in the handbook.” No doubt the notion that they could actually derive the equation from first principles would have been completely foreign to them.
Recently I was shocked by a nasty review of my book, Math Toolkit for RealTime Programming. It turns out that the reviewer needed a specific algorithmthe log functionand turned to that section for it. He seems to have done something wrong in coding it, and it didn't work. So he blamed me.
The lesson here is an important one. He didn't get it. He didn't even want to get it. He just wanted a canned function that he could use without thinking about it.
Let's be clear: we all need canned results, sometimes. We don't always have the luxury of analyzing the problem in four colors and shades of gray. But that reader ignored my caveat in the book's preface, that the book was not intended for use as a handbook. Neither is this column. If that's what you need, there are lots of books and columns that might help. This isn't one of them. This column, like my book, is intended for those who want to learn, and learn in depth.
The RTMI algorithm is an excellent case in point. If you just want to use the algorithm, the code is there for the using at Embedded.com (ftp://ftp.embedded.com/pub/2002/05crens.txt). On the other hand, if you want to understand itI mean, really understand it, I'm your man. At this point, I've finally progressed past my level of understanding of 1968, and gone to a whole new level. The experience has been almost epiphanous. It also strikes me as a wonderful object lesson in the difference between merely getting it, and Getting It. In the rest of this column, I'm going to explain. I hope you enjoy the trip.
Background
Let y =ƒ(x ) be some function of x. We're looking for the value of x where the curve crosses the xaxis; that is, the “root” where:
ƒ(x )=0
(1)
Presumably, ƒ(x ) is a nasty function, or else we'd solve for the root analytically.
To narrow down the root, we start with three points equally spaced in x. Call them P_{0} , P_{1} , and P_{2} . In practice, we calculate x_{1} by bisecting the other two, and we compute the yvalues of each point by evaluating ƒ(x ) for that x. We require starting values that straddle the root. That means that the signs of y_{0} and y_{2} must be different, which in turn means that:
y_{0} y_{2} <>
(2)
We can go even further. We can require that the root is between x_{0} and x_{1} , meaning that:
y_{0} y_{1} <>
(3)
Of course, we have no idea, in advance, what sign y_{1} is going to have. However, if it turns out that Equation 3 isn't satisfied, we can always force it to be, simply by swapping the points P_{0} and P_{2} . This might make the point counting go from right to left, but who cares?
Our starting geometry, then, looks something like Figure 1.
Figure 1: The starting configuration
The RTMI/TWBRF algorithm fits a horizontallyopening parabola through the three points. The point where this parabola cuts the xaxis represents our best guess, based on the current data, for the root of the function.
We'll get to the math soon enough, but for now I want you to just grasp the concept with both hands: We're given two points straddling the root; we get a midpoint by bisection, we fit a parabola through the three points, and we takeor at least hope to takethe point where the parabola crosses the xaxis, as our root.
Figure 2: A poor fit
I say “hope to take” because the geometry may be ridiculous. We might have a situation like Figure 2, where the fitting parabola crosses the xaxis far away from the true root. As a matter of fact, Figure 2 is precisely the firststep fit for the function,
(4)
with an initial range, 0..6. In such a case, we'd be a lot better off simply bisecting again.
It's clear that we need some kind of criterion to decide whether or not we're wasting our time using the parabolic fit. The writers of RTMI knew this, of course, and included a test for reasonableness. I knew it too, these many years ago, and came up with my own criterion.
Our criteria were different. The reason why they were different has occupied entirely too much of my time since. Oh, it's not as though I stayed awake nights worrying about it. Well, maybe some nights, especially lately. But not many.
Looking at Figure 2, it stands to reason that we should not accept an estimated root that lies outside the original range. Heck, we don't even know if ƒ(x ) is defined out there. So a naive requirement (which is precisely the one I started with), would require that the root be between x_{0} and x_{2} . If we know that x_{2} lies to the right of x_{0} (remember, it may not), we could write two conditions,
x >x_{0} ; x x_{2}
(5)
A more general version, insensitive to the numbering rules, is:
(x –x_{0} ) (x –x_{2} )<>
(6)
There's no point allowing the case where the root, x, is equal to one of the bounding points. In fact, that path leads us to two successive, wrong estimates, which was the problem with John Lull's case.
More to the point, Equation 6 isn't restrictive enough. We've already said and enforced it by relabellingthat we're going to maintain the root between x_{0} and x_{1} . Therefore, at the very least, Equation 6 should read:
(x –x_{0} ) (x –x_{1} )<>
(7)
The deep math
To fit the parabola, we assume a secondorder function in a very special form:
x (y )=A +B (y –y_{0} )[1+C (y –y_{1} )]
(8)
This form is chosen only to make the solution easier.
Requiring that this equation fit at the three points gives us three equations:
x_{0} =A
x_{1} =A +B (y_{1} –y_{0} )
x_{2} =A +B (y_{2} –y_{0} )[1+C (y_{2} –y_{1} )])
(9)
Solving them for A, B, and C gives:
A =x_{0}
(10)
(11)and:
(12)
Remember, the root is the place where the fitted parabola has a y value of zero. Setting this to zero in Equation 8 gives:
x =x_{1} +By_{1} (1Cy_{2} )
(13)
Now that we have a formula for x, we can write down the conditions that satisfy Equation 7. After a little math, we find those to be:
(1Cy_{1} )(1Cy_{2} )>0
(14)
or, even more concisely:
uv >0
(15)
where:
u =y_{3} (y_{3} –y_{1} ) 2y_{2} (y_{2} –y_{1} )
(16)
and:
v =y_{3} (y_{3} –y_{2} )y_{1} (y_{2} –y_{1} )
(17)
(Don't be dismayed if you see that my previous definitions seem to have swapped u and v. It doesn't really matter, since we only care about the product. This notation is consistent with my old, original analysis.)
The flaw that Jack built
Now we get to the raw, exposed, evil heart of the analysis. All that's gone before tells us that to be certain that the estimated root is valid, we must evaluate u and v, and test their product. IBM didn't do that in RTMI. They tested only the sign of u.
When I was analyzing this algorithm, lo those 30odd years ago, I was at a loss to explain their reasoning. The RTMI algorithm worked just fine, but according to my own analysis, it shouldn't have. What followed in that analysis was a longwinded, armwaving rationalization why IBM was right to do what they did. What resulted was one of my fourcolor graphs, which I'm going to repeat for you here. The conclusion may have been flawed, but the graph is really neat.
We can divide the situation into four possible cases:
1. u <0,>v <>
2. u <0,>v >0 Forbidden
3. u >0, v <>
4. u >0, v >0 Okay
We can understand what's going on by examining the places where both u and v change signs. Of course, a function has to pass through zero to change signs, so the equations u =0 and v =0 define the curves that separate forbidden from acceptable regions.
Take another look at Equations 16 and 17, and you'll see that if we scale every value of y by some factor, we won't change the sign. We might change the value, but not the sign. Only the ratios of the y' s matters. Accordingly, we can simplify the equations a lot by introducing two new parameters:
(18)
Except for a factor y_{1} (which doesn't affect the signs), we can redefine u and v to be:
u =γ (γ 1)2β (β 1)
(19)
v =γ (γ –β )(β 1)
(20)
Setting these two equations to zero gives us the equations for the curves separating regions in β –γ space.
The canonical equation for conic sections Anytime you see an equation that's quadratic in two variables, be assured that it represents one of the conic sections. The standard, or canonical, form for an ellipse is:
where a and b are called the semimajor and semiminor axes of the ellipse. By convention, we assume that a is the larger distance, even though this may require us to relabel x and y. The corresponding equation for a hyperbola has one subtle difference: The sign of the second term:
We can think of the circle as a special case of the ellipse, for which a and b are equal. The one conic section that doesn't fit this picture is the parabola, which we can think of as the degenerative case of a hyperbola or ellipse (it's really just between the two). Its equation is enough different so that I won't bother showing it here. The equations for conic sections can be made to look strikingly different by rotating the reference axes, or adding offsets to the center. The most general form is: Ax^{2} +Bxy +Cy^{2} +Dx +Ey +F =0 Because of such transformations, it's not always obvious which figure a given equation represents. But be assured that any form of Equation i3 represents some conic section, and there are definite rules on how to transform such an equation into canonical form. In this column, I've used these rules to recognize the form of Equations 19 and 20.

Because both equations are quadratic in β and γ , they both describe curves that are conic sectionsin this case, hyperbolas (see the sidebar for more details). The reason you probably don't see this is that the hyperbolas are translated and, in one case, rotated from their “standard” forms. Applying the rules for how to get them back into canonical form, the transformed equation for u = 0 turns out to be:
(21)
This is the equation of a hyperbola opening along the baxis, centered at the point {1/2, 1/2}, and having semimajor and semiminor axes given by the square root of the inverse of the coefficients. That is to say:
(22)
The angle of the asymptotes is given by:
(23)or:
α=35.36^{º}
(24)
The two lobes of the hyperbola are plotted in Figure 3. Again, let me explain what you're looking at. The curves represent the places where u changes sign. It's not apparent from the figure, but u is positive in the middle region; negative above and below the hyperbolas.
Figure 3: Switching curves for u
We can repeat the process for v, but I hope by now you'll allow me to skip some of the details. We still get a hyperbola for the curve v = 0. The difference is that this one is centered at β = 2, γ = 1, and is rotated 22.5^{º} , so that one asymptote is vertical, the other inclined at 45^{º} . Figure 4 shows the curves. Again, v is positive in the region between curves; negative in the parts enclosed by curves.
Figure 4 Switching curves for v
Figure 5: Combined switching curves
Things begin to get really interesting when we draw both sets of curves together. That's shown in Figure 5. Remember, each curve represents a switching curve, where one of the functions goes through zero. Where two curves intersect, the sign of the function uv goes through four changes. Here, the acceptable areas are colored blue; forbidden red. You can see that things get really interesting in the area around β = 1, γ = 1/2. Here the lobes of two different hyperbolas overlap, just barely, leaving a very complicated crossover region.
Fortunately, however, that's all academic. Remember the definitions of β and γ , given in Equation 18. Recall also the inequality conditions, which require that both y_{1} and y_{2} lie on opposite sides of the xaxis, from y_{0} . That means that we are limited by our rules to only negative values of β and γ . That'll be the lower left quadrant of Figure 5, which is a whole lot more peaceful.
Even so, however, there's still that issue of how IBM managed to ignore one of the switching functions. Of our four possible cases, three are visible, even in the more limited region. They are:
1. u <0,>v <>
2. u <0,>v >0 Forbidden
4. u>0 , v >0 Okay
However, as you can see, the three cases include both positive and negative values for u and v . How did IBM, in their implementation, manage to test only one? Because they tested only u, it appeared that they were ignoring Case 1, which is represented by the (rather large) lobe inside the lower hyperbolic curve. Even though it seems to be a nice place to live, IBM's algorithm excludes it. They use only the triangular area above the top curve.
And this is where my arm waving began, so long ago. I won't go into the particular arm motions, since they were wrong anyhow. Suffice it to say that I got a proof that satisfied me at the time. But it was wrong.
When I reprised the analysis for Embedded Systems Programming in 1994, I had trouble accepting the armwaving. At least my baloney meter was running, albeit not on all cylinders. However, pushed for time, I managed to invent an even more frantically armwaving argument, which I foisted on the readers of that year.
That argument centered around the fact that C, the measure of curvature in the fitting parabola, would trend to zero as we homed in on a solution. A C trending to zero would make us wish to be near the γ axis. That claim turns out to be false also (though the wish is real).
So, what's really going on? Well, I finally got it. After all these years, I've got it.
The answer is simple. Look at the values the y' s must have to fall into the Case 1 region. A point firmly in the region is β =10, γ =5. Setting y_{0} =1 gives us the geometry shown in Figure 6.
Figure 6: A Case 1 configuration
Whoops! y_{1} is lower than y_{2} . That makes it seem likely that y_{0} is on one leg of the fitted parabola, while y_{1} and y_{2} are on another. As you can see, this causes the fitting algorithm to generate a curve that's a horrible fit to the actual data. Even though the estimated root that we get might be reasonably close to the true root, as we see in this figure, we get it more by good luck than anything else. And, to be honest, we'd probably be a lot better off simply to bisect the lefthand region again, and try for a better fit.
How can we avoid poor geometries like that shown in Figure 6? Clearly, it's by not allowing y_{1} to be further from the xaxis than y_{2} . Looking back at the definitions of β and γ , that means that we must require b to be the numerically smaller of the two. The switching curve is:
β = γ
(25)
which is shown as the dashed line in Figure 5. Anything below this line is verboten, which of course rules out all of that tempting Case 1 area. What remains is only the small, triangular Case 4 area near the baxis. This may seem small, but it's still where all the good configurations live. Anywhere else, we're better off bisecting.
Bottom line: IBM's designers threw out Case 1 completely, because it produces bad results. Now, after all this time, I finally understand.
Getting better
So, where are we? Simply stated, we finally, fully understand why the algorithm works the way it does. There are no more mysteries lurking within Figure 4. Rest assured that the code posted at ftp://ftp.embedded.com/pub/2002/05crens.txt works nicely.
In June, though, I said I saw ways to improve the algorithm. Was I wrong there also? Nope, I wasn't. You see, even though excluding Case 1 leads to better fits, it doesn't preclude all the bad ones, entirely. I had seen a few in my rather exhaustive testing, before. Using a tool like Mathcad gives me a chance to graph every single curve fit and see how well it works. That's an advantage neither IBM nor I had 30 years ago. And, sure enough, I saw some fits that didn't look all that great.
I don't need to show you yet another figure; Figure 6 will suffice, even though we know it to be forbidden, now. The point of the figure, and the reason the fit is poor, is because the vertex of the parabola lies in the yaxis range of interest. Even though Case 1 is ruled out, I found other cases where this was also true. Because the estimated root met all the conditions, it was accepted. However, the fit was poor, and I could see the algorithm hunting a bit, before it got a parabola it really liked. In every case, the hunting occurred when the vertex was in the range of the fit.
To avoid this, I'm adding yet another condition, that the vertex of the fitted parabola be outside the range of y_{0} ..y_{2} . With that added condition, only parabolic estimates that are really good, not just in range, will be used. I've worked out the math for the new condition, but it's too long to include here.
I'm also adding a test to catch John Lull's semipathological case, and others like it. If the range of β and γ are ridiculously large or small, I force a bisection.
Finally, I've addressed the issue of the stop criterion. Some readersincluding Johnhave suggested that we should be testing the range of both the x – and y values. As it turns out, this isn't really necessary, because the fit always becomes closer and closer to a straight line, as we get near convergence. However, I have found it useful to change the value of the convergence criterion, as that slope changes during the iteration. That saves us from pathological functions with weird kinks in them, near the root.
I think I won't bore you with the grungy details next time. Now that the fog has cleared away, I'm confident enough in the algorithm to change it off line. I'll just present the final version of the algorithm by posting it on the web site. I'll briefly discuss the changes made next time, then move on back to ztransforms and the Rosetta Stone equation.
Jack Crenshaw is a senior software engineer at SpectrumAstro and the author of Math Toolkit for RealTime Programming, from CMP Books. He holds a PhD in physics from Auburn University. Email him at .
References
 Crenshaw, J.W., “Programmer's ToolboxRoots of a Function,” Embedded Systems Programming, August 1994, pp. 1323.
 Crenshaw, J.W., “Programmer's ToolboxNone Other than e,” Embedded Systems Programming, September 1994, pp. 1324.
 Crenshaw, J.W., “Programmer's ToolboxRoots by Iteration,” Embedded Systems Programming, October 1994, pp. 1325.