Logo Search packages:      
Sourcecode: csound version File versions

cmath.c

/*
    cmath.c:

    Copyright (C) 1994 Paris Smaragdis, John ffitch

    This file is part of Csound.

    The Csound Library is free software; you can redistribute it
    and/or modify it under the terms of the GNU Lesser General Public
    License as published by the Free Software Foundation; either
    version 2.1 of the License, or (at your option) any later version.

    Csound is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU Lesser General Public License for more details.

    You should have received a copy of the GNU Lesser General Public
    License along with Csound; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
    02111-1307 USA
*/

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

#include "csoundCore.h"
#include "cmath.h"
#include <math.h>

int ipow(CSOUND *csound, POW *p)        /*      Power for i-rate */
{
    MYFLT in = *p->in;
    MYFLT powerOf = *p->powerOf;
    if (UNLIKELY(in == FL(0.0) && powerOf == FL(0.0)))
      return csound->PerfError(csound, Str("NaN in pow\n"));
    else
      *p->sr = POWER(in, powerOf) / *p->norm;
    return OK;
}

int apow(CSOUND *csound, POW *p)        /* Power routine for a-rate  */
{
    int   n, nsmps = csound->ksmps;
    MYFLT *in = p->in, *out = p->sr;
    MYFLT powerOf = *p->powerOf;

    if (UNLIKELY(powerOf == FL(0.0))) {
      MYFLT yy = FL(1.0) / *p->norm;
      for (n = 0; n < nsmps; n++) {
        MYFLT xx = in[n];
        if (UNLIKELY(xx == FL(0.0))) {
          return csound->PerfError(csound, Str("NaN in pow\n"));
        }
        else
          out[n] = yy;
      }
    }
    else {
      for (n = 0; n < nsmps; n++)
        out[n] = POWER(in[n], powerOf) / *p->norm;
    }
    return OK;
}

int seedrand(CSOUND *csound, PRAND *p)
{
    uint32_t  seedVal = (uint32_t)0;

    if (*p->out > FL(0.0))
      seedVal = (uint32_t)((double)*p->out + 0.5);
    else if (!seedVal) {
      seedVal = (uint32_t)csound->GetRandomSeedFromTime();
      csound->Message(csound, Str("Seeding from current time %u\n"),
                              (unsigned int)seedVal);
    }
    else
      csound->Message(csound, Str("Seeding with %u\n"), (unsigned int)seedVal);
    csound->SeedRandMT(&(csound->randState_), NULL, seedVal);
    csound->holdrand = (int)(seedVal & (uint32_t) 0x7FFFFFFF);
    while (seedVal >= (uint32_t)0x7FFFFFFE)
      seedVal -= (uint32_t)0x7FFFFFFE;
    csound->randSeed1 = ((int)seedVal + 1);

    return OK;
}

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

#define UInt32toFlt(x) ((double)(x) * (1.0 / 4294967295.03125))

#define unirand(c) ((MYFLT) UInt32toFlt(csoundRandMT(&((c)->randState_))))

static inline MYFLT unifrand(CSOUND *csound, MYFLT range)
{
    return (range * unirand(csound));
}

/* linear distribution routine */

static inline MYFLT linrand(CSOUND *csound, MYFLT range)
{
    uint32_t  r1, r2;

    r1 = csoundRandMT(&(csound->randState_));
    r2 = csoundRandMT(&(csound->randState_));

    return ((MYFLT)UInt32toFlt(r1 < r2 ? r1 : r2) * range);
}

/* triangle distribution routine */

static inline MYFLT trirand(CSOUND *csound, MYFLT range)
{
    uint64_t  r1;

    r1 = (uint64_t)csoundRandMT(&(csound->randState_));
    r1 += (uint64_t)csoundRandMT(&(csound->randState_));

    return ((MYFLT) ((double)((int64_t)r1 - (int64_t)0xFFFFFFFFU)
                     * (1.0 / 4294967295.03125)) * range);
}

/* exponential distribution routine */

