More Effective Thinking Design Patterns
in the Exceptional C++ Programming Language

These notes share new tricks for creative expression in the C++ programming language.
Your feedback and contributions are welcome!
Mail: colas.schretter [at] gmail.com.

I.Foreword
II.Preliminary Readings
III.The Singulator Template
IV.Variable Column Types Matrix
V.Variable Row Types Matrix
VI.Evaluation of Expressions
VII.Command Line Parsing
VIII.Uniform Pseudo-Random Numbers
IX.Gaussian Pseudo-Random Numbers
X.Uniform Pseudo-Random Points on Sphere

I. Foreword

The aim of these notes is to remember some of the thoughts I had while discussing with friends about how to program in C++. Maybe it was worth explaining and sharing some of them instead of filling a dusty drawer with them.

The title "More Effective Thinking Design Patterns in the Exceptional C++ Programming Language" has been first coined by my friend Utku Salihoglu, in an ultimate attempt to express our grateful acknowledgment to the books that steeped our learning curve.

Each article is built upon the literate programming paradigm introduced by Don. Knuth. The presentation emphasizes explanations while an implementation is systematically provided on a gray background. Each program compiles without any dependency.

The page should have been larger, but I'm now quite busy with other projects. Some of the solutions can be considered as free-style figures, others as handy code snippets. I hope that a part of these notes could be at least useful for someone.

II. Preliminary Readings

The C++ Programming Language is the reference book from the author of C++, Bjarne Stroustrup. While somewhat incomplete, this book provides guidelines for C++ programming and introduces the STL.

Thinking in C++ is an introductory book to C++, written with humor. The author, Bruce Eckel shows great pedagogic skills. This book is an ideal first reading to dive in C++.

Design Patterns, by Erich Gamma et al. is the true revelation in software engineering from the nineties. You cannot pretend to understand what Object Orientation is without reading this catalog of software building-blocks. Indeed, modern OOD rely solely on Design Patterns.

Effective C++ by Scott Meyer is a collection of advices and hints. The book open the curtain to unveil how creative, complicated and enjoying C++ is. Based on the success of the first book, More Effective C++ and Effective STL complete the collection.

Exceptional C++ by Herb Sutter is a representative book of the C++ In-Depth series books, edited by Bjarne Stroustrup himself. The book challenges the reader with programming puzzles and focus on the, often misunderstood, exception mechanism of C++.

III. The Singulator Template

One day, a friend asked my preferred naming to get the unique instance of a singleton object of class Foo. Basically, there is two mainstream schools: singleton.instance() and singleton.foo(). We came to the conclusion that none of these ways really fit since we can not directly use the singleton object anyway.

A singleton type definition ensures that only one instance of this type can be created in a program. Therefore, the usual trick is to declare the constructors private and delegate the creation of a hidden instance to a public method.

I propose an alternative singularization template that ensures unique instantiation of the templatized type. This universal template definition allows to singularize any non-singleton type.

Let the definition of a class Foo:

class Foo {
  int value;

public:
  Foo(): value(0) {}

  void set_value(const int i) {
    value = i;
  }

  int get_value() const {
    return value;
  }
};

An unique instance of the class can be singularized by a template specialization.

template<typename Type>
class Singulator {
  static Type* instance;

public:
  Singulator() {
    if(!instance)
      instance = new Type();
  }

  Type* operator -> () {
    return instance;
  }
};

template<>
Foo* Singulator<Foo>::instance = (Foo*)0;

Several singularized objects can be created while sharing the same underlying instance as shown in the following example:

#include <iostream>
using namespace std;

int main() {
  Singulator<Foo> bar;
  Singulator<Foo> baz;

  bar->set_value(3);
  cout << baz->get_value() << endl; // print "3"

  return EXIT_SUCCESS;
}

IV. Variable Column Types Matrix

A spreadsheet editor typically presents information as a table where each column contains a different type (floating point value, characters string, integer value...). It would be nice to create such a matrix in C++ and manipulate the rows with regular [] accessors.

For that purpose, I define a collection of matrix elements; each of them inherits from a common interface:

#include <stdexcept>
using namespace std;

