/// Josh's 3D Math Library
/// A C++20 Library for 3D Math, Computer Graphics, and Scientific Computing.
/// Developed and Maintained by Josh O'Leary @ Redacted Software.
/// Special Thanks to William Tomasine II and Maxine Hayes.
/// (c) 2024 Redacted Software
/// This work is dedicated to the public domain.

/// @file J3ML.h
/// @desc Core mathematical and utility functions, concrete types, and math constants.
/// @edit 2024-07-05

#pragma once

#include <cstdint>
#include <cmath>
#include <string>
#include <cassert>
#include <vector>

/// This set of functions may be set to use lookup tables or SIMD operations.
/// If no options are set, they will default to using standard library implementation.
#undef USE_LOOKUP_TABLES /// Pre-computed lookup tables.
#undef USE_SSE           /// Streaming SIMD Extensions (x86)
#undef USE_NEON          /// ARM Vector Processing
#undef USE_AVX           /// Advanced Vector Extensions (x86)

/// TODO: Implement lookup tables.
/// TODO: Implement constexpr Trigonometric LUT generators that are parameterized (samples, samples-per-period, etc.)

#ifdef USE_LOOKUP_TABLES

#define LUT_SAMPLES 1024

#pragma region Trigonometric Lookup Tables

// Formula: sin(2*pi*t/T)
/** Generated using Dr LUT - Free Lookup Table Generator
  * https://github.com/ppelikan/drlut
  **/
