randomNumbers_Mseries.cpp

Go to the documentation of this file.

#include "randomNumbers_Mseries.h"

using std::string;

//====================================================================
void RandomNumbers_Mseries::set_readfile(string filename)
{
  vout.general(m_vl, "RandomNumbers_Mseries: read from file = %s\n", filename.c_str());

  std::ifstream infile;
  infile.open(filename.c_str());

  if (!infile) {
    vout.crucial(m_vl, "RandomNumbers_Mseries: unable to open input file.\n");
    return;
  }

  infile >> jr >> kr;

  for (int i = 0; i < Np; i++) {
    infile >> w[i];
  }

  infile.close();
}


//====================================================================
void RandomNumbers_Mseries::initset(const int ndelay)
{
  vout.general(m_vl, " Random number init: ndelay =%10d\n", ndelay);

  int IP = Np;
  int IQ = Nq;

  int      MACRM = 40;
  int      MACRI = 1;
  unsigned IB[IP];

  unsigned IX = MACRI;
  for (int i = 0; i < IP; ++i) {
    IX    = IX * 69069;
    IB[i] = IX >> 31;
  }

  int JR = 0;
  int KR = IP - IQ;
  for (int j = 0; j < IP; ++j) {
    unsigned iwork = 0;

    for (int i = 0; i < 32; ++i) {
      iwork  = iwork * 2 + IB[JR];
      IB[JR] = IB[JR] ^ IB[KR];
      ++JR;
      if (JR == IP) JR = 0;
      ++KR;
      if (KR == IP) KR = 0;
    }

    w[j] = iwork >> 1;
  }

  jr = JR;
  kr = KR;

  delay3(ndelay);
}


//====================================================================
void RandomNumbers_Mseries::delay3(const int ndelay)
{
  int IP = Np;
  int IQ = Nq;

  int LAMBDA = ndelay;

  int MACRM = 40;
  int MACRI = 1;

  unsigned IWK[2 * IP - 1];
  unsigned C[2 * IP];
  unsigned IB[IP + 32];

  int MU = MACRM;

  for (int i = 0; i < IP; ++i) {
    IWK[i] = w[i];
  }

  for (int i = IP; i < 2 * IP - 1; ++i) {
    IWK[i] = IWK[i - IP] ^ IWK[i - IQ];
  }

  for (int i = 0; i < MU; ++i) {
    IB[i] = 0;
  }
  int M  = LAMBDA;
  int NB = MU - 1;

continued220:
  if (M <= IP - 1) goto continued300;
  ++NB;
  IB[NB] = M % 2;
  M      = M / 2;
  goto continued220;

continued300:
  for (int i = 0; i < IP; ++i) {
    C[i] = 0;
  }

  C[M] = 1;
  for (int j = NB; j >= 0; --j) {
    for (int i = IP - 1; i >= 0; --i) {
      C[2 * i + IB[j]]     = C[i];
      C[2 * i + 1 - IB[j]] = 0;
    }
    for (int i = 2 * IP - 1; i >= IP; --i) {
      C[i - IP] = C[i - IP] ^ C[i];
      C[i - IQ] = C[i - IQ] ^ C[i];
    }
  }

  for (int j = 0; j < IP; ++j) {
    unsigned iwork = 0;
    for (int i = 0; i < IP; ++i) {
      iwork = iwork ^ (C[i] * IWK[j + i]);
    }
    w[j] = iwork;
  }
}


//====================================================================
void RandomNumbers_Mseries::writefile(const string filename)
{
  vout.general(m_vl, "RandomNumbers_Mseries: write down to file = %s\n", filename.c_str());

  std::ofstream outfile;
  outfile.open(filename.c_str());

  if (!outfile) {
    vout.crucial(m_vl, "RandomNumbers_Mseries: unable to open output file.\n");
    return;
  }

  outfile << " " << jr << " " << kr << std::endl;

  for (int i = 0; i < Np; i++) {
    outfile << " " << w[i] << std::endl;
  }

  outfile.close();
}


//====================================================================
//============================================================END=====