Statistics for MySQL  0.9
sqlrand.cc
Go to the documentation of this file.
00001 /* sqlrand.cc (Random numbers) */
00002 
00003 /***********************************************************************
00004 *  This code is part of Statistics for MySQL.
00005 *
00006 *  Copyright (C) 2011 Heinrich Schuchardt (xypron.glpk@gmx.de)
00007 *
00008 *  Licensed under the Apache License, Version 2.0 (the "License");
00009 *  you may not use this file except in compliance with the License.
00010 *  You may obtain a copy of the License at
00011 *
00012 *      http://www.apache.org/licenses/LICENSE-2.0
00013 *
00014 *  Unless required by applicable law or agreed to in writing, software
00015 *  distributed under the License is distributed on an "AS IS" BASIS,
00016 *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
00017 *  See the License for the specific language governing permissions and
00018 *  limitations under the License.
00019 ***********************************************************************/
00020 
00027 #include "sqlrand.h"
00028 #include "sqlthread.h"
00029 #include <stdlib.h>
00030 #include <time.h>
00031 #include <math.h>
00032 
00033 #include <stdio.h>
00034 
00035 namespace sqlstat{
00036 
00037   static pthread_mutex_t THR_LOCK_mt;
00038   static pthread_mutex_t THR_LOCK_gaussian;
00039   static bool mt_initialized = false;
00040   static int mt_index;
00041   static unsigned long mt_buffer[MT_LEN];
00042   static double nextGaussian;
00043   static bool haveNextGaussian = false;
00044 
00045   void MersenneTwister::deinit() {
00046     pthread_mutex_destroy(&THR_LOCK_gaussian);
00047     pthread_mutex_destroy(&THR_LOCK_mt);
00048     mt_initialized = false;
00049   }
00050   
00051   void MersenneTwister::init() {
00052     int i;
00053     if (mt_initialized) {
00054       return;
00055     }
00056     pthread_mutex_init(&THR_LOCK_mt, NULL);
00057     pthread_mutex_init(&THR_LOCK_gaussian, NULL);
00058     pthread_mutex_lock(&THR_LOCK_mt);
00059     mt_initialized = true;
00060     srand((unsigned int) time(NULL));
00061     for (i = 0; i < MT_LEN; i++) {
00062       mt_buffer[i] =  (unsigned long) rand() << 16;
00063       mt_buffer[i] += (unsigned long) rand() & 0xFFFF;
00064     }
00065     mt_index = MT_LEN;
00066     pthread_mutex_unlock(&THR_LOCK_mt);
00067   }
00068     
00069   unsigned long MersenneTwister::random() {
00070     unsigned long s;
00071     unsigned long r;
00072     int i;
00073     
00074     init();
00075     pthread_mutex_lock(&THR_LOCK_mt);
00076     
00077     if (mt_index >= MT_LEN) {
00078       mt_index = 0;
00079       i = 0;
00080       for (; i < MT_IB; i++) {
00081         s = TWIST(mt_buffer, i, i+1);
00082         mt_buffer[i] = mt_buffer[i + MT_IA] ^ (s >> 1) ^ MAGIC(s);
00083       }
00084       for (; i < MT_LEN-1; i++) {
00085         s = TWIST(mt_buffer, i, i+1);
00086         mt_buffer[i] = mt_buffer[i - MT_IB] ^ (s >> 1) ^ MAGIC(s);
00087       }
00088       
00089       s = TWIST(mt_buffer, MT_LEN-1, 0);
00090       mt_buffer[MT_LEN-1] = mt_buffer[MT_IA-1] ^ (s >> 1) ^ MAGIC(s);
00091     }
00092     mt_index++;
00093     r = mt_buffer[mt_index];
00094     pthread_mutex_unlock(&THR_LOCK_mt);
00095     /* Tempering transform to compensate for the reduced dimensionality of 
00096      * equidistribution */
00097     r ^= (r >> 11);
00098     r ^= (r << 7) & 0x9D2C5680;
00099     r ^= (r << 15) & 0xEFC60000;
00100     r ^= (r >> 18);
00101     return r;
00102   }
00103 
00104   double MersenneTwister::nextDouble() {
00105     return (double)
00106       ((((unsigned long long) random() & 0x00000000FFFFFFC0L) << 21) 
00107       ^ ((unsigned long long) random() >> 5)) 
00108       / (double) 0x0020000000000000L;
00109     }
00110 
00111   double MersenneTwister::gaussian() {
00112     double a1, a2, q, ret;
00113     pthread_mutex_lock(&THR_LOCK_gaussian);
00114     if (haveNextGaussian) {
00115       haveNextGaussian = false;
00116       ret = nextGaussian;
00117     } else {
00118       do {
00119       a1 = 2. * nextDouble() - 1.;
00120       a2 = 2. * nextDouble() - 1.;
00121       q = a1 * a1 + a2 * a2;
00122       } while ( q <= 0 || q > 1 );
00123       q = sqrt(-2 * log(q) / q );
00124       nextGaussian = a1 * q;
00125       haveNextGaussian = true;
00126       ret = a2 * q;
00127     }
00128     pthread_mutex_unlock(&THR_LOCK_gaussian);
00129     return ret;
00130   }
00131   
00132 }
 All Classes Files Functions Variables Typedefs Defines