// Formula: sin(2*pi*t/T)
const uint8_t u8_sin_lut[1024] = {
127,128,129,129,130,131,132,132,133,134,135,136,136,
137,138,139,139,140,141,142,143,143,144,145,146,146,
147,148,149,149,150,151,152,153,153,154,155,156,156,
157,158,159,159,160,161,162,162,163,164,165,165,166,
167,168,168,169,170,171,171,172,173,173,174,175,176,
176,177,178,178,179,180,181,181,182,183,183,184,185,
185,186,187,188,188,189,190,190,191,192,192,193,194,
194,195,196,196,197,198,198,199,199,200,201,201,202,
203,203,204,205,205,206,206,207,208,208,209,209,210,
211,211,212,212,213,213,214,215,215,216,216,217,217,
218,218,219,220,220,221,221,222,222,223,223,224,224,
225,225,226,226,227,227,228,228,229,229,229,230,230,
231,231,232,232,233,233,233,234,234,235,235,236,236,
236,237,237,238,238,238,239,239,239,240,240,240,241,
241,241,242,242,242,243,243,243,244,244,244,245,245,
245,245,246,246,246,247,247,247,247,248,248,248,248,
249,249,249,249,249,250,250,250,250,250,251,251,251,
251,251,251,252,252,252,252,252,252,252,253,253,253,
253,253,253,253,253,253,253,253,254,254,254,254,254,
254,254,254,254,254,254,254,254,254,254,254,254,254,
254,254,254,254,254,254,254,254,254,254,254,253,253,
253,253,253,253,253,253,253,253,253,252,252,252,252,
252,252,252,251,251,251,251,251,251,250,250,250,250,
250,249,249,249,249,249,248,248,248,248,247,247,247,
247,246,246,246,245,245,245,245,244,244,244,243,243,
243,242,242,242,241,241,241,240,240,240,239,239,239,
238,238,238,237,237,236,236,236,235,235,234,234,233,
233,233,232,232,231,231,230,230,229,229,229,228,228,
227,227,226,226,225,225,224,224,223,223,222,222,221,
221,220,220,219,218,218,217,217,216,216,215,215,214,
213,213,212,212,211,211,210,209,209,208,208,207,206,
206,205,205,204,203,203,202,201,201,200,199,199,198,
198,197,196,196,195,194,194,193,192,192,191,190,190,
189,188,188,187,186,185,185,184,183,183,182,181,181,
180,179,178,178,177,176,176,175,174,173,173,172,171,
171,170,169,168,168,167,166,165,165,164,163,162,162,
161,160,159,159,158,157,156,156,155,154,153,153,152,
151,150,149,149,148,147,146,146,145,144,143,143,142,
141,140,139,139,138,137,136,136,135,134,133,132,132,
131,130,129,129,128,127,126,125,125,124,123,122,122,
121,120,119,118,118,117,116,115,115,114,113,112,111,
111,110,109,108,108,107,106,105,105,104,103,102,101,
101,100, 99, 98, 98, 97, 96, 95, 95, 94, 93, 92, 92,
 91, 90, 89, 89, 88, 87, 86, 86, 85, 84, 83, 83, 82,
 81, 81, 80, 79, 78, 78, 77, 76, 76, 75, 74, 73, 73,
 72, 71, 71, 70, 69, 69, 68, 67, 66, 66, 65, 64, 64,
 63, 62, 62, 61, 60, 60, 59, 58, 58, 57, 56, 56, 55,
 55, 54, 53, 53, 52, 51, 51, 50, 49, 49, 48, 48, 47,
 46, 46, 45, 45, 44, 43, 43, 42, 42, 41, 41, 40, 39,
 39, 38, 38, 37, 37, 36, 36, 35, 34, 34, 33, 33, 32,
 32, 31, 31, 30, 30, 29, 29, 28, 28, 27, 27, 26, 26,
 25, 25, 25, 24, 24, 23, 23, 22, 22, 21, 21, 21, 20,
 20, 19, 19, 18, 18, 18, 17, 17, 16, 16, 16, 15, 15,
 15, 14, 14, 14, 13, 13, 13, 12, 12, 12, 11, 11, 11,
 10, 10, 10,  9,  9,  9,  9,  8,  8,  8,  7,  7,  7,
  7,  6,  6,  6,  6,  5,  5,  5,  5,  5,  4,  4,  4,
  4,  4,  3,  3,  3,  3,  3,  3,  2,  2,  2,  2,  2,
  2,  2,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,
  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
  0,  0,  0,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,
  1,  2,  2,  2,  2,  2,  2,  2,  3,  3,  3,  3,  3,
  3,  4,  4,  4,  4,  4,  5,  5,  5,  5,  5,  6,  6,
  6,  6,  7,  7,  7,  7,  8,  8,  8,  9,  9,  9,  9,
 10, 10, 10, 11, 11, 11, 12, 12, 12, 13, 13, 13, 14,
 14, 14, 15, 15, 15, 16, 16, 16, 17, 17, 18, 18, 18,
 19, 19, 20, 20, 21, 21, 21, 22, 22, 23, 23, 24, 24,
 25, 25, 25, 26, 26, 27, 27, 28, 28, 29, 29, 30, 30,
 31, 31, 32, 32, 33, 33, 34, 34, 35, 36, 36, 37, 37,
 38, 38, 39, 39, 40, 41, 41, 42, 42, 43, 43, 44, 45,
 45, 46, 46, 47, 48, 48, 49, 49, 50, 51, 51, 52, 53,
 53, 54, 55, 55, 56, 56, 57, 58, 58, 59, 60, 60, 61,
 62, 62, 63, 64, 64, 65, 66, 66, 67, 68, 69, 69, 70,
 71, 71, 72, 73, 73, 74, 75, 76, 76, 77, 78, 78, 79,
 80, 81, 81, 82, 83, 83, 84, 85, 86, 86, 87, 88, 89,
 89, 90, 91, 92, 92, 93, 94, 95, 95, 96, 97, 98, 98,
 99,100,101,101,102,103,104,105,105,106,107,108,108,
109,110,111,111,112,113,114,115,115,116,117,118,118,
119,120,121,122,122,123,124,125,125,126 };


#pragma endregion
#endif

#include <J3ML/Algorithm/Reinterpret.hpp>

/// Swaps two elements in-place without copying their data.
template <typename T>
void Swap(T &a, T &b)
{
    T temp = std::move(a);
    a = std::move(b);
    b = std::move(temp);
}

namespace J3ML {
    /// Clean symbolic names for integers of specific size.
    namespace SizedIntegralTypes {
        using u8 = uint8_t;
        using u16 = uint16_t;
        using u32 = uint32_t;
        using u64 = uint64_t;
        using s8 = int8_t;
        using s16 = int16_t;
        using s32 = int32_t;
        using s64 = int64_t;
    }
    //using namespace SizedIntegralTypes; // Bring into J3ML namespace.

    namespace SizedFloatTypes { // TODO: Use C++23 <stdfloat>
        using f16  = float;
        using f32 = float;
        using f64 = double;
        using f128 = long double;
    }
    //using namespace SizedFloatTypes;    // Bring into J3ML namespace.
}

using namespace J3ML::SizedIntegralTypes;
using namespace J3ML::SizedFloatTypes;

namespace J3ML::BitTwiddling {
    /// Parses a string of form "011101010" to a u32
    u32 BinaryStringToValue(const char* s);

    /// Returns the number of 1's set in the given value.
    //inline int CountBitsSet(u32 value);
}

