Gibbs sampler for 1d Ising model

This depends on Matsumoto's Mersenne Twister.

gibbs1d.h:
#ifndef __GIBBS1D_H__
#define __GIBBS1D_H__

#include <vector>

/**
 * Produces a sequence of variates with P(x) propto exp(-beta * H(nodes))
 * as its stationary distribution, where H(nodes) = the number of consecutive
 * pairs with equal state.
 */
class Gibbs1d {
  
public:

  /**
   * beta -- larger values imply larger affinities for forming clusters.
   * numNodes -- number of points along the line.
   * q -- number of states for each point.
   */
  Gibbs1d(double beta, int numNodes, int q = 2);

  /**
   * Update the state of each node by Gibbs sampling.
   */
  void iterate();

  /**
   * Returns a pointer to a numNodes length array giving
   * the current state of each node.
   */
  const int* current() const;

private:
  //Hide the copy ctor and assignment operator
  Gibbs1d(const Gibbs1d&);
  Gibbs1d& operator=(const Gibbs1d&);

  const double beta;
  const int numNodes;
  const int q;
  std::vector<int> states;
  const double probNotEq1;
  const double probNotEq2;
  const double probEq;
  const double prob1Out;
};

#endif

gibbs1d.cpp:
#include <algorithm>
#include <cmath>
#include "gibbs1d.h"
#include "mtrandom.h"
using namespace std;

/**
 * beta -- larger values imply larger affinities for forming clusters.
 * numNodes -- number of points along the line.
 * q -- number of states for each point.
 */
Gibbs1d::Gibbs1d(double beta, int numNodes, int q)
  : beta(beta),
    numNodes(numNodes),
    q(q),
    states(numNodes + 2),
    probNotEq1(exp(beta) / (2.0 * exp(beta) + q - 2)),
    probNotEq2(2.0 * probNotEq1),
    probEq(exp(2.0 * beta) / (exp(2.0 * beta) + q - 1)),
    prob1Out(exp(beta) / (exp(beta) + q - 1))
{
  //Add sentinel values to the ends to avoid special casing
  //the boundaries.
  states[0] = states[numNodes + 1] = -1;
}

/**
 * Returns an integer in the interval [0, q-1]
 * != exclude, assigning equal probability to each
 * allowed value.
 * Assumes exclude is in [0, q-1].
 */
static int uniformNotEq(int q, int exclude)
{
  int r = (int) ((q - 1) * genrand_real2());
  if (r >= exclude) {
    return r + 1;
  }
  else {
    return r;
  }
}

/**
 * Returns an integer in the interval [0, q-1]
 * != exclude1 or exclude2, assigning equal probability to each
 * allowed value.
 * Assumes exclude1 and exclude2 are in [0, q-1]
 * and that exclude1 != exclude2.
 */
static int uniformNotEq(int q, int exclude1, int exclude2)
{
  if (exclude1 > exclude2) {
    swap(exclude1, exclude2);
  }

  int r = (int) ((q - 2) * genrand_real2());
  if (r >= exclude1) {
    r = r + 1;
  }
  
  if (r >= exclude2) {
    r = r + 1;
  }

  return r;
}

/**
 * Returns the given state with the given probability.
 * Otherwise, chooses one of the other q - 1 states
 * uniformly.
 */
static int withProbElseUniformNotEq(double prob, int state, int q)
{
  const double u = genrand_real2();
  if (u <= prob) {
    return state;
  }
  else {
    return uniformNotEq(q, state);
  }
}

void Gibbs1d::iterate()
{
  for (int i = 1; i <= numNodes; ++i) {
    const int left = states[i - 1];
    const int right = states[i + 1];

    if (left >= 0) {
      if (right >= 0) {
 if (left != right) {
   const double u = genrand_real2();
   if (u <= probNotEq1) {
     states[i] = left;
   }
   else if (u <= probNotEq2) {
     states[i] = right;
   }
   else {
     states[i] = uniformNotEq(q, left, right);
   }
 }
 else { // left == right
   states[i] = withProbElseUniformNotEq(probEq, left, q);
 }
      }
      else { //right < 0
 states[i] = withProbElseUniformNotEq(prob1Out, left, q);
      }
    }
    else {
      if (right >= 0) {
 states[i] = withProbElseUniformNotEq(prob1Out, right, q);
      }
      else {
 states[i] = (int) (q * genrand_real2());
      }
    }
  }
}

/**
 * Returns a pointer to a numNodes length array giving
 * the current state of each node.
 */
const int* Gibbs1d::current() const
{
  return &states[1];
}

I think there must be a much more elegant implementation of this. For one thing, the use of "sentinel" values isn't providing much of a payoff here. And the code for computing the probabilities is just very ugly. What's a better alternative to the nested if-statements?

//I really hate how Blogger cuts off the edges of wide elements. I suppose I can edit the style to fix this.