static MYFLT exprand(CSOUND *csound, MYFLT l)
{
    uint32_t  r1;

    if (l < FL(0.0)) return (FL(0.0)); /* for safety */

    do {
      r1 = csoundRandMT(&(csound->randState_));
    } while (!r1);

    return -((MYFLT)log(UInt32toFlt(r1)) * l);
}

/* bilateral exponential distribution routine */

static MYFLT biexprand(CSOUND *csound, MYFLT l)
{
    int32_t r1;

    if (l < FL(0.0)) return (FL(0.0)); /* For safety */

    do {
      r1 = (int32_t)csoundRandMT(&(csound->randState_));
    } while (!r1);

    if (r1 < (int32_t)0) {
      return -(LOG(-(r1) * (FL(1.0) / FL(2147483648.0))) * l);
    }
    return (LOG(r1 * (FL(1.0) / FL(2147483648.0))) * l);
}

/* gaussian distribution routine */

static MYFLT gaussrand(CSOUND *csound, MYFLT s)
{
    int64_t   r1 = -((int64_t)0xFFFFFFFFU * 6);
    int       n = 12;
    double    x;

    do {
      r1 += (int64_t)csoundRandMT(&(csound->randState_));
    } while (--n);
    x = (double)r1;
    return (MYFLT)(x * ((double)s * (1.0 / (3.83 * 4294967295.03125))));
}

/* cauchy distribution routine */

static MYFLT cauchrand(CSOUND *csound, MYFLT a)
{
    uint32_t  r1;
    MYFLT     x;

    do {
      r1 = csoundRandMT(&(csound->randState_)); /* Limit range artificially */
    } while (r1 > (uint32_t)2143188560U && r1 < (uint32_t)2151778735U);
    x = TAN((MYFLT)r1 * (PI_F / FL(4294967295.0))) * (FL(1.0) / FL(318.3));
    return (x * a);
}

/* positive cauchy distribution routine */

static MYFLT pcauchrand(CSOUND *csound, MYFLT a)
{
    uint32_t  r1;
    MYFLT     x;

    do {
      r1 = csoundRandMT(&(csound->randState_));
    } while (r1 > (uint32_t)4286377121U);      /* Limit range artificially */
    x = TAN((MYFLT)r1 * (PI_F * FL(0.5) / FL(4294967295.0))) * (FL(1.0) / FL(318.3));
    return (x * a);
}

/* beta distribution routine */

