linrdist.cpp

#include <fstream>
#include "linrdist.h"

int static const RANTYPE=2;
int SIMSTART;

#define IA 16807
#define IM 2147483647
#define AM (1.0/IM)
#define IQ 127773
#define IR 2836
#define NTAB 32
#define NDIV (1+(IM-1)/NTAB)
#define EPS 1.2e-7
#define RNMX (1.0-EPS)


#define IM1 2147483563
#define IM2 2147483399
#define AMran2 (1.0/IM1)
#define IMM1 (IM1-1)
#define IA1 40014
#define IA2 40692
#define IQ1 53668
#define IQ2 52774
#define IR1 12211
#define IR2 3791
#define NTAB 32
#define NDIVran2 (1+IMM1/NTAB)
//#define EPS 1.2e-7
//#define RNMX (1.0-EPS)


#define MBIG 1000000000
#define MSEED 161803398
#define MZ 0
#define FAC (1.0/MBIG)


double ran1(long *idum)
{
  int j;
  long k;
  static long iy=0;
  static long iv[NTAB];
  double temp;

  if (*idum <= 0 || !iy) {
    if (-(*idum) < 1) *idum=1;
    else *idum = -(*idum);
    for (j=NTAB+7;j>=0;j--) {
      k=(*idum)/IQ;
      *idum=IA*(*idum-k*IQ)-IR*k;
      if (*idum < 0) *idum += IM;
      if (j < NTAB) iv[j] = *idum;
    }
    iy=iv[0];
  }
  k=(*idum)/IQ;
  *idum=IA*(*idum-k*IQ)-IR*k;
  if (*idum < 0) *idum += IM;
  j=iy/NDIV;
  iy=iv[j];
  iv[j] = *idum;
  if ((temp=AM*iy) > RNMX) return RNMX;
  else return temp;
}


double ran2(long *idum)
{
  int j;
  long k;
  static long idum2=123456789;
  static long iy=0;
  static long iv[NTAB];
  double temp;

  if (*idum <= 0) {
    if (-(*idum) < 1) *idum=1;
    else *idum = -(*idum);
    idum2=(*idum);
    for (j=NTAB+7;j>=0;j--) {
      k=(*idum)/IQ1;
      *idum=IA1*(*idum-k*IQ1)-k*IR1;
      if (*idum < 0) *idum += IM1;
      if (j < NTAB) iv[j] = *idum;
    }
    iy=iv[0];
  }
  k=(*idum)/IQ1;
  *idum=IA1*(*idum-k*IQ1)-k*IR1;
  if (*idum < 0) *idum += IM1;
  k=idum2/IQ2;
  idum2=IA2*(idum2-k*IQ2)-k*IR2;
  if (idum2 < 0) idum2 += IM2;
  j=iy/NDIVran2;
  iy=iv[j]-idum2;
  iv[j] = *idum;
  if (iy < 1) iy += IMM1;
  if ((temp=AMran2*iy) > RNMX) return RNMX;
  else return temp;
}



double ran3(long *idum)
{
  static int inext,inextp;
  static long ma[56];
  static int iff=0;
  long mj,mk;
  int i,ii,k;

  if (*idum < 0 || iff == 0) {
    iff=1;
    mj=MSEED-(*idum < 0 ? -*idum : *idum);
    mj %= MBIG;
    ma[55]=mj;
    mk=1;
    for (i=1;i<=54;i++) {
      ii=(21*i) % 55;
      ma[ii]=mk;
      mk=mj-mk;
      if (mk < MZ) mk += MBIG;
      mj=ma[ii];
    }
    for (k=1;k<=4;k++)
      for (i=1;i<=55;i++) {
        ma[i] -= ma[1+(i+30) % 55];
        if (ma[i] < MZ) ma[i] += MBIG;
      }
    inext=0;
    inextp=31;
    *idum=1;
  }
  if (++inext == 56) inext=1;
  if (++inextp == 56) inextp=1;
  mj=ma[inext]-ma[inextp];
  if (mj < MZ) mj += MBIG;
  ma[inext]=mj;
  return mj*FAC;
}


#undef IA
#undef IM
#undef AM
#undef IQ
#undef IR
#undef NTAB
#undef NDIV
#undef EPS
#undef RNMX

