# pow(*,1./7.)

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?

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.