static MYFLT betarand(CSOUND *csound, MYFLT range, MYFLT a, MYFLT b)
{
    double  r1, r2;

    if (a <= FL(0.0) || b <= FL(0.0))
      return FL(0.0);

    do {
      uint32_t  tmp;
      do {
        tmp = csoundRandMT(&(csound->randState_));
      } while (!tmp);
      r1 = pow(UInt32toFlt(tmp), 1.0 / (double)a);
      do {
        tmp = csoundRandMT(&(csound->randState_));
      } while (!tmp);
      r2 = r1 + pow(UInt32toFlt(tmp), 1.0 / (double)b);
    } while (r2 > 1.0);

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

/* weibull distribution routine */

static MYFLT weibrand(CSOUND *csound, MYFLT s, MYFLT t)
{
    uint32_t  r1;
    double    r2;

    if (t <= FL(0.0)) return FL(0.0);

    do {
      r1 = csoundRandMT(&(csound->randState_));
    } while (!r1 || r1 == (uint32_t)0xFFFFFFFFU);

    r2 = 1.0 - ((double)r1 * (1.0 / 4294967295.0));

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

/* Poisson distribution routine */

static MYFLT poissrand(CSOUND *csound, MYFLT l)
{
    MYFLT r1, r2, r3;

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

    r1 = unirand(csound);
    r2 = EXP(-l);
    r3 = FL(0.0);

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

    return (r3);
}

 /* ------------------------------------------------------------------------ */

int auniform(CSOUND *csound, PRAND *p)  /* Uniform distribution */
{
    MYFLT   *out = p->out, *endp = p->out + csound->ksmps;
    double  scale = (double)*p->arg1 * (1.0 / 4294967295.03125);

    do {
      *(out++) = (MYFLT)((double)csoundRandMT(&(csound->randState_)) * scale);
    } while (out < endp);
    return OK;
}

int ikuniform(CSOUND *csound, PRAND *p)
{
    *p->out = unifrand(csound, *p->arg1);
    return OK;
}

int alinear(CSOUND *csound, PRAND *p)   /* Linear random functions      */
{
    int   n, nsmps = csound->ksmps;
    MYFLT *out = p->out;
    MYFLT arg1 = *p->arg1;

    for (n = 0; n < nsmps; n++)
      out[n] = linrand(csound, arg1);
    return OK;
}

int iklinear(CSOUND *csound, PRAND *p)
{
    *p->out = linrand(csound, *p->arg1);
    return OK;
}

int atrian(CSOUND *csound, PRAND *p)    /* Triangle random functions  */
{
    int   n, nsmps = csound->ksmps;
    MYFLT *out = p->out;
    MYFLT arg1 = *p->arg1;

    for (n = 0; n < nsmps; n++)
      out[n] = trirand(csound, arg1);
    return OK;
}

int iktrian(CSOUND *csound, PRAND *p)
{
    *p->out = trirand(csound, *p->arg1);
    return OK;
}

int aexp(CSOUND *csound, PRAND *p)      /* Exponential random functions */
{
    int   n, nsmps = csound->ksmps;
    MYFLT *out = p->out;
    MYFLT arg1 = *p->arg1;

    for (n = 0; n < nsmps; n++)
      out[n] = exprand(csound, arg1);
    return OK;
}

int ikexp(CSOUND *csound, PRAND *p)
{
    *p->out = exprand(csound, *p->arg1);
    return OK;
}

int abiexp(CSOUND *csound, PRAND *p)    /* Bilateral exponential rand */
{                                       /* functions */
    int   n, nsmps = csound->ksmps;
    MYFLT *out = p->out;
    MYFLT arg1 = *p->arg1;

    for (n = 0; n < nsmps; n++)
      out[n] = biexprand(csound, arg1);
    return OK;
}

int ikbiexp(CSOUND *csound, PRAND *p)
{
    *p->out = biexprand(csound, *p->arg1);
    return OK;
}

int agaus(CSOUND *csound, PRAND *p)     /* Gaussian random functions */
{
    int   n, nsmps = csound->ksmps;
    MYFLT *out = p->out;
    MYFLT arg1 = *p->arg1;

    for (n = 0; n < nsmps; n++)
      out[n] = gaussrand(csound, arg1);
    return OK;
}

int ikgaus(CSOUND *csound, PRAND *p)
{
    *p->out = gaussrand(csound, *p->arg1);
    return OK;
}

int acauchy(CSOUND *csound, PRAND *p)   /* Cauchy random functions */
{
    int   n, nsmps = csound->ksmps;
    MYFLT *out = p->out;
    MYFLT arg1 = *p->arg1;

    for (n = 0; n < nsmps; n++)
      out[n] = cauchrand(csound, arg1);
    return OK;
}

int ikcauchy(CSOUND *csound, PRAND *p)
{
    *p->out = cauchrand(csound, *p->arg1);
    return OK;
}

int apcauchy(CSOUND *csound, PRAND *p)  /* +ve Cauchy random functions */
{
    int   n, nsmps = csound->ksmps;
    MYFLT *out = p->out;
    MYFLT arg1 = *p->arg1;

    for (n = 0; n < nsmps; n++)
      out[n] = pcauchrand(csound, arg1);
    return OK;
}

int ikpcauchy(CSOUND *csound, PRAND *p)
{
    *p->out = pcauchrand(csound, *p->arg1);
    return OK;
}

int abeta(CSOUND *csound, PRAND *p)     /* Beta random functions   */
{
    int   n, nsmps = csound->ksmps;
    MYFLT *out = p->out;
    MYFLT arg1 = *p->arg1;
    MYFLT arg2 = *p->arg2;
    MYFLT arg3 = *p->arg3;

    for (n = 0; n < nsmps; n++)
      out[n] = betarand(csound, arg1, arg2, arg3);
    return OK;
}

int ikbeta(CSOUND *csound, PRAND *p)
{
    *p->out = betarand(csound, *p->arg1, *p->arg2, *p->arg3);
    return OK;
}

int aweib(CSOUND *csound, PRAND *p)     /* Weibull randon functions */
{
    int   n, nsmps = csound->ksmps;
    MYFLT *out = p->out;
    MYFLT arg1 = *p->arg1;
    MYFLT arg2 = *p->arg2;

    for (n = 0; n < nsmps; n++)
      out[n] = weibrand(csound, arg1, arg2);
    return OK;
}

int ikweib(CSOUND *csound, PRAND *p)
{
    *p->out = weibrand(csound, *p->arg1, *p->arg2);
    return OK;
}

int apoiss(CSOUND *csound, PRAND *p)    /*      Poisson random funcions */
{
    int   n, nsmps = csound->ksmps;
    MYFLT *out = p->out;
    MYFLT arg1 = *p->arg1;

    for (n = 0; n < nsmps; n++)
      out[n] = poissrand(csound, arg1);
    return OK;
}

int ikpoiss(CSOUND *csound, PRAND *p)
{
    *p->out = poissrand(csound, *p->arg1);
    return OK;
}

 /* ------------------------------------------------------------------------ */

int gen21_rand(FGDATA *ff, FUNC *ftp)
{
    CSOUND  *csound = ff->csound;
    int32    i, n;
    MYFLT   *ft;
    MYFLT   scale;
    int     nargs = ff->e.pcnt - 4;

    ft = ftp->ftable;
    scale = (nargs > 1 ? ff->e.p[6] : FL(1.0));
    n = ff->flen;
    if (ff->guardreq)
      n++;
    switch ((int) ff->e.p[5]) {
    case 1:                     /* Uniform distribution */
      for (i = 0 ; i < n ; i++)
        ft[i] = unifrand(csound, scale);
      break;
    case 2:                     /* Linear distribution */
      for (i = 0 ; i < n ; i++)
        ft[i] = linrand(csound, scale);
      break;
    case 3:                     /* Triangular about 0.5 */
      for (i = 0 ; i < n ; i++)
        ft[i] = trirand(csound, scale);
      break;
    case 4:                     /* Exponential */
      for (i = 0 ; i < n ; i++)
        ft[i] = exprand(csound, scale);
      break;
    case 5:                     /* Bilateral exponential */
      for (i = 0 ; i < n ; i++)
        ft[i] = biexprand(csound, scale);
      break;
    case 6:                     /* Gaussian distribution */
      for (i = 0 ; i < n ; i++)
        ft[i] = gaussrand(csound, scale);
      break;
    case 7:                     /* Cauchy distribution */
      for (i = 0 ; i < n ; i++)
        ft[i] = cauchrand(csound, scale);
      break;
    case 8:                     /* Positive Cauchy */
      for (i = 0 ; i < n ; i++)
        ft[i] = pcauchrand(csound, scale);
      break;
    case 9:                     /* Beta distribution */
      if (nargs < 3) {
        return -1;
      }
      for (i = 0 ; i < n ; i++)
        ft[i] = betarand(csound, scale, (MYFLT) ff->e.p[7], (MYFLT) ff->e.p[8]);
      break;
    case 10:                    /* Weibull Distribution */
      if (nargs < 2) {
        return -1;
      }
      for (i = 0 ; i < n ; i++)
        ft[i] = weibrand(csound, scale, (MYFLT) ff->e.p[7]);
      break;
    case 11:                    /* Poisson Distribution */
      for (i = 0 ; i < n ; i++)
        ft[i] = poissrand(csound, scale);
      break;
    default:
      return -2;
    }
    return OK;
}


Generated by  Doxygen 1.6.0   Back to index