Math.hpp

Go to the documentation of this file.

//
//  CeresEngine - A game development framework
//
//  Created by Rogiel Sulzbach.
//  Copyright (c) 2018-2022 Rogiel Sulzbach. All rights reserved.
//
 
#pragma once
 
#include "FloatingPoint.hpp"
#include "Matrix.hpp"
#include "Quaternion.hpp"
#include "Vector.hpp"
 
#include "CeresEngine/Foundation/Container/Array.hpp"
 
CE_DISABLE_WARNING_PUSH
CE_DISABLE_WARNING_DEPRECATED_VOLATILE
CE_DISABLE_WARNING_DEPRECATED_DECLARATIONS
#include <glm/common.hpp>
#include <glm/exponential.hpp>
#include <glm/geometric.hpp>
#include <glm/gtc/constants.hpp>
#include <glm/gtc/matrix_inverse.hpp>
#include <glm/gtc/matrix_transform.hpp>
#include <glm/trigonometric.hpp>
 
#include <glm/vector_relational.hpp>
 
#define GLM_ENABLE_EXPERIMENTAL
#include <glm/ext/matrix_clip_space.hpp>
#include <glm/gtc/matrix_transform.hpp>
#include <glm/gtc/type_ptr.hpp>
#include <glm/gtx/matrix_decompose.hpp>
#include <glm/gtx/matrix_major_storage.hpp>
#include <glm/gtx/matrix_transform_2d.hpp>
#include <glm/gtx/norm.hpp>
#include <glm/gtx/normal.hpp>
 
#include <glm/gtx/rotate_vector.hpp>
#include <glm/gtx/string_cast.hpp>
CE_DISABLE_WARNING_POP
 
#include <algorithm>
#include <cmath>
 
#define IMPORT_FROM_GLM(Symbol)          \
    namespace CeresEngine::inline Math { \
        using glm::Symbol;               \
    }
 
// Common functions
IMPORT_FROM_GLM(abs)
IMPORT_FROM_GLM(sign)
IMPORT_FROM_GLM(floor)
IMPORT_FROM_GLM(trunc)
IMPORT_FROM_GLM(round)
IMPORT_FROM_GLM(roundEven)
IMPORT_FROM_GLM(ceil)
IMPORT_FROM_GLM(mod)
IMPORT_FROM_GLM(min)
IMPORT_FROM_GLM(max)
IMPORT_FROM_GLM(clamp)
IMPORT_FROM_GLM(mix)
IMPORT_FROM_GLM(step)
IMPORT_FROM_GLM(smoothstep)
IMPORT_FROM_GLM(isnan)
IMPORT_FROM_GLM(isinf)
IMPORT_FROM_GLM(fma) //< a///b + c.
IMPORT_FROM_GLM(frexp) //< significand///exp(2, exponent)
IMPORT_FROM_GLM(ldexp) //< significand///exp(2, exponent)
 
// Exponential functions
IMPORT_FROM_GLM(pow)
IMPORT_FROM_GLM(exp)
IMPORT_FROM_GLM(log)
IMPORT_FROM_GLM(exp2)
IMPORT_FROM_GLM(log2)
IMPORT_FROM_GLM(sqrt)
IMPORT_FROM_GLM(inversesqrt)
 
// Trigonometric functions
IMPORT_FROM_GLM(radians)
IMPORT_FROM_GLM(degrees)
IMPORT_FROM_GLM(sin)
IMPORT_FROM_GLM(cos)
IMPORT_FROM_GLM(tan)
IMPORT_FROM_GLM(asin)
IMPORT_FROM_GLM(acos)
IMPORT_FROM_GLM(atan)
IMPORT_FROM_GLM(sinh)
IMPORT_FROM_GLM(cosh)
IMPORT_FROM_GLM(tanh)
IMPORT_FROM_GLM(asinh)
IMPORT_FROM_GLM(acosh)
IMPORT_FROM_GLM(atanh)
 
// Geometric functions
IMPORT_FROM_GLM(length)
IMPORT_FROM_GLM(length2)
IMPORT_FROM_GLM(distance)
IMPORT_FROM_GLM(dot)
IMPORT_FROM_GLM(cross)
IMPORT_FROM_GLM(normalize)
IMPORT_FROM_GLM(faceforward)
IMPORT_FROM_GLM(reflect)
IMPORT_FROM_GLM(refract)
 
