Some simple constexpr math function (C++11)

ConstMath.hxx

/* This algorithm is based on a very popular algorithm among C programmers, used to calculate power of
   a number. I used a tail recursivity, as I wanted it to be C++11 compliant, and thus I had no
   unconstrained constexpr, thus no loops.
   The algorithm is based on this simple observation :
   for instance, we want to compute 6^7. We have that 6^7 = 6*(6^3*6^3) = 6*(6*(6^2)*6*(6^2))
   We can see that there's some pattern in the expression. To see things more clearly, we can rewrite
   the expression as follow:
   6^7 = 6*6^2*(6^2)*(6^2)
   Actually, the algorithm is doing the opposite order of the decomposition. That make sens as
   we are starting from the base. To summarize the algorithm way of functionning, I would say the following :
   In every case, do base*base
   If the exponent is odd, do result*base
   We can see that we get an extra "base" multiplication when the exponent is odd. This is caused by the flooring
   of the exponent, if we do not add this extra base multiplication, we would miss one multiplication
   every odd exponent value.
   So, enough chit chat, let follow this algorithm step by step with our previous example :
   1 : base = 6, exponent = 7, result = 1
       return powImpl(6*6, 7/2, 1*6)
   2 : base = 6^2, exponent = 3, result = 6 // See, we dropped one exponent, thus the 1*6 in the previous step
       return powImpl(6^2*6^2, exponent = 3/2, 6*6^2)
	   // If you look more closely, above we have exaclty the decomposition we wanted
   3 : base = 6^4, exponent = 1, result = 6^3 // Now you should really see how that's going :)
       return powImpl(6^4, exponent = 1/2, 6^3*6^4)
   4 : base = 6^8, exponent= 0, result = 6^7
       return result, so return 6^7

   I could dwelve more into theory for this algorithm, but I'm sure that if you're interested enough, you'll find
   on the internet mountains of proofs made by people way more clever and competent than me :)
   Anyway, I hope that this explanation helped you to understand this algorithm if you didn't now it beforehand */

template<typename T,
		typename std::enable_if<std::is_integral<T>::value>::type* = nullptr>
constexpr T powImpl(T base, size_t exponent, T result)
{
	return exponent > 0 ? powImpl(base*base, exponent/2, exponent&1 ? result*base : result) : result;
}

template<typename T,
		 typename std::enable_if<std::is_integral<T>::value>::type* = nullptr>
constexpr T pow(T base, size_t exponent)
{
	return exponent > 1 ? powImpl(base, exponent, T(1)) : base;
}

template<typename T,
		 typename std::enable_if<std::is_signed<T>::value>::type* = nullptr>
constexpr T signByte(T num)
{
	constexpr size_t shift = sizeof(T); // Sign byte is always the last;
	return num && (size_t(1) << (shift - 1));
}

template<typename T,
		 typename std::enable_if<!std::is_signed<T>::value>::type* = nullptr>
constexpr T signByte(T num)
{
	return 0;
}

template<typename T,
         typename std::enable_if<std::is_floating_point<T>::value>::type* = nullptr>
constexpr int64 floor(T num)
{
	//constexpr size_t sign = signByte(num);
	//ceil_impl(num);
	return int64(num + pow(2L, (sizeof(max_type) / 2) * 8)) - pow(2L, (sizeof(max_type) / 2) * 8);
}

/* Here, we use the fact that floor(-num) = -ceil(num)
   Indeed, if we have -5.5, floor(-5.5) = -ceil(5.5) = -6
   For those who need something more "formal", here is a simple assertion
   floor rounds toward infinity for negatives, and toward 0 for positives, while ceil
   rounds toward 0 for negatives, and toward infinity for positives.
   In fact, the minus sign invert the two propositions. Thus, -floor rounds toward 0 for negatives
   and toward infinity for positives. We thus proved that -floor = ceil.
   More with mathematic formulas ? (at this point I find it almost useless :) )
   Let floor:R->N be the function that satisfy that for all x€R, it exists n€N, and
   floor(x) <= x < floor(x) + 1
   and let ceil:R->N be the function that satisfy that for all x€R, it exists n€N, and
   ceil(x) - 1 < x <= ceil(x)
   thus, we have
   -floor(x) >= x > -floor(x) - 1
<=>-floor(x) - 1 < x <= -floor(x)

    I don't think we need any more proof :) (excuse me for my bad proof writing, english is not my "math" language)
	We can now write the floor function. If we take the ceil implementation, we get
	-(int64(-num+pow(2, sizeof(int64)) - pow(2, sizeof(int64))))
	= -int(-num+pow) + pow
	= pow - int(-num + pow)
	= pow - int(pow - num)
	 */