class base_element {
public:
  base_element() {}

  virtual ~base_element() {}

  virtual char operator = (char) {
    throw runtime_error("assignation is not allowed");
  }

  virtual int operator = (int) {
    throw runtime_error("assignation is not allowed");
  }

  virtual float operator = (float) {
    throw runtime_error("assignation is not allowed");
  }

  virtual string operator = (string) {
    throw runtime_error("assignation is not allowed");
  }

  virtual void print(ostream&) const {}

  friend ostream& operator << (ostream& os, const base_element& e) {
    e.print(os);
    return os;
  }
};

The default behavior of assignations is triggering a run-time error. Each derived class will implement one of the possible assignation operators, depending on the type of the template attribute:

template<typename T>
class element: public base_element {
  T value;

public:
  T operator = (T rhs) {
    return value = rhs;
  }

  void print(ostream& os) const {
    os << value;
  }
};

I declare a row containing different types, the number and the type of elements are given by a vector of type identifiers.

#include <vector>
#include <string>

class vct_row {
private:
  vector<base_element*> row;

public:
  vct_row(const vector<string> column_types);

  ~vct_row();

  base_element& operator [] (const unsigned int index) {
    return *row[index];
  }
};

At construction, a new matrix element is created according to a string containing the corresponding type identifier. At destruction, the allocated objects are independently freed. The actual definition of the constructor and destructor follow.

#include <typeinfo>

vct_row::vct_row(const vector<string> column_types) {
  typedef vector<string>::const_iterator iter;
  for(iter i = column_types.begin(); i != column_types.end(); ++i) {
    if(*i == typeid(char).raw_name())
      row.push_back(new element<char>);
    else if(*i == typeid(int).raw_name())
      row.push_back(new element<int>);
    else if(*i == typeid(float).raw_name())
      row.push_back(new element<float>);
    else if(*i == typeid(string).raw_name())
      row.push_back(new element<string>);
    else
      throw runtime_error("column type '" + *i + "' is not allowed");
  }
}

vct_row::~vct_row() {
  typedef vector<base_element*>::const_iterator iter;
  for(iter i = row.begin(); i != row.end(); ++i)
    delete *i;
}

To create a matrix, one can define a regular array, or better, a vector of rows:

class vct_matrix {
private:
  vector<vct_row*> rows;

public:
  vct_matrix(const int nbr_rows, const vector<string> column_types);

  ~vct_matrix();

  vct_row& operator [] (const unsigned int index) {
    return *rows[index];
  }
};

vct_matrix::vct_matrix(const int nbr_rows, const vector<string> column_types) {
  for(int i = 0; i != nbr_rows; ++i)
    rows.push_back(new vct_row(column_types));
}

vct_matrix::~vct_matrix() {
  typedef vector<vct_row*>::const_iterator iter;
  for(iter i = rows.begin(); i != rows.end(); ++i)
    delete *i;
}

Finally, the matrix can be defined at run-time and used as a regular array. If a value of a different type than the matrix element is assigned, an exception will be thrown.

#include <iostream>

int main() {
  // list the type names for each column of the matrix
  vector<string> column_types;

  // fill the list of type names dynamically at execution
  column_types.push_back(typeid(char).raw_name());
  column_types.push_back(typeid(float).raw_name());
  column_types.push_back(typeid(string).raw_name());

  // matrix containing five rows
  vct_matrix matrix(5, column_types);

  // assign some values
  matrix[1][0] = 'a';
  matrix[0][1] = 3.14f;
  matrix[4][2] = string("hello world");

  // print some values
  cout << matrix[1][0] << endl;
  cout << matrix[0][1] << endl;
  cout << matrix[4][2] << endl;

  return EXIT_SUCCESS;
}

V. Variable Row Types Matrix

If you liked the previous definition of the variable column types matrix, it's time to address the transpose problem: defining a matrix where each row can contain a different type.

Implementing a solution for the transpose problem is not straightforward because the size of each row can be different, depending on the size in byte that takes each type. I propose to alleviate this problem by storing only pointers that allow indirect access to values. Indeed, every pointer has identical size regardless of the pointed type.

