how to handle efficiently different probability distribution functions?

Question

I have a particle system that as usually creates new particles, updates them and destroys...

In the emitter module there is for loop that resets particles:

foreach p in particles
    p.position = rand()
    p.velocity = rand()

usually when using C's rand() function we get uniform distibution, but what when I would like to use some other distribution (gaussian for instance) ?

how to change that code so that it will handle several (or at least two) different ways of generating new particles' parameters?

Of course you can create some object: like RandomGenerator, and use some virtual function calls and handle those different behaviours. But this piece of code should be very fast (when updating thousands of particles), so using virtual functions is not good I think.

or maybe I should not care and simply write:

foreach p in particles
    p.position = useGaussian ? gausRand() : UniRand()
    p.velocity = useGaussian ? gausRand() : UniRand()

we can narrow number of different distribution and use only two or three of them...

please notice that my example is very simple, but in real code you have several particle parametrs configurations.

I would like to get some general advice on that subject.

Answer 1

While @gavinb's answer is a perfectly valid way to go, I'd suggest to avoid reinventing the wheel and go with the standard facilities: if you have c++11 support, use std::normal_distribution and its relatives (see, e.g. C++ TR1: how to use the normal_distribution?). Otherwise, use boost::random library.

Since these are header-only (the boost version at least), there are no polymorphic calls involved, so you don't have to worry about them. Which, of course, does not obviate the utmost relevance of @Oli Charlesworth's advice.

EDIT: if the overhead due to polymorphic calls is non-negligible, you can always template your functions on an enumerated type of distributions and specialize them as necessary.

In a nutshell, it's as simple as this:

#include<iostream>

// template on an int selector
template<int N> void foo(){ std::cout<<"42\n"; }
template<> void foo<1>() {std::cout<<"1\n";}

//now use an enum
enum  distr_types {UNIF, NORMAL, UNKNOWN}; 
template<distr_types T> void bar() {std::cout<<"fourty two\n";}
template<> void bar<UNIF>() {std::cout<<"UNIF\n";}
template<> void bar<NORMAL>(){std::cout<<"NORMAL\n";}

int main(){
  foo<3>();
  foo<1>();

  bar<UNIF>();
  bar<NORMAL>();
  bar<UNKNOWN>();
}

But if you find yourself doing things of this sort, it's worth having a look at one of good C++ books.

Answer 2

Usually when using C's rand() function we get uniform distibution, but what when I would like to use some other distribution (gaussian for instance) ?

The Box Muller Transform is a very clever algorithm which uses trig functions to generate random numbers from a Gaussian distribution using a uniform distribution as input (ie. using rand()). You specify the mean and standard deviation, and invoke this function to generate new variates. The only drawback is that it is more expensive than simply calling rand, as it also calls sin() and cos() (though only every second invocation).

how to change that code so that it will handle several (or at least two) different ways of generating new particles' parameters?

I suggest you start with the RandomGenerator and virtual methods approach. It will be the easiest to maintain. Start with the simplest approach and profile before you try to optimise things.

Given the computational complexity of generating random numbers, the cost of producing a variate will far outweigh the overhead of a virtual method call versus a static function call.

If this really, really isn't fast enough you could always generate a pool of random numbers, generating more in the background as required.