template<typename T,
	     typename std::enable_if<std::is_floating_point<T>::value>::type* = nullptr>
constexpr int64 ceil(T num)
{
	return pow(2L, (sizeof(max_type) / 2) * 8) - int64(pow(2L, (sizeof(max_type) / 2) * 8) - num);
}

template<typename T,
		 typename std::enable_if<std::is_arithmetic<T>::value && std::is_signed<T>::value>::type* = nullptr>
constexpr T abs(T num)
{
	return signByte(num) ? -num : num;
}

template<typename T,
		 typename std::enable_if<std::is_arithmetic<T>::value>::type* = nullptr>
constexpr T clip(T num, T min, T max)
{
	return std::max(min, std::min(num, max));
}

// Some demo examples of two type traits.
// The standards ones are called std::is_arithmetic and std::is_signed.
// This is just for as an exercise purpose, you should use the standards ones ;)
// (Like I did in the code above)
template<typename T>
struct is_number : std::conditional<std::is_integral<T>::value || std::is_floating_point<T>::value, std::true_type, std::false_type>::type
{};

template<typename T, bool = std::is_number<T>::value>
struct is_signed: std::integral_constant<bool, T(-1) < T(0)>
{};

template<typename T>
struct is_signed<T, false>: std::false_type
{};

Any clue as to what max_type is supposed to be?

Any clue as to what max_type is supposed to be?

Hi. I'm fairly sure that what I meant to write was the biggest of the two types, so, something along those lines :

template<class T1, class T2>
using max_type = std::conditional_t<sizeof(T1) >= sizeof(T2), T1, T2>;

Although upon re-reading this admittedly old gist, I do believe it to be incorrect. Indeed, the trick is to make the following branchless :

template<typename T,
         typename std::enable_if<std::is_floating_point<T>::value>::type* = nullptr>
constexpr int64 floor(T num)
{
    int64_t res{};
    res = static_cast<int64_t>(num);
    if(res > num) --res;
    
    return res;
}

Hence, as we know that the floating point is signed, we shift it to a positive value (it could be already positive) using a shift of max_val / 2, we floor, and then we shift back to where it was, by doing the inverse shift. This allow us to always do the flooring in the positive space, avoiding the branch. But here is where it gets dicey : what is max_val.
If max_val is the max val for the integer, well, we expose ourselves to floating point overflow, and wrong results. Consider the following :

template<typename T,
         typename std::enable_if<std::is_floating_point<T>::value>::type* = nullptr>
constexpr int64_t floor(T num)
{
    using max_type = max_type<int64_t, T>;
    return int64_t(pow(2L, (sizeof(max_type) / 2) * 8) + num) - pow(2L, (sizeof(max_type) / 2) * 8);
}

int main()
{
    std::cout << floor(double(6.4)) << std::endl;
    std::cout << floor(float(6.4)) << std::endl;
    
    return 0;
}

First result will be 6, as expected, but the second will be 0... well, what obviously makes the most sense is to use the size of our input floating point type T, and assume that it will always be <= 8 bytes (reasonable assumption on most platforms).

template<typename T,
         typename std::enable_if<std::is_floating_point<T>::value>::type* = nullptr>
constexpr int64_t floor(T num)
{
    return int64_t(pow(2L, (sizeof(T) / 2) * 8) + num) - pow(2L, (T) / 2) * 8);
}

But here's the issue, the shift works well if we assume that floating points behave like integers, but they do absolutely NOT. Let's remember that the first integer that can not be represented exactly by a 32-bit floating point number is 16777217, which is quite a bit less than what 32-bit signed integers are capable of representing. Now, it's easy to get a value that will cause our function to fail, something close to this integer, let's say 16777215.

int main()
{
    std::cout << float(16777215) << std::endl;
    std::cout << floor(float(16777215)) << std::endl;
    
    return 0;
}

Yeah, despite being representable by a 32-bit floating point number, the floor of 16777215 is 16777216, not great. There's not real way to avoid that without delving into inline assembly or just using branches. Unless you're willing to deal with the aforementioned danger, replacing max_type by T will work, but you're always running the risk of getting into trouble with the above corner case. In production I would simply use a branch, ESPECIALLY if you're gonna evaluate them at compile time whenever you can!

Hope this helps.