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
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;
|
|
}
|