#ifndef FIXED_POINT_HPP_
#define FIXED_POINT_HPP_

#include <Serialize/Serialize.hpp>
#include <Utilities/CircularInt.hpp>
#include <Utilities/Extern.hpp>
#include <Utilities/Utilities.hpp>
#include <cstddef>
#include <ios>
#include <sstream>

template<typename TInner, size_t NFractionBits>
class fixed
{
private:
    TInner repr;

public:
    static constexpr size_t fractionBits = NFractionBits;
    static constexpr size_t combinedBits = sizeof(TInner) * 8;
    static constexpr size_t integerBits = combinedBits - fractionBits;

    fixed() : repr(0) {}

    template<typename TSource>
    explicit constexpr fixed(TSource value)
    {
        if constexpr (std::is_floating_point_v<TSource>)
        {
            repr = static_cast<TInner>(value * (1ULL << fractionBits));
        }
        else if constexpr (std::is_integral_v<TSource>)
        {
            repr = static_cast<TInner>(value) << fractionBits;
        }
        else
        {
            static_assert(!std::is_same_v<TSource, TSource>, "Unsupported type for fixed constructor.");
        }
    }

    static constexpr fixed fromRaw(TInner rawValue)
    {
        fixed result;
        result.repr = rawValue;
        return result;
    }

    operator double() const
    {
        return static_cast<double>(repr) / (1ULL << fractionBits);
    }

    operator TInner() const
    {
        if (repr >= TInner(0))
        {
            return repr >> fractionBits;
        }
        else
        {
            return (repr + TInner((1LL << fractionBits) - 1)) >> fractionBits;
        }
    }

    fixed operator+(const fixed& other) const
    {
        return fixed::fromRaw(repr + other.repr);
    }

    fixed operator-(const fixed& other) const
    {
        return fixed::fromRaw(repr - other.repr);
    }

    fixed operator-() const
    {
        return fixed::fromRaw(-repr);
    }

    fixed operator*(const fixed& other) const
    {
        // Compute in parts such that we don't overflow in any individual multiplication, then sum up the parts.
        // For partD (the product of fractional parts), we round it towards the nearest representable value.
        // Furthermore, since partD is the product of the *fractional parts*, it's always a positive number between
        // 0 and 1, and the raw multiplication gives us only the fractional bits. Therefore there is no sign bit,
        // and it is implicitly positive, so we must convert to unsigned before shifting.
        auto partA = fromRaw(((repr >> fractionBits) * (other.repr >> fractionBits)) << fractionBits);
        auto partB = fromRaw((repr >> fractionBits) * other.fract().repr);
        auto partC = fromRaw((other.repr >> fractionBits) * fract().repr);
        auto partD = fromRaw(typename TInner::UnsignedT(other.fract().repr * fract().repr + (TInner(1) << (fractionBits - 1))) >> fractionBits);

        return partA + partB + partC + partD;
    }

    fixed operator/(const fixed& divisor) const
    {
        if (repr == TInner(0))
        {
            return fixed(0);
        }

        if (divisor.repr == TInner(0))
        {
            crash(formatString("Math error: division of $+ by zero", *this));
        }

        // Turn this division problem into a division of positive numbers.
        if (repr < TInner(0))
        {
            return -((-*this)/divisor);
        }

        if (divisor.repr < TInner(0))
        {
            return -((*this)/(-divisor));
        }

        // Now we have a positive division problem. We are trying to find x = a/b.
        // It must satisfy the condition that x * b == a.
        // 
        // But there's also a 'convergence test'. For any "y":
        //   y * b < a -> y < x
        //
        // So for any candidate, we can test whether it's smaller or larger than the
        // result. We compute the largest possible power-of-2, "k", which can be multiplied
        // by "b" without overflowing. If k < x, it means the result cannot be represented.
        //
        // Otherwise we let result = 0. We start with "k", and iterate over consecutive smaller
        // powers of 2, for each one doing a test:
        //
        // If (result + k) > x, do nothing; else, add "k" to "result".
        //
        // At the end, we will have the largest possible value that is smaller than or equal
        // to "x". If it is equal then we have an exact answer; otherwise:
        //   result < x < (result + step)
        //
        // So then we round to the nearest representable value (result or (result + step)).

        // Returns positive values if y is *too big*.
        auto getRelativeDistanceFromResult = [this, &divisor](fixed y)
        {
            return y * divisor - *this;
        };

        fixed step = fromRaw(TInner(1));

        fixed k(0);
        if (divisor <= fixed(1))
        {
            k = fromRaw(TInner(1) << (combinedBits - 2));
        }
        else
        {
            auto overflowBit = TInner(1) << (combinedBits - 1);
            auto overflowBoundary = TInner((uint64_t)overflowBit / (uint64_t)divisor.repr);

            for (size_t i = 0; i < integerBits; i++)
            {
                auto mask = TInner(1) << i;
                if (overflowBoundary & mask)
                {
                    k.repr = mask << fractionBits;
                }
            }

            if (getRelativeDistanceFromResult(fromRaw(overflowBoundary << fractionBits) - step) < fixed(0))
            {
                // Even the largest representable value without overflow is too small, so we cannot represent the result.
                crash(formatString("Math error: $+ / $+ is too big to represent", *this, divisor));
            }
        }

        fixed result(0);

        while (k.repr != TInner(0))
        {
            if ((result + k) * divisor >= fixed(0) && getRelativeDistanceFromResult(result + k) <= fixed(0))
            {
                result += k;
            }

            k.repr >>= 1;
        }

        if (getRelativeDistanceFromResult(result + step).abs() < getRelativeDistanceFromResult(result).abs())
        {
            result += step;
        }

        return result;
    }

