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.

Computing orbits (1)

I want to simulate the orbits of planets and moons in our solar-system. They will form a carousel of gravitational bodies that my little spaceships will lift-off from, orbit around, and slingshot past, on their imaginary journeys of adventure and derring-do.

The wikipedia article on orbit modeling is a nice introduction to some of the concepts involved. Newton’s laws form the basis of computing orbits – every body attracts every other body in proportion to the product of their masses and inverse to the square of the distance between them.

If you imagine ttmp1240_thumbwo bodies (say the Sun and Earth) orbiting in isolation, then an analytic solution (a directly computable equation) to the orbit is possible. It’s a conic section (circle, ellipse or parabola) called a Kepler orbit.

A Kepler orbit, being a nice geometric figure, can be described rather simply with a small set of parameters.

Of course, the Sun and the Earth aren’t the only two bodies involved. Though the Sun’s mass dwarfs all others, planets and moons also pull on each other causing deviations – called perturbations – from the Keplerian orbit, which are evident when you measure carefully enough. There are no analytical solutions to the general problem that considers more than two bodies. Orbit prediction (also called propagation) is done using a variety of methods with the complexity dependent on the level of accuracy needed.

The most direct way that comes to mind – which is what I have used in previous code – is formally called an N-body integrator. Here, if you have N bodies in the simulation you simply compute the force F between every pair of bodies, vector sum the resultant force for each body, find the accelerations, and then use numerical integration to update the velocities and positions. I’ll call this the N-body method. The basic principle is straightforward, but there are important details to do with the numerical integration itself.

NASA’s HORIZONS system (pdf introduction) is the goto place for ephemeris data (ephemeris means “diary” in Latin and Greek and is the historical name for tables of positions of astronomical objects).  The whole HORIZONS system has a quaint 90s scientific computer system feel to it. There are three access methods – via web, via email and via telnet. Yes, telnet. I haven’t heard that word in a decade. I was pretty surprised when I found it still existed on my mac.

The web interface is the nod to “modernity”. It allows you to experiment with queries and familiarize yourself with the system.  The web interface is hard to scale, and I started to look for a Python based programatic interface to the HORIZONS system. My first hit was callhorizons. Callhorizons seems to have been created with a lot of love but it returns data in a somewhat specialized format. Scanning through the code and poking a little around the batch web interface it turns out it’s quite simple to create a programatic interface to query the system. I’ve put some barebones Python code to query the HORIZONS system here.

So, another method to determine orbits of planets for the simulation would be to get ephemeris data (I had to get used to that term too) from NASA’s Horizons system and store it as a look-up table. In this case the planetary orbits are completely predetermined and given with a high degree of fidelity. I’ll call this the look-up table method.

After I started doing my research I came across the paper: Keplerian Elements for Approximate Positions of the Major Planets. Here the author takes high accuracy data over a certain time period (epoch) from NASA and fits Keplerian orbits to the data – Keplerian orbits are ellipses and are completely defined by a small number of parameters (called elements). One neat thing the author does is to include delta terms that change the Keplerian elements with time – our ellipses precess round slowly – and this more or less captures the perturbation of the planetary orbits due to the presence of other planets. I’ll call this the “keplerian fit” approach.

My inclination is to start with the N-body simulation approach which is what we’ll explore in the next post.