// Vector relational functions
IMPORT_FROM_GLM(lessThan);
IMPORT_FROM_GLM(lessThanEqual);
IMPORT_FROM_GLM(greaterThan);
IMPORT_FROM_GLM(greaterThanEqual);
IMPORT_FROM_GLM(equal);
IMPORT_FROM_GLM(notEqual);
IMPORT_FROM_GLM(any);
IMPORT_FROM_GLM(all);
IMPORT_FROM_GLM(not_)
 
// Matrix operations
IMPORT_FROM_GLM(inverse)
IMPORT_FROM_GLM(transpose)
IMPORT_FROM_GLM(inverseTranspose)
IMPORT_FROM_GLM(identity)
 
// Matrix transformations functions
IMPORT_FROM_GLM(translate)
IMPORT_FROM_GLM(rotate)
IMPORT_FROM_GLM(scale)
 
// Projections and transform operations functions
IMPORT_FROM_GLM(unProject)
IMPORT_FROM_GLM(pickMatrix)
IMPORT_FROM_GLM(lookAt)
 
// String casting functions
// IMPORT_FROM_GLM(to_string)
IMPORT_FROM_GLM(decompose)
 
#undef IMPORT_FROM_GLM
 
namespace CeresEngine::inline Math {
 
    static constexpr double PI = 3.14159265358979323846;
 
    template<typename T> [[nodiscard]] TMatrix4<T> inline constexpr flipProjectionY(TMatrix4<T> proj) noexcept {
        // Flip Y in clipspace. X = -1, Y = -1 is topLeft in Vulkan.
        proj[1][1] *= T(-1);
        return proj;
    }
 
    template<typename T> [[nodiscard]] inline constexpr TMatrix4<T> frustum(T left, T right, T bottom, T top, T zNear, T zFar) noexcept {
        return flipProjectionY(glm::frustum(left, right, bottom, top, zNear, zFar));
    }
 
    template<typename T> [[nodiscard]] inline constexpr TMatrix4<T> perspective(T fovy, T aspect, T zNear, T zFar) noexcept {
        return flipProjectionY(glm::perspective(fovy, aspect, zNear, zFar));
    }
 
    template<typename T> [[nodiscard]] inline constexpr TMatrix4<T> perspectiveFov(T fov, T width, T height, T zNear, T zFar) noexcept {
        return flipProjectionY(glm::perspectiveFov(fov, width, height, zNear, zFar));
    }
 
    template<typename T> [[nodiscard]] inline constexpr TMatrix4<T> infinitePerspective(T fovy, T aspect, T zNear) noexcept {
        return flipProjectionY(glm::infinitePerspective(fovy, aspect, zNear));
    }
 
    template<typename T> [[nodiscard]] inline constexpr TMatrix4<T> tweakedInfinitePerspective(T fovy, T aspect, T zNear) noexcept {
        return flipProjectionY(glm::tweakedInfinitePerspective(fovy, aspect, zNear));
    }
 
    template<typename T> [[nodiscard]] inline constexpr TMatrix4<T> tweakedInfinitePerspective(T fovy, T aspect, T zNear, T ep) noexcept {
        return flipProjectionY(glm::tweakedInfinitePerspective(fovy, aspect, zNear, ep));
    }
 
    template<typename T> [[nodiscard]] inline constexpr TMatrix4<T> ortho(T left, T right, T bottom, T top, T zNear, T zFar) {
        return flipProjectionY(glm::ortho(left, right, bottom, top, zNear, zFar));
    }
 
    // ---------------------------------------------------------------------------------------------
 
    template<typename T> inline double smoothStep(T val1, T val2, T t) {
        t = clamp((t - val1) / (val2 - val1), 0.0, 1.0);
        return t * t * (3.0 - 2.0 * t);
    }
 
    template<typename T> inline T quintic(T val) { return val * val * val * (val * (val * 6.0 - 15.0) + 10.0); }
 
    template<typename T, typename U> inline auto cubic(T val1, T val2, T val3, T val4, U f) {
        auto t = (val4 - val3) - (val1 - val2);
        return f * f * f * t + f * f * ((val1 - val2) - t) + f * (val3 - val1) + val2;
    }
 
    // ---------------------------------------------------------------------------------------------
 
