How to initialise a floating point array at compile time?

Question

I have found two good approaches to initialise integral arrays at compile times here and here.

Unfortunately, neither can be converted to initialise a float array straightforward; I find that I am not fit enough in template metaprogramming to solve this 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, as far as the availability of information goes, this array could be declared constexpr.

However, for the first approach linked, the generator function f would have to look like this:

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

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

Now, the first idea that springs to mind would be to store the float in an int variable - nothing happens to the array indices after their calculation, so pretending that they were of another type than they are (as long as byte-length is equal) is perfectly fine.

But see:

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

is not a valid constexpr, as 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)));
}

does not work either; even though one might think that reinterpret_cast does exactly what is needed here (namely nothing), it apparently only works on pointers.

Following the second approach, the generator function would look slightly different:

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

With what is essentially the same problem: That enum cannot be of type float and the type cannot be masked as int.

---

Now, even though I have only approached the problem on the path of "pretend the float is an int", I do not actually like that path (aside from it not working). I would much prefer a way that actually handled the float as float (and would just as well handle a double as double), but I see no way to get around the type restriction imposed.

Sadly, there are many questions about this issue, which always refer to integral types, swamping the search for this specialised issue. Similarly, questions about masking one type as the other typically do not consider the restrictions of a constexpr or template parameter environment.
See [1][2][3] and [4][5] etc.

Answer 1

Assuming your actual goal is to have a concise way to initialize an array of floating point numbers and it isn't necessarily spelled 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 significance of the type of array is that std::array<T, N> can be copied - unlike T[N]. If it can be copied, you can obtain the content of the array from a function call, e.g.:

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

How does that help us? Well, when we can call a function taking an integer as an argument, we can use std::make_index_sequence<N> to give use a compile-time sequence of std::size_t from 0 to N-1. If we have that, we can initialize an array easily with a formula based on the index 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 the function used to initialize the array elements is actually a constexpr expression, this approach can generate a constexpr. The function const_sin() which is there just for demonstration purpose does that but it, obviously, doesn't compute a reasonable approximation of sin(x).

The comments indicate that the answer so far doesn't quite explain what's going on. So, let's break it down into digestible parts:

1. The goal is to produce a constexpr array filled with suitable sequence of values. However, the size of the array should be easily changeable by adjusting just the array size N. That is, conceptually, the objective is to create

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

   Where f() is a suitable function producing a constexpr. For example, f() could be defined as

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

   However, typing in the calls to f(0), f(1), etc. would need to change with every change of N. So, essentially the same as the above declaration should be done but without extra typing.

2. The first step towards the solution is to replace float[N] by std:array<float, N>: built-in arrays cannot be copied while std::array<float, N> can be copied. That is, the initialization could be delegated to to a function parameterized by N. That is, we'd 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>();
   

3. Within the function we can't simply loop over the array because the non-const subscript operator isn't a constexpr. Instead, the array needs to be initialized upon creation. If we had a parameter pack std::size_t... I which represented the sequence 0, 1, .., N-1 we could just do

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

   as the expansion would effectively become equivalent to typing

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

   So the question becomes: how to get such a parameter pack? I don't think it can be obtained directly in the function but it can be obtained by calling another function with a suitable parameter.

4. The using alias std::make_index_sequence<N> is an alias for the type std::index_sequence<0, 1, .., N-1>. The details of the implementation are a bit arcane but std::make_index_sequence<N>, std::index_sequence<...>, and friends are part of C++14 (they were proposed by N3493 based on, e.g., on this answer from me). That is, all we need to do is call an auxiliary function with a parameter of type std::index_sequence<...> and get the parameter pack 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 to this function is only used to have the parameter pack std::size_t... I be deduced.

Answer 2

Here's a working example that generates a table of sin values, and that you can easily adapt to logarithm tables by passing a different 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";
}

Live Example.

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

Also note that you need a constexpr function object (both pi and std::sin are not constexpr) to initialize the array truly at compile-time rather than at program initialization.

Answer 3

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

1. std::array is a little broken in that the operator[] is not constexpr in the case of a mutable constexpr std::array (I believe this will be fixed in the next release of the standard).

2. the functions in std::math are not marked constexpr!

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

I decided to see if it could be done with a constexpr cache with interpolation to get speed without loss of accuracy.

An advantage of making the cache constexpr is that the calculation of sin(x) for a value known at compile-time is that the sin is pre-computed and simply exists in the code as an immediate register load! In the worst case of a runtime argument, it's merely a constant array lookup followed by w weighted average.

This code will need to be compiled with -fconstexpr-steps=2000000 on clang, or the 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 4

I am just keeping this answer for documentation. As the comments say, I was mislead by gcc being permissive. It fails, when f(42) is used e.g. as a template parameter like this:

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

sorry, this was not a solution

Separate the calculation of your float and the conversion to an 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)));
}