Physics engine: use double or single precision?

Question 1

This is all going to depend on your application. You pretty much already understand the tradeoffs between the two:

Single-precision

Less accurate
Faster computations even on todays hardware. Take up less memory and operations are faster. Get more out of cache optimizations, etc.

Double-precision

More accurate
Slower computations.

Typically in graphics applications the precision for floats is plenty given the number of pixels on the screen and scaling of the scene. In scientific settings or smaller scale simulation you may need the extra precision. It also may depend on your hardware. For instance, I coded a physically based simulation for rigid bodies on a netbook and switching to float gained on average 10-15 FPS which almost doubled the FPS at that point in my implementation.

My recommendation is that if this is an educational activity use floats and target the graphics application. If you find in your studies and timing and personal experience you need double-precision then head in that direction.

Question 2

Surely the general rule is correctness first and performance second? That means using doubles unless you can convince yourself that you'll get fidelity required using floats.

The thing to look at is the effective size of one bit the coordinate system relative to the smallest size you intend to model.

For example, if you use earth coordinates, 100 degrees works around to around 1E7 metres. An IEEE 754 float has only 23 bits of precision, so that gives a relative precision of only about 1E-7. Hence the coordinate is only accurate to around 1 meter. This may or may not be sufficient for the problem.

I have learnt from experience to always use doubles for the physics and physical modelling calculations, but concede that cannot be a universal requirement.

It does not of course follow that the rendering should be using double; you may well want that as a float.

Question 3

I used typedef in a common header and went with float as my default.

typedef real_t float;

I would not recommend using templates for this because it causes huge design problem as you try to use polymorphic/virtual function.

Why floats does work

The floats worked pretty fine for me for 3 reasons:

First, almost every physical simulation would involve adding some noise to the forces and torques to be realistic. This random noises are usually far larger in magnitude than precision of floats.

Second, having limited precision is actually beneficial on many instances. Consider that almost all of the classical mechanics for rigid body doesn't apply in real world because there is no such thing as perfect rigid body. So when you apply force to less than perfect rigid body you don't gets perfect acceleration to the 7th digit.

Third, many simulations are for short duration so the accumulated errors remain small enough. Using double precision doesn't change this automatically. Creating long running simulations that matches the real world is extremely difficult and would be very specialized project.

When floats don't work

Here are the situation where I had to consider using double.

Latitude and longitudes should be double. Floats simply doesn't have good enough resolution for most purposes for these quantities.
Computing integral of very small quantities over time. For example, Gaussian Markov process is good way to represent random walks in sensor bias. However the values will typically be very small and accumulates. Errors in calculation could be much bigger in floats than doubles.
Specialized simulations that goes beyond usual classical mechanics of linear and rotational motions of rigid body. For example, if you do things with protein molecules, crystal growth, micro-gravity physics etc then you probably want to use double.

When doubles don't work

There are actually times when higher precision in double hurts, although its rare. An example from What every computer scientists should know...: if you have some quantity that is converging to 1 over time. You take its log and do something if result is 0. When using double, you might never get to 1 because rounding might not happen but with floats it might.

Another example: You need to use special code to compare real values. These code often has default rounding to epsilon which for float is fairly reasonable 1E-6 but for double its 1E-15. If you are not careful, this can give lot of surprises.

Performance

Here's another surprise: On modern x86 hardware there is little difference between raw performance of float vs double. The memory alignment, caching etc almost overwhelmingly dominates more than floating point types. On my machine a simple summation test of 100M random numbers with floats took 22 sec and with double it takes 25 sec. So floats are 12% faster indeed but I still think its too low to abandon double just for performance. However if you use SSE instructions or GPUs or embedded/mobile hardware like Arduino then floats would be much more faster and that can most certainly be driving factor.

A Physics engine that does nothing but linear and rotational motions of rigid body can run at 2000Hz on today's desktop-grade hardware on single thread. You can trivially parallelize this to many cores. Lot of simple low end simulations require just 50Hz. At 100Hz things starts to get pretty smooth. If you have things like PID controllers, you might have to go up to 500Hz. But even at that worse-case rate, you can still simulate plenty of objects with good enough desktop.

In summary, don't let performance be your driving factor unless you actually measure it.

What to do

A rule of thumb is to use as much precision as you need to get your code work. For simple physics engine for rigid body, float are often good enough. However you want to be able to change your mind without revamping your code. So the best approach is to use typedef as mentioned at the start and make sure you have your code working for float as well as double. Then measure often and chose the type as your project evolves.

Another vital thing in your case: Keep physics engine religiously separated from rendering system. Output from physics engine could be either double or float and should be typecasted to whatever rendering system needs.

Question 4

Here's the short answer.

Q. Why does OpenGL use float rather than double? 
A. Because most of the time you don't need the precision and doubles 
are twice the size.

Another thing to consider is that you shouldn't use doubles everywhere, just as some things may take require using a double as opposed to a float. For example, if you are drawing a circle by drawing squares by looping through the angles, there can only be so many squares shown on the screen. They will overlap, and in this case, doubles would be pointless. However if you're doing arbitrary floating point arithmetic, you may need the extra precision if you're trying to accurately represent the Mandelbrot series (although that totally depends on your algorithm.)

Either way, in the end, you will need to usually cast back to float if you intend to use those values in drawing.

Question 5

Single prec operations are faster and the data uses less memory less network bandwidth. So you only use double if you gain something in exchange for slower ops and more mem and bandwidth required. There are certainly applications of rigid body physics where the extra precision would be worth it, such as in manipulating lat\lon where single precision only gives you meter accuracy but is this your case?

Since it's educational purpose, maybe you want to educate yourself in the use of high precision physics algorithms where the extra accuracy would matter but lots of rigid body phys involves processes that can only be approximately quantified such as friction between 2 solids, collision reaction after detection etc, that extra precision wont matter you just get more precise approximate behavior :)