namespace J3ML::Math {
    enum class Quadrant { I, II, III, IV };

    // Zero technically isn't a sign, but zero also isn't positive, or negative, so bite me.
    enum class Sign { ZERO, POSITIVE, NEGATIVE};

}

namespace J3ML::Math::Constants { // TODO: Consider double precision for these.
        /// sqrt(2pi) ^ -1
        constexpr float RecipSqrt2Pi = 0.3989422804014326779399460599343818684758586311649346576659258296706579258993018385012523339073069364;
        /// pi - https://www.mathsisfun.com/numbers/pi.html
        constexpr float Pi = 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679;
        constexpr float TwoPi = Pi*2.0;
        constexpr float PiOverTwo = Pi/2.0;
        constexpr float ThreePiOverTwo = 3.0*Pi/2.0;
        /// e - https://www.mathsisfun.com/numbers/e-eulers-number.html
        constexpr float EulersNumber = 2.7182818284590452353602874713526624977572470936999595749669676277240766303535475945713821785251664274;
        /// 2pi - The ratio of a circle's circumferecne to its radius, and the number of radians in one turn.
        constexpr float Tau = 6.28318530717958647692;
        /// sqrt(2)
        constexpr float PythagorasConstant = 1.41421356237309504880;
        /// sqrt(3)
        constexpr float TheodorusConstant = 1.73205080756887729352;
        /// Golden Ratio
        constexpr float Phi = 1.61803398874989484820;
        /// ln 2
        constexpr float NaturalLog2 = 0.6931471805599453094172321214581765680755001343602552541206800094933936219696947156058633269964186875;
        /// ln 10
        constexpr float NaturalLog10 = 2.3025850929940456840179914546843642076011014886287729760333279009675726096773524802359972050895982983;
        constexpr float Infinity = INFINITY;
        constexpr float NegativeInfinity = -INFINITY;
        constexpr float NotANumber = NAN;
}

namespace J3ML::Math {
    using namespace Constants; // Bring into J3ML::Math namespace.
}

namespace J3ML::Math::Functions {
         // TODO: Implement "Wrappers" for most standard math functions.
         // We want to later-on implement lookup tables and SSE as conditional macros.

    /// Clamps the given input value to the range [min, max].
    /** @see Clamp01(), Min(), Max(). */
    template<typename T>
    T Clamp(const T &val, const T &floor, const T &ceil)
    {
        assert(floor <= ceil);
        return val <= ceil ? (val >= floor ? val : floor) : ceil;
    }
    /// Clamps the given input value to the range [0, 1].
    /** @see Clamp(), Min(), Max(). */
    template<typename T>
    T Clamp01(const T &val) { return Clamp(val, T(0), T(1)); }

    /// Computes the smaller of the two values.
    /** @see Clamp(), Clamp01(), Max() */
    template <typename T>
    T Min(const T& a, const T& b) {
        return a <= b ? a : b;
    }

    /// Computes the larger of two values.
    /** @see Clamp(), Clamp01(), Max() */
    template <typename T>
    T Max(const T& a, const T& b) {
        return a >= b ? a : b;
    }

    /// Computes the smallest in an arbitrary list of values.
    /** @see Clamp(), Clamp01(), Max() */
    template <typename T>
    T Min(const std::initializer_list<T>& list) {
        T minimum = list[0];

        for (T entry : list) {
            if (entry <= minimum)
                minimum = entry;
        }
        return minimum;
    }

    /// Computes the largest in an arbitrary list of values.
    /** @see Clamp(), Clamp01(), Max() */
    template <typename T>
    T Max(const std::initializer_list<T>& list) {
        T maximum = list[0];
        for (T entry : list) {
            if (entry >= maximum)
                maximum = entry;
        }
        return maximum;
    }

    /** @return True if a > b. */
    template <typename T>
    bool GreaterThan(const T& a, const T& b) {
        return a > b;
    }

    /** @return True if a < b. */
    template <typename T>
    bool LessThan(const T& a, const T& b) {
        return a < b;
    }

    /** @return The absolute value of a. */
    template <typename T>
    T Abs(const T& a) {
        return a >= 0 ? a : -a;
    }

    template<> inline float Abs(const float& a) {
#ifdef USE_SSE
#else
        return a >= 0 ? a : -a;
#endif
    }

    /// @return True if a and b are equal, using operator ==().
    template <typename T>
    bool EqualExact(const T& a, const T& b) {
        return a == b;
    }

