lecuyeruniformrng.cpp

/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
 
/*
 Copyright (C) 2000, 2001, 2002, 2003 RiskMap srl
 
 This file is part of QuantLib, a free-software/open-source library
 for financial quantitative analysts and developers - http://quantlib.org/
 
 QuantLib is free software: you can redistribute it and/or modify it
 under the terms of the QuantLib license.  You should have received a
 copy of the license along with this program; if not, please email
 <quantlib-dev@lists.sf.net>. The license is also available online at
 <http://quantlib.org/license.shtml>.
 
 This program 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 license for more details.
*/
 
#include <ql/math/randomnumbers/lecuyeruniformrng.hpp>
#include <ql/math/randomnumbers/seedgenerator.hpp>
 
namespace QuantLib {
 
    const long LecuyerUniformRng::m1 = 2147483563L;
    const long LecuyerUniformRng::a1 = 40014L;
    const long LecuyerUniformRng::q1 = 53668L;
    const long LecuyerUniformRng::r1 = 12211L;
 
    const long LecuyerUniformRng::m2 = 2147483399L;
    const long LecuyerUniformRng::a2 = 40692L;
    const long LecuyerUniformRng::q2 = 52774L;
    const long LecuyerUniformRng::r2 = 3791L;
 
    const int LecuyerUniformRng::bufferSize = 32;
 
    // int(1+m1/bufferSize) = int(1+(m1-1)/bufferSize)
    const long LecuyerUniformRng::bufferNormalizer = 67108862L;
 
    const long double LecuyerUniformRng::maxRandom = 1.0-QL_EPSILON;
 
    LecuyerUniformRng::LecuyerUniformRng(long seed)
    : buffer(LecuyerUniformRng::bufferSize) {
        // Need to prevent seed=0, so use seed=0 to have a "random" seed
        temp2 = temp1 = (seed != 0 ? seed : SeedGenerator::instance().get());
        // Load the shuffle table (after 8 warm-ups)
        for (int j=bufferSize+7; j>=0; j--) {
            long k = temp1/q1;
            temp1 = a1*(temp1-k*q1)-k*r1;
            if (temp1 < 0)
                temp1 += m1;
            if (j < bufferSize)
                buffer[j] = temp1;
        }
        y = buffer[0];
    }
 
    LecuyerUniformRng::sample_type LecuyerUniformRng::next() const {
        long k = temp1/q1;
        // Compute temp1=(a1*temp1) % m1
        // without overflows (Schrage's method)
        temp1 = a1*(temp1-k*q1)-k*r1;
        if (temp1 < 0)
            temp1 += m1;
        k = temp2/q2;
        // Compute temp2=(a2*temp2) % m2
        // without overflows (Schrage's method)
        temp2 = a2*(temp2-k*q2)-k*r2;
        if (temp2 < 0)
            temp2 += m2;
        // Will be in the range 0..bufferSize-1
        int j = y/bufferNormalizer;
        // Here temp1 is shuffled, temp1 and temp2 are
        // combined to generate output
        y = buffer[j]-temp2;
        buffer[j] = temp1;
        if (y < 1)
            y += m1-1;
        double result = y/double(m1);
        // users don't expect endpoint values
        if (result > maxRandom)
            result = (double) maxRandom;
        return {result, 1.0};
    }
 
}