    template<typename T, glm::length_t D>
    inline bool approxEquals(const TVector<D, T>& a, const TVector<D, T>& b, T tolerance = std::numeric_limits<T>::epsilon()) {
        TVector<D, bool> result;
        for(glm::length_t i = 0; i < D; i++) {
            result[i] = approxEquals(a[i], b[i], tolerance);
        }
        return all(result);
    }
 
    template<typename T, glm::length_t D> inline bool approxEquals(const TVector<D, T>& a, T b, T tolerance = std::numeric_limits<T>::epsilon()) {
        return approxEquals(a, TVector<D, T>(b), tolerance);
    }
 
    template<typename T> inline bool approxEquals(const TQuaternion<T>& a, const TQuaternion<T>& b, T tolerance = std::numeric_limits<T>::epsilon()) {
        return approxEquals(a.x, b.x, tolerance) && approxEquals(a.y, b.y, tolerance) && approxEquals(a.z, b.z, tolerance) && approxEquals(a.w, b.w, tolerance);
    }
 
    // ---------------------------------------------------------------------------------------------
 
    template<typename T> inline T lerp(T t, T min, T max) { return (1.0 - t) * min + t * max; }
 
    template<typename T> inline T invLerp(T val, T min, T max) { return clamp((val - min) / std::max(max - min, 0.0001), T(0), T(1)); }
 
    template<typename A, typename B> inline std::common_type_t<A, B> gcd(const A& a, const B& b) { return (b == 0) ? a : gcd(b, a % b); }
 
    template<typename A, typename B> inline std::common_type_t<A, B> lcm(const A& a, const B& b) { return (a * b) / gcd(a, b); }
 
    // ---------------------------------------------------------------------------------------------
 
    template<typename T> inline int solveLinear(T A, T B, T* roots) {
        if(!approxEquals(A, T(0))) {
            roots[0] = -B / A;
            return 1;
        }
 
        roots[0] = T(0);
        return 1;
    }
 
    template<typename T> inline int solveQuadratic(T A, T B, T C, T* roots) {
        if(!approxEquals(A, T(0))) {
            T p = B / (2 * A);
            T q = C / A;
            T D = p * p - q;
 
            if(!approxEquals(D, T(0))) {
                if(D < T(0)) {
                    return 0;
                }
 
                T sqrtD = sqrt(D);
                roots[0] = sqrtD - p;
                roots[1] = -sqrtD - p;
 
                return 2;
            }
            roots[0] = -p;
            roots[1] = -p;
 
            return 1;
        }
        return solveLinear(B, C, roots);
    }
 
    template<typename T> inline int solveCubic(T A, T B, T C, T D, T* roots) {
        static constexpr T THIRD = (1 / T(3));
 
        T invA = 1 / A;
        A = B * invA;
        B = C * invA;
        C = D * invA;
 
        T sqA = A * A;
        T p = THIRD * (-THIRD * sqA + B);
        T q = (T(0.5)) * ((2 / T(27)) * A * sqA - THIRD * A * B + C);
 
        T cbp = p * p * p;
        D = q * q + cbp;
 
        int numRoots = 0;
        if(!approxEquals(D, T(0))) {
            if(D < 0.0) {
                T phi = THIRD * ::acos(-q / sqrt(-cbp));
                T t = 2 * sqrt(-p);
 
                roots[0] = t * cos(phi);
                roots[1] = -t * cos(phi + PI * THIRD);
                roots[2] = -t * cos(phi - PI * THIRD);
 
                numRoots = 3;
            } else {
                T sqrtD = sqrt(D);
                T u = cbrt(sqrtD + fabs(q));
 
                if(q > T(0)) {
                    roots[0] = -u + p / u;
                } else {
                    roots[0] = u - p / u;
                }
 
                numRoots = 1;
            }
        } else {
            if(!approxEquals(q, T(0))) {
                T u = cbrt(-q);
                roots[0] = 2 * u;
                roots[1] = -u;
 
                numRoots = 2;
            } else {
                roots[0] = 0.0;
                numRoots = 1;
            }
        }
 
        T sub = THIRD * A;
        for(int i = 0; i < numRoots; i++) {
            roots[i] -= sub;
        }
 
        return numRoots;
    }
 
