![]() |
Statistics for MySQL
0.9
|
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 }