00001 #ifndef COMMONS_RAND_H 00002 #define COMMONS_RAND_H 00003 00004 #include <cassert> 00005 #include <cstdlib> // random, RAND_MAX 00006 00007 namespace commons 00008 { 00015 class posix_rand { 00016 private: 00017 unsigned long randx; 00018 public: 00019 inline posix_rand(long s = 0) : randx(s) {} 00020 inline void seed(long s) { randx = s; } 00021 inline long abs(long x) { return x & 0x7fffffff; } 00022 inline long draw() { return randx = randx * 1103515245 + 12345; } 00023 inline long operator()() { return abs(draw()); } 00024 }; 00025 00030 inline int 00031 randint(int max = RAND_MAX) 00032 { 00033 // My flow of thought: 00034 // 00035 // It seems simplest to just mod. But what if not each bit was 00036 // identically random? What if the randomness was only a uniform spread 00037 // over RAND_MAX, but all the integers are even (and we only ended up 00038 // selecting the lowest bit because max = 2)? 00039 // 00040 // We want to solve for x: 00041 // 00042 // random() ./ RAND_MAX = x ./ max 00043 // x = max * random() ./ RAND_MAX 00044 // 00045 // But the multiplication could overflow if RAND_MAX is near the max of 00046 // long int (random()'s return type), so re-write: 00047 // 00048 // random() / (RAND_MAX / max) 00049 // 00050 // Now nothing overflows, but information is lost. Is that OK? No! 00051 // If RAND_MAX = 5 and max = 3 and random() = 4, then 00052 // 00053 // 4 / ((5 / 3 = 1) + 1) = 4 which is >= 3! 00054 // 00055 // Makes sense, since: 00056 // 00057 // RAND_MAX / max <= RAND_MAX ./ max 00058 // 00059 // so 00060 // 00061 // random() / quotient >= random() / real_quotient 00062 // 00063 // More "spatious" example: 00064 // 00065 // (random() = 20000) / ( ( (RAND_MAX = 32768) / (max = 16637) ) = 1 ) = 00066 // 20000 which is > max. 00067 // 00068 // Let's tweak that formula by incrementing the denominator so that it's 00069 // impossible for the quotient to exceed max. 00070 // 00071 // random() / (RAND_MAX / max + 1) 00072 // 00073 // How does this fare? 00074 // 00075 // 4 / ((5 / 2 = 2) + 1 = 3) = 1, which is < 2 00076 // 4 / ((5 / 5 = 1) + 1 = 2) = 2, which is < 5, but too much less 00077 // 32767 / ((32768 / 32760 = 1) + 1 = 2) = 16633, 00078 // which is < 32760, but too much less 00079 // 00080 // Think of the problem in terms of pixel line drawing algorithms; we can 00081 // only take steps up at fixed periods, so it's impossible to choose a 00082 // period that will span enough of the space up to max (in the worst case, 00083 // you'll span nearly half the space). 00084 // 00085 // 0 / ((5 / 4 = 1) + 1 = 2) = 0, which is < 4 00086 // 1 / ((5 / 4 = 1) + 1 = 2) = 0, which is < 4 00087 // 2 / ((5 / 4 = 1) + 1 = 2) = 1, which is < 4 00088 // 3 / ((5 / 4 = 1) + 1 = 2) = 1, which is < 4 00089 // 4 / ((5 / 4 = 1) + 1 = 2) = 2, which is < 4 00090 // 4 / ((5 / 5 = 1) + 1 = 2) = 2, which is < 4 00091 // 00092 // random() / (RAND_MAX / (max - 1)) doesn't work either (it's even worse): 00093 // 00094 // 4 / (5 / (3 - 1 = 2) = 2) = 2, which is < 3 00095 // 4 / (5 / (5 - 1 = 4) = 1) = 4, which is < 5 00096 // 4 / (5 / (2 - 1 = 1) = 5) = 0, which is < 2, but too much less 00097 // 5 / (6 / (3 - 1 = 2) = 3) = 1, which is < 3, but too much less 00098 // 00099 // OK, this path isn't taking us anywhere. Let's go back to the original 00100 // formula. Sure, it exceeds max. But what if we just wrap it back around 00101 // when it does? We can do so using % max. Success! 00102 00103 return static_cast<int>( random() / ( RAND_MAX / max ) % max ); 00104 } 00105 00110 inline int 00111 randint(int min, int max) 00112 { 00113 return max > min + 1 ? min + randint(max - min) : min; 00114 } 00115 } 00116 00117 #endif
1.5.6