#undef IM1
#undef IM2
#undef AMran2
#undef IMM1
#undef IA1
#undef IA2
#undef IQ1
#undef IQ2
#undef IR1
#undef IR2
#undef NTAB
#undef NDIVran2
#undef EPS
#undef RNMX

#undef MBIG
#undef MSEED
#undef MZ
#undef FAC

double runif(long *idum)
{
  double ran;
  
  switch(RANTYPE)
    {
    case 1: ran = ran1(idum); break;
    case 2: ran = ran2(idum); break;
    case 3: ran = ran3(idum); break;
    default: std::cout << "Wrong input in runif (linrdist.h)"; break;
    }

  return(ran);
}


double rexp(long *idum)
{
  /* q[k-1] = sum(alog(2.0)**k/k!) k=1,..,n, */
  /* The highest n (here 8) is determined by q[n-1] = 1.0 */
  /* within standard precision */
  static double q[] =
    {
      0.6931471805599453,
      0.9333736875190459,
      0.9888777961838675,
      0.9984959252914960,
      0.9998292811061389,
      0.9999833164100727,
      0.9999985691438767,
      0.9999998906925558,
      0.9999999924734159,
      0.9999999995283275,
      0.9999999999728814,
      0.9999999999985598,
      0.9999999999999289,
      0.9999999999999968,
      0.9999999999999999,
      1.0000000000000000
    };
  double a, u, ustar, umin;
  int i;

  a = 0.0;
  u = runif(idum);
  for (;;) {
    u = u + u;
    if (u > 1.0)
      break;
    a = a + q[0];
  }
  u = u - 1.0;

  if (u <= q[0])
    return a + u;

  i = 0;
  ustar = runif(idum);
  umin = ustar;
  do {
    ustar = runif(idum);
    if (ustar < umin)
      umin = ustar;
    i = i + 1;
  }
  while (u > q[i]);
  return a + umin * q[0];
}


static double a[32] =
  {
    0.0000000, 0.03917609, 0.07841241, 0.1177699,
    0.1573107, 0.19709910, 0.23720210, 0.2776904,
    0.3186394, 0.36012990, 0.40225010, 0.4450965,
    0.4887764, 0.53340970, 0.57913220, 0.6260990,
    0.6744898, 0.72451440, 0.77642180, 0.8305109,
    0.8871466, 0.94678180, 1.00999000, 1.0775160,
    1.1503490, 1.22985900, 1.31801100, 1.4177970,
    1.5341210, 1.67594000, 1.86273200, 2.1538750
  };

static double d[31] =
  {
    0.0000000, 0.0000000, 0.0000000, 0.0000000,
    0.0000000, 0.2636843, 0.2425085, 0.2255674,
    0.2116342, 0.1999243, 0.1899108, 0.1812252,
    0.1736014, 0.1668419, 0.1607967, 0.1553497,
    0.1504094, 0.1459026, 0.1417700, 0.1379632,
    0.1344418, 0.1311722, 0.1281260, 0.1252791,
    0.1226109, 0.1201036, 0.1177417, 0.1155119,
    0.1134023, 0.1114027, 0.1095039
  };

static double t[31] =
  {
    7.673828e-4, 0.002306870, 0.003860618, 0.005438454,
    0.007050699, 0.008708396, 0.010423570, 0.012209530,
    0.014081250, 0.016055790, 0.018152900, 0.020395730,
    0.022811770, 0.025434070, 0.028302960, 0.031468220,
    0.034992330, 0.038954830, 0.043458780, 0.048640350,
    0.054683340, 0.061842220, 0.070479830, 0.081131950,
    0.094624440, 0.112300100, 0.136498000, 0.171688600,
    0.227624100, 0.330498000, 0.584703100
  };

static double h[31] =
  {
    0.03920617, 0.03932705, 0.03950999, 0.03975703,
    0.04007093, 0.04045533, 0.04091481, 0.04145507,
    0.04208311, 0.04280748, 0.04363863, 0.04458932,
    0.04567523, 0.04691571, 0.04833487, 0.04996298,
    0.05183859, 0.05401138, 0.05654656, 0.05953130,
    0.06308489, 0.06737503, 0.07264544, 0.07926471,
    0.08781922, 0.09930398, 0.11555990, 0.14043440,
    0.18361420, 0.27900160, 0.70104740
  };

