sqlrand.cc

Go to the documentation of this file.

/* sqlrand.cc (Random numbers) */

/***********************************************************************
*  This code is part of Statistics for MySQL.
*
*  Copyright (C) 2011 Heinrich Schuchardt (xypron.glpk@gmx.de)
*
*  Licensed under the Apache License, Version 2.0 (the "License");
*  you may not use this file except in compliance with the License.
*  You may obtain a copy of the License at
*
*      http://www.apache.org/licenses/LICENSE-2.0
*
*  Unless required by applicable law or agreed to in writing, software
*  distributed under the License is distributed on an "AS IS" BASIS,
*  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
*  See the License for the specific language governing permissions and
*  limitations under the License.
***********************************************************************/

#include "sqlrand.h"
#include "sqlthread.h"
#include <stdlib.h>
#include <time.h>
#include <math.h>

#include <stdio.h>

namespace sqlstat {

static pthread_mutex_t THR_LOCK_mt;
static pthread_mutex_t THR_LOCK_gaussian;
static bool mt_initialized = false;
static int mt_index;
static unsigned long mt_buffer[MT_LEN];
static double nextGaussian;
static bool haveNextGaussian = false;

void MersenneTwister::deinit() {
  pthread_mutex_destroy(&THR_LOCK_gaussian);
  pthread_mutex_destroy(&THR_LOCK_mt);
  mt_initialized = false;
}

void MersenneTwister::init() {
  int i;
  if (mt_initialized) {
    return;
  }
  pthread_mutex_init(&THR_LOCK_mt, NULL);
  pthread_mutex_init(&THR_LOCK_gaussian, NULL);
  pthread_mutex_lock(&THR_LOCK_mt);
  mt_initialized = true;
  srand((unsigned int) time(NULL));
  for (i = 0; i < MT_LEN; i++) {
    mt_buffer[i] =  (unsigned long) rand() << 16;
    mt_buffer[i] += (unsigned long) rand() & 0xFFFF;
  }
  mt_index = MT_LEN;
  pthread_mutex_unlock(&THR_LOCK_mt);
}

unsigned long MersenneTwister::random() {
  unsigned long s;
  unsigned long r;
  int i;

  init();
  pthread_mutex_lock(&THR_LOCK_mt);

  if (mt_index >= MT_LEN) {
    mt_index = 0;
    i = 0;
    for (; i < MT_IB; i++) {
      s = TWIST(mt_buffer, i, i+1);
      mt_buffer[i] = mt_buffer[i + MT_IA] ^ (s >> 1) ^ MAGIC(s);
    }
    for (; i < MT_LEN-1; i++) {
      s = TWIST(mt_buffer, i, i+1);
      mt_buffer[i] = mt_buffer[i - MT_IB] ^ (s >> 1) ^ MAGIC(s);
    }

    s = TWIST(mt_buffer, MT_LEN-1, 0);
    mt_buffer[MT_LEN-1] = mt_buffer[MT_IA-1] ^ (s >> 1) ^ MAGIC(s);
  }
  mt_index++;
  r = mt_buffer[mt_index];
  pthread_mutex_unlock(&THR_LOCK_mt);
  /* Tempering transform to compensate for the reduced dimensionality of
   * equidistribution */
  r ^= (r >> 11);
  r ^= (r << 7) & 0x9D2C5680;
  r ^= (r << 15) & 0xEFC60000;
  r ^= (r >> 18);
  return r;
}

double MersenneTwister::nextDouble() {
  return (double)
         ((((unsigned long long) random() & 0x00000000FFFFFFC0L) << 21)
          ^ ((unsigned long long) random() >> 5))
         / (double) 0x0020000000000000L;
}

double MersenneTwister::gaussian() {
  double a1, a2, q, ret;
  pthread_mutex_lock(&THR_LOCK_gaussian);
  if (haveNextGaussian) {
    haveNextGaussian = false;
    ret = nextGaussian;
  } else {
    do {
      a1 = 2. * nextDouble() - 1.;
      a2 = 2. * nextDouble() - 1.;
      q = a1 * a1 + a2 * a2;
    } while ( q <= 0 || q > 1 );
    q = sqrt(-2 * log(q) / q );
    nextGaussian = a1 * q;
    haveNextGaussian = true;
    ret = a2 * q;
  }
  pthread_mutex_unlock(&THR_LOCK_gaussian);
  return ret;
}

}