module std.math.traits;

import std.traits;

enum RealFormat
{
	ieeeHalf,
	ieeeSingle,
	ieeeDouble,
	ieeeExtended,   // x87 80-bit real
	ieeeExtended53, // x87 real rounded to precision of double.
	ibmExtended,    // IBM 128-bit extended
	ieeeQuadruple,
}

template floatTraits(T)
{
	import std.traits : Unqual;

	// EXPMASK is a ushort mask to select the exponent portion (without sign)
	// EXPSHIFT is the number of bits the exponent is left-shifted by in its ushort
	// EXPBIAS is the exponent bias - 1 (exp == EXPBIAS yields ×2^-1).
	// EXPPOS_SHORT is the index of the exponent when represented as a ushort array.
	// SIGNPOS_BYTE is the index of the sign when represented as a ubyte array.
	// RECIP_EPSILON is the value such that (smallest_subnormal) * RECIP_EPSILON == T.min_normal
	enum Unqual!T RECIP_EPSILON = (1/T.epsilon);
	static if (T.mant_dig == 24)
	{
		// Single precision float
		enum ushort EXPMASK = 0x7F80;
		enum ushort EXPSHIFT = 7;
		enum ushort EXPBIAS = 0x3F00;
		enum uint EXPMASK_INT = 0x7F80_0000;
		enum uint MANTISSAMASK_INT = 0x007F_FFFF;
		enum realFormat = RealFormat.ieeeSingle;
		version (LittleEndian)
		{
			enum EXPPOS_SHORT = 1;
			enum SIGNPOS_BYTE = 3;
		}
		else
		{
			enum EXPPOS_SHORT = 0;
			enum SIGNPOS_BYTE = 0;
		}
	}
	else static if (T.mant_dig == 53)
	{
		static if (T.sizeof == 8)
		{
			// Double precision float, or real == double
			enum ushort EXPMASK = 0x7FF0;
			enum ushort EXPSHIFT = 4;
			enum ushort EXPBIAS = 0x3FE0;
			enum uint EXPMASK_INT = 0x7FF0_0000;
			enum uint MANTISSAMASK_INT = 0x000F_FFFF; // for the MSB only
			enum ulong MANTISSAMASK_LONG = 0x000F_FFFF_FFFF_FFFF;
			enum realFormat = RealFormat.ieeeDouble;
			version (LittleEndian)
			{
				enum EXPPOS_SHORT = 3;
				enum SIGNPOS_BYTE = 7;
			}
			else
			{
				enum EXPPOS_SHORT = 0;
				enum SIGNPOS_BYTE = 0;
			}
		}
		else static if (T.sizeof == 12)
		{
			// Intel extended real80 rounded to double
			enum ushort EXPMASK = 0x7FFF;
			enum ushort EXPSHIFT = 0;
			enum ushort EXPBIAS = 0x3FFE;
			enum realFormat = RealFormat.ieeeExtended53;
			version (LittleEndian)
			{
				enum EXPPOS_SHORT = 4;
				enum SIGNPOS_BYTE = 9;
			}
			else
			{
				enum EXPPOS_SHORT = 0;
				enum SIGNPOS_BYTE = 0;
			}
		}
		else
			static assert(false, "No traits support for " ~ T.stringof);
	}
	else static if (T.mant_dig == 64)
	{
		// Intel extended real80
		enum ushort EXPMASK = 0x7FFF;
		enum ushort EXPSHIFT = 0;
		enum ushort EXPBIAS = 0x3FFE;
		enum realFormat = RealFormat.ieeeExtended;
		version (LittleEndian)
		{
			enum EXPPOS_SHORT = 4;
			enum SIGNPOS_BYTE = 9;
		}
		else
		{
			enum EXPPOS_SHORT = 0;
			enum SIGNPOS_BYTE = 0;
		}
	}
	else static if (T.mant_dig == 113)
	{
		// Quadruple precision float
		enum ushort EXPMASK = 0x7FFF;
		enum ushort EXPSHIFT = 0;
		enum ushort EXPBIAS = 0x3FFE;
		enum realFormat = RealFormat.ieeeQuadruple;
		version (LittleEndian)
		{
			enum EXPPOS_SHORT = 7;
			enum SIGNPOS_BYTE = 15;
		}
		else
		{
			enum EXPPOS_SHORT = 0;
			enum SIGNPOS_BYTE = 0;
		}
	}
	else static if (T.mant_dig == 106)
	{
		// IBM Extended doubledouble
		enum ushort EXPMASK = 0x7FF0;
		enum ushort EXPSHIFT = 4;
		enum realFormat = RealFormat.ibmExtended;

		// For IBM doubledouble the larger magnitude double comes first.
		// It's really a double[2] and arrays don't index differently
		// between little and big-endian targets.
		enum DOUBLEPAIR_MSB = 0;
		enum DOUBLEPAIR_LSB = 1;

		// The exponent/sign byte is for most significant part.
		version (LittleEndian)
		{
			enum EXPPOS_SHORT = 3;
			enum SIGNPOS_BYTE = 7;
		}
		else
		{
			enum EXPPOS_SHORT = 0;
			enum SIGNPOS_BYTE = 0;
		}
	}
	else
		static assert(false, "No traits support for " ~ T.stringof);
}

