/*	Math functions for Csound coded by Paris Smaragdis 1994         */
/*	Berklee College of Music Csound development team                */

#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;
}

void seedrand(PRAND *p)
{
    if ((unsigned int)*p->out == 0) {
      printf("Seeding from current time\n");
      srand((unsigned int)time(NULL));
    }
    else {
      printf("Seeding with %.3f\n", *p->out);
      srand((unsigned int)*p->out);
    }
}

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 * * * * * */


#define unirand()	(float)((double)rand()/(double)RAND_MAX)

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);
}