1

I have written the following code:

#include<array>
#include<type_traits>

namespace math{
    namespace detail{

        template<std::size_t... Is> struct seq{};

        template<std::size_t N, std::size_t... Is>
        struct gen_seq : gen_seq<N-1, N-1, Is...>{};

        template<std::size_t... Is>
        struct gen_seq<0, Is...> : seq<Is...>{};

        template<class T,std::size_t N>
        struct sin_coeffs{
            using array_type = std::array<T,N>;
            constexpr static inline T coeff(std::size_t n){
                return power(-1, n-1) * inverse((T)factorial((2 * n)-1));
            }
            template<std::size_t...NS>
            constexpr static array_type _coeffs(seq<NS...>){
                return {{coeff(NS)...}};
            }
            constexpr static array_type coeffs=_coeffs(gen_seq<N>{});
        };
    }

    template<class T,std::size_t N = max_factorial, class dcy = std::decay_t<T>>
    constexpr std::enable_if_t<std::is_floating_point<dcy>::value,dcy> sin(T x) noexcept{
        constexpr std::array<dcy,N>& coeffs = detail::sin_coeffs<dcy,N>::coeffs;
        const dcy x_2 = x*x;
        dcy pow = x;
        dcy result = 0;
        for(std::size_t i=0;i<N;++i){
            result += coeffs[i] * pow;
            pow*=x_2;
        }
        return result;
    }
}

example main:

int main()
{
    constexpr double d = math::sin(0.0);
}

The code is designed to be a constexpr sin function that uses a constexpr array that holds the coefficients to do the needed calculations.

All of the functions not listed are present in a separate header and compile with no problem.

I am trying to use the "indicies" trick to fill the array using a constexpr function based off this answer to another question.

I am compiling with GCC 5.3.1 with the flags

--std=c++1z -pthread -g -O3 -MMD -MP -Wall -pedantic

The compiler emits no errors to my code but stalls when compiling.

I have let the compilation run for a few minutes but it does not do anything.

I have tested all of the functions used in the code and they all compile just fine independently of this section.

int main()
{
    math::detail::sin_coeffs<double,20>::coeffs[0];
}

This code snippet also reproduces the issue which leads me to believe that it is not related to the sin function itself but the sin_coeffs struct.

EDIT: Here are the other functions as requested:

#include <type_traits>
namespace math{

    template<class T,class dcy = std::decay_t<T>>
    constexpr inline std::enable_if_t<std::is_floating_point<T>::value,dcy> inverse(T value){
        return (value == 0) ? 0.0 : 1.0 / value;
    }
    template <class T>
    constexpr inline std::decay_t<T> sign(T value) {
        return value < 0 ? -1 : 1;
    }
    template <typename T>
    constexpr inline std::decay_t<T> abs(T value) {
        return value * sign(value);
    }
    template<class T>
    constexpr inline std::decay_t<T> power(T const& base, std::size_t const& pow){
        if(pow==0){return 1;}
        else if(pow == 1){return base;}
        else{
            T result = base;
            for(std::size_t i=1;i<pow;++i){
                result*=base;
            }
            return result;
        }
    }
    constexpr std::intmax_t factorial(std::intmax_t const& n){
        if(n==0){return 1;}
        std::intmax_t result = n;
        for(std::intmax_t i=n-1;i>0;--i){
            result *=i;
        }
        return result;
    }
    constexpr static std::size_t max_factorial = 20;//replace with calculated version later
}
Community
  • 1
  • 1