Hence, I modify the previous base_element class by adding a pointer, accessible only by derived classes.

#include <stdexcept>
using namespace std;

class base_element {
protected:
  void* value;

public:
  base_element() {}

  virtual ~base_element() {}

  virtual char operator = (char) {
    throw runtime_error("assignation is not allowed");
  }

  virtual int operator = (int) {
    throw runtime_error("assignation is not allowed");
  }

  virtual float operator = (float) {
    throw runtime_error("assignation is not allowed");
  }

  virtual string operator = (string) {
    throw runtime_error("assignation is not allowed");
  }

  virtual void print(ostream&) const {}

  friend ostream& operator << (ostream& os, const base_element& e) {
    e.print(os);
    return os;
  }
};

As before, the matrix elements redefine one of the assignations; moreover, they handle implicitly the indirection to access values.

template<typename T>
class element: public base_element {
public:
  element() {
    value = new T;
  }

  ~element() {
    delete (T*)value;
  }

  T operator = (T rhs) {
    return *(T*)value = rhs;
  }

  void print(ostream& os) const {
    os << *(T*)value;
  }
};

A collection of type identifiers defines the number and the types of each row of the matrix. Furthermore, the number of columns must also be specified. The matrix allocates and retains a pointer to a regular array of elements for each row.

#include <vector>
#include <string>

class vrt_matrix {
private:
  vector<base_element*> rows;

public:
  vrt_matrix(const int nbr_columns, const vector<string> row_types);

  ~vrt_matrix();

  base_element* operator [] (const unsigned int index) {
    return rows[index];
  }
};

At construction, a new array of elements is created according to a string containing the corresponding type identifier of the row. The length of the arrays is the number of columns of the matrix. At destruction, the allocated arrays are independently freed. The actual definition of the constructor and destructor follow.

#include <typeinfo>

vrt_matrix::vrt_matrix(const int nbr_columns, const vector<string> row_types) {
  typedef vector<string>::const_iterator iter;
  for(iter i = row_types.begin(); i != row_types.end(); ++i) {
    if(*i == typeid(char).raw_name())
      rows.push_back(new element<char>[nbr_columns]);
    else if(*i == typeid(int).raw_name())
      rows.push_back(new element<int>[nbr_columns]);
    else if(*i == typeid(float).raw_name())
      rows.push_back(new element<float>[nbr_columns]);
    else if(*i == typeid(string).raw_name())
      rows.push_back(new element<string>[nbr_columns]);
    else
      throw runtime_error("row type '" + *i + "' is not allowed");
  }
}

vrt_matrix::~vrt_matrix() {
  typedef vector<base_element*>::const_iterator iter;
  for(iter i = rows.begin(); i != rows.end(); ++i)
    delete[] *i;
}

Finally, the variable row type matrix can be used in a very similar way than the previous variable column types matrix, as shown in the following example:

#include <iostream>

int main() {
  // list the type names for each row of the matrix
  vector<string> row_types;

  // fill the list of type names dynamically at execution
  row_types.push_back(typeid(char).raw_name());
  row_types.push_back(typeid(float).raw_name());
  row_types.push_back(typeid(string).raw_name());

  // matrix containing five columns
  vrt_matrix matrix(5, row_types);

  // assign some values
  matrix[0][1] = 'a';
  matrix[1][0] = 3.14f;
  matrix[2][4] = string("hello world");

  // print some values
  cout << matrix[0][1] << endl;
  cout << matrix[1][0] << endl;
  cout << matrix[2][4] << endl;

  return EXIT_SUCCESS;
}

Based on these two examples, it should be easy to implement a fully general matrix, allowing definition of arbitrary types, element-by-element.

VI. Evaluation of Expressions

Bjarne Stroustrup once said: "Example of design that meet most of the criteria for "goodness" (easy to understand, flexible, efficient) is a recursive descent parser, which is traditional procedural code."

As every computer scientist, I wrote a parser to evaluate infix expressions. Since there are quite some computer scientists around, you can imagine the number of expression parser that is already available. Most of them are bloated code. My contribution enters in the lightweight category of parsers in about 200 lines of code.