    template<typename T> inline int solveQuartic(T A, T B, T C, T D, T E, T* roots) {
        T invA = 1 / A;
        A = B * invA;
        B = C * invA;
        C = D * invA;
        D = E * invA;
 
        T sqA = A * A;
        T p = -(3 / T(8)) * sqA + B;
        T q = (1 / T(8)) * sqA * A - (T)0.5 * A * B + C;
        T r = -(3 / T(256)) * sqA * sqA + (1 / (T)16) * sqA * B - (1 / (T)4) * A * C + D;
 
        int numRoots = 0;
        if(!approxEquals(r, T(0))) {
            T cubicA = 1;
            T cubicB = -T(0.5) * p;
            T cubicC = -r;
            T cubicD = T(0.5) * r * p - (1 / (T)8) * q * q;
 
            solveCubic(cubicA, cubicB, cubicC, cubicD, roots);
            T z = roots[0];
 
            T u = z * z - r;
            T v = 2 * z - p;
 
            if(approxEquals(u, T(0))) {
                u = 0;
            } else if(u > 0) {
                u = sqrt(u);
            } else {
                return 0;
            }
 
            if(approxEquals(v, T(0))) {
                v = 0;
            } else if(v > 0) {
                v = sqrt(v);
            } else {
                return 0;
            }
 
            T quadraticA = 1;
            T quadraticB = q < 0 ? -v : v;
            T quadraticC = z - u;
 
            numRoots = solveQuadratic(quadraticA, quadraticB, quadraticC, roots);
 
            quadraticA = 1;
            quadraticB = q < 0 ? v : -v;
            quadraticC = z + u;
 
            numRoots += solveQuadratic(quadraticA, quadraticB, quadraticC, roots + numRoots);
        } else {
            numRoots = solveCubic(q, p, T(0), T(1), roots);
            roots[numRoots++] = 0;
        }
 
        T sub = (1 / T(4)) * A;
        for(int i = 0; i < numRoots; i++) {
            roots[i] -= sub;
        }
 
        return numRoots;
    }
 
    // ---------------------------------------------------------------------------------------------
 
    template<class T> inline T cubicHermite(const double t, const T& pointA, const T& pointB, const T& tangentA, const T& tangentB) {
        const double t2 = t * t;
        const double t3 = t2 * t;
 
        double a = 2 * t3 - 3 * t2 + 1;
        double b = t3 - 2 * t2 + t;
        double c = -2 * t3 + 3 * t2;
        double d = t3 - t2;
 
        return T(a * pointA + b * tangentA + c * pointB + d * tangentB);
    }
 
    template<class T> inline T cubicHermiteD1(const double t, const T& pointA, const T& pointB, const T& tangentA, const T& tangentB) {
        const double t2 = t * t;
 
        double a = 6 * t2 - 6 * t;
        double b = 3 * t2 - 4 * t + 1;
        double c = -6 * t2 + 6 * t;
        double d = 3 * t2 - 2 * t;
 
        return T(a * pointA + b * tangentA + c * pointB + d * tangentB);
    }
 
    template<class T> inline Array<T, 4> cubicHermiteCoefficients(const T& pointA, const T& pointB, const T& tangentA, const T& tangentB) {
        T diff = pointA - pointB;
        return {2 * diff + tangentA + tangentB, -3 * diff - 2 * tangentA - tangentB, tangentA, pointA};
    }
 
    template<class T, typename L>
    inline Array<T, 4> cubicHermiteCoefficients(const T& pointA, const T& pointB, const T& tangentA, const T& tangentB, L length) {
        L length2 = length * length;
        L invLength2 = 1.0 / length2;
        L invLength3 = 1.0 / (length2 * length);
 
        T scaledTangentA = T(tangentA * length);
        T scaledTangentB = T(tangentB * length);
 
        T diff = pointA - pointB;
 
        return {T((2.0 * diff + scaledTangentA + scaledTangentB) * invLength3), T((-3.0 * diff - 2.0 * scaledTangentA - scaledTangentB) * invLength2), tangentA,
            pointA};
 
        //        double length2 = length * length;
        //        double invLength2 = 1.0 / length2;
        //        double invLength3 = 1.0 / (length2 * length);
 
        //        double invLength2 = 1.0;
        //        double invLength3 = 1.0;
        //
        //        T scaledTangentA = T(tangentA * length);
        //        T scaledTangentB = T(tangentB * length);
        //
        //        T diff = pointA - pointB;
        //
        //        return {
        //                T((2.0 * diff + scaledTangentA + scaledTangentB) *
        //                invLength3), T((-3.0 * diff - 2.0 * scaledTangentA -
        //                scaledTangentB) * invLength2), tangentA, pointA
        //        };
    }
 