#define repeat for(;;)

double rnorm1(long *idum)
{
  double s, u, w, y, ustar, aa, tt;
  int i;

  u = runif(idum);
  s = 0.0;
  if (u > 0.5)
    s = 1.0;
  u = u + u - s;
  u *= 32.0;
  i = (int) u;
  if (i == 32)
    i = 31;
  if (i != 0) {
    ustar = u - i;
    aa = a[i - 1];
    while (ustar <= t[i - 1]) {
      u = runif(idum);
      w = u * (a[i] - aa);
      tt = (w * 0.5 + aa) * w;
      repeat {
        if (ustar > tt)
          goto deliver;
        u = runif(idum);
        if (ustar < u)
          break;
        tt = u;
        ustar = runif(idum);
      }
      ustar = runif(idum);
    }
    w = (ustar - t[i - 1]) * h[i - 1];
  }
  else {
    i = 6;
    aa = a[31];
    repeat {
      u = u + u;
      if (u >= 1.0)
        break;
      aa = aa + d[i - 1];
      i = i + 1;
    }
    u = u - 1.0;
    repeat {
      w = u * d[i - 1];
      tt = (w * 0.5 + aa) * w;
      repeat {
        ustar = runif(idum);
        if (ustar > tt)
          goto jump;
        u = runif(idum);
        if (ustar < u)
          break;
        tt = u;
      }
      u = runif(idum);
    }
  jump:;
  }

 deliver:
  y = aa + w;
  return (s == 1.0) ? -y : y;

}

double rnorm2(long *idum)
{
  static int iset = 0;
  static double gset;
  double fac,rsq,v1,v2;

  if (SIMSTART == 1) {
    iset = 0;
    SIMSTART = 0;
  }

  if (iset == 0) {
    do {
      v1 = 2.0*runif(idum) - 1.0;
      v2 = 2.0*runif(idum) - 1.0;
      rsq = v1*v1 + v2*v2;
    } while (rsq >= 1.0 || rsq == 0.0);
    fac = sqrt(-2.0*log(rsq)/rsq);
    gset = v1*fac;
    iset = 1;
    return v2*fac;
  }
  else {
    iset = 0;
    return gset;
  }
}

double ranexp(long *idum)
{

  double dum;

  do
    dum = runif(idum);
  while (dum == 0.0);
  return -log(dum);
}

double rangamma(double a, long *idum) //a = shape parameter
{
  double e, u1, u2, x=0;
  
  e = exp(1.0);

  if (a < 1) {
    do {
      u1 = runif(idum);
      u2 = runif(idum);
      if (u1 < e/(a+e)) x = -log((a+e)*u1/e);
    } while (u1 < 1);
  }
  else if (a > 1) {

  }
  else  x = ranexp(idum);
  return x;
}
      

static double a1 = 0.3333333;
static double a2 = -0.250003;
static double a3 = 0.2000062;
static double a4 = -0.1662921;
static double a5 = 0.1423657;
static double a6 = -0.1367177;
static double a7 = 0.1233795;
static double e1 = 1.0;
static double e2 = 0.4999897;
static double e3 = 0.166829;
static double e4 = 0.0407753;
static double e5 = 0.010293;
static double q1 = 0.04166669;
static double q2 = 0.02083148;
static double q3 = 0.00801191;
static double q4 = 0.00144121;
static double q5 = -7.388e-5;
static double q6 = 2.4511e-4;
static double q7 = 2.424e-4;
static double sqrt32 = 5.656854;

static double aa = 0.;
static double aaa = 0.;

#define repeat for(;;)

