ATTENTION: This algorithm presented here is cool, and teaches a lot, but is not necessarily much faster than the alternative of using (float)rand()/RAND_MAX; Test it to make sure it can really help you in your code.

Generating random numbers to use in your program is one of those things that has many ways to be done. And usually you are in the middle of something when you suddenly have to use a random number generator, and it’s always a wall you hit in the middle of your work, because you suddenly find yourself thinking “but is this the best way? What should I do?” while you wanted to keep working on your stuff. Usually functions are also not perfectly fit to what you need, and that makes it worse to decide what code to use.

The best reference that I know to start learning about how pseudo-random number generators work is Knuth’s TAOCP. He has an entire chapter (3) dedicated to the subject, and it’s a great read. He explains how one of the most common generators work, the one based on modular arithmetic and that create a sequence of random integers inside a very large range, limited by the precision of your register.

In Python, Matlab and other high-level languages it’s usually easy to find good floating-point generators. But in C, unless you are using some fancy scientific library, you only have access to the standard rand() function, which returns a number between 0 and RAND_MAX. So when you need a floating-point number the first solution most people end up using is:

(float)rand()/RAND_MAX;

This converts two large integers to float, and then divides them to produce a number between 0.0 and 1.0, that you can now convert to, for example, -1.0 and 1.0 by multiplying by two and subtracting by one. And of course, before that you also need a srandom(time(NULL)); to set the seed of the generator with the current time to (hopefully) make the sequence change every time you run the program again.

That is not cool because it involves a division, and this is usually a very slow operation. The best thing would be to have a more direct a way to generate this random number, by messing around with the bits of the floating-point value.

A floating point value in the IEEE754 standard has a bit to represent the signal, some bits for the exponent, and the last bits are the mantissa. The exponent bits determine the range covered by the mantissa bits. For example, the number 1.0 is 0x3f800000. The first bit is 0, the followed by ‘01111111’, which means 0 because you subtract 127 from this number to find out the exponent. The rest are the “decimal places” after a supposed “1.”. The number 1.5 is given by 0x3fc00000, the first of the mantissa bits is turned on to represent the 0.5 added to the 1.0 base. if we increment the exponent to ‘10000000’, this base is now 2.0, so 0x40000000 represent 2.0, and 0x40004000 is 3.0. The same mantissa of 1.5, but multiplied by 2 because we incremented the exponent.

To generate the random floating-point numbers directly, my idea (that I’m pretty sure is not original at all) was to replace these mantissa bits with random values, and leave a proper content on the exponent and signal.

If we used the 1.0-2.0 range, we would then pick up the 0x3f800000 value, and replace the last 23 bits of the mantissa with a random number, produced by the integer pseudo-generator. This would directly translate a random number in the 0-8388607 interval (2^23-1) to a floating point between in the range [1.0,2.0). A subtraction wou make this a number between 0.0 and 1.0, and then we can perform the transformation to -1.0 to 1.0.

These linear transformations can be optimized if we instead generate the number already in a range with size 2.0. That means making the exponent the one for the 2.0 base. So what we do is pick up the number 0x40000000, and replace again the last 23 bits with random bits. After that we only neet to subtract 3.0 from the number to generate something between -1.0 and 1.0.

You can’t run away from this subtractions because crossing the 0 causes a lot of trouble. We are sweeping all the range of possible exponents from 1 to the infinitesimal… So you need this operation to calculate the proper exponents for the numbers close to 0.

So here is the code.

static inline double myrand(){
  /* ******************************************************************
     Generates a random floating-point number between -1.0 and
     1.0. Not sure if this is the best precision possible, but is
     based on bitwise operations, and I believe it should be faster
     than using (float)rand()/RAND_MAX.

     We pick up the number generated by rand, and leave just the bits
     over the mantissa region. The signal and exponent part are
     replaced to generate a number between 2.0 and 3.9999999. We
     then subtract 3.0 to bring the range down to -1.0 and 0.9999999.

     Coded in 2010-09-24 by Nicolau Werneck <nwerneck@gmail.com>.
     *************************************************************** */
  static union {
    unsigned  int i;
    float f;
  } myrand;
  myrand.i = rand() & 0x007fffff | 0x40000000;
  return myrand.f-3.0;
};