    template<typename T, typename Integrand> inline T rombergIntegration(T a, T b, int order, Integrand&& integrand) {
        T h[order + 1];
        T r[order + 1][order + 1];
 
        for(int i = 1; i < order + 1; ++i) {
            h[i] = (b - a) / pow(2, i - 1);
        }
 
        r[1][1] = h[1] / 2 * (integrand(a) + integrand(b));
 
        for(int i = 2; i < order + 1; ++i) {
            T coeff = 0;
            for(int k = 1; k <= pow(2, i - 2); ++k) {
                coeff += integrand(a + (2 * k - 1) * h[i]);
            }
 
            r[i][1] = 0.5 * (r[i - 1][1] + h[i - 1] * coeff);
        }
 
        for(int i = 2; i < order + 1; ++i) {
            for(int j = 2; j <= i; ++j) {
                r[i][j] = r[i][j - 1] + (r[i][j - 1] - r[i - 1][j - 1]) / (pow(4, j - 1) - 1);
            }
        }
 
        return r[order][order];
    }
 
    template<typename T, typename Integrand> inline T gaussianQuadrature(T a, T b, T* roots, T* coefficients, Integrand&& integrand) {
        const T half = T(0.5);
        const T radius = half * (b - a);
        const T center = half * (b + a);
        T res = T(0);
 
        for(int i = 0; i < sizeof(roots) / sizeof(*roots); ++i) {
            res += coefficients[i] * integrand(radius * roots[i] + center);
        }
 
        res *= radius;
 
        return res;
    }
 
    template<typename T> inline T haltonSequence(int index, int base) {
        T output = T(0.0);
        T invBase = T(1.0) / base;
        T frac = invBase;
        while(index > 0) {
            output += (index % base) * frac;
            index /= base;
            frac *= invBase;
        }
 
        return output;
    }
 
    // Function substitutes
    template<typename T> T angle(const TVector2<T>& v1, const TVector2<T>& v2) {
        using std::acos;
        using std::numeric_limits;
        using std::sqrt;
 
        const T len = sqrt(dot(v1, v1) * dot(v2, v2));
        if(len <= std::numeric_limits<T>::epsilon())
            return T{0};
        return acos(clamp(dot(v1, v2) / len, T{-1}, T{1}));
    }
 
    template<typename T> T angle(const TVector3<T>& v1, const TVector3<T>& v2) {
        using std::acos;
        using std::numeric_limits;
        using std::sqrt;
 
        const T len = sqrt(dot(v1, v1) * dot(v2, v2));
        if(len <= std::numeric_limits<T>::epsilon())
            return T{0};
        return acos(clamp(dot(v1, v2) / len, T{-1}, T{1}));
    }
 
    template<typename T> TVector2<T> cross(const TVector2<T>& v) { return TVector2<T>(-v.y, v.x); }
 
    template<typename T> TVector3<T> transform(const TQuaternion<T>& q, const TVector3<T>& v) {
        const TVector3<T> temp = T{2.0} * cross(TVector3<T>(q.x, q.y, q.z), v);
        return v + q.w * temp + cross(TVector3<T>(q.x, q.y, q.z), temp);
    }
 
    template<typename T> TQuaternion<T> qrotate(const T& angle, const TVector3<T>& axis) {
        const T a = angle / T{2.0};
        return TQuaternion<T>{cos(a), sin(a) * axis};
    }
 
    template<typename T> TVector3<T> normal(const TVector3<T>& p1, const TVector3<T>& p2, const TVector3<T>& p3) { return normalize(cross(p2 - p1, p3 - p1)); }
 
    template<typename T, typename TI, typename TS>
    TVector3<T> project(const TVector3<T>& v, const TMatrix4<T>& modelViewProj, const TVector2<TI>& viewportOrigin, const TVector2<TS>& viewportSize) {
        TVector4<T> in = modelViewProj * TVector4<T>{v, static_cast<T>(1)};
 
        in[0] /= in[3];
        in[1] /= in[3];
        in[2] /= in[3];
 
        const T half = static_cast<T>(0.5);
 
        in[0] = in[0] * half + half;
        in[1] = in[1] * half + half;
        in[2] = in[2] * half + half;
 
        in[0] = in[0] * static_cast<T>(viewportSize[0]) + static_cast<T>(viewportOrigin[0]);
        in[1] = in[1] * static_cast<T>(viewportSize[1]) + static_cast<T>(viewportOrigin[1]);
 
        return TVector3<T>{in};
    }
 