double rgamma(double a, double scale, long *idum)
{
  static double b, c, d, e, p, q, r, s, t, u, v, w, x;
  static double q0, s2, si;
  double ret_val;

  if (a < 1.0) {
    /* alternate method for parameters a below 1 */
    /* 0.36787944117144232159 = exp(-1) */
    aa = 0.0;
    b = 1.0 + 0.36787944117144232159 * a;
    repeat {
      p = b * runif(idum);
      if (p >= 1.0) {
        ret_val = -log((b - p) / a);
        if (rexp(idum) >= (1.0 - a) * log(ret_val))
          break;
      } else {
        ret_val = exp(log(p) / a);
        if (rexp(idum) >= ret_val)
          break;
      }
    }
    return scale * ret_val;
  }
  /* Step 1: Recalculations of s2, s, d if a has changed */
  if (a != aa) {
    aa = a;
    s2 = a - 0.5;
    s = sqrt(s2);
    d = sqrt32 - s * 12.0;
  }
  /* Step 2: t = standard normal deviate, */
  /* x = (s,1/2)-normal deviate. */
  /* immediate acceptance (i) */

  t = rnorm2(idum);
  x = s + 0.5 * t;
  ret_val = x * x;
  if (t >= 0.0)
    return scale * ret_val;

  /* Step 3: u = 0,1 - uniform sample. squeeze acceptance (s) */
  u = runif(idum);
  if (d * u <= t * t * t) {
    return scale * ret_val;
  }
  /* Step 4: recalculations of q0, b, si, c if necessary */

  if (a != aaa) {
    aaa = a;
    r = 1.0 / a;
    q0 = ((((((q7 * r + q6) * r + q5) * r + q4)
            * r + q3) * r + q2) * r + q1) * r;

    /* Approximation depending on size of parameter a */
    /* The constants in the expressions for b, si and */
    /* c were established by numerical experiments */

    if (a <= 3.686) {
      b = 0.463 + s + 0.178 * s2;
      si = 1.235;
      c = 0.195 / s - 0.079 + 0.16 * s;
    } else if (a <= 13.022) {
      b = 1.654 + 0.0076 * s2;
      si = 1.68 / s + 0.275;
      c = 0.062 / s + 0.024;
    } else {
      b = 1.77;
      si = 0.75;
      c = 0.1515 / s;
    }
  }
  /* Step 5: no quotient test if x not positive */

  if (x > 0.0) {
    /* Step 6: calculation of v and quotient q */
    v = t / (s + s);
    if (fabs(v) <= 0.25)
      q = q0 + 0.5 * t * t * ((((((a7 * v + a6)
                                  * v + a5) * v + a4) * v + a3)
                               * v + a2) * v + a1) * v;
    else
      q = q0 - s * t + 0.25 * t * t + (s2 + s2)
        * log(1.0 + v);


    /* Step 7: quotient acceptance (q) */

    if (log(1.0 - u) <= q)
      return scale * ret_val;
  }
  /* Step 8: e = standard exponential deviate */
  /* u= 0,1 -uniform deviate */
  /* t=(b,si)-double exponential (laplace) sample */

  repeat {
    e = rexp(idum);
    u = runif(idum);
    u = u + u - 1.0;
    if (u < 0.0)
      t = b - si * e;
    else
      t = b + si * e;
    /* Step  9:  rejection if t < tau(1) = -0.71874483771719 */
    if (t >= -0.71874483771719) {
      /* Step 10:  calculation of v and quotient q */
      v = t / (s + s);
      if (fabs(v) <= 0.25)
        q = q0 + 0.5 * t * t * ((((((a7 * v + a6)
                                    * v + a5) * v + a4) * v + a3)
                                 * v + a2) * v + a1) * v;
      else
        q = q0 - s * t + 0.25 * t * t + (s2 + s2)
          * log(1.0 + v);
      /* Step 11:  hat acceptance (h) */
      /* (if q not positive go to step 8) */
      if (q > 0.0) {
        if (q <= 0.5)
          w = ((((e5 * q + e4) * q + e3)
                * q + e2) * q + e1) * q;
        else
          w = exp(q) - 1.0;
        /* if t is rejected */
        /* sample again at step 8 */
        if (c * fabs(u) <= w * exp(e - 0.5 * t * t))
          break;
      }
    }
  }
  x = s + 0.5 * t;
  return scale* x*x;

}


#ifdef TEST

#include <iostream>
#include <cstdlib>

int main() {
  try {

    
    
    long idum = 12345213;
    double ran1(long *idum);

    for (unsigned i=0; i<100; i++) {
      double val = ran1(&idum);
      std::cerr << val << std::endl;
    }
      

    

    std::cerr << "All tests in linrdist.cpp were successful.\n";
  } catch(std::string const& err) {
    std::cerr << err << std::endl;
    return EXIT_FAILURE;
  }
  return EXIT_SUCCESS;
}


#endif//TEST

---