I’m not sure this is the most perfect algorithm and implementation possible, but I guess it’s pretty good. I took the int/float union thing from the fast reciprocal-of-square-root code example at wikipedia, that I mentioned in the previous post. The two bitwise operations just zero all the non-mantissa bits first, and then set the bits we need there (a single one, BTW). If you need to generate values in another range, like -0.5 and 0.5, you should adjust this part to set the proper exponent. Notice that it’s only easy to make ranges that are powers of two, anything different requires more complicated code. And you will always need a final addition for e.g. ranges that start at 0.

There is probably some smart optimization we can still make to both codes. Here is the output generated by the compiler for this code of mine:

  4006e8:	e8 cb fe ff ff       	callq  4005b8 <rand@plt>
  4006ed:	25 ff ff 7f 00       	and    $0x7fffff,%eax
  4006f2:	0d 00 00 00 40       	or     $0x40000000,%eax
  4006f7:	89 c3                	mov    %eax,%ebx
  4006f9:	89 5d ec             	mov    %ebx,-0x14(%rbp)
  4006fc:	f3 0f 10 45 ec       	movss  -0x14(%rbp),%xmm0
  400701:	f3 0f 10 0d a3 06 00 	movss  0x6a3(%rip),%xmm1        # 400dac <_IO_stdin_used+0x6c>
  400708:	00
  400709:	f3 0f 5c c1          	subss  %xmm1,%xmm0

The first line is the call to rand(). Then it performs the two bitwise operations, and we see a step that I would like to eliminate: saving the resulting value to a memory variable that we end up not using for anything else. The value is saved, then loaded again, and moved to an MMX register (this is x86 code, as you should have noticed by now). Then another value gets loaded from memory, that is most certainly the 3.0 in our code. The two are subtracted, and the result is returned afterwards.

I hope this can help someone. There are many blogs out there that just tell you how to use the method based on converting from integer to floating-point, and then using a division. It’s not much hard to do this when you know how floating-point representation works.

Division in more complex processors can be quite fast, actually, so this code may not perform much better. You should always test in your applications if using the floating-point method is really much better. And you should always pay attention if the real bottlenecks in your code aren’t related for example to memory access.

My final inline assembly macro to produce a random floating point number between -1 and 1 is here:

#define NIC_SSE_RAND2(A) asm ("movd %[rand],%%xmm0;" "addss %%xmm0,%[out];" :[out] "=x" ((A)) : "x" (-3.0f) , [rand] "r" ((rand() & 0x007fffff | 0x40000000)) : "%xmm0");

And the assembly looked very nice (Notice how it’s quite similar to the other one, only a few differences). It was approximately 30% faster than using division in my tests.

  400b58:	44 8b 65 e4          	mov    -0x1c(%rbp),%r12d
  400b5c:	e8 57 fa ff ff       	callq  4005b8 <rand@plt>
  400b61:	25 ff ff 7f 00       	and    $0x7fffff,%eax
  400b66:	0d 00 00 00 40       	or     $0x40000000,%eax
  400b6b:	f3 0f 10 0d 35 07 00 	movss  0x735(%rip),%xmm1        # 4012a8 <_IO_stdin_used+0x98>
  400b72:	00
  400b73:	66 0f 6e c0          	movd   %eax,%xmm0
  400b77:	f3 0f 58 c8          	addss  %xmm0,%xmm1
  400b7b:	48 8b 85 70 ff ff ff 	mov    -0x90(%rbp),%rax
  400b82:	49 63 d4             	movslq %r12d,%rdx
  400b85:	f3 0f 11 0c 90       	movss  %xmm1,(%rax,%rdx,4)