    /** Compares the two values for equality up to a small epsilon. */
    bool Equal(float a, float b, float epsilon = 1e-3f);

    /** Compares the two values for equality up to a small epsilon. */
    bool Equal(double a, double b, float epsilon = 1e-3f);

    /** Compares the two values for equality, allowing the given amount of absolute error. */
    bool EqualAbs(float a, float b, float epsilon = 1e-3f);

    /// Computes the relative error of the two variables.
    float RelativeError(float a, float b);

    template <typename T> bool IsFinite(T) { return true;}

    template<> inline bool IsFinite(float f) { return (ReinterpretAs<u32>(f) << 1) < 0xFF000000u; }
    template<> inline bool IsFinite(double d) { return (ReinterpretAs<u64>(d) << 1) < 0xFFE0000000000000ULL; }
    template<> inline bool IsFinite(u32 i) { return (i << 1) < 0xFF000000u; }
    template<> inline bool IsFinite(u64 i) { return (i << 1) < 0xFFE0000000000000ULL;}

    inline bool IsNotANumber(float f) { return (ReinterpretAs<u32>(f) << 1) > 0xFF000000u; }
    inline bool IsNotANumber(double d) { return (ReinterpretAs<u64>(d) << 1) > 0xFFE0000000000000ULL; }

    inline bool IsInfinite(float f) { return (ReinterpretAs<u32>(f) << 1) == 0xFF000000u; }
    inline bool IsInfinite(double d) { return (ReinterpretAs<u64>(d) << 1) == 0xFFE0000000000000ULL; }


    namespace Trigonometric {
        Sign SignOfSin(float radians);
        Sign SignOfCos(float radians);
        Sign SignOfTan(float radians);

        Quadrant QuadrantOf(float radians);


        float Radians(float deg); /// Converts the given amount of degrees into radians.
        float Degrees(float rad); /// Converts the given amount of radians into degrees.

        float Sin(float x); /// Computes the sine of x, in radians.
        float Cos(float x); /// Computes the cosine of x, in radians.
        float Tan(float x); /// Computes the tangent of x, in radians.

        /// Simultaneously computes both sine and cosine of x, in radians.
        /// This yields a small performance increase over computing them separately.
        /// @see Sin(), Cos().
        void SinCos(float x, float& outSin, float& outCos);

        float Asin(float x); /// Computes the inverse sine of x, in radians.
        float Acos(float x); /// Computes the inverse cosine of x, in radians.
        float Atan(float x); /// Computes the inverse tangent of x, in radians.
        float Atan2(float y, float x); /// Computes the signed (principal value) inverse tangent of y/x, in radians.
        float Sinh(float x); /// Computes the hyperbolic sine of x, in radians.
        float Cosh(float x); /// Computes the hyperbolic cosine of x, in radians.
        float Tanh(float x); /// Computes the hyperbolic tangent of x, in radians.
    }


    using namespace Trigonometric;

    bool IsPow2(u32 number); /// Returns true if the given number is a power of 2.
    bool IsPow2(u64 number); /// Returns true if the given number is a power of 2.

    float PowInt(float base, int exponent); /// Raises the given base to an integral exponent.
    float Pow(float base, float exponent); /// Raises the given base to a float exponent.
    float Exp(float exp); /// Returns e (the constant 2.71828...) raised to the given power.
    float Log(float base, float value); /// Computes a logarithm of the given value in the specified base.
    float Log2(float value); /// Computes a logarithm in base-2.
    float Ln(float value); /// Computes a logarithm in the natural base (using e as the base).
    float Log10(float value); /// Computes a logarithm in base-10;

    float Ceil(float f);    /// Returns f rounded up to the next integer, as float.
    int CeilInt(float f);   /// Returns f rounded up to the next integer, as integer.
    float Floor(float f);   /// Returns f rounded down to the previous integer, as float.
    int FloorInt(float f);  /// Returns f rounded down to the previous integer, as integer.
    float Round(float f);   /// Returns f rounded to the nearest integer, as float.
    int RoundInt(float f);  /// Returns f rounded to the nearest integer, as int.


    float Round(float f, float decimalPlaces); /// Returns f rounded to the given decimal places.


    float Sign(float f); /// Returns -1 or 1 depending on the sign of f.
    ///
    float SignOrZero(float f, float epsilon = 1e-8f);

    /// Formats larger numbers into shortened 'Truncated' string representations.
    /// 2241 -> 2.2k, 55421 -> 55.4k, 1000000 -> 1.0M
    std::string Truncate(float input);

    float RecipFast(float x);


    struct NumberRange
    {
        float LowerBound;
        float UpperBound;
    };


    float NormalizeToRange(float input, float fromLower, float fromUpper, float toLower, float toUpper);
    float NormalizeToRange(float input, const NumberRange& from, const NumberRange& to);
    // auto rotation_normalized = NormalizeToRange(inp, {0, 360}, {-1, 1});

    /// Linearly interpolates between a and b.
    /** @param t A value between [0,1].
        @param a The first endpoint to lerp between.
        @param b The second endpoint to lerp between.
        @return This function computes a + t*(b-a). That is, if t==0, this function returns a. If t==1, this function returns b.
            Otherwise, the returned value linearly moves from a to b as t ranges from 0 to 1.
        @see LerpMod(), InvLerp(), Step(), SmoothStep(), PingPongMod(), Mod(), ModPos(), Frac(). */
    float Lerp(float a, float b, float t);

    /// Linearly interpolates from a to b, under the modulus mod.
    /** This function takes evaluates a and b in the range [0, mod] and takes the shorter path to reach from a to b.
        @see Lerp(), InvLerp(), Step(), SmoothStep(), PingPongMod(), Mod(), ModPos(), Frac(). */
    float LerpMod(float a, float b, float mod, float t);

    /// Computes the lerp factor a and b have to be Lerp()ed to get x.
    /// /** @see Lerp(), LerpMod(), Step(), SmoothStep(), PingPongMod(), Mod(), ModPos(), Frac(). */
    float InverseLerp(float a, float b, float x);
    /// See http://msdn.microsoft.com/en-us/library/bb509665(v=VS.85).aspx
    float Step(float y, float x);
    /// See http://msdn.microsoft.com/en-us/library/bb509658(v=vs.85).aspx
    float Ramp(float min, float max, float x);
    /// Limits x to the range [0, mod], but instead of wrapping around from mod to 0, the result will move back
    /// from mod to 0 as x goes from mod to 2*mod.
    /** @see Lerp(), LerpMod(), InvLerp(), Step(), SmoothStep(), Mod(), ModPos(), Frac(). */
    float PingPongMod(float x, float mod);

    /// Computes a floating-point modulus.
    /// This function returns a value in the range ]-mod, mod[.
    /** @see Lerp(), LerpMod(), InvLerp(), Step(), SmoothStep(), PingPongMod(), ModPos(), Frac(). */
    float Mod(float x, float mod);
    /// Computes a floating-point modulus using an integer as the modulus.
    float Mod(float x, int mod);

    /// Computes a floating-point modulus, but restricts the output to the range [0, mod[.
    /** @see Lerp(), LerpMod(), InvLerp(), Step(), SmoothStep(), PingPongMod(), Mod(), Frac(). */
    float ModPos(float x, float mod);
    /// Computes a floating-point modulus, but restricts the output to the range [0, mod[.
    float ModPos(float x, int mod);
    /// Returns the fractional part of x.
    /** @see Lerp(), LerpMod(), InvLerp(), Step(), SmoothStep(), PingPongMod(), Mod(), ModPos(). */
    float Frac(float x);
    /// Returns the square root of x.
    float Sqrt(float x);
    /// Computes a fast approximation of the square root of x.
    float FastSqrt(float x);
    /// Returns 1/Sqrt(x). (The reciprocal of the square root of x)
    float RSqrt(float x);
    /// SSE implementation of reciprocal square root.
    float FastRSqrt(float x);
    /// Returns 1/x, the reciprocal of x.
    float Recip(float x);
    /// Returns 1/x, the reciprocal of x, using a fast approximation (SSE rcp instruction).
    float RecipFast(float x);
    /// Carmack (Quake) implementation of inverse (reciprocal) square root.
    /// Relies on funky type-casting hacks to avoid floating-point division and sqrt, which is very slow on legacy hardware.
    /// This technique is superseded by modern processors having built-in support, but is included for its historical significance.
    /// https://en.wikipedia.org/wiki/Fast_inverse_square_root
    float QRSqrt(float x);
}

namespace J3ML::Math::Functions::Interpolation
{
    inline float SmoothStart(float t);
}


namespace J3ML::Math {
    using namespace Functions;
}

namespace J3ML::Math::Types {


    struct Radians { // TODO: Fill in with relevant members.
        float value;
        float operator()() const { return value; }
    };

    struct Degrees { // TODO: Fill in with relevant members.
        float value;
        float operator()() const { return value; }
    };


}