stochastic.h

Go to the documentation of this file.

// SPDX-License-Identifier: Apache-2.0
// Copyright Contributors to the OpenQMC Project.
 
 
#pragma once
 
#include "lookup.h"
#include "pcg.h"
#include "unused.h"
 
#include <cassert>
#include <cstdint>
 
namespace oqmc
{
 
inline void stochasticPmjInit(int nsamples, std::uint32_t table[][4])
{
    constexpr auto maxIndexSize = 0x10000; // 2^16 index upper limit.
    OQMC_MAYBE_UNUSED(maxIndexSize);
 
    assert(nsamples >= 1);
    assert(nsamples <= maxIndexSize);
 
    // clang-format off
    constexpr std::uint16_t pmjXors[2][16] = {
        {
        0b0000000000000000,
        0b0000000000000000,
        0b0000000000000010,
        0b0000000000000110,
        0b0000000000000110,
        0b0000000000001110,
        0b0000000000110110,
        0b0000000001001110,
        0b0000000000010110,
        0b0000000000101110,
        0b0000001001110110,
        0b0000011011001110,
        0b0000011100010110,
        0b0000110000101110,
        0b0011000001110110,
        0b0100000011001110,
        },
 
        {
        0b0000000000000000,
        0b0000000000000001,
        0b0000000000000011,
        0b0000000000000011,
        0b0000000000000111,
        0b0000000000011011,
        0b0000000000100111,
        0b0000000000001011,
        0b0000000000010111,
        0b0000000100111011,
        0b0000001101100111,
        0b0000001110001011,
        0b0000011000010111,
        0b0001100000111011,
        0b0010000001100111,
        0b0000000010001011,
        },
    };
    // clang-format on
 
    const auto buffer = new std::uint32_t[nsamples][2];
 
    auto state = pcg::init();
 
    for(int k = 0; k < 2; ++k)
    {
        buffer[0][k] = pcg::rng(state);
    }
 
    for(int prevLen = 1, logN = 0; prevLen < nsamples; prevLen *= 2, ++logN)
    {
        for(int i1 = 0, i2 = prevLen; i1 < prevLen && i2 < nsamples; ++i1, ++i2)
        {
            for(int k = 0; k < 2; ++k)
            {
                const auto swapBit = 0x80000000u >> logN;
                const auto bitMask = swapBit - 1;
 
                const auto j = i1 ^ pmjXors[k][logN];
 
                const auto prevStratum = buffer[j][k] & ~bitMask;
                const auto nextStratum = prevStratum ^ swapBit;
 
                buffer[i2][k] = nextStratum | (pcg::rng(state) & bitMask);
            }
        }
    }
 
    for(int i = 0; i < nsamples; ++i)
    {
        shuffledScrambledLookup<2, 2>(i, pcg::hash(0), buffer, &table[i][0]);
        shuffledScrambledLookup<2, 2>(i, pcg::hash(1), buffer, &table[i][2]);
    }
 
    delete[] buffer;
}
 
} // namespace oqmc