The following program evaluate mathematical expressions such as for example: "sin(pi/2) + max(1.2, atan(cos(e)), 5 / (8 * 2)) + sqrt(0.1)".

I will assume that you are already familiar with recursive descent parsing and LL grammars. The book dragon book: Compilers: Principles, Techniques, and Tools by Alfred Aho, Ravi Sethi and Jeffrey Ullman could have been included in the preliminary readings section.

First I define a very simple symbol table that holds some constants (e and pi). This class can be considered as a placeholder for a more realistic implementation.

#include <stdexcept>
#include <string>
using namespace std;

// values of identifiers are collected in a shared table
struct Symbols {
  static bool find(const string& name) {
    return name == "pi" || name == "e";
  }

  static float value(const string& name) {
    if(name == "pi")
      return 3.1415926535f;
    else if(name == "e")
      return 2.7182818284f;
    else
      throw runtime_error("undefined identifier '" + name + "'");
  }
};

The following snippet is a universal function evaluator. Given the name of the function and a collection of parameters, the right function is evaluated. The features of the parser are strongly dependent on the functions supported here.

#include <sstream>
#include <vector>
#include <cmath>

// evaluate a function, given its name and a list of parameters
float function(const string name, const vector<float> parameters) {
  // detect arbitrary arity functions
  if(name == "min" && parameters.size() > 0) {
    float minimum = parameters[0];
    for(unsigned int i = 1; i < parameters.size(); ++i)
      if(parameters[i] < minimum)
        minimum = parameters[i];
    return minimum;
  }

  if(name == "max" && parameters.size() > 0) {
    float maximum = parameters[0];
    for(unsigned int i = 1; i < parameters.size(); ++i)
      if(parameters[i] > maximum)
        maximum = parameters[i];
    return maximum;
  }

  if(name == "avg" && parameters.size() > 0) {
    float average = 0;
    for(unsigned int i = 0; i < parameters.size(); ++i)
      average += parameters[i];
    return average / parameters.size();
  }

  // detect unary functions
  if(parameters.size() == 1) {
    if(name == "abs")
      return fabs(parameters[0]);
    else if(name == "floor")
      return floor(parameters[0]);
    else if(name == "ceil")
      return ceil(parameters[0]);
    else if(name == "trunc")
      return int(parameters[0]);
    else if(name == "sqrt")
      return sqrt(parameters[0]);
    else if(name == "sin")
      return sin(parameters[0]);
    else if(name == "cos")
      return cos(parameters[0]);
    else if(name == "tan")
      return tan(parameters[0]);
    else if(name == "asin")
      return asin(parameters[0]);
    else if(name == "acos")
      return acos(parameters[0]);
    else if(name == "atan")
      return atan(parameters[0]);
  }

  // detect binary functions
  if(parameters.size() == 2) {
    if(name == "pow")
      return pow(parameters[0], parameters[1]);
  }

  // unknown function name or mismatch with the number of parameters
  ostringstream nbr_parameters;
  nbr_parameters << parameters.size();
  throw runtime_error("undefined function '" + name + "' taking "
    + nbr_parameters.str() + " parameter(s)");
}

If the syntax of the input expression is not recognized, a parse error exception will be generated by a call to the following function:

// something wrong happens?
void parse_error(istream& is) {
  string context;
  getline(is, context);
  throw runtime_error("parse error before '" + context + "'");
}

The "lexer" reads, char by char, the input and cut the string in tokens. Each token is an identifier, a function name, a separator or a constant numeric value. One token is read ahead. Based on the next token, we can decide which grammar rule to apply.

// look-ahead token
string next_token;

// advance through the stream of chars
void consume(istream& is) {
  // skip leading spaces
  while(isspace(is.peek()))
    is.get();

  // return an non-initialized token if we reached the end of the stream
  if(is.peek() == EOF) return;

  // put at least one char in the next token
  next_token = is.get();

  // separator chars are tokens on themselves
  if(!isalnum(next_token[0])) return;

  // extend further the token if it is not a separator char
  while(isalnum(is.peek()) || is.peek() == '_' || is.peek() == '.')
    next_token += is.get();
}

