include/plugin_interface/SC_RGen.h

   1 /*
   2         SuperCollider real time audio synthesis system
   3     Copyright (c) 2002 James McCartney. All rights reserved.
   4         http://www.audiosynth.com
   5
   6     This program is free software; you can redistribute it and/or modify
   7     it under the terms of the GNU General Public License as published by
   8     the Free Software Foundation; either version 2 of the License, or
   9     (at your option) any later version.
  10
  11     This program is distributed in the hope that it will be useful,
  12     but WITHOUT ANY WARRANTY; without even the implied warranty of
  13     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14     GNU General Public License for more details.
  15
  16     You should have received a copy of the GNU General Public License
  17     along with this program; if not, write to the Free Software
  18     Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301  USA
  19 */
  20
  21 //----------------------------------------------------------------------------//
  22 // Ran088: L'Ecuyer's 1996 three-component Tausworthe generator "taus88"
  23 //----------------------------------------------------------------------------//
  24 //
  25 // Returns an integer random number uniformly distributed within [0,4294967295]
  26 //
  27 // The period length is approximately 2^88 (which is 3*10^26).
  28 // This generator is very fast and passes all standard statistical tests.
  29 //
  30 // Reference:
  31 //   (1) P. L'Ecuyer, Maximally equidistributed combined Tausworthe generators,
  32 //       Mathematics of Computation, 65, 203-213 (1996), see Figure 4.
  33 //   (2) recommended in:
  34 //       P. L'Ecuyer, Random number generation, chapter 4 of the
  35 //       Handbook on Simulation, Ed. Jerry Banks, Wiley, 1997.
  36 //
  37 //----------------------------------------------------------------------------//
  38
  39 //----------------------------------------------------------------------------//
  40 // I chose this random number generator for the following reasons:
  41 //              fast.
  42 //              easier and faster to seed than other high quality rng's such as Mersenne Twister.
  43 //              the internal state is only 12 bytes.
  44 //              the period is long enough for music/audio.
  45 //              possible to code in altivec in future if needed.
  46 // - James McCartney
  47 //----------------------------------------------------------------------------//
  48
  49 #ifndef _SC_RGen_
  50 #define _SC_RGen_
  51
  52 #include "SC_Endian.h"
  53 #include "SC_Types.h"
  54 #include "SC_BoundsMacros.h"
  55 #include "Hash.h"
  56 #include <math.h>
  57
  58 struct RGen
  59 {
  60         void init(uint32 seed);
  61
  62         uint32 trand();
  63
  64         int32 irand(int32 scale);
  65         int32 irand2(int32 scale);
  66         int32 ilinrand(int32 scale);
  67         int32 ibilinrand(int32 scale);
  68
  69         float fcoin();
  70         float frand();
  71         float frand2();
  72         float frand0();
  73         float frand8();
  74         double drand();
  75         double drand2(double scale);
  76         double linrand(double scale);
  77         double bilinrand(double scale);
  78         double exprandrng(double lo, double hi);
  79         double exprand(double scale);
  80         double biexprand(double scale);
  81         double sum3rand(double scale);
  82
  83         uint32 s1, s2, s3;              // random generator state
  84 };
  85
  86 inline void RGen::init(uint32 seed)
  87 {
  88         // humans tend to use small seeds - mess up the bits
  89         seed = (uint32)Hash((int)seed);
  90
  91         // initialize seeds using the given seed value taking care of
  92         // the requirements. The constants below are arbitrary otherwise
  93         s1 = 1243598713U ^ seed; if (s1 <  2) s1 = 1243598713U;
  94         s2 = 3093459404U ^ seed; if (s2 <  8) s2 = 3093459404U;
  95         s3 = 1821928721U ^ seed; if (s3 < 16) s3 = 1821928721U;
  96 }
  97
  98 inline uint32 trand( uint32& s1, uint32& s2, uint32& s3 )
  99 {
 100         // This function is provided for speed in inner loops where the
 101         // state variables are loaded into registers.
 102         // Thus updating the instance variables can
 103         // be postponed until the end of the loop.
 104         s1 = ((s1 &  (uint32)-2) << 12) ^ (((s1 << 13) ^  s1) >> 19);
 105         s2 = ((s2 &  (uint32)-8) <<  4) ^ (((s2 <<  2) ^  s2) >> 25);
 106         s3 = ((s3 & (uint32)-16) << 17) ^ (((s3 <<  3) ^  s3) >> 11);
 107         return s1 ^ s2 ^ s3;
 108 }
 109
 110 inline uint32 RGen::trand()
 111 {
 112         // generate a random 32 bit number
 113         return ::trand(s1, s2, s3);
 114 }
 115
 116 inline double RGen::drand()
 117 {
 118         // return a double from 0.0 to 0.999...
 119 #if BYTE_ORDER == BIG_ENDIAN
 120         union { struct { uint32 hi, lo; } i; double f; } du;
 121 #else
 122         union { struct { uint32 lo, hi; } i; double f; } du;
 123 #endif
 124         du.i.hi = 0x41300000;
 125         du.i.lo = trand();
 126         return du.f - 1048576.;
 127 }
 128
 129 inline float RGen::frand()
 130 {
 131         // return a float from 0.0 to 0.999...
 132         union { uint32 i; float f; } u;         // union for floating point conversion of result
 133         u.i = 0x3F800000 | (trand() >> 9);
 134         return u.f - 1.f;
 135 }
 136
 137 inline float RGen::frand0()
 138 {
 139         // return a float from +1.0 to +1.999...
 140         union { uint32 i; float f; } u;         // union for floating point conversion of result
 141         u.i = 0x3F800000 | (trand() >> 9);
 142         return u.f;
 143 }
 144
 145 inline float RGen::frand2()
 146 {
 147         // return a float from -1.0 to +0.999...
 148         union { uint32 i; float f; } u;         // union for floating point conversion of result
 149         u.i = 0x40000000 | (trand() >> 9);
 150         return u.f - 3.f;
 151 }
 152
 153 inline float RGen::frand8()
 154 {
 155         // return a float from -0.125 to +0.124999...
 156         union { uint32 i; float f; } u;         // union for floating point conversion of result
 157         u.i = 0x3E800000 | (trand() >> 9);
 158         return u.f - 0.375f;
 159 }
 160
 161 inline float RGen::fcoin()
 162 {
 163         // only return one of the two values -1.0 or +1.0
 164         union { uint32 i; float f; } u;         // union for floating point conversion of result
 165         u.i = 0x3F800000 | (0x80000000 & trand());
 166         return u.f;
 167 }
 168
 169 inline int32 RGen::irand(int32 scale)
 170 {
 171         // return an int from 0 to scale - 1
 172         return (int32)floor(scale * drand());
 173 }
 174
 175 inline int32 RGen::irand2(int32 scale)
 176 {
 177         // return a int from -scale to +scale
 178         return (int32)floor((2. * scale + 1.) * drand() - scale);
 179 }
 180
 181 inline int32 RGen::ilinrand(int32 scale)
 182 {
 183         int32 a = irand(scale);
 184         int32 b = irand(scale);
 185         return sc_min(a,b);
 186 }
 187
 188 inline double RGen::linrand(double scale)
 189 {
 190         double a = drand();
 191         double b = drand();
 192         return sc_min(a,b) * scale;
 193 }
 194
 195 inline int32 RGen::ibilinrand(int32 scale)
 196 {
 197         int32 a = irand(scale);
 198         int32 b = irand(scale);
 199         return a - b;
 200 }
 201
 202 inline double RGen::bilinrand(double scale)
 203 {
 204         double a = drand();
 205         double b = drand();
 206         return (a - b) * scale;
 207 }
 208
 209 inline double RGen::exprandrng(double lo, double hi)
 210 {
 211         return lo * exp(log(hi / lo) * drand());
 212 }
 213
 214 inline double RGen::exprand(double scale)
 215 {
 216         double z;
 217         while ((z = drand()) == 0.0) {}
 218         return -log(z) * scale;
 219 }
 220
 221 inline double RGen::biexprand(double scale)
 222 {
 223         double z;
 224         while ((z = drand2(1.)) == 0.0 || z == -1.0) {}
 225         if (z > 0.0) z = log(z);
 226         else z = -log(-z);
 227         return z * scale;
 228 }
 229
 230 inline double RGen::sum3rand(double scale)
 231 {
 232         // larry polansky's poor man's gaussian generator
 233         return (drand() + drand() + drand() - 1.5) * 0.666666667 * scale;
 234 }
 235
 236 inline double drand( uint32& s1, uint32& s2, uint32& s3 )
 237 {
 238         union { struct { uint32 hi, lo; } i; double f; } u;
 239         u.i.hi = 0x41300000;
 240         u.i.lo = trand(s1,s2,s3);
 241         return u.f - 1048576.;
 242 }
 243
 244 inline float frand( uint32& s1, uint32& s2, uint32& s3 )
 245 {
 246         // return a float from 0.0 to 0.999...
 247         union { uint32 i; float f; } u;
 248         u.i = 0x3F800000 | (trand(s1,s2,s3) >> 9);
 249         return u.f - 1.f;
 250 }
 251
 252 inline float frand0( uint32& s1, uint32& s2, uint32& s3 )
 253 {
 254         // return a float from +1.0 to +1.999...
 255         union { uint32 i; float f; } u;
 256         u.i = 0x3F800000 | (trand(s1,s2,s3) >> 9);
 257         return u.f;
 258 }
 259
 260 inline float frand2( uint32& s1, uint32& s2, uint32& s3 )
 261 {
 262         // return a float from -1.0 to +0.999...
 263         union { uint32 i; float f; } u;
 264         u.i = 0x40000000 | (trand(s1,s2,s3) >> 9);
 265         return u.f - 3.f;
 266 }
 267
 268 inline float frand8( uint32& s1, uint32& s2, uint32& s3 )
 269 {
 270         // return a float from -0.125 to +0.124999...
 271         union { uint32 i; float f; } u;
 272         u.i = 0x3E800000 | (trand(s1,s2,s3) >> 9);
 273         return u.f - 0.375f;
 274 }
 275
 276 inline float fcoin( uint32& s1, uint32& s2, uint32& s3 )
 277 {
 278         // only return one of the two values -1.0 or +1.0
 279         union { uint32 i; float f; } u;
 280         u.i = 0x3F800000 | (0x80000000 & trand(s1,s2,s3));
 281         return u.f;
 282 }
 283
 284 #endif
 285