bool isNaN(X)(X x) @nogc @trusted pure nothrow
if (isFloatingPoint!(X))
{
    version (all)
    {
        return x != x;
    }
    else
    {
        /*
        Code kept for historical context. At least on Intel, the simple test
        x != x uses one dedicated instruction (ucomiss/ucomisd) that runs in one
        cycle. Code for 80- and 128-bits is larger but still smaller than the
        integrals-based solutions below. Future revisions may enable the code
        below conditionally depending on hardware.
        */
        alias F = floatTraits!(X);
        static if (F.realFormat == RealFormat.ieeeSingle)
        {
            const uint p = *cast(uint *)&x;
            // Sign bit (MSB) is irrelevant so mask it out.
            // Next 8 bits should be all set.
            // At least one bit among the least significant 23 bits should be set.
            return (p & 0x7FFF_FFFF) > 0x7F80_0000;
        }
        else static if (F.realFormat == RealFormat.ieeeDouble)
        {
            const ulong  p = *cast(ulong *)&x;
            // Sign bit (MSB) is irrelevant so mask it out.
            // Next 11 bits should be all set.
            // At least one bit among the least significant 52 bits should be set.
            return (p & 0x7FFF_FFFF_FFFF_FFFF) > 0x7FF0_0000_0000_0000;
        }
        else static if (F.realFormat == RealFormat.ieeeExtended ||
                        F.realFormat == RealFormat.ieeeExtended53)
        {
            const ushort e = F.EXPMASK & (cast(ushort *)&x)[F.EXPPOS_SHORT];
            const ulong ps = *cast(ulong *)&x;
            return e == F.EXPMASK &&
                ps & 0x7FFF_FFFF_FFFF_FFFF; // not infinity
        }
        else static if (F.realFormat == RealFormat.ieeeQuadruple)
        {
            const ushort e = F.EXPMASK & (cast(ushort *)&x)[F.EXPPOS_SHORT];
            const ulong psLsb = (cast(ulong *)&x)[MANTISSA_LSB];
            const ulong psMsb = (cast(ulong *)&x)[MANTISSA_MSB];
            return e == F.EXPMASK &&
                (psLsb | (psMsb& 0x0000_FFFF_FFFF_FFFF)) != 0;
        }
        else
        {
            return x != x;
        }
    }
}

bool isInfinity(X)(X x) @nogc @trusted pure nothrow
if (isFloatingPoint!(X))
{
    import std.math.traits : floatTraits, RealFormat, MANTISSA_MSB, MANTISSA_LSB;

    alias F = floatTraits!(X);
    static if (F.realFormat == RealFormat.ieeeSingle)
    {
        return ((*cast(uint *)&x) & 0x7FFF_FFFF) == 0x7F80_0000;
    }
    else static if (F.realFormat == RealFormat.ieeeDouble)
    {
        return ((*cast(ulong *)&x) & 0x7FFF_FFFF_FFFF_FFFF)
            == 0x7FF0_0000_0000_0000;
    }
    else static if (F.realFormat == RealFormat.ieeeExtended ||
                    F.realFormat == RealFormat.ieeeExtended53)
    {
        const ushort e = cast(ushort)(F.EXPMASK & (cast(ushort *)&x)[F.EXPPOS_SHORT]);
        const ulong ps = *cast(ulong *)&x;

        // On Motorola 68K, infinity can have hidden bit = 1 or 0. On x86, it is always 1.
        return e == F.EXPMASK && (ps & 0x7FFF_FFFF_FFFF_FFFF) == 0;
    }
    else static if (F.realFormat == RealFormat.ieeeQuadruple)
    {
        const long psLsb = (cast(long *)&x)[MANTISSA_LSB];
        const long psMsb = (cast(long *)&x)[MANTISSA_MSB];
        return (psLsb == 0)
            && (psMsb & 0x7FFF_FFFF_FFFF_FFFF) == 0x7FFF_0000_0000_0000;
    }
    else
    {
        return (x < -X.max) || (X.max < x);
    }
}

version (LittleEndian)
{
    enum MANTISSA_LSB = 0;
    enum MANTISSA_MSB = 1;
}
else
{
    enum MANTISSA_LSB = 1;
    enum MANTISSA_MSB = 0;
}