Alex Zywicki
  • 2,263
  • 1
  • 19
  • 34
  • 1
    Did you try with a recent [Clang](http://clang.llvm.org/) (3.7 or 3.8) or a [GCC6](http://gcc.gnu.org/gcc-6) compiler? – Basile Starynkevitch Jun 28 '16 at 06:57
  • Same error with GCC6 and I dont have Clang right now but will install and try when i can – Alex Zywicki Jun 28 '16 at 06:58
  • 3
    Could you provide at least the declarations of the missing function (without body) so that we could compile? – Holt Jun 28 '16 at 07:02
  • 2
    To echo @Holt, without `power` etc., the code you've shown is not _that_ useful. See [sscce.org](http://sscce.org/). – ildjarn Jun 28 '16 at 07:10
  • 1
    I am not 100% sure, but since the terminal condition for `gen_seq` is `N = 0`, you probably have a `0` in your seq, thus you are going to call `coeff(0)` which will call `power(-1, ((size_t)0) - 1)` which is `power(-1, 18446744073709551615)` on my architecture, which will obviously never compile. – Holt Jun 28 '16 at 07:15
  • I downloaded clang which is actually giving me errors to work from now. I will look into them ASAP – Alex Zywicki Jun 28 '16 at 07:28
  • 1
    @AlexZywicki See my (now edited) answer for indications on *why* your code did not compile and how to fix everything. – Holt Jun 28 '16 at 08:06

1 Answers1

6

How to fix your code?

  1. You are missing a const qualifier here:
constexpr const std::array<dcy,N>& coeffs = /* ... */ ;
          ^^^^^
  1. Your gen_seq generates values from 0 to N - 1, but regarding your coeff function, you want values from 1 to N. You can fix this by either:
  • Changing the base template to generate values from 1 to N (see the bottom of this answer for detailed explanation):
template<std::size_t N, std::size_t... Is>
struct gen_seq: gen_seq<N-1, N, Is...> {};
//                           ^--- N instead of N - 1
  • Changing the way you call coeff (thanks to @T.C., see comments):
template<std::size_t...NS>
constexpr static array_type _coeffs(seq<NS...>){
    return {{coeff(NS + 1)...}};
//                   ^^^^
}
  1. Your should not use power to compute -1 ** n, use a ternary condition:
constexpr static inline T coeff(std::size_t n){
    return (n%2 ? 1 : -1) * inverse((T)factorial((2 * n)-1));
}

With this, I can compute the coeffs:

auto arr = math::detail::sin_coeffs<double, 10>::coeffs;
for (auto x: arr) {
    std::cout << x << " ";
}

Output:

1 -0.166667 0.00833333 -0.000198413 2.75573e-06 -2.50521e-08 1.6059e-10 -7.64716e-13 ...

As far as I know, these are the correct coefficients (1, -1/3!, 1/5!, ...). Note that I had to use N = 10, or I would have overflow std::intmax_t (on my architecture) - Clang warns you at compile time if you overflow std::intmax_t with factorial (what a nice compiler!).

With the two above modifications, your code works fine (except for max_factorial value, but you can tune this as you like).

See code on rextester: http://rextester.com/CRR35028

What was the main problem in your code?

Your gen_seq<N> generated a sequence from 0 to N - 1, so you were calling coeff(0), which was calling power(-1, static_cast<size_t>(0) - 1), which was actually power(-1, 18446744073709551615) (on my architecture), which cannot compiles. Adding a trivial case in power "fix" the compilation (showing that this was an issue, but not resolving the real issue):

template<class T>
constexpr inline std::decay_t<T> power(T const& base, std::size_t const& pow) {
    if (pow == static_cast<size_t>(0) - 1) { // Stupid test
      return 1;
    }
    /* ... */
}

Also, your max_factorial value was also probably too big. After the correction to power, I could not compile for max_factorial > 11 (I probably have 32 bits std::intmax_t so you can probably go above this, but I think that 20 is too big in all cases).

Also, for future problem, clang seems to provide better information on why it does not compile:

  • It has a limit on the number of iterations, so I was able to find out that power was doing an (almost) infinite loop.
  • It gives a warning for the return value of factorial which overflow intmax_t.

How does the gen_seq trick work?

Your seq struct is basically a structure holder for a std::size_t... (because you cannot store this directly in a variable).

The gen_seq is a "recursive" structure that build gen_seq<N, ...> using gen_seq<N - 1, ...>. How does it work:

gen_seq<3>: gen_seq<2, 2>
gen_seq<2, 2>: gen_seq<1, 1, 2>
gen_seq<1, 1, 2>: gen_seq<0, 0, 1, 2>
gen_seq<0, 0, 1, 2>: seq<0, 1, 2> // Specialized case

As you can see, since you inherit gen_seq<N - 1, N - 1, ...>, the last inheritance from seq has the 0 value, which you do not want. Since what is "send" to seq is the variadic std::size_t..., you want to avoid having 0, so you change to gen_seq<N - 1, N, ...>:

gen_seq<3>: gen_seq<2, 3>
gen_seq<2, 3>: gen_seq<1, 2, 3>
gen_seq<1, 2, 3>: gen_seq<0, 1, 2, 3>
gen_seq<0, 1, 2, 3>: seq<1, 2, 3> // Specialized case

Now, when you call _coeffs(gen_seq<N>{}), you allow the compiler to deduce the template argument NS of _coeffs. From this, you can use NS in _coeffs to do pack expansion:

{{coeff(NS)...}};
Holt
  • 36,600
  • 7
  • 92
  • 139
  • I will get these changes in place as soon as I can. As a side not I plan to implement "max_factorial" as a function that does a compile time check to detect at which N the factorial overflows on a given system hopefully making if much more portable. – Alex Zywicki Jun 28 '16 at 12:26
  • 1
    @AlexZywicki You could do a factorial on `double` instead of `std::intmax_t`, this would allow you to go to much bigger factorial (you can probably go up to `factorial(170)` with double). You could also try to do the "factorial" on the inverse of the number, e.g. `1/1 * 1/2 * 1/3 * 1/4 * ...` instead of `1/(1 * 2 * 3 * ...)`, you may get better accuracy (did not test thought, I am not 100% sure). – Holt Jun 28 '16 at 12:30
  • I will give that a try! Im interested to see if that would yield the correct results. That would allow me to increase the accuracy of my Taylor series significantly – Alex Zywicki Jun 28 '16 at 13:33
  • @AlexZywicki I did some testing - Using `double` for factorial, I can go up to `factorial(50)` using clang, which do improve the results I get using `intmax_t` which can only go up to `factorial(10)`. Using inverse factorial does not yield better results... Actually, with `factorial(50)`, the difference between `math::sin(x)` and `std::sin(x)` is about `1e-16` which is very close to the `double` epsilon (~2.0e-16) so I am pretty sure you cannot get better approximation. – Holt Jun 28 '16 at 14:08
  • Using what range of X? I usually apply a transform to X before the calculations that will keep X closer to zero where the approximation is more accurate but still yeilds the correct results – Alex Zywicki Jun 28 '16 at 15:16
  • @AlexZywicki I tried between -pi up to +pi with a step of 0.1, "worst" values where near +/- pi (~5e-15) as expected. – Holt Jun 28 '16 at 15:17
  • Let us [continue this discussion in chat](http://chat.stackoverflow.com/rooms/115861/discussion-between-alex-zywicki-and-holt). – Alex Zywicki Jun 28 '16 at 17:04
  • 1
    I would argue that `gen_seq` is correctly written, but incorrectly used. – T.C. Jun 29 '16 at 04:24
  • @T.C. You mean keeping a sequence from `0` to `N - 1` and changing the formula in coeff to `(n%2 ? -1 : 1) * inverse((T)factorial((2 * n)+1))`? – Holt Jun 29 '16 at 05:46
  • 1
    Or do `{{coeff(NS + 1) ...}}` instead. – T.C. Jun 29 '16 at 05:55
  • @T.C. Oh... I did not think about that ;) – Holt Jun 29 '16 at 06:00