    template<typename T> TMatrix4<T> ortho2D(const T& left, const T& right, const T& bottom, const T& top) {
        return ortho(left, right, bottom, top, T{-1}, T{1});
    }
 
    template<typename T> TVector3<T> slerp(const TVector3<T>& v1, const TVector3<T>& v2, const T& a) {
        using std::sin;
        const T theta = angle(v1, v2);
        const T sine = sin(theta);
        return sin((T{1} - a) * theta) / sine * v1 + sin(a * theta) / sine * v2;
    }
 
    template<typename T> TMatrix4<T> rotate(const TVector3<T>& angle) {
        const T sy = sin(angle[2]);
        const T cy = cos(angle[2]);
        const T sp = sin(angle[1]);
        const T cp = cos(angle[1]);
        const T sr = sin(angle[0]);
        const T cr = cos(angle[0]);
 
        const T data[16] = {cp * cy, sr * sp * cy + cr * -sy, cr * sp * cy + -sr * -sy, T{0}, cp * sy, sr * sp * sy + cr * cy, cr * sp * sy + -sr * cy, T{0},
            -sp, sr * cp, cr * cp, T{0}, T{0}, T{0}, T{0}, T{1}};
        return glm::rowMajor4(glm::make_mat4(data));
    }
 
    template<typename T> TMatrix3<T> rotate(const T& angle) {
        const T s = sin(angle);
        const T c = cos(angle);
 
        const T data[9] = {c, -s, T{0}, s, c, T{0}, T{0}, T{0}, T{1}};
 
        return glm::rowMajor3(glm::make_mat3(data));
    }
 
    template<typename T> TVector2<T> transform(const TMatrix3<T>& m, const TVector2<T>& v) { return TVector2<T>(m * TVector3<T>(v, 1.0)); }
 
    namespace detail {
 
        template<typename VecT, typename T> VecT bezierImpl(const VecT* p, const int n, T t1, T t2, int stride = 1) {
            if(n == 1)
                return *p;
            if(n == 2)
                return t1 * p[0] + t2 * p[stride];
            return t1 * bezierImpl(p, n - 1, t1, t2, stride) + t2 * bezierImpl(p + stride, n - 1, t1, t2, stride);
        }
 
    } // namespace detail
 
    template<int D, typename T> TVector2<T> bezier(const TVector2<T> (&p)[D], T t) {
        static_assert(D > 0, "At least one control point needed.");
        return detail::bezierImpl(&p[0], D, static_cast<T>(1) - t, t);
    }
 
    template<int D, typename T> TVector3<T> bezier(const TVector3<T> (&p)[D], T t) {
        static_assert(D > 0, "At least one control point needed.");
        return detail::bezierImpl(&p[0], D, static_cast<T>(1) - t, t);
    }
 
    template<int D0, int D1, typename T> TVector3<T> bezier2(const TVector3<T> (&p)[D1][D0], const TVector2<T>& t) {
        static_assert(D0 > 0, "At least one control point needed.");
        static_assert(D1 > 0, "At least one control point needed.");
 
        TVector3<T> temp[D1];
        for(int i = 0; i < D1; ++i) {
            temp[i] = bezier(p[i], t[0]);
        }
        return bezier(temp, t[1]);
    }
 
    namespace detail {
 
        template<int O, int D, typename VecT, typename T> struct bezierDerivativeImpl {
            static VecT calc(const VecT (&p)[D], T t) {
                VecT temp[D - 1];
                for(int i = 0; i < D - 1; ++i) {
                    temp[i] = static_cast<T>(D - 1) * (p[i + 1] - p[i]);
                }
                return bezierDerivativeImpl<O - 1, D - 1, VecT, T>::calc(temp, t);
            }
        };
 
        template<int D, typename VecT, typename T> struct bezierDerivativeImpl<0, D, VecT, T> {
            static VecT calc(const VecT (&p)[D], T t) { return bezier(p, t); }
        };
 
        template<typename VecT, typename T> struct bezierDerivativeImpl<0, 1, VecT, T> {
            static VecT calc(const VecT (&p)[1], T t) { return bezier(p, t); }
        };
 
