A couple weeks ago, Mike Acton, head of Core Technologies at Insomniac (where I’m now working, btw) posted this reaming of some standard C++ math code. I think it’s an excellent read if you’re interested in deeply optimizing your math code, in particular for consoles, and in particular for the PS3.
Of course, posting this at all is a bit uncomfortable from where I’m sitting. I think it’s safe to say that if you were to compare our library code to what he is, er, analyzing, you’d find some significant similarities. However, our libraries were mainly designed to teach mathematical concepts rather than optimization concepts. As such, you end up with classes that represent mathematical entities, e.g. Planes, Spheres, Lines, etc. and simple code that emphasizes clarity. That said — as Mike’s slides aptly demonstrate — if you’re trying to squeeze every cycle out of the machine there are other ways to organize your data and code, and in fairness we didn’t address that much at all. In retrospect this was an unfortunate oversight. I’d say more emphasis on optimization would be a good addition to the Third Edition, but I’m not sure what my family would think of another book…
As part of the book review (see previous article), they have images from the Indexed weblog.  If you like graphical math jokes, you might find them amusing.
For those who are not aware, there is a new chapter in the book about Random Numbers. Because it is only a 50 page chapter in a much larger work, I don’t go into as much detail as I would have liked about probability. So I recommend a couple of longer works dedicated to the subject. One is Grinstead and Snell’s Introduction to Probability, which I had the honor of being taught from (in a pre-release edition) by John Kemeny. The other is Larry Gonick’s Cartoon Guide to Statistics, which is a surprisingly informative little book. However, after seeing a review of The Drunkard’s Walk: How Randomness Rules Our Lives, by Leonard Mlodinow, I may have a third. I recommend at least reading the review, and possibly checking out the book.
And if you’re into podcasts (and who isn’t, really? But I jest), the BBC’s In Our Time podcast has a episode about Probability. I haven’t had a chance to listen to it yet, but if it’s like their other programmes, I’m sure it will be quite interesting and informative.
Now we come to the final part: merging Newton-Raphson and bisection together, and choosing the correct one depending on the situation.
(more…)
Due to preparing for a week of vacation, and then going on a week of vacation, and then recovering from a week of vacation, this blog has fallen from the wayside a bit. But I’m back, and I’ll try to wrap up this topic in a couple of postings.
When last we left our intrepid heroes, we were discussing Newton-Raphson root finding, how it’s used for reparameterizing a curve, and when it fails (namely when the speed on the curve approaches zero). The solution hinted at was the bisection method.
(more…)
I’m still optimizing memory and speed so nothing exciting to report under the lighting topic. Maybe later this week or next week.
So let’s look at a bit of code from the book, shall we?
(more…)
One of the questions asked during the tutorial was whether matrices or quaternions were more efficient for angular rigid body dynamics. My response was that a matrix would be more efficient, because a quaternion requires a conversion to a matrix to compute the transformed inertial tensor matrix. If we analyze the operation count, this certainly appears to be true:
(more…)
Having promised to post about center-of-mass calculations, and subsequently having put it off for a week, I’m going to punt.
First, you can read the seminal Brian Mirtich paper on the subject. Then Dave Eberly responds with a more efficient method. Both rely computing solid integrals across a polytope of constant density by treating them as surface integrals across the polytope’s faces. In the end they end up with total mass (assuming a density of 1), center of mass, and the inertial tensor matrix for the polytope.
If those haven’t scared you off, Jonathan Blow has a more intuitive approach, noting that the solid integral is basically a volume calculation. By breaking the object into tetrahedrons, we can compute the total solid integral as a weighted sum of tetrahedral solid integrals. The tetrahedrons can be computed by selecting a single point — the origin or the centroid will work — and using that point together with the points on each triangular face. As he notes, it’s likely that the Eberly and Mirtich approaches are more efficient, but not as easy to understand from a geometric standpoint. Blow doesn’t provide an implementation, but you can grab a similar one from Stan Melax here.
Finally, Erin Catto pointed out to me after the tutorial that you can compute the center of mass by computing the centroids of the tetrahedrons and then performing a weighted sum of all the centroids, where the weight is the volume of a particular centroid’s tetrahedron divided by the total volume of the polytope. Blow also covers this in his paper as well. To get an intuitive sense of this, I recommend messing around with some origin-centered triangles (e.g. (1,0),(-.5,-.5),(0,-.5),(.5,-.5)) and using areas instead of volumes.
William Brennan points out that on page 84, the equation
\mathbf{w} = \mathbf{u}(\mathbf{v}\otimes\mathbf{w}))
is not correct. There are two ways to rewrite this. In one the intended order is correct, but is missing the transpose operator to indicate that it’s a row vector:
\mathbf{w} = \mathbf{u}^T(\mathbf{v}\otimes\mathbf{w}))
Alternatively, you can generate the same result with a column vector by doing:
\mathbf{w} = (\mathbf{w}\otimes\mathbf{v})\mathbf{u})
In Chapter 3, we simplify the Rodrigues formula (page 123) and general planar reflections (page 128) by using the tensor product. Since the arguments of the tensor product are the same in these cases, the ordering doesn’t matter. However, the ordering of the arguments needs to be reversed in the shear matrices on page 132, so
and
Other tensor product arguments may need to be reversed elsewhere in the text, though I can’t find any at this time.
There is a bug in the code that gets the axis and angle information from a rotation matrix. Where the angle is 180 degrees (the last else statement) it computes the reciprocal value incorrectly. The original code is:
The fixed code should be:
This should be corrected in both IvMatrix33 and IvMatrix44.
This is both a correction, and an improvement to the creation of a quaternion from two vectors. Or as it is commonly known, Shortest RotationArc.
The current code doesn't handle the case where the two vectors are opposing; there's an if statement to detect it, and a comment, but no recovery code. The problem is that there are an infinite number of orthogonal rotation axes -- we just need to pick one. The solution is to check to see which axis the start vector (it could be the finish vector, it doesn't matter) is generally pointing along. If it's pointing generally along the 
-axis, then we take the cross product with the 
-axis to get our rotation axis. Otherwise we take the cross product with the -axis.
(more...)
Ah, our first errata. How bittersweet.
Eric Mumau sent me this one. On page 499, equation 10.11 is wrong. It should read:
\mathbf{p}+2(\mathbf{v} \bullet \mathbf{p})\mathbf{v} + 2w(\mathbf{v} \times \mathbf{p}))
Correspondingly, the following equation should read:
 = (2 \cos^2 (\theta/2)-1)\mathbf{p} + 2(\hat{\mathbf{r}}\sin(\theta/2)\bullet\mathbf{p})\hat{\mathbf{r}}\sin(\theta/2) + 2\cos(\theta/2)(\hat{\mathbf{r}}\sin(\theta/2)\times\mathbf{p}))
The rest of the page follows as before.
This affects the code which starts at the bottom of page 499. At the top of 500, the first three lines should be rewritten as:
C++:
-
-
float vMult = 2.0f*(x*vector.x + y*vector.y + z*vector.z);
-
float crossMult = 2.0f*w;
-
float pMult = crossMult*w - 1.0f;
-
Note that the code on the CD is not the same as that given in the text, and in fact will actually work. However, it is slightly less efficient than what is given here.