lecuyer.h

使用R语言的马尔科夫链蒙特卡洛模拟（MCMC）源代码程序。

        for (int i = 0; i < 6; ++i)
          Cg[i] = Bg[i] = Ig[i];
      }

      /*! \brief Reset the current substream.
       *
       *
       * This method resets the stream to the first state of its
       * current substream.
       *
       * \see ResetStartStream()
       * \see ResetNextSubstream()
       * \see SetSeed(unsigned long seed[6])
       * 
       */
      void
      ResetStartSubstream ()
      {
        for (int i = 0; i < 6; ++i)
          Cg[i] = Bg[i];
      }

      /*! \brief Jump to the next substream.
       *
       * This method resets the stream to the first state of its next
       * substream.
       *
       * \see ResetStartStream()
       * \see ResetStartSubstream()
       * \see SetSeed(unsigned long seed[6])
       * 
       */
      void
      ResetNextSubstream ()
      {
        MatVecModM(A1p76, Bg, Bg, m1);
        MatVecModM(A2p76, &Bg[3], &Bg[3], m2);
        for (int i = 0; i < 6; ++i)
          Cg[i] = Bg[i];
      }

      /*! \brief Set the package seed.
       *
       *  This method sets the overall package seed.  The default
       *  initial seed is (12345, 12345, 12345, 12345, 12345, 12345).
       *  The package seed is the seed used to initialize the first
       *  constructed random number stream in a given program.
       *
       *  \param seed An array of six integers to seed the package.
       *  The first three values cannot all equal 0 and must all be
       *  less than 4294967087 while the second trio of integers must
       *  all be less than 4294944443 and not all 0.
       *
       * \see SetSeed(unsigned long seed[6])
       *
       * \throw scythe_randseed_error (Level 0)
       */
      static void
      SetPackageSeed (unsigned long seed[6])
      {
         if (CheckSeed (seed)) return;
         for (int i = 0; i < 6; ++i)
            nextSeed[i] = seed[i];
      }

      /*! \brief Set the stream seed.
       *
       *  This method sets the stream seed which is used to initialize
       *  the state of the given stream.
       *
       *  \warning This method sets the stream seed in isolation and
       *  does not coordinate with any other streams.  Therefore,
       *  using this method without care can result in multiple
       *  streams that overlap in the course of their runs.
       *
       *  \param seed An array of six integers to seed the stream.
       *  The first three values cannot all equal 0 and must all be
       *  less than 4294967087 while the second trio of integers must
       *  all be less than 4294944443 and not all 0.
       *
       * \see SetPackageSeed(unsigned long seed[6])
       * \see ResetStartStream()
       * \see ResetStartSubstream()
       * \see ResetNextSubstream()
       *
       * \throw scythe_randseed_error (Level 0)
       */
      void
      SetSeed (unsigned long seed[6])
      {
        if (CheckSeed (seed)) return;
          for (int i = 0; i < 6; ++i)
            Cg[i] = Bg[i] = Ig[i] = seed[i];
      }

      // XXX: get the cases formula working!
      /*! \brief Advances the state of the stream.
       *
       * This method advances the input \f$n\f$ steps, using the rule:
       * \f[
       * n =
       * \begin{cases}
       *  2^e + c \quad if~e > 0, \\
       *  -2^{-e} + c \quad if~e < 0, \\
       *  c \quad if~e = 0.
       * \end{cases}
       * \f]
       * 
       * \param e This parameter controls state advancement.
       * \param c This parameter also controls state advancement.
       *
       * \see GetState()
       * \see ResetStartStream()
       * \see ResetStartSubstream()
       * \see ResetNextSubstream()
       */
      void
      AdvanceState (long e, long c)
      {
        double B1[3][3], C1[3][3], B2[3][3], C2[3][3];

        if (e > 0) {
          MatTwoPowModM (A1p0, B1, m1, e);
          MatTwoPowModM (A2p0, B2, m2, e);
        } else if (e < 0) {
          MatTwoPowModM (InvA1, B1, m1, -e);
          MatTwoPowModM (InvA2, B2, m2, -e);
        }

        if (c >= 0) {
          MatPowModM (A1p0, C1, m1, c);
          MatPowModM (A2p0, C2, m2, c);
        } else {
          MatPowModM (InvA1, C1, m1, -c);
          MatPowModM (InvA2, C2, m2, -c);
        }

        if (e) {
          MatMatModM (B1, C1, C1, m1);
          MatMatModM (B2, C2, C2, m2);
        }

        MatVecModM (C1, Cg, Cg, m1);
        MatVecModM (C2, &Cg[3], &Cg[3], m2);
      }

      /*! \brief Get the current state.
       *
       * This method places the current state of the stream, as
       * represented by six integers, into the array argument.  This
       * is useful for saving and restoring streams across program
       * runs.
       *
       * \param seed An array of six integers that will hold the state values on return.
       *
       * \see AdvanceState()
       */
      void
      GetState (unsigned long seed[6]) const
      {
        for (int i = 0; i < 6; ++i)
          seed[i] = static_cast<unsigned long> (Cg[i]);
      }

      /*! \brief Toggle generator precision.
       *
       * This method sets the precision level of the given stream.  By
       * default, streams generate random numbers with 32 bit
       * resolution.  If the user invokes this method with \a incp =
       * true, then the stream will begin to generate variates with
       * greater precision (53 bits on machines following the IEEE 754
       * standard).  Calling this method again with \a incp = false
       * will return the precision of generated numbers to 32 bits.
       *
       * \param incp A boolean value where true implies high (most
       * likely 53 bit) precision and false implies low (32 bit)
       * precision.
       *
       * \see SetAntithetic(bool)
       */
      void
      IncreasedPrecis (bool incp)
      {
        incPrec = incp;
      }

      /*! \brief Toggle the orientation of generated random numbers.
       *
       * This methods causes the given stream to generate antithetic
       * (1 - U, where U is the default number generated) when called
       * with \a a = true.  Calling this method with \a a = false will
       * return generated numbers to their default orientation.
       *
       * \param a A boolean value that selects regular or antithetic
       * variates.
       *
       * \see IncreasedPrecis(bool)
       */
      void
      SetAntithetic (bool a)
      {
        anti = a;
      }

      /*! \brief Generate a random uniform variate on (0, 1).
       *
       * This routine returns a random double precision floating point
       * number from the uniform distribution on the interval (0,
       * 1).  This method overloads the pure virtual method of the
       * same name in the rng base class.
       *
       * \see runif(unsigned int, unsigned int)
       * \see RandInt(long, long)
       * \see rng
       */
      double
      runif ()
      {
        if (incPrec)
          return U01d();
        else
          return U01();
      }

      /* We have to override the overloaded form of runif because
       * overloading the no-arg runif() hides the base class
       * definition; C++ stops looking once it finds the above.
       */
      /*! \brief Generate a Matrix of random uniform variates.
       *
       * This routine returns a Matrix of double precision random
       * uniform variates. on the interval (0, 1).  This method
       * overloads the virtual method of the same name in the rng base
       * class.
       *
       * This is the general template version of this method and
       * is called through explicit template instantiation.
       *
       * \param rows The number of rows in the returned Matrix.
       * \param cols The number of columns in the returned Matrix.
       * 
       * \see runif()
       * \see rng
       *
       * \note We are forced to override this overloaded method
       * because the 1-arg version of runif() hides the base class's
       * definition of this method from the compiler, although it
       * probably should not.
       */
      template <matrix_order O, matrix_style S>
      Matrix<double,O,S> runif(unsigned int rows, unsigned int cols)
      {
        return rng<lecuyer>::runif<O,S>(rows,cols);
      }

      /*! \brief Generate a Matrix of random uniform variates.
       *
       * This routine returns a Matrix of double precision random
       * uniform variates on the interval (0, 1).  This method
       * overloads the virtual method of the same name in the rng base
       * class.
       *
       * This is the default template version of this method and
       * is called through implicit template instantiation.
       *
       * \param rows The number of rows in the returned Matrix.
       * \param cols The number of columns in the returned Matrix.
       * 
       * \see runif()
       * \see rng
       *
       * \note We are forced to override this overloaded method
       * because the 1-arg version of runif() hides the base class's
       * definition of this method from the compiler, although it
       * probably should not.
       */
      Matrix<double,Col,Concrete> runif(unsigned int rows,
                                        unsigned int cols)
      {
        return rng<lecuyer>::runif<Col,Concrete>(rows, cols);
      }

      /*! \brief Generate the next random integer.
       *
       * This method generates a random integer from the discrete
       * uniform distribution on the interval [\a low, \a high].
       *
       * \param low The lower bound of the interval to evaluate.
       * \param high the upper bound of the interval to evaluate.
       *
       * \see runif()
       */
      long
      RandInt (long low, long high)
      {
        return low + static_cast<long> ((high - low + 1) * runif ());
      }

    protected:
      // Generate the next random number.
      //
      double
      U01 ()
      {
        long k;
        double p1, p2, u;

        /* Component 1 */
        p1 = a12 * Cg[1] - a13n * Cg[0];
        k = static_cast<long> (p1 / m1);
        p1 -= k * m1;
        if (p1 < 0.0) p1 += m1;
        Cg[0] = Cg[1]; Cg[1] = Cg[2]; Cg[2] = p1;

        /* Component 2 */
        p2 = a21 * Cg[5] - a23n * Cg[3];
        k = static_cast<long> (p2 / m2);
        p2 -= k * m2;
        if (p2 < 0.0) p2 += m2;
        Cg[3] = Cg[4]; Cg[4] = Cg[5]; Cg[5] = p2;

        /* Combination */
        u = ((p1 > p2) ? (p1 - p2) * norm : (p1 - p2 + m1) * norm);

        return (anti == false) ? u : (1 - u);
      }

      // Generate the next random number with extended (53 bits) precision.
      double 
      U01d ()
      {
        double u;
        u = U01();
        if (anti) {
          // Don't forget that U01() returns 1 - u in the antithetic case
          u += (U01() - 1.0) * fact;
          return (u < 0.0) ? u + 1.0 : u;
        } else {
          u += U01() * fact;
          return (u < 1.0) ? u : (u - 1.0);
        }
      }


      // Public members of the class start here
      
      // The default seed of the package; will be the seed of the first
      // declared RngStream, unless SetPackageSeed is called.
      static double nextSeed[6];

      /* Instance variables */
      double Cg[6], Bg[6], Ig[6];


      bool anti, incPrec;


      std::string streamname_;

  };

}

#endif /* SCYTHE_LECUYER_H */