An expression is recursively defined as a terminal and a sub-expression infixed by a binary operator. Multiplications and divisions have conventionally higher precedence over the additions and subtractions. Hence, an additional parameter keeps track of the precedence level that allows deciding the associativity.

// forward declaration
float parse_terminal(istream& is);

// evaluate the expression with recursive descent parsing
float parse_expression(istream& is, const int precedence = 0) {
  float result = parse_terminal(is);

  while(true) {
    switch(next_token[0]) {
      case '+':
        if(precedence > 0) return result;
        consume(is);
        result += parse_expression(is, 1);
        break;
      case '-':
        if(precedence > 0) return result;
        consume(is);
        result -= parse_expression(is, 1);
        break;
      case '*':
        if(precedence > 1) return result;
        consume(is);
        result *= parse_expression(is, 2);
        break;
      case '/':
        if(precedence > 1) return result;
        consume(is);
        result /= parse_expression(is, 2);
        break;
      default:
        return result;
    }
  }
}

The main features supported by LL grammars are defined in the parsing of terminals, when an expression can not be split in sub-expressions.

// evaluate leafs of the syntax tree
float parse_terminal(istream& is) {
  // just some braces to override precedences?
  if(next_token == "(") {
    consume(is);
    const float result = parse_expression(is);

    // check that braces are well-paired
    if(next_token != ")") parse_error(is);
    consume(is);

    return result;
  }

  // is the token the negation unary operator?
  if(next_token == "-") {
    consume(is);
    return -parse_expression(is, 1);
  }

  // the token holds either a function name, an identifier or a constant
  const string token = next_token;
  consume(is);

  // is the token a function name?
  if(next_token == "(") {
    consume(is);

    // evaluate all parameters of the function
    vector<float> parameters;
    while(next_token != ")") {
      parameters.push_back(parse_expression(is));
      if(next_token == ",") {
        consume(is); // skip the argument separator
        if(next_token == ")") parse_error(is); // at least one parameter
      }
    };

    // consume the closing brace
    consume(is);

    return function(token, parameters);
  }

  // then, probably an identifier?
  if(Symbols::find(token))
    return Symbols::value(token);

  // else, it must contain a constant value!
  float value;
  if(!(istringstream(token) >> value)) parse_error(is);
  return value;
}

Finally, a function encapsulates a string with an input stream and launches the parsing. After evaluation, we check that there is no remaining character left in the input stream.

// parse and evaluate an expression given as a string
float evaluate(const string& expression) {
  istringstream is(expression);
  consume(is); // look-ahead the next token
  const float result = parse_expression(is);
  if(is.peek() != EOF) parse_error(is);
  return result;
}

Finally, the main function shows how to convert the above to an universal calculator that evaluate any expression passed on the command line.

#include <iostream>

int main(int argc, char** argv) {
  try {
    // eat all parameters to form an unique string
    string expression = "";
    for(int i = 1; i < argc; ++i)
      expression += argv[i];

    // evaluate the expression and output the result
    cout << evaluate(expression) << endl;
  } catch(exception& e) {
    cerr << e.what() << endl;
    return EXIT_FAILURE;
  }

  return EXIT_SUCCESS;
}

VII. Command Line Parsing

Anyone needs a straightforward way to extract values from command lines. The system passes arguments as a list of char*. This convention is a legacy from C and also fits well to be parsed by procedural code.

The following program extract values from command lines such as for example: "./program -epsilon 0.01 -answer 42 -message "hello world"":

First, I define a function that seeks for a specific word in the list of possible argument:

#include <stdexcept>
#include <string>
using namespace std;

// common part of the scanning of the command line
int find_argument(const string name, const int argc, const char** argv) {
  // seek for the argument name
  int index = 1;
  while(index < argc && string(argv[index]) != name) ++index;

  // test if the search is successful
  if(index == argc)
    throw runtime_error("unspecified argument '" + name + "'");

  // check that a value follows
  if(index == argc - 1)
    throw runtime_error("missing value after argument '" + name + "'");

  return index + 1;
}