        template<int O, typename VecT, typename T> struct bezierDerivativeImpl<O, 1, VecT, T> {
            static VecT calc(const VecT (&)[1], T) { return VecT{static_cast<T>(0)}; }
        };
 
    } // namespace detail
 
    template<int O, int D, typename T> TVector2<T> bezierDerivative(const TVector2<T> (&p)[D], T t) {
        static_assert(O > 0, "The derivative order must be at least one.");
        static_assert(D > 0, "At least one control point needed.");
        return detail::bezierDerivativeImpl<O, D, TVector2<T>, T>::calc(p, t);
    }
 
    template<int O, int D, typename T> TVector3<T> bezierDerivative(const TVector3<T> (&p)[D], T t) {
        static_assert(O > 0, "The derivative order must be at least one.");
        static_assert(D > 0, "At least one control point needed.");
        return detail::bezierDerivativeImpl<O, D, TVector3<T>, T>::calc(p, t);
    }
 
    template<int O, int D0, int D1, typename T> TMatrix2x3<T> bezier2Jacobian(const TVector3<T> (&p)[D1][D0], const TVector2<T>& t) {
        static_assert(O > 0, "Order of the Jacobian must be at least one.");
        static_assert(D0 > 0, "At least one control point needed.");
        static_assert(D1 > 0, "At least one control point needed.");
 
        TVector3<T> temp0[D0];
        for(int i = 0; i < D0; ++i) {
            temp0[i] = detail::bezierImpl(&p[0][i], D1, static_cast<T>(1) - t[1], t[1], D0);
        }
 
        TVector3<T> temp1[D1];
        for(int i = 0; i < D1; ++i) {
            temp1[i] = bezier(p[i], t[0]);
        }
 
        return TMatrix2x3<T>{bezierDerivative<O>(temp0, t[0]), bezierDerivative<O>(temp1, t[1])};
    }
 
} // namespace CeresEngine::inline Math
 
namespace CeresEngine::inline Math {
 
    namespace Constant {
 
#define IMPORT_GLM_CONSTANT_NAMED(Symbol, Name) template<typename T> inline const T Name = glm::Symbol<T>()
#define IMPORT_GLM_CONSTANT(Symbol) IMPORT_GLM_CONSTANT_NAMED(Symbol, Symbol)
 
        IMPORT_GLM_CONSTANT(epsilon);
 
        IMPORT_GLM_CONSTANT_NAMED(zero, _0);
 
        IMPORT_GLM_CONSTANT_NAMED(one, _1);
 
        IMPORT_GLM_CONSTANT(pi);
 
        IMPORT_GLM_CONSTANT_NAMED(two_pi, _2pi);
 
        IMPORT_GLM_CONSTANT_NAMED(root_pi, sqrtPi);
 
        IMPORT_GLM_CONSTANT_NAMED(half_pi, pi_2);
 
        IMPORT_GLM_CONSTANT_NAMED(three_over_two_pi, _3pi_2);
 
        IMPORT_GLM_CONSTANT_NAMED(quarter_pi, pi_4);
 
        IMPORT_GLM_CONSTANT_NAMED(one_over_pi, _1_pi);
 
        IMPORT_GLM_CONSTANT_NAMED(one_over_two_pi, _1_2pi);
 
        IMPORT_GLM_CONSTANT_NAMED(two_over_pi, _2_pi);
 
        IMPORT_GLM_CONSTANT_NAMED(four_over_pi, _4_pi);
 
        IMPORT_GLM_CONSTANT_NAMED(two_over_root_pi, _2_sqrtpi);
 
        IMPORT_GLM_CONSTANT_NAMED(one_over_root_two, _1_sqrt2);
 
        IMPORT_GLM_CONSTANT_NAMED(root_half_pi, sqrt_pi_2);
 
        IMPORT_GLM_CONSTANT_NAMED(root_two_pi, sqrt_2pi);
 
        IMPORT_GLM_CONSTANT(e);
 
        IMPORT_GLM_CONSTANT(euler);
 
        IMPORT_GLM_CONSTANT_NAMED(root_two, sqrt2);
 
        IMPORT_GLM_CONSTANT_NAMED(root_three, sqrt3);
 
        IMPORT_GLM_CONSTANT_NAMED(root_five, sqrt5);
 
        IMPORT_GLM_CONSTANT_NAMED(golden_ratio, goldenRatio);
 
    } // namespace Constant
 
    using Constant::pi;
} // namespace CeresEngine::inline Math