如果你想在编译时初始化一个浮点数组，有几个问题需要克服：
std::array 有点破损，因为在可变 constexpr std::array 的情况下，operator[] 不是 constexpr（我相信这将在标准的下一个版本中得到修复）。
std::math 中的函数未标记为 constexpr！
我最近遇到了类似的问题域。我想创建一个准确但快速的sin(x) 版本。
我决定看看是否可以使用带有插值的constexpr 缓存来获得速度而不损失准确性。
使高速缓存 constexpr 的一个优点是 sin(x) 在编译时已知的值的计算是 sin 是预先计算的，并且只是作为立即寄存器加载存在于代码中！在运行时参数的最坏情况下，它只是一个常量数组查找，然后是 w 加权平均值。
此代码需要在 clang 上使用 -fconstexpr-steps=2000000 编译，或在 windows 中使用等效项。
享受：
#include <iostream>
#include <cmath>
#include <utility>
#include <cassert>
#include <string>
#include <vector>
namespace cpputil {
// a fully constexpr version of array that allows incomplete
// construction
template<size_t N, class T>
struct array
{
// public constructor defers to internal one for
// conditional handling of missing arguments
constexpr array(std::initializer_list<T> list)
: array(list, std::make_index_sequence<N>())
{
}
constexpr T& operator[](size_t i) noexcept {
assert(i < N);
return _data[i];
}
constexpr const T& operator[](size_t i) const noexcept {
assert(i < N);
return _data[i];
}
constexpr T& at(size_t i) noexcept {
assert(i < N);
return _data[i];
}
constexpr const T& at(size_t i) const noexcept {
assert(i < N);
return _data[i];
}
constexpr T* begin() {
return std::addressof(_data[0]);
}
constexpr const T* begin() const {
return std::addressof(_data[0]);
}
constexpr T* end() {
// todo: maybe use std::addressof and disable compiler warnings
// about array bounds that result
return &_data[N];
}
constexpr const T* end() const {
return &_data[N];
}
constexpr size_t size() const {
return N;
}
private:
T _data[N];
private:
// construct each element from the initialiser list if present
// if not, default-construct
template<size_t...Is>
constexpr array(std::initializer_list<T> list, std::integer_sequence<size_t, Is...>)
: _data {
(
Is >= list.size()
?
T()
:
std::move(*(std::next(list.begin(), Is)))
)...
}
{
}
};
// convenience printer
template<size_t N, class T>
inline std::ostream& operator<<(std::ostream& os, const array<N, T>& a)
{
os << "[";
auto sep = " ";
for (const auto& i : a) {
os << sep << i;
sep = ", ";
}
return os << " ]";
}
}
namespace trig
{
constexpr double pi() {
return M_PI;
}
template<class T>
auto constexpr to_radians(T degs) {
return degs / 180.0 * pi();
}
// compile-time computation of a factorial
constexpr double factorial(size_t x) {
double result = 1.0;
for (int i = 2 ; i <= x ; ++i)
result *= double(i);
return result;
}
// compile-time replacement for std::pow
constexpr double power(double x, size_t n)
{
double result = 1;
while (n--) {
result *= x;
}
return result;
}
// compute a term in a taylor expansion at compile time
constexpr double taylor_term(double x, size_t i)
{
int powers = 1 + (2 * i);
double top = power(x, powers);
double bottom = factorial(powers);
auto term = top / bottom;
if (i % 2 == 1)
term = -term;
return term;
}
// compute the sin(x) using `terms` terms in the taylor expansion
constexpr double taylor_expansion(double x, size_t terms)
{
auto result = x;
for (int term = 1 ; term < terms ; ++term)
{
result += taylor_term(x, term);
}
return result;
}
// compute our interpolatable table as a constexpr
template<size_t N = 1024>
struct sin_table : cpputil::array<N, double>
{
static constexpr size_t cache_size = N;
static constexpr double step_size = (pi() / 2) / cache_size;
static constexpr double _360 = pi() * 2;
static constexpr double _270 = pi() * 1.5;
static constexpr double _180 = pi();
static constexpr double _90 = pi() / 2;
constexpr sin_table()
: cpputil::array<N, double>({})
{
for(int slot = 0 ; slot < cache_size ; ++slot)
{
double val = trig::taylor_expansion(step_size * slot, 20);
(*this)[slot] = val;
}
}
double checked_interp_fw(double rads) const {
size_t slot0 = size_t(rads / step_size);
auto v0 = (slot0 >= this->size()) ? 1.0 : (*this)[slot0];
size_t slot1 = slot0 + 1;
auto v1 = (slot1 >= this->size()) ? 1.0 : (*this)[slot1];
auto ratio = (rads - (slot0 * step_size)) / step_size;
return (v1 * ratio) + (v0 * (1.0-ratio));
}
double interpolate(double rads) const
{
rads = std::fmod(rads, _360);
if (rads < 0)
rads = std::fmod(_360 - rads, _360);
if (rads < _90) {
return checked_interp_fw(rads);
}
else if (rads < _180) {
return checked_interp_fw(_90 - (rads - _90));
}
else if (rads < _270) {
return -checked_interp_fw(rads - _180);
}
else {
return -checked_interp_fw(_90 - (rads - _270));
}
}
};
double sine(double x)
{
if (x < 0) {
return -sine(-x);
}
else {
constexpr sin_table<> table;
return table.interpolate(x);
}
}
}
void check(float degs) {
using namespace std;
cout << "checking : " << degs << endl;
auto mysin = trig::sine(trig::to_radians(degs));
auto stdsin = std::sin(trig::to_radians(degs));
auto error = stdsin - mysin;
cout << "mine=" << mysin << ", std=" << stdsin << ", error=" << error << endl;
cout << endl;
}
auto main() -> int
{
check(0.5);
check(30);
check(45.4);
check(90);
check(151);
check(180);
check(195);
check(89.5);
check(91);
check(270);
check(305);
check(360);
return 0;
}
预期输出：
checking : 0.5
mine=0.00872653, std=0.00872654, error=2.15177e-09
checking : 30
mine=0.5, std=0.5, error=1.30766e-07
checking : 45.4
mine=0.712026, std=0.712026, error=2.07233e-07
checking : 90
mine=1, std=1, error=0
checking : 151
mine=0.48481, std=0.48481, error=2.42041e-08
checking : 180
mine=-0, std=1.22465e-16, error=1.22465e-16
checking : 195
mine=-0.258819, std=-0.258819, error=-6.76265e-08
checking : 89.5
mine=0.999962, std=0.999962, error=2.5215e-07
checking : 91
mine=0.999847, std=0.999848, error=2.76519e-07
checking : 270
mine=-1, std=-1, error=0
checking : 305
mine=-0.819152, std=-0.819152, error=-1.66545e-07
checking : 360
mine=0, std=-2.44929e-16, error=-2.44929e-16