The following template function extracts the value, given an argument name and its type.

#include <sstream>

template<typename Type>
Type argument(const string name, const int argc, const char** argv) {
  Type value;
  if(!(istringstream(argv[find_argument(name, argc, argv)]) >> value))
    throw runtime_error("invalid value for argument '" + name + "'");
  return value;
}

A template specialization overrides the special case for string arguments.

template<>
string argument(const string name, const int argc, const char** argv) {
  return string(argv[find_argument(name, argc, argv)]);
}

Here is a little main function that extracts three different types of parameters. Only one line of code is used to extract value and the above code implicitly handles errors.

#include <iostream>

int main(const int argc, const char** argv) {
  try {
    // extract the value of some arguments
    const string message = argument<string>("-message", argc, argv);
    cout << "message = '" << message << "'" << endl;

    const float epsilon = argument<float>("-epsilon", argc, argv);
    cout << "epsilon = " << epsilon << endl;

    const int answer = argument<int>("-answer", argc, argv);
    cout << "answer = " << answer << endl;
  } catch(exception& e) {
    cerr << e.what() << endl;
    return EXIT_FAILURE;
  }

  return EXIT_SUCCESS;
}

VIII. Uniform Pseudo-Random Numbers

Generation of uniform pseudo-random numbers is more and more needed, mainly for Monte-Carlo simulations of physical processes. The built-in implementation of the C language is known to have very a short limit cycle. So whenever an application relies on random numbers, a custom-made implementation is often needed.

This snippet implements the famous Mersenne twister pseudo-random number generator by Takuji Nishimura and Makoto Matsumoto. The length of the limit cycle is 2¹⁹⁹³⁷-1 and the implementation is very fast while only using integer arithmetic operations. The class generates 624 integers at once and caches the values. Some helper functions allow to return random floating point and integer values in any given a range.

Let the definition of the class Random, the implementation of methods will follow. The class can be used by calling the template methods next() or by considering an instance of Random as a function object (functor) to toss floating point numbers in [0,1).

class Random {
  unsigned int state[624];
  int index;

  // Mersenne twisting
  unsigned int twist(const unsigned int u, const unsigned int v) const;

  // generate 624 integers at once
  void update_state();

public:
  Random(const unsigned int s = 5489);

  void seed(const unsigned int s = 5489);

  // return a number in [0,1) for float and in [0,2^32) for unsigned int
  template<typename Type>
  Type next();

  // return a number in [min,max]
  template<typename Type>
  Type next(const Type min, const Type max);

  // Random instances are function objects (functor)
  float operator() ();
};

The implementation of the algorithm is given by the update_state() method that call the twist() method to combine two 32 bits unsigned integers. I refer the reader to the original paper of Takuji Nishimura and Makoto Matsumoto to understand what happens here.

unsigned int Random::twist(const unsigned int u, const unsigned int v) const {
  if(v & 1)
    return (((u & 0x80000000) | (v & 0x7FFFFFFF)) >> 1) ^ 2567483615;
  else
   return (((u & 0x80000000) | (v & 0x7FFFFFFF)) >> 1);
}

void Random::update_state() {
  for(int i = 0; i < 624 - 397; ++i)
    state[i] = state[i + 397] ^ twist(state[i], state[i + 1]);

  for(int i = 624 - 397; i < 623; ++i)
    state[i] = state[i + 397 - 624] ^ twist(state[i], state[i + 1]);

  state[623] = state[396] ^ twist(state[623], state[0]);
  index = 0; // reset the internal counter
}

At the creation of the Random object, the generator must be seed. I use here the method of Don. Knuth, see The Art of Computer Programming, Vol 2. 3rd Ed. P.106.

Random::Random(const unsigned int s) {
  seed(s);
}

void Random::seed(const unsigned int s) {
  state[0] = s;

  for(int i = 1; i < 624; ++i)
    state[i] = i + (state[i - 1] ^ (state[i - 1] >> 30)) * 1812433253;

  update_state();
}

A final step to ensure an uniform distribution for high order of dimensions (up to 624), is to temper the numbers according to the following procedure, also explained in the original reference.

