This patch file will alter three files: entry.c cmath.c cmath.h to add a 31 bit internal Pseudo Random Number Generator for the "x-class noise" opcodes. It also adds ibunirand, kbunirand and abunirand which are bipolar uniformly distributed noise / random number generators. More details at http://www.firstpr.com.au/csound/ Robin Whittle 6 July 1998 --- entry.c-orig Thu Apr 23 01:02:00 1998 +++ entry.c Mon Jul 6 15:47:51 1998 @@ -194,6 +194,8 @@ void pvinterpset(void*), pvinterp(void*); void nlfiltset(void*), nlfilt(void*), Xsynth(void*), Xsynthset(void*); void auniform(void*), ikuniform(void*); +void abuniform(void*), ikbuniform(void*); + void tblsetw(void*), ktablew(void*), tablew(void*), itablew(void*); void tblsetwkt(void*), ktablewkt(void*), tablewkt(void*); void tableng(void*), itableng(void*), tablegpw(void*), itablegpw(void*); @@ -662,6 +664,11 @@ { "iunirand", S(PRAND), 1, "i", "k", ikuniform, NULL, NULL }, { "kunirand", S(PRAND), 2, "k", "k", NULL, ikuniform, NULL}, { "aunirand", S(PRAND), 4, "a", "k", NULL, NULL, auniform }, + +{ "ibunirand", S(PRAND),1, "i", "k", ikbuniform, NULL, NULL }, +{ "kbunirand", S(PRAND),2, "k", "k", NULL, ikbuniform, NULL}, +{ "abunirand", S(PRAND),4, "a", "k", NULL, NULL, abuniform }, + { "diskin",S(SOUNDINEW),5, "mmmm", "Skooo", newsndinset,NULL, soundinew}, { "ion", S(OUT_ON), 1, "", "iii", iout_on, NULL, NULL }, { "ioff", S(OUT_ON), 1, "", "iii", iout_off, NULL, NULL }, --- cmath.c-orig Sat Feb 28 02:51:00 1998 +++ cmath.c Mon Jul 6 15:55:22 1998 @@ -1,6 +1,24 @@ /* Math functions for Csound coded by Paris Smaragdis 1994 */ /* Berklee College of Music Csound development team */ +/* Changes 6 July 1998 - Robin Whittle http://www.firstpr.com.au + * + * 31 bit Park Miller "Minimal Standard" Pseudo Random Number Generator + * now used for all these noise sources, with a seed setting function and a + * raw 31 bit function available for other Csound functions too. + * + * This removes the dependency on the library function rand() which could + * be as simple as a 16 bit PRNG, which would be lousy for audio purposes. + * + * Changes made to seedrand(). + * A new function randint31() is created. + * Macro "unirand()" is now defined in terms of randint31() instead of + * the library call rand(). + * New functions for bipolar uniform noise -1 to +1. Three new lines + * for this xbunirand family need to be added to entry.c, as well as + * declaring the new function names: + */ + #include "cs.h" #include "cmath.h" #include @@ -30,18 +48,151 @@ *p->sr = (float)pow(*p->in, *p->pow) / *p->norm; } +/*========================================================================= + * + * seedrand() and the PRNG seed. + * + * Driven by opcode "seed". Takes a float input parm and seeds the + * 31 bit state of the PRNG as follows: + * + * Input = 0 Use time of day to generate seed. + * + * Input != 0 Make it positive, multiply by 1,000, mod it by + * (2^31 - 1) turn it into an integer using floor, + * so it can never be 2^31, make it 1 if the result + * was 0, and use this as the seed. + * The PRNG will put out 0s if its state is 0! + * + * The seed is a 32 bit integer. Initialise it with an + * auspicious number - pi to the 17th power. + * + * The seed opcode uses the same data structure PRAND as the generators + * - its sole parameter, an input, is in the first position, which is + * called "out". + */ + +#define CS_RANDINT_MAX 2147483647L /* 2**31 - 1 = 0x7FFFFFFF */ + + /* Add one in case y2k voodoo causes time() + * to return 0. + */ +static long seed31 = 282844563L; + void seedrand(PRAND *p) { - if ((unsigned int)*p->out == 0) { - printf("Seeding from current time\n"); - srand((unsigned int)time(NULL)); + /* Make sure that the time of day function is returning + * non-zero. + */ + if ( (*p->out == 0) && (time(NULL) != 0) ) + { + printf("0 -> Setting ramdom seed from current time and date. "); + + /* Since the time changes rather slowly, and the + * first generated randint31() result depends + * on it and therefore changes rather slowly + * call randint31() twice the mash the seed + * well before an opcode can use it. + */ + + seed31 = (long) time(NULL); + seed31 = randint31(); + seed31 = randint31(); } - else { - printf("Seeding with %.3f\n", *p->out); - srand((unsigned int)*p->out); + else + { + seed31 = (long) floor( fmod( (double)(1000 * fabs(*p->out)), (double)2147483647 ) ); + printf("Seeding with %.3f ", *p->out); } + + if (seed31 == 0) seed31 = 1; + printf("Seed = %10lu \n", seed31) ; } + +/*========================================================================= + * + * randint31() + * + * 31 bit Park Miller PRNG using Linus Schrage's method for doing it all + * with 32 bit variables. + * + * Code adapted from Ray Garder's public domain code of July 1997 at: + * http://www.snippets.org/RG_RAND.C Thanks! + * + * Based on "Random Number Generators: Good Ones Are Hard to Find", + * S.K. Park and K.W. Miller, Communications of the ACM 31:10 (Oct 1988), + * and "Two Fast Implementations of the 'Minimal Standard' Random + * Number Generator", David G. Carta, Comm. ACM 33, 1 (Jan 1990), p. 87-88 + * + * Linear congruential generator: f(z) = (16807 * z) mod (2 ** 31 - 1) + * + * Uses L. Schrage's method to avoid overflow problems. + */ + +#define ria 16807 /* multiplier */ +#define rim 2147483647L /* 2**31 - 1 */ +#define riq 127773L /* m div a */ +#define rir 2836 /* m mod a */ + +long randint31() +{ + unsigned long rilo, rihi; + + rilo = ria * (long)(seed31 & 0xFFFF); + rihi = ria * (long)((unsigned long)seed31 >> 16); + rilo += (rihi & 0x7FFF) << 16; + if (rilo > rim) + { + rilo &= rim; + ++rilo; + } + rilo += rihi >> 15; + if (rilo > rim) + { + rilo &= rim; + ++rilo; + } + return ( seed31 = (long)rilo ); +} + + +/*========================================================================= + * + * unirand() macro to call the above PRNG function. + */ + +#define unirand() (float)((double)randint31()/(double)CS_RANDINT_MAX) + + +/*========================================================================= + * + * Funcions for xbunirand opcodes in entry.c. + * Based on auinform() and ikuniform() below. + * unirand was 0 to +1. bunirand is -1 to +1 + */ + +void abuniform(PRAND *p) /* Uniform bipolar distribution */ +{ + long n = ksmps; + float *out = p->out; + float arg1 = *p->arg1; + + do *out++ = arg1 * (unifrand(2) - 1); + while( --n); +} + +void ikbuniform(PRAND *p) +{ + float *out = p->out; + + *out = *p->arg1 * (unifrand(2) - 1); +} + + +/* Back to Paris' code. - - - - - - - - - - - - - - - - - - - - - - - - - - */ + + + void auniform(PRAND *p) /* Uniform distribution */ { long n = ksmps; @@ -223,8 +374,13 @@ /* * * * * * RANDOM NUMBER GENERATORS * * * * * */ -#define unirand() (float)((double)rand()/(double)RAND_MAX) - +/* This was the old macro which called the library function rand(). + * + * #define unirand() (float)((double)rand()/(double)RAND_MAX) + * + * See above for the new unirand() macro. + */ + float unifrand(float range) { return range*unirand(); --- cmath.h-orig Mon Jul 22 17:50:11 1996 +++ cmath.h Mon Jul 6 00:41:11 1998 @@ -1,3 +1,12 @@ +/* cmath.h + * + * Originally by Paris Smaragdis 1994. + * + * rantint31() added by Robin Whittle 6 July 1998. + */ + +long randint31(void); + float besseli(double); float unifrand(float); float linrand(float);