I was reading through some code (this file from the REBOUND N-body simulator) when I came across this function:

// Machine independent implementation of pow(*,1./7.)
static double sqrt7(double a){
    double x = 1.;
    for (int k=0; k<20;k++){  // A smaller number should be ok too.
        double x6 = x*x*x*x*x*x;
        x += (a/x6-x)/7.;
    return x;

and I got very curious as to how it worked, and why one would need it and mostly about whether it would generalize to arbitrary values of base and exponent.

We see that we have an initial guess (x) which we then increment using the curious factor \frac{1}{7}(\frac{a}{x^6} - x) = \frac{a - x^7}{7x^6}. We can see that when x approaches a^{\frac{1}{7}} this factor vanishes and x stops moving.

But why the other components?

Untitled drawing

If we create a function f = x^7 - a we see that it goes to zero when x = a^{\frac{1}{7}} (The point marked out in the sketch. Now, we know that the differential of the curve (f' = 7x^6 = \frac{\Delta f}{\Delta x}) gives us the slope of the curve at x. At x we are \Delta f = x^7 - a away from our desired value. So, by moving by \Delta x = \frac{\Delta f}{f'} we get closer to our desired point in a principled way – we move faster when the slope is greater, and slower when it’s smaller (and close to the desired location)

Does this look familiar? Of course it does! It’s gradient descent which they taught you when you were a baby, and now you’ve forgotten it, because who does things like that any more? Formally, it’s called the nth root algorithm and the wikipedia page suggests that you probably only need 6 or so iterations to get 64 bit accuracy.

I got curious about how the pow function was implemented in the C++ standard library. The internet tells me that – it depends. Thrashing about a bit in the cmath code for gcc we see different prototypes for different types of input base and exponent and it looks like the actual function may be implemented in Fortran.

Stackoverflow pointed me to an opensource Apple implementation. The vast bulk of the function is devoted to dealing with edge cases (It’s impressive to read, if a little messy) and there is an interesting algorithm (which uses a lookup table included in the source) which I couldn’t figure out.

On this topic, if your exponent is a positive integer, Exponentiation by squaring is a fun method to compute A^n.


Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google+ photo

You are commenting using your Google+ account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )


Connecting to %s