2022/11/02: following the correction of a bug in the rational implementation (cf this article), I've checked that the results introduced in the present article were not influenced by that bug and stay the same.
In a previous article I've explained why and how one would like to implement rational numbers. A set of unit tests helped me verify that the implementation I gave was valid, then I wanted a real use case to qualify it and prove there was really a possible advantage to use them instead of double.
Around the same time I happened to have looked at a video by Physics Simulation about the n-body problem (check also that video from PBS space time). "In physics, the n-body problem is the problem of predicting the motion of a group of celestial objects interacting with each other gravitationnally" (cf Wikipedia). This is exactly what I needed. The implementation is very easy and covers all the operations I have implemented. The problem can be set up easily to ensure a mix of arbitrarily small and large values during the calculation. The system is chaotic, hence numerical imprecision shows up clearly.
To solve the differential equations the method varies depending on how accurate you want to calculate the trajectory. The Euler method is the simplest one but not very accurate. The Runge-Kutta method is more accurate but more complex. Here I don't care about the accuracy of the calculated trajectory, rather I want to measure the accuracy of the calculation itself, whatever it is actually calculating. Euler's method is then good enough.
How is it possible to measure the accuracy of the calculation ? There exists no such a thing as a truth table to which I could compare the result of my calculation. And even if there was, there would still remain the question of how much of the measured inaccuracy comes from the measurement itself. The solution I used is to compare the calculation "to itself" in the most straightful way possible. As initial condition, I set 4 bodies with positions and velocities exactly symmetrical relatively to the origin (see videos below). Take aside the symmetry, the calculation for each body is exactly the same as the calculation for each other body. Hence, as long as the calculation is perfectly accurate, the positions of a pair of opposite bodies stays perfectly symmetric. Calculating the accuracy can then be reduced to a substraction of the coordinate along one arbitrarily choosen axis.
The problem of checking accuracy for one implementation is cleared. What about the comparison between the two implementations ? To be fair I need to ensure the initial conditions are the same in both implementation. Due to the chaotic nature of the n-body problem, the slightest difference in the initial conditions may produce misleading results. There is a way to ensure these initial conditions are perfectly equals: both double and CapyRatio implementation ensure that (within a range) integer values are exactly encoded. Choosing initial acceleration, velocity and position having integer values ensure then both simulations start in the exact same state and the comparison is fair.
For the implementation I'm using LibCapy of course, the code is available at the end of this article. It is well commented and not complicated so I skip the details here. I ran the calculation for both versions simultaneously, output the inaccuacy measurement and used LibreOffice Calc to plot graphs. I also recorded a video of the bodies trajectory because it's cool to look at and it helps understand what's happening.
The plot for the inaccuracy is as follow (click to enlarge):
and the video (only the end of the simulation):
The double implementation is in red, and the CapyRatio implementation is in blue. We see clearly the inaccuracy building up during the first 6 revolutions, until it reaches a point where the trajectory completely diverges. Both implementations follow the same pattern, suffer from an accuracy of the same order of magnitude, and diverge at around the same time. The only difference is at the beginning of the simulation: the inaccuracy for CapyRatio appears right from the beginning, while those of double stays equal to 0.0 for until around t=28.
Why this difference ? Well if you know about how double are implemented, or have read my article about random floating-point numbers, you already know that their imprecision grows proportionnaly with the value they represent. At the beginning of the simulation, the position is set to 100.0. In double implementation the value right below 100.0 is 99.99999999999998578915, a gap of the order of 1e-13, which matches the first recorded values of the recorded inaccuracy. CapyRatio has a constant precision, hence the apparition of the inaccuracy right from the beginning for that implementation.
On the previous graph, I measured the lost of symmetry based on the position, but I could have used the acceleration instead, which equally should stay exactly symmetric. At the start of the simulation the acceleration is null, so the inaccuracy should be measurable for both implementation right from the beginning. Indeed I obtained the following plot.
So, what I think is happening is as follow: at the beginning, the relatively large gap kind of dampens the inaccuracy of position while the acceleration and velocity have low values which also help to mitigate the inaccuracy (because the precision is better for these values). Then the simulation moves forward, the inaccuracy finally build up in the acceleration, then velocity, then position, until the inaccuracy in position (relative to its value) get large enough to be measured. Meanwhile CapyRatio measures the inaccuracy correctly right from the beginning for the position, meanwhile the inaccuracy for acceleration and velocity build up in the same manner.
Under the initial conditions of this first simulation both implementation ended up with similar results. Somehow it is disappointing as I was expecting to show that CapyRatio is better suited for that kind of simulation. However, the previous paragraphs remind me that the point is, the strength of rationals compare to double shows up the better the wider the range of values involved. In that first simulation I used initial values in the order of 1e2. Scaling up everything by a factor of 1.000.000 and running the exact same simulation gives the following results:
The double implementation diverges in the same fashion, but quicker, as expected. The CapyRatio implementation stays put... better than expected ! Why is that ? First, we know that the inaccuracy is not proportional to the value. So, while the values manipulated are much biggers, the inaccuracy stays the same, hence the relative inaccuracy is much smaller, helping to stabilize the simulation. Still, it doesn't explain why the inaccuracies, even if relatively smaller, do not accumulate over time.
I found the answer in the particular shape of the plot for rationals. The inaccuracy is measured in absolute value, those troughs are actually when the plot jumps from positive to negative (or the opposite). Indeed, the inaccuracy happens as well in both direction, sometime underestimating, sometime overestimating. So, my guess is that there is also some kind of self cancellation that's occuring here. As long as the relative inaccuracy stays low enough, the errors cancel out and the simulation stays stable. At some point the circumstances lead to inaccuracy leaning more in one direction and from then the divergence starts. In that second simulation the rationals implementation seems to be forever stable, but we see also variations in the inaccuracy. Very probably, if I had left the simulation run forever, at some point it would diverge too.
In conclusion, what's the takeaway of these simulations. First, it definitely qualifies the implementation of rationals I gave in the previous article. Second, it clearly shows that, under circumstances, the usage of rationals rather than double can be justified and at least should be considered. But, don't forget that rationals are awfully slow !
The code of the simulation:
Compile with:
A quick side note. For this simulation I needed to calculate the square root of rationals, which I haven't yet done at the time of my previous article and added later. If you wonder how I did it, it's a dichotomic search using the Heron's method.