How to initialize a floating point array at compile time?

Question

I found two good approaches for initializing integrated arrays at compile time here and here .

Unfortunately, none of them can be converted to initialize the float array. I found that I am not well suited for metaprogramming templates to solve this problem through trial and error.

First let me declare a use case:

constexpr unsigned int SineLength = 360u; constexpr unsigned int ArrayLength = SineLength+(SineLength/4u); constexpr double PI = 3.1415926535; float array[ArrayLength]; void fillArray(unsigned int length) { for(unsigned int i = 0u; i < length; ++i) array[i] = sin(double(i)*PI/180.*360./double(SineLength)); } 

As you can see how accessible the information is, this array can be declared constexpr .

However, for the first coupled approach, the generator function f should look something like this:

 constexpr float f(unsigned int i) { return sin(double(i)*PI/180.*360./double(SineLength)); } 

And that means a template argument of type float needed. What is not allowed.

Now the first idea that comes to mind was to store a float in an int variable - nothing happens to the array indices after they are calculated, so pretending to be of a different type than they are (as long as the length is equal) excellent.

But look:

 constexpr int f(unsigned int i) { float output = sin(double(i)*PI/180.*360./double(SineLength)); return *(int*)&output; } 

constexpr is not valid since it contains more than the return statement.

 constexpr int f(unsigned int i) { return reinterpret_cast<int>(sin(double(i)*PI/180.*360./double(SineLength))); } 

also does not work; even if you might think that reinterpret_cast does exactly what is needed (namely nothing), it seems to work only with pointers.

Following the second approach, the generator function would look a little different:

 template<size_t index> struct f { enum : float{ value = sin(double(index)*PI/180.*360./double(SineLength)) }; }; 

With the same problem: this enumeration cannot be of type float , and the type cannot be disguised as int .

---

Now, although I only approached the problem on the path "pretend that float is int ", I really don't like this path (except it doesn't work). I would prefer the way that float actually treated as a float (and it would have treated double as a double as well), but I see no way around the type restriction introduced.

Unfortunately, there are many questions on this problem that always relate to integral types, waving the search for this specialized problem. Similarly, questions about disguising one type as another usually do not take into account the constraints of the constexpr environment or template. See [1] [2] [3] and [4] [5] etc.

Answer 1

Assuming your actual goal is to briefly initialize an array of floating point numbers, and this is not necessarily written float array[N] or double array[N] , but rather std::array<float, N> array or std::array<double, N> array , this can be done.

The value of the array type is that std::array<T, N> can be copied - unlike T[N] . If it can be copied, you can get the contents of the array from a function call, for example:

 constexpr std::array<float, ArrayLength> array = fillArray<N>(); 

How does this help us? Well, when we can call a function that takes an integer as an argument, we can use std::make_index_sequence<N> to use the compile time sequence std::size_t from 0 to N-1 . If we have this, we can easily initialize the array with an index-based formula, like this:

 constexpr double const_sin(double x) { return x * 3.1; } // dummy... template <std::size_t... I> constexpr std::array<float, sizeof...(I)> fillArray(std::index_sequence<I...>) { return std::array<float, sizeof...(I)>{ const_sin(double(I)*M_PI/180.*360./double(SineLength))... }; } template <std::size_t N> constexpr std::array<float, N> fillArray() { return fillArray(std::make_index_sequence<N>{}); } 

Assuming that the function used to initialize the elements of the array is actually a constexpr expression, this approach can generate a constexpr . The const_sin() function, which exists only for a demonstration purpose, does this, but it obviously does not calculate a reasonable approximation of sin(x) .

Comments show that the answer does not yet fully explain what is happening. So let me break it down to digestible parts:

