HeteQuadrature/algoim/organizer/halton_algoim.hpp

/**
 * @file halton.hpp
 * @author Wang Zhicheng (modified from John Burkardt's)
 * @brief Halton sequence generator
 * @attention
 * 1. This code is modified from John Burkardt's code (https://people.sc.fsu.edu/~jburkardt/cpp_src/halton/halton.html) \n
 * 2. The original code is under the MIT license. \n
 * 3. The modified code is under the same license as the original code. \n
 * 4. Reference: \n
 *    John Halton, \n
 *    On the efficiency of certain quasi-random sequences of points \n
 *    in evaluating multi-dimensional integrals, \n
 *    Numerische Mathematik, \n
 *    Volume 2, pages 84-90, 1960. \n \n
 *    Milton Abramowitz, Irene Stegun, \n
 *    Handbook of Mathematical Functions, \n
 *    US Department of Commerce, 1964, \n
 *    ISBN: 0-486-61272-4, \n
 *    LC: QA47.A34. \n \n
 *    Daniel Zwillinger, \n
 *    CRC Standard Mathematical Tables and Formulae, \n
 *    30th Edition, \n
 *    CRC Press, 1996, pages 95-98.
 */
#pragma once

#include <cstdint>
#include <iostream>
#include <vector>

#include "uvector.hpp"

static constexpr size_t  g_prime_max = 1600;
static constexpr int32_t g_npvec[g_prime_max]{
#include "default_halton_prime.inl"
};

/**
 * @brief PRIME returns any of the first PRIME_MAX prime numbers.
 * @note
    PRIME_MAX is 1600, and the largest prime stored is 13499. \n
    Thanks to Bart Vandewoestyne for pointing out a typo, 18 February 2005.
 * @param[in] n
 *  the index of the desired prime number. In general, is should be true that 0 <= N <= PRIME_MAX. \n
 *  N = -1 returns PRIME_MAX, the index of the largest prime available. \n
 *  N = 0 is legal, returning PRIME = 1.
 * @return the N-th prime.  If N is out of range, PRIME is returned as -1.
 */
static int32_t prime(size_t n) noexcept
{
    if (n > g_prime_max) std::cerr << "Fatal: PRIME - Unexpected input value of n = " << n << std::endl;
    return g_npvec[std::max(n - 1, size_t{0})];
}

static constexpr int32_t prime_guaranteed(size_t n) noexcept { return g_npvec[n - 1]; }

static constexpr int32_t* prime_guaranteed_ptr(size_t n) noexcept { return (int32_t*)&g_npvec[n - 1]; }

/**
 * @brief HALTON computes an element of a Halton sequence.
 * @param[in] i the index of the element to compute (i >= 0).
 * @param[in] m the spatial dimension.
 * @return the element of the sequence with index I.
 */
template <size_t m, typename _Tp>
auto halton(size_t i) noexcept
{
    // Map<Matrix<value_type, m, Dynamic>> r(&*d_first);
    // Map<Vector<int32_t, m>>             prime(prime_guaranteed_ptr(1));
    // Vector<value_type, m>               prime_inv_keep, prime_inv;
    // Vector<int32_t, m>                  t = Vector<int32_t, m>::Constant(static_cast<int32_t>(i));
    // r.setZero();
    // for (auto i = 0; i < m; ++i) {
    //     prime_inv_keep[i] = 1.0 / prime_guaranteed(i + 1);
    //     prime_inv[i]      = 1.0 / prime_guaranteed(i + 1);
    // }

    algoim::uvector<_Tp, m>     res{};
    algoim::uvector<_Tp, m>     prime_inv_keep, prime_inv;
    algoim::uvector<int32_t, m> t = static_cast<int32_t>(i);
    for (auto i = 0; i < m; ++i) {
        prime_inv_keep[i] = 1.0 / prime_guaranteed(i + 1);
        prime_inv[i]      = prime_inv_keep[i];
    }

    while (algoim::sum(t) > 0) {
        for (auto i = 0; i < m; ++i) {
            res[i]       += prime_inv[i] * (t[i] % prime_guaranteed(i + 1));
            prime_inv[i] *= prime_inv_keep[i];
            t[i]          = static_cast<int32_t>(t[i] * prime_inv_keep[i]);
        }
    }
}

/**
 * @brief HALTON_SEQUENCE computes elements in [i_first, i_last) of a Halton sequence.
 * @param[in] i_first the indices of the first element of the sequence. (i_first >= 0).
 * @param[in] i_last the indices of the last + 1 element of the sequence. (i_last >= 0).
 * @param[in] m the spatial dimension.
 * @return the elements of the sequence with indices [i_first, i_last).
 */
template <size_t m, typename _Tp>
auto halton_sequence(size_t i_first, size_t i_last) noexcept
{
    const auto n = i_last - i_first;

    std::vector<algoim::uvector<_Tp, m>> res(n);
    algoim::uvector<_Tp, m>              prime_inv, prime_inv_keep;
    algoim::uvector<int32_t, m>          t{};
    for (auto i = 0; i < m; ++i) { prime_inv_keep[i] = 1.0 / prime_guaranteed(i + 1); }
    for (auto i = 0; i < n; ++i) {
        t         = static_cast<int32_t>(i + i_first);
        prime_inv = prime_inv_keep;
        while (algoim::sum(t) > 0) {
            for (auto j = 0; j < m; ++j) {
                res[i][j]    += prime_inv[j] * (t[j] % prime_guaranteed(j + 1));
                prime_inv[j] *= prime_inv_keep[j];
                t[j]          = static_cast<int32_t>(t[j] * prime_inv_keep[j]);
            }
        }
    }

    return res;
}