zufall.cc

模糊聚类分析的源程序！

/********************************************************/
/*   filename: zufall.cc                                */
/*                                                      */
/********************************************************/
/* programmed by: Oliver Wagner                         */
/* last change:  (XXX: not updated)			*/
/********************************************************/
#include "zufall.h"

unsigned int a[RANDOM_r], x_i[RANDOM_r], c;
time_t timer;

void 
init_rand ()
{
  int i;

  srand ((unsigned) time (&timer));

  for (i = 0; i < RANDOM_r; i++) {
    a[i] = rand () % USHRT_MAX;
    x_i[i] = rand () % USHRT_MAX;
  }
  c = rand () % USHRT_MAX;
}

int 
grand ()
{
  int i, j, Result;
  static random_i = 0;

  x_i[random_i] = 0;
  for (i = random_i + 1, j = 0; i != random_i; i = (i + 1) % RANDOM_r)
    x_i[random_i] += a[j++] * x_i[i];

  x_i[random_i] %= RANDOM_m;
  Result = x_i[random_i];
  random_i = (random_i + 1) % RANDOM_r;

  return (Result);
}

double 
nrand ()
{
  double u1, u2, u3, z /* Ergebnis */ ;

  u1 = ((double) grand ()) / (double) USHRT_MAX;
  if (u1 < 0.8638) {
    u1 = ((double) grand ()) / (double) USHRT_MAX;
    u2 = ((double) grand ()) / (double) USHRT_MAX;
    u3 = ((double) grand ()) / (double) USHRT_MAX;
    z = 2 * (u1 + u2 + u3 - 3 / 2);			// XXX: ???
  } else if (u1 < 0.9745) {
    u1 = ((double) grand ()) / (double) USHRT_MAX;
    u2 = ((double) grand ()) / (double) USHRT_MAX;
    z = (u1 + u2 - 1) * 3 / 2;			// XXX: ???
  } else if (u1 < 0.9973002039) {
    do {
      u1 = ((double) grand ()) / (double) USHRT_MAX;
      u2 = ((double) grand ()) / (double) USHRT_MAX;
      u1 = 6 * u1 - 3;
      u2 = 0.358 * u2;
    } while (u2 <= u1);
    z = u1;
  } else {
    do {
      u1 = ((double) grand ()) / (double) USHRT_MAX;
      u2 = ((double) grand ()) / (double) USHRT_MAX;
    } while (u2 <= 3 / sqrt (9 - 2 * log (u1)));
    z = sqrt (9 - 2 * log (u1));
  }

  return (z);
}

void 
nrand (double &z1, double &z2)
{
  double u1, u2, v1, v2, s;

/* Polarmethode */
  do {
    u1 = ((double) grand ()) / (double) USHRT_MAX;
    u2 = ((double) grand ()) / (double) USHRT_MAX;
    v1 = u1 + u1 - 1;
    v2 = u2 + u2 - 1;
    s = v1 * v1 + v2 * v2;
  } while (s >= 1);
  z1 = sqrt (-2 * log (s) / s);
  z2 = v2 * z1;
  z1 = v1 * z1;
}