template<>
unsigned int Random::next() {
  if(index == 624)
    update_state();

  // tempering of the number
  unsigned int r = state[index++];
  r ^= (r >> 11);
  r ^= (r << 7) & 2636928640;
  r ^= (r << 15) & 4022730752;
  r ^= (r >> 18);
  return r;
}

The generation of integer values can be used to generate floating points values by dividing the random integer by it's maximal value 2³².

template<>
float Random::next() {
  return next<unsigned int>() * (1.0f / 4294967296); // 1.0 / (2^32)
}

I propose to rely on a simple rejection sampling method to generate integer numbers in a given range [min,max].

template<>
unsigned int Random::next(const unsigned int min, const unsigned int max) {
  unsigned int mask = max - min;

  // fill the lower bits with 1
  mask |= mask >> 1;
  mask |= mask >> 2;
  mask |= mask >> 4;
  mask |= mask >> 8;
  mask |= mask >> 16;

  // rejection sampling: toss numbers until one is found in [0,max-min]
  unsigned int r;
  do r = next<unsigned int>() & mask; while(r > max - min);

  // return a number in the appropriate range
  return min + r;
}

The generation of floating points values in a given range is similar to the previous template specialization for float.

template<>
float Random::next(const float min, const float max) {
  return min + next<unsigned int>() * (max - min) * (1.0f / 4294967296);
}

Finally, the functor definition:

float Random::operator() () {
  return next<float>();
}

The following main program tests the class Random by evaluating the uniformity of a sequence of numbers.

#include <iostream>
using namespace std;

int main() {
  Random random;
  random.seed(42);

  const int nbr_samples = 1000000;
  const int nbr_bins = 100;
  unsigned int bins[nbr_bins] = {}; // initialize values to zero

  // assign randomly samples to bins
  for(int i = 0; i < nbr_samples; ++i)
    ++bins[random.next<unsigned int>(0, nbr_bins - 1)];

  // retain the upper and lower bounds (worst case)
  unsigned int minimum = nbr_samples, maximum = 0;
  for(int i = 0; i < nbr_bins; ++i)
    if(bins[i] < minimum)
      minimum = bins[i];
    else if(bins[i] > maximum)
      maximum = bins[i];

  // this uniformity factor is very conservative
  const float uniformity = 1 - float(maximum - minimum) / nbr_samples;
  cout << "uniformity = " << 100 * uniformity << '%' << endl;

  return EXIT_SUCCESS;
}

IX. Gaussian Pseudo-Random Numbers

This note provides a solution to pick pseudo-random numbers according to a Gaussian distribution, given a uniform generator. First, let a placeholder function to generate uniform pseudo-random numbers in [0,1):

#include <cstdlib>
using namespace std;

// placeholder function
float uniform() {
  return float(rand()) / RAND_MAX; // [0,1)
}

The most efficient and numerically stable way I found to produce a Gaussian distribution is to rely on the polar form of the Box-Muller transform. This technique leads to the following code:

#include <cmath>

float gaussian(float mean = 0, float std_deviation = 1) {
  static float t = 0;
  float x1, x2, r;

  // reuse previous calculations
  if(t) {
    const float tmp = t;
    t = 0;
    return mean + std_deviation * tmp;
  }

  // pick randomly a point inside the unit disk
  do {
    x1 = 2 * uniform() - 1;
    x2 = 2 * uniform() - 1;
    r = x1 * x1 + x2 * x2;
  } while(r >= 1);

  // Box-Muller transform
  r = sqrt(-2.0 * log(r) / r);

  // save for next call
  t = r * x2;

  // only use one of the coordinates of a bivariate distribution
  return mean + std_deviation * r * x1;
}

The main program tests the distribution of the numbers. They should have a standard deviation (and variance) of one around the mean of zero.

#include <iostream>

int main() {
  const int nbr_samples = 100000;
  float samples[nbr_samples];

  float mean = 0;
  for(int i = 0; i < nbr_samples; ++i)
    mean += samples[i] = gaussian();
  mean /= nbr_samples;

  float variance = 0;
  for(int i = 0; i < nbr_samples; ++i)
    variance += (samples[i] - mean) * (samples[i] - mean);
  variance /= nbr_samples;

  cout << "mean = " << mean << ", std dev. = " << sqrt(variance) << endl;

  return EXIT_SUCCESS;
}

