You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
 
 

124 lines
4.7 KiB

/**
* @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;
}