Lemire's almost divisionless random int in a range, extended with some fancy 32/64 extensible arithmetic

divless.cpp

#include <inttypes.h>
#include <stdio.h>
#include <stdlib.h>
#include <random>

typedef uint32_t (*randfnc32)(void);

/*
  To generate a random integer within a range, we start
  with a random fractional number (in [0,1)) and multiply it
  by the range.  The integer part of the result is the output,
  but we check the fractional part and reject the result if
  it biases the result (some integers have one more random value 
  than others).
 */
struct FixedPoint{
  uint64_t val;

  FixedPoint() : val(0) {};
  void setFraction(uint32_t fraction) {val = fraction;};
  void addFraction(uint32_t fraction) {val += fraction;};
  FixedPoint& operator *= (uint32_t x) { val *= x; return *this; };

  uint32_t fraction() { return (uint32_t) val; };
  uint32_t floor() { return val>>32; }
};

/*
  RangeGenerator generates random integers in a range without bias.

  It uses fixed point arithmetic to generate a value in [0,1) and multiplies it by range to get
  a number x in [0, range).  It returns floor(x) unless x is an "extra" value which is defined as follows:
 
  Since x is a multiple of range, 

  1. The fractional part of x is the first possible fraction corresponding to floor(x).  

  2. This first possible fraction is small enough for there to be room for one more value at the end.

  
 */
template <randfnc32 Random32>
class RangeGenerator
{
private:
  uint32_t range;
  uint32_t _extraValueThreshold;

  /*
    fraction values occur every
    (range) positions, so a fraction less than (range) must be the first
   */
  bool isFirstFraction(uint32_t f) { return f < range; };

  /*
    Calculate 2^32 mod range (a 33-bit operation)
    
    (2^32 - range) mod range is the same thing and fits in 32 bits

    This is the most expensive operation (40+ cycles), so we 
    lazy-evaluate and save the result.

   */
  uint64_t extraValueThreshold() 
  { 
    return _extraValueThreshold == UINT_MAX ? 
      _extraValueThreshold = -range % range : _extraValueThreshold;
  }

  /*
    Although we know fraction values occur at a regular interval, their alignment 
    with the start of the fraction range [0,2^32) may vary.

   If the first value is small enough (less than (2^32 mod range)), 
   we know there is one extra value, so we discard it to 
   maintain unbiased-ness.
   */
  bool isExtraValue(uint32_t f) { return f < extraValueThreshold(); };
  bool isRejectedValue(uint32_t f) { return isFirstFraction(f) && isExtraValue(f); }

public:
  RangeGenerator(uint32_t _range): range(_range), _extraValueThreshold(UINT_MAX)
  {
    if(_range <= 1)
      throw std::invalid_argument("range must be greater than 1");
  };

  uint32_t generate()
  {
    FixedPoint x;
    do {
      x.setFraction(Random32());
      x *= range;
    } while(isRejectedValue(x.fraction()));

    return x.floor();
  }
};

// Similar idea but using 32.32 extending to 32.64 fixed point only as needed
/*
  With the RangeGenerator above, when range approaches 2^32, a large number of 
  generated values are less than range, meaning they have to be checked against 
  the expensive modulus.  To lessen that case, we expand the fraction value range to 2^64.

  This version uses a 64-bit fraction which is only fully calculated if the most significant 32 bits 
  don't preclude rejection.

  Rejection requires the 64-bit fraction to be less than (range), meaning its upper 32 bits must be 0.
  So we check these conditions in order before rejecting:

      1. upper 32 bits of fraction could become 0 if the product is extended by 32 more random bits.
      2. the upper 32 bits of the fraction do in fact become 0.
      3. the 64-bit fraction is less than 2^64 mod range.

  To become zero after the 32 bit extension, the upper 32 bits must be zero or within
  (range) below zero, since extending another 32 bits will only add up to (range-1) 
  carried from the right.  Negation, even for unsigned, gives us this distance.

 */
template <randfnc32 Random32>
class RangeGeneratorExtended 
{
private:
  uint32_t range;

  uint32_t extraValueThreshold() { return - (uint64_t) range % range; };  // 2^64 mod range
  bool isFirstFraction(uint32_t f) { return f < range; };
  bool isExtraValue(uint32_t f) { return f < extraValueThreshold(); };
  bool isRejectedValue(uint32_t f) { return isFirstFraction(f) && isExtraValue(f); }

public:

  RangeGeneratorExtended(uint32_t _range): range(_range)
  {}

  uint32_t generate()
  {
    FixedPoint x, addend; 

    do {
      x.setFraction(Random32());
      x *= range;

      if( - x.fraction() < range )
	{
	  addend.setFraction(Random32());
	  addend *= range;

	  x.addFraction(addend.floor());
	}
    } while(x.fraction() == 0 && isRejectedValue(addend.fraction()));

    return x.floor();
  }
      
};

uint32_t rnd32(void) {
  return random()<<1 ^ random();
}
 
int main(){

  RangeGenerator<rnd32> rg(10);
  RangeGeneratorExtended<rnd32> rge(10);

  for(int i=0; i < 50; i++)
    printf("val:  %d %d\n", rg.generate(), rge.generate());

}