X. Uniform Pseudo-Random Points on Sphere

Monte-Carlo simulations of the photons transport require to pick, very quickly, random points on the surface of a unit sphere. I tried three ways to randomly generate unit direction vectors in three dimensions.

First, let a placeholder function to generate uniform pseudo-random numbers in [0,1):

#include <cstdlib>
using namespace std;

// placeholder function
float uniform() {
  return float(rand()) / RAND_MAX; // [0,1)
}

I also define a basic Vector structure that implements two methods: length2() that returns the squared length of the vector and normalize() that normalizes the vector.

#include <cmath>

struct Vector {
  float x, y, z;

  Vector() {}
  Vector(const float x, const float y, const float z): x(x), y(y), z(z) {}

  float length2() const {
    return x*x + y*y + z*z;
  }

  Vector& normalize() {
    const float inverse = 1 / sqrt(length2());
    x *= inverse, y *= inverse, z *= inverse;
    return *this;
  }
};

The description of the three methods is given along with an implementation:

Naive method: A point is picked randomly in the unit cube with a uniform random generator. If this point lies inside the unit sphere, it is projected on the surface of the sphere by normalizing its coordinates. If the point do not lies inside the unit sphere, then we iterate this procedure. Since the volume of the unit sphere is PI/6, the average number of iterations is equal to 1.909.

inline Vector naive() {
  Vector direction;

  do {
    direction.x = 2 * uniform() - 1,
    direction.y = 2 * uniform() - 1,
    direction.z = 2 * uniform() - 1;
  } while(direction.length2() > 1);

  return direction.normalize();
}

Trig method: Any point on the surface of a unit sphere can be represented by two coordinates. The trig method pick randomly the z in [-1,+1) coordinate to select a slice of the sphere then a point is uniformly selected on this circle by randomly picking an azimuthal angle phi in [0, 2PI). Note that the random generator is only called twice, and no trial-error loop is needed.

inline Vector trig() {
  const float two_pi = 6.2831853071795864770;
  const float phi = two_pi * uniform();
  const float cos_theta = 2 * uniform() - 1;
  const float sin_theta = sqrt(1 - cos_theta * cos_theta);

  return Vector(cos_theta, sin_theta * cos(phi), sin_theta * sin(phi));
}

Hybrid method: A point is picked randomly in the unit disk with a trial-error loop. Since the area of the unit disk is PI/4, the average number of iterations is equal to 1.2732. Hence, the probability of rejection is lower than in three dimensions. The third coordinate is deduced by using the fact that the coordinates of a point on the unit sphere ensure x + y + z = 1. Finally the x and y coordinates are scaled according to the value of z.

inline Vector hybrid() {
  Vector direction;
  float s;

  do {
    direction.x = 2 * uniform() - 1,
    direction.y = 2 * uniform() - 1,
    s = direction.x * direction.x + direction.y * direction.y;
  } while(s > 1);

  const float r = 2 * sqrt(1 - s);

  direction.x *= r;
  direction.y *= r;
  direction.z = 2 * s - 1;

  return direction;
}

When I do not switch on optimization flags of my compiler, then the timings confirm the results of Ken Sloan, referenced in the FAQ of comp.graphics.algorithms. However, with optimization, the naive method is 20% faster than the trig method. In any cases, the hybrid method is the winner.

Similar tests have been conducted for the two-dimensional case. Using polar coordinates, picking a point on a circle only requires one random toss and two calls to sin/cos trigonometric functions. I observed that the performances of the naive rejection method also outperform the trigonometric calculations in two dimensions.

Here finally a little example that uses my favorite method:

#include <iostream>

int main() {
  Vector v = hybrid();
  cout << '<' << v.x << ',' << v.y << ',' << v.z << '>' << endl;

  return EXIT_SUCCESS;
}

More Effective Thinking Design Patternsin the Exceptional C++ Programming Language

Contents