• The goal is to create a constexpr array filled with a suitable sequence of values. However, the size of the array should be easily replaced by adjusting the size of array N That is, conceptually, the goal is to create

   constexpr float array[N] = { f(0), f(1), ..., f(N-1) }; 

  Where f() is a suitable function that creates a constexpr . For example, f() can be defined as

   constexpr float f(int i) { return const_sin(double(i) * M_PI / 180.0 * 360.0 / double(Length); } 

  However, with every change in N when entering calls f(0) , f(1) , etc. need to change. Thus, it is essentially the same as the above expression, but without additional customization.

• The first step to solving is to replace float[N] with std:array<float, N> : inline arrays cannot be copied, and std::array<float, N> can be copied. That is, initialization can be delegated to a function parameterized by N That is, we will use

   template <std::size_t N> constexpr std::array<float, N> fillArray() { // some magic explained below goes here } constexpr std::array<float, N> array = fillArray<N>(); 

• Inside the function, we cannot simply iterate over the array due to the fact that the index operator is not const not constexpr . Instead, the array should be initialized at creation. If we had a package of parameters std::size_t... I , which represented the sequence 0, 1, .., N-1 , we could just do

   std::array<float, N>{ f(I)... }; 

  since decomposition will actually become equivalent to type

   std::array<float, N>{ f(0), f(1), .., f(N-1) }; 

  So the question is: how to get such a package of parameters? I do not think that it can be obtained directly in the function, but it can be obtained by calling another function with the appropriate parameter.

• The alias for using std::make_index_sequence<N> is an alias for the type std::index_sequence<0, 1, .., N-1> . The implementation details are a little cryptic, but std::make_index_sequence<N> , std::index_sequence<...> , and friends are part of C ++ 14 (they were proposed by N3493 based on, for example, this answer from me ). That is, all we need to do is call a helper function with a parameter of type std::index_sequence<...> and get a package of parameters from there:

   template <std::size_t...I> constexpr std::array<float, sizeof...(I)> fillArray(std::index_sequence<I...>) { return std::array<float, sizeof...(I)>{ f(I)... }; } template <std::size_t N> constexpr std::array<float, N> fillArray() { return fillArray(std::make_index_sequence<N>{}); } 

  The [unnamed] parameter for this function is used only to output the parameter package std::size_t... I

Answer 2

There are a few issues to overcome if you want to initialize a floating point array at compile time:

• std::array little broken by the fact that operator[] not constexpr in the case of mutable constexpr std :: array (I believe that this will be fixed in the next version of the standard).

• functions in std :: math are not marked constexpr!

I recently had a similar problem area. I wanted to create an accurate but quick version of sin(x) .

I decided to see if this could be done with the constexpr interpolation cache to get speed without losing accuracy.

The advantage of creating a constexpr cache is that evaluating sin(x) for the value known at compile time is that sin precomputed and simply exists in the code as directly loading the register! In the worst case, the runtime argument is just a constant search for the array, and then w is the weighted average.

This code needs to be compiled with -fconstexpr-steps=2000000 on clang or equivalent in windows.

enjoy:

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

expected output:

 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 

Answer 3

Here is a working example that generates a table of sin values, and that you can easily adapt to logarithmic tables by passing another function object

 #include <array> // array #include <cmath> // sin #include <cstddef> // size_t #include <utility> // index_sequence, make_index_sequence #include <iostream> namespace detail { template<class Function, std::size_t... Indices> constexpr auto make_array_helper(Function f, std::index_sequence<Indices...>) -> std::array<decltype(f(std::size_t{})), sizeof...(Indices)> { return {{ f(Indices)... }}; } } // namespace detail template<std::size_t N, class Function> constexpr auto make_array(Function f) { return detail::make_array_helper(f, std::make_index_sequence<N>{}); } static auto const pi = std::acos(-1); static auto const make_sin = [](int x) { return std::sin(pi * x / 180.0); }; static auto const sin_table = make_array<360>(make_sin); int main() { for (auto elem : sin_table) std::cout << elem << "\n"; } 

Living example.

Note that you need to use -stdlib=libc++ , because libstdc++ has a rather inefficient index_sequence implementation.

Also note that you need a constexpr function constexpr (both pi and std::sin not constexpr ) to initialize the array really at compile time, and not during program initialization.

Answer 4

I just keep this answer for documentation. As the comments say, I was misled when gcc was permissive. It does not work when f(42) , for example. as a template parameter like this:

 std::array<int, f(42)> asdf; 

sorry this is not a solution

Separate the calculation of your float and the conversion to int in two different constexpr functions:

 constexpr int floatAsInt(float float_val) { return *(int*)&float_val; } constexpr int f(unsigned int i) { return floatAsInt(sin(double(i)*PI/180.*360./double(SineLength))); }