Moment of Inertia by Random Iterations

I am in the middle of writing a program (in C) that will calculate the moment of inertia for a 3D rigid-body ("object"). This object is non-symmetrical and not of any standard shape (i.e. sphere, cube, etc.). I have the exact x,y,z-coordinates of every point in or on this object (actually, it is a protein structure in PDB format. So, I have the x,y,z-coordinates of every atom in the structure). The program (or algorithm) assumes that the object's center of mass (CM) is already at the origin.

There are various approaches and algorithms out there. However, it is important, for me, the be able to code my own version. I have to calculate the moment of inertia for several thousand structures.

I have settled on finding the eigenvalues of a 3x3 matrix (the inertia tensor matrix). It is quite simple for me to sit down with pencil-and-paper and calculate these values (or use my TI-86 or Mathematica, etc.), however, I need to teach the computer to do this for me. It appears that this is not as simple as I had originally thought.

My approach is use random iterations that will converge upon the approximate eigenvalues for a given 3x3 matrix.

Below are the steps I have taken so far:

First, assume the following 3x3 matrix:
| I11 I12 I13 |
| I21 I22 I23 |
| I31 I32 I33 |

1. For each atom/point in the structure/object, calculate and sum the products, etc. of the matrix (i.e. I11 += y^2 + z^2, etc.).
2. Calculate the determinate, cofactor, and trace for the resultant inertia tensor matrix (using the sum of the values from step #1).
3. Iterate randomly through, say one million calculations (any large value should work here), converging upon an approximate eigenvalue whose difference is less than some pre-defined value. Use the following equation during the iterations:

x = [ ( (2*r^3) - (t*r^2) + d ) / ( (3*r^2) - (2*t*r) + c ) ],

where r = random value at given iteration,
t = trace, d = determinate, and c = cofactor

The value of x will be calculated again and again until it converges upon an eigenvalue (i.e. its difference is less than a given value).

Okay. So, I have coded the process above and it works. It is extremely fast and its approximate eigenvalues (the three of them) are extremely close to the actual values (as calculated by Mathematica).

I know why each step is done and I understand the overall process (I should. I coded it.). Here is what I don't understand, and why I am posting here: Can anyone tell me what that equation above is for or does? I found that equation on a website (can't remember which). The website presented it as "the one to use" but gave no reason why. A simple plot of the iterative values generated by this equation (during a sample run of the algorithm) produces a standard convergence x,y-plot. The equation appears to be an integral of some sort (maybe a harmonic series?). I just can't seem to figure out how it was integrated or why.

Any help with the above will be much appreciated. Also, I apologise for the mess that my post is. I hope the above makes sense.

Thank you!