    fixed& operator+=(const fixed& other)
    {
        repr += other.repr;
        return *this;
    }

    fixed& operator-=(const fixed& other)
    {
        repr -= other.repr;
        return *this;
    }

    fixed& operator*=(const fixed& other)
    {
        return (*this = *this * other);
    }

    fixed& operator/=(const fixed& other)
    {
        return (*this = *this / other);
    }

    bool operator==(const fixed& other) const
    {
        return repr == other.repr;
    }

    bool operator!=(const fixed& other) const
    {
        return repr != other.repr;
    }

    bool operator<(const fixed& other) const
    {
        return repr < other.repr;
    }

    bool operator<=(const fixed& other) const
    {
        return repr <= other.repr;
    }

    bool operator>(const fixed& other) const
    {
        return repr > other.repr;
    }

    bool operator>=(const fixed& other) const
    {
        return repr >= other.repr;
    }

    fixed operator<<(int bits) const
    {
        return fixed::fromRaw(repr << bits);
    }

    fixed operator>>(int bits) const
    {
        return fixed::fromRaw(repr >> bits);
    }

    fixed& operator<<=(int bits)
    {
        repr <<= bits;
        return *this;
    }

    fixed& operator>>=(int bits)
    {
        repr >>= bits;
        return *this;
    }

    TInner raw() const
    {
        return repr;
    }

    fixed fract() const
    {
        return fixed::fromRaw(repr & ((1ULL << fractionBits) - 1));
    }

    fixed floor() const
    {
        return fixed::fromRaw(repr & ~((1ULL << fractionBits) - 1));
    }

    fixed ceil() const
    {
        return fixed::fromRaw((repr + (TInner) ((1ULL << fractionBits) - 1)) & ~((1ULL << fractionBits) - 1));
    }

    fixed abs() const
    {
        if (repr < TInner(0))
        {
            return fixed::fromRaw(-repr);
        }

        return *this;
    }

    /**
     * @brief Returns the exact number represented by this `fixed` as a decimal string.
     */
    std::string exactDecimal() const
    {
        std::stringstream ss;
        
        if (*this < fixed(0))
        {
            ss << "-";
        }

        TInner intPart = abs();
        ss << intPart;

        auto rest = abs().fract();
        if (rest == fixed(0))
        {
            return ss.str();
        }

        ss << ".";

        while (rest != fixed(0))
        {
            rest *= fixed(10);
            ss << rest.operator TInner();
            rest = rest.fract();
        }

        return ss.str();
    }

    TInner asBits() const
    {
        return repr;
    }

    friend struct SerialCodec<fixed<TInner, NFractionBits>>;

    friend std::ostream& operator<<(std::ostream& os, const fixed& f)
    {
        std::streamsize precision = os.precision();
        auto base = os.flags() & std::ios_base::basefield;
        int radix = 10; // default

        switch (base)
        {
            case std::ios_base::hex: radix = 16; break;
            case std::ios_base::oct: radix = 8; break;
            default: radix = 10; break;
        }

        auto transformed = f.abs();
        if (f < fixed(0))
        {
            os << "-";
        }

        os << transformed.operator TInner() << ".";

        auto fraction = transformed.fract();
        for (size_t i = 0; i < precision; i++)
        {
            fraction *= fixed(radix);
            os << fraction.operator TInner();
            fraction = fraction.fract();
        }
        
        return os;
    }
};

template<typename TInner, size_t NFractionBits>
DERIVE_STRUCT_SERIALIZER_PARTIAL(SINGLE_ARG(fixed<TInner, NFractionBits>), self.repr);

using wfixed = fixed<cint64, 32>;
using wfvec3 = glm::tvec3<wfixed>;

template<typename TInner, size_t NFractionBits, int NComps>
glm::vec<NComps, fixed<TInner, NFractionBits>> operator<<(
    const glm::vec<NComps, fixed<TInner, NFractionBits>> &f,
    glm::vec<NComps, int> offsets
)
{
    glm::vec<NComps, fixed<TInner, NFractionBits>> result = f;
    for (int i = 0; i < NComps; i++)
    {
        result[i] <<= offsets[i];
    }

    return result;
}

// Overload some GLM functions.
namespace glm
{
    template <typename TInner, size_t NFractionBits, length_t L, qualifier Q>
    inline vec<L, fixed<TInner, NFractionBits>, Q> floor(const vec<L, fixed<TInner, NFractionBits>, Q>& v)
    {
        vec<L, fixed<TInner, NFractionBits>, Q> result(v);
        for (length_t i = 0; i < L; ++i)
        {
            result[i] = v[i].floor();
        }

        return result;
    }

    template <typename TInner, size_t NFractionBits, length_t L, qualifier Q>
    inline vec<L, fixed<TInner, NFractionBits>, Q> fract(const vec<L, fixed<TInner, NFractionBits>, Q>& v)
    {
        vec<L, fixed<TInner, NFractionBits>, Q> result(v);
        for (length_t i = 0; i < L; ++i)
        {
            result[i] = v[i].fract();
        }

        return result;
    }
}

#endif