1 /******************************************************************************* 2 3 With gratitude to Dr Jurgen A Doornik. See his paper entitled 4 "Conversion of high-period random numbers to floating point" 5 6 Copyright: 7 Copyright (c) 2008 Tango contributors. 8 Some parts copyright (c) 2009-2016 dunnhumby Germany GmbH. 9 All rights reserved. 10 11 License: 12 Tango Dual License: 3-Clause BSD License / Academic Free License v3.0. 13 See LICENSE_TANGO.txt for details. 14 15 Version: Initial release: May 2008 16 17 Authors: Various 18 19 *******************************************************************************/ 20 21 module ocean.math.random.Kiss; 22 23 import core.sys.posix.sys.time : timeval, gettimeofday; 24 25 /****************************************************************************** 26 27 KISS (from George Marsaglia) 28 29 The idea is to use simple, fast, individually promising 30 generators to get a composite that will be fast, easy to code 31 have a very long period and pass all the tests put to it. 32 The three components of KISS are 33 --- 34 x(n)=a*x(n-1)+1 mod 2^32 35 y(n)=y(n-1)(I+L^13)(I+R^17)(I+L^5), 36 z(n)=2*z(n-1)+z(n-2) +carry mod 2^32 37 --- 38 39 The y's are a shift register sequence on 32bit binary vectors 40 period 2^32-1; The z's are a simple multiply-with-carry sequence 41 with period 2^63+2^32-1. The period of KISS is thus 42 --- 43 2^32*(2^32-1)*(2^63+2^32-1) > 2^127 44 --- 45 46 Note that this should be passed by reference, unless you really 47 intend to provide a local copy to a callee 48 49 ******************************************************************************/ 50 51 struct Kiss 52 { 53 /// 54 public alias natural toInt; 55 /// 56 public alias fraction toReal; 57 58 private uint kiss_k; 59 private uint kiss_m; 60 private uint kiss_x = 1; 61 private uint kiss_y = 2; 62 private uint kiss_z = 4; 63 private uint kiss_w = 8; 64 private uint kiss_carry = 0; 65 66 private enum double M_RAN_INVM32 = 2.32830643653869628906e-010, 67 M_RAN_INVM52 = 2.22044604925031308085e-016; 68 69 /********************************************************************** 70 71 A global, shared instance, seeded via startup time 72 73 **********************************************************************/ 74 75 public static Kiss instance; 76 77 static this () 78 { 79 instance.seed; 80 } 81 82 /********************************************************************** 83 84 Creates and seeds a new generator with the current time 85 86 **********************************************************************/ 87 88 static Kiss opCall () 89 { 90 Kiss rand; 91 rand.seed; 92 return rand; 93 } 94 95 /********************************************************************** 96 97 Seed the generator with current time 98 99 **********************************************************************/ 100 101 void seed () 102 { 103 ulong s; 104 105 timeval tv; 106 107 gettimeofday (&tv, null); 108 s = tv.tv_usec; 109 110 return seed (cast(uint) s); 111 } 112 113 /********************************************************************** 114 115 Seed the generator with a provided value 116 117 **********************************************************************/ 118 119 void seed (uint seed) 120 { 121 kiss_x = seed | 1; 122 kiss_y = seed | 2; 123 kiss_z = seed | 4; 124 kiss_w = seed | 8; 125 kiss_carry = 0; 126 } 127 128 /********************************************************************** 129 130 Returns X such that 0 <= X <= uint.max 131 132 **********************************************************************/ 133 134 uint natural () 135 { 136 kiss_x = kiss_x * 69069 + 1; 137 kiss_y ^= kiss_y << 13; 138 kiss_y ^= kiss_y >> 17; 139 kiss_y ^= kiss_y << 5; 140 kiss_k = (kiss_z >> 2) + (kiss_w >> 3) + (kiss_carry >> 2); 141 kiss_m = kiss_w + kiss_w + kiss_z + kiss_carry; 142 kiss_z = kiss_w; 143 kiss_w = kiss_m; 144 kiss_carry = kiss_k >> 30; 145 return kiss_x + kiss_y + kiss_w; 146 } 147 148 /********************************************************************** 149 150 Returns X such that 0 <= X < max 151 152 Note that max is exclusive, making it compatible with 153 array indexing 154 155 **********************************************************************/ 156 157 uint natural (uint max) 158 { 159 return natural % max; 160 } 161 162 /********************************************************************** 163 164 Returns X such that min <= X < max 165 166 Note that max is exclusive, making it compatible with 167 array indexing 168 169 **********************************************************************/ 170 171 uint natural (uint min, uint max) 172 { 173 return (natural % (max-min)) + min; 174 } 175 176 /********************************************************************** 177 178 Returns a value in the range [0, 1$(RPAREN) using 32 bits 179 of precision (with thanks to Dr Jurgen A Doornik) 180 181 **********************************************************************/ 182 183 double fraction () 184 { 185 return ((cast(int) natural) * M_RAN_INVM32 + (0.5 + M_RAN_INVM32 / 2)); 186 } 187 188 /********************************************************************** 189 190 Returns a value in the range [0, 1$(RPAREN) using 52 bits 191 of precision (with thanks to Dr Jurgen A Doornik) 192 193 **********************************************************************/ 194 195 double fractionEx () 196 { 197 return ((cast(int) natural) * M_RAN_INVM32 + (0.5 + M_RAN_INVM52 / 2) + 198 ((cast(int) natural) & 0x000FFFFF) * M_RAN_INVM52); 199 } 200 } 201 202 203 204 /****************************************************************************** 205 206 207 ******************************************************************************/ 208 209 debug (Kiss) 210 { 211 import ocean.io.Stdout; 212 import ocean.time.StopWatch; 213 214 void main() 215 { 216 auto dbl = Kiss(); 217 auto count = 100_000_000; 218 StopWatch w; 219 220 w.start; 221 double v1; 222 for (int i=count; --i;) 223 v1 = dbl.fraction; 224 Stdout.formatln ("{} fraction, {}/s, {:f10}", count, count/w.stop, v1); 225 226 w.start; 227 for (int i=count; --i;) 228 v1 = dbl.fractionEx; 229 Stdout.formatln ("{} fractionEx, {}/s, {:f10}", count, count/w.stop, v1); 230 231 for (int i=count; --i;) 232 { 233 auto v = dbl.fraction; 234 if (v <= 0.0 || v >= 1.0) 235 { 236 Stdout.formatln ("fraction {:f10}", v); 237 break; 238 } 239 v = dbl.fractionEx; 240 if (v <= 0.0 || v >= 1.0) 241 { 242 Stdout.formatln ("fractionEx {:f10}", v); 243 break; 244 } 245 } 246 } 247 }