/*	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 <math.h>
#include <time.h>   /* for random seeding from the time - mkc July 1997 */

#ifndef RAND_MAX
#define RAND_MAX        (32767)
#endif

void ipow(POW *p)		/*      Power for i-rate */
{
    *p->sr = (float)pow(*p->in, *p->pow);	/*      sophisticated code!     */
}

void apow(POW *p)		/*      Power routine for a-rate  */
{
    long n = ksmps;
    float *in = p->in, *out = p->sr;

    do {
	*out++ = (float)pow(*in++, *p->pow) / *p->norm;
    } while (--n);
}

void kpow(POW *p)		/*      Power routine for k-rate        */
{
    *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)
{
			/* 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 
    {
	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;
    float *out = p->out;
    float arg1 = *p->arg1;

    do *out++ = unifrand(arg1);
    while( --n);
}

void ikuniform(PRAND *p)
{
    float *out = p->out;

    *out = unifrand(*p->arg1);
}

void alinear(PRAND *p)		/* Linear random functions      */
{
    long n = ksmps;
    float *out = p->out;

    do {
	*out++ = linrand(*p->arg1);
    } while (--n);
}

void iklinear(PRAND *p)
{
    float *out = p->out;
    *out++ = linrand(*p->arg1);
}

void atrian(PRAND *p)		/*      Triangle random functions       */
{
    long n = ksmps;
    float *out = p->out;

    do {
	*out++ = trirand(*p->arg1);
    } while (--n);
}

void iktrian(PRAND *p)
{
    float *out = p->out;
    *out++ = trirand(*p->arg1);
}

void aexp(PRAND *p)		/*      Exponential random functions    */
{
    long n = ksmps;
    float *out = p->out;

    do {
	*out++ = exprand(*p->arg1);
    } while (--n);
}

void ikexp(PRAND *p)
{
    float *out = p->out;
    *out++ = exprand(*p->arg1);
}

void abiexp(PRAND *p)		/*      Bilateral exponential rand. functions */
{
    long n = ksmps;
    float *out = p->out;

    do {
	*out++ = biexprand(*p->arg1);
    } while (--n);
}

void ikbiexp(PRAND *p)
{
    float *out = p->out;
    *out++ = biexprand(*p->arg1);
}

void agaus(PRAND *p)		/*      Gaussian random functions       */
{
    long n = ksmps;
    float *out = p->out;

    do {
	*out++ = gaussrand(*p->arg1);
    } while (--n);
}

void ikgaus(PRAND *p)
{
    float *out = p->out;
    *out++ = gaussrand(*p->arg1);
}

void acauchy(PRAND *p)		/*      Cauchy random functions */
{
    long n = ksmps;
    float *out = p->out;

    do {
	*out++ = cauchrand(*p->arg1);
    } while (--n);
}

void ikcauchy(PRAND *p)
{
    float *out = p->out;
    *out++ = cauchrand(*p->arg1);
}

void apcauchy(PRAND *p)		/*      Positive Cauchy random functions */
{
    long n = ksmps;
    float *out = p->out;

    do {
	*out++ = pcauchrand(*p->arg1);
    } while (--n);
}

void ikpcauchy(PRAND *p)
{
    float *out = p->out;
    *out++ = pcauchrand(*p->arg1);
}

void abeta(PRAND *p)		/*      Beta random functions   */
{
    long n = ksmps;
    float *out = p->out;

    do {
	*out++ = betarand(*p->arg1, *p->arg2, *p->arg3);
    } while (--n);
}

void ikbeta(PRAND *p)
{
    float *out = p->out;
    *out++ = betarand(*p->arg1, *p->arg2, *p->arg3);
}

void aweib(PRAND *p)		/*      Weibull randon functions        */
{
    long n = ksmps;
    float *out = p->out;

    do {
	*out++ = weibrand(*p->arg1, *p->arg2);
    } while (--n);
}

void ikweib(PRAND *p)
{
    float *out = p->out;
    *out++ = weibrand(*p->arg1, *p->arg2);
}

void apoiss(PRAND *p)		/*      Poisson random funcions */
{
    long n = ksmps;
    float *out = p->out;

    do {
	*out++ = poissrand(*p->arg1);
    } while (--n);
}

void ikpoiss(PRAND *p)
{
    float *out = p->out;
    *out++ = poissrand(*p->arg1);
}


/* * * * * * RANDOM NUMBER GENERATORS * * * * * */


/* 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();
}

float linrand(float range)	/*      linear distribution routine     */
{
    float r1, r2;

    r1 = unirand();
    r2 = unirand();

    if (r1 > r2)
	r1 = r2;

    return (r1 * range);
}

float trirand(float range)	/*      triangle distribution routine   */
{
    float r1, r2;

    r1 = unirand();
    r2 = unirand();

    return (((r1 + r2) - 1.0f) * range);
}
		
float exprand(float l)		/*      exponential distribution routine */
{
    float r1;

    if (l < 0.0f) return (0.0f);

    do {
	r1 = unirand();
    } while (r1 == 0.0f);

    return (-(float)log(r1 *l));
}

float biexprand(float l)	/* bilateral exponential distribution routine */
{
    float r1;

    if (l < 0.0f) return (0.0f);

    do {
	r1 = 2.0f * unirand();
    } while (r1 == 0.0f || r1 == 2.0f); 

    if (r1 > 1.0f)     {
	r1 = 2.0f - r1;
	return (-(float)log(r1 * l));
    }
    return ((float)log(r1 * l));
}

float gaussrand(float s)	/*      gaussian distribution routine   */
{
    float r1 = 0.0f;
    int n = 12;
    s /= 3.83f;

    do {
	r1 += unirand();
    } while (--n);

    return (s * (r1 - 6.0f));
}

float cauchrand(float a)	/*      cauchy distribution routine     */
{
    float r1;
    a /= 318.3f;

    do {
      do {
	r1 = unirand();
      } while (r1 == 0.5f);

      r1 = a * (float)tan(pi*(double)r1);
    } while (r1>1.0f || r1<-1.0f); /* Limit range artificially */
    return r1;
}

float pcauchrand(float a)	/*      positive cauchy distribution routine */
{
    float r1;
    a /= 318.3f;

    do {
      do {
	r1 = unirand();
      } while (r1 == 1);

      r1 = a * (float)tan( pi * 0.5 * r1);
    } while (r1>1.0f);
    return r1;
}

float betarand(float range, float a, float b) /* beta distribution routine  */
{
    float r1, r2;

    if (a < 0.0f || b < 0.0f ) return (0.0f); 

    do {
	do {
	    r1 = unirand();
	} while (r1 == 0.0f);
	
	do {
	    r2 = unirand();
	} while (r2 == 0.0f);
	
	r1 = (float)pow(r1, 1.0 / (double)a);
	r2 = (float)pow(r2, 1.0 / (double)b);
    } while ((r1 + r2) > 1.0f);

    return ((r1 / (r1 +r2)) * range);
}

float weibrand(float s, float t) /*      weibull distribution routine    */
{
    float r1, r2;

    if (t < 0.0f ) return (0.0f);

    do {
	r1 = unirand();
    } while (r1 == 0.0f || r1 == 1.0f);

    r2 = 1.0f /  (1.0f - r1);

    return (s * (float)pow (log((double)r2),  (1.0 /(double)t)));
}

float poissrand(float l)	/*      poisson distribution routine    */
{
    float r1, r2, r3;

    if (l < 0.0f ) return (0.0f);

    r1 = unirand();
    r2 = (float)exp(-l);
    r3 = 0.0f;

    while (r1 >= r2) {
	r3++;
	r1 *= unirand();
    }

    return (r3);
}