// Compiler implementation of the D programming language
// Copyright (c) 1999-2011 by Digital Mars
// All Rights Reserved
// written by Rainer Schuetze
// http://www.digitalmars.com
// License for redistribution is by either the Artistic License
// in artistic.txt, or the GNU General Public License in gnu.txt.
// See the included readme.txt for details.

// 80 bit floating point value implementation for Microsoft compiler

#if _MSC_VER
#include "longdouble.h"

#include "assert.h"

#include <float.h>
#include <stdio.h>
#include <string.h>

extern "C"
{
    // implemented in ldfpu.asm for _WIN64
    int ld_initfpu(int bits, int mask);
    void ld_expl(longdouble* ld, int exp);
    longdouble ld_add(longdouble ld1, longdouble ld2);
    longdouble ld_sub(longdouble ld1, longdouble ld2);
    longdouble ld_mul(longdouble ld1, longdouble ld2);
    longdouble ld_div(longdouble ld1, longdouble ld2);
    longdouble ld_mod(longdouble ld1, longdouble ld2);
    bool ld_cmpb(longdouble ld1, longdouble ld2);
    bool ld_cmpbe(longdouble ld1, longdouble ld2);
    bool ld_cmpa(longdouble ld1, longdouble ld2);
    bool ld_cmpae(longdouble ld1, longdouble ld2);
    bool ld_cmpe(longdouble ld1, longdouble ld2);
    bool ld_cmpne(longdouble ld1, longdouble ld2);
    longdouble ld_sqrt(longdouble ld1);
    longdouble ld_sin(longdouble ld1);
    longdouble ld_cos(longdouble ld1);
    longdouble ld_tan(longdouble ld1);
}

bool initFPU()
{
#ifdef _WIN64
//    int old_cw = ld_initfpu(_RC_NEAR);
    int old_cw = ld_initfpu(0x300 /*_PC_64  | _RC_NEAR*/, // #defines NOT identical to CPU FPU control word!
                            0xF00 /*_MCW_PC | _MCW_RC*/);
#else
    int old_cw = _control87(_MCW_EM | _PC_64  | _RC_NEAR,
                            _MCW_EM | _MCW_PC | _MCW_RC);
#endif
    return true;
}
static bool doInitFPU = initFPU();

#ifndef _WIN64
extern "C"
{

double ld_read(const longdouble* pthis)
{
    double res;
    __asm
    {
        mov eax, pthis
        fld tbyte ptr [eax]
        fstp res
    }
    return res;
}
long long ld_readll(const longdouble* pthis)
{
#if 1
    return ld_readull(pthis);
#elif defined _WIN64
    return ld_readll(this);
#else
    longdouble* pthis = this;
    long long res;
    __asm
    {
        mov eax, pthis
        fld tbyte ptr [eax]
        fistp qword ptr res
    }
    return res;
#endif
}

unsigned long long ld_readull(const longdouble* pthis)
{
#if 1
    // somehow the FPU does not respect the CHOP mode of the rounding control
    // in 64-bit mode
    // so we roll our own conversion (it also allows the usual C wrap-around
    // instead of the "invalid value" created by the FPU)
    int expo = pthis->exponent - 0x3fff;
    unsigned long long u;
    if(expo < 0 || expo > 127)
        return 0;
    if(expo < 64)
        u = pthis->mantissa >> (63 - expo);
    else
        u = pthis->mantissa << (expo - 63);
    if(pthis->sign)
        u = ~u + 1;
    return u;
#else
    longdouble* pthis = this;
    long long res; // cannot use unsigned, VC will not generate "fistp qword"
    longdouble twoPow63 = { 1ULL << 63, 0x3fff + 63, 0 };
    __asm
    {
        mov eax, pthis
        fld tbyte ptr [eax]
        fld tbyte ptr twoPow63
        fsubp ST(1),ST(0)  // move it into signed range

        lea eax, res
        fistp qword ptr [eax]
    }
    res ^= (1LL << 63);
    return res;
#endif
}

void ld_set(longdouble* pthis, double d)
{
    __asm
    {
        mov eax, pthis
        fld d
        fstp tbyte ptr [eax]
    }
}
void ld_setll(longdouble* pthis, long long d)
{
    __asm
    {
        fild qword ptr d
        mov eax, pthis
        fstp tbyte ptr [eax]
    }
}
void ld_setull(longdouble* pthis, unsigned long long d)
{
    d ^= (1LL << 63);
    longdouble twoPow63 = { 1ULL << 63, 0x3fff + 63, 0 };
    __asm
    {
        fild qword ptr d
        fld tbyte ptr twoPow63
        faddp ST(1),ST(0)
        mov eax, pthis
        fstp tbyte ptr [eax]
    }
}

} // extern "C"
#endif // !_WIN64

longdouble ldexpl(longdouble ld, int exp)
{
#ifdef _WIN64
    ld_expl(&ld, exp);
#else
    __asm
    {
        fild    dword ptr exp
        fld     tbyte ptr ld
        fscale                  // ST(0) = ST(0) * (2**ST(1))
        fstp    ST(1)
        fstp    tbyte ptr ld
    }
#endif
    return ld;
}

///////////////////////////////////////////////////////////////////////
longdouble operator+(longdouble ld1, longdouble ld2)
{
#ifdef _WIN64
    return ld_add(ld1, ld2);
#else
    longdouble res;
    __asm
    {
        fld tbyte ptr ld1
        fld tbyte ptr ld2
        fadd
        fstp tbyte ptr res;
    }
    return res;
#endif
}

longdouble operator-(longdouble ld1, longdouble ld2)
{
#ifdef _WIN64
    return ld_sub(ld1, ld2);
#else
    longdouble res;
    __asm
    {
        fld tbyte ptr ld1
        fld tbyte ptr ld2
        fsub
        fstp tbyte ptr res;
    }
    return res;
#endif
}

longdouble operator*(longdouble ld1, longdouble ld2)
{
#ifdef _WIN64
    return ld_mul(ld1, ld2);
#else
    longdouble res;
    __asm
    {
        fld tbyte ptr ld1
        fld tbyte ptr ld2
        fmul
        fstp tbyte ptr res;
    }
    return res;
#endif
}

longdouble operator/(longdouble ld1, longdouble ld2)
{
#ifdef _WIN64
    return ld_div(ld1, ld2);
#else
    longdouble res;
    __asm
    {
        fld tbyte ptr ld1
        fld tbyte ptr ld2
        fdiv
        fstp tbyte ptr res;
    }
    return res;
#endif
}

bool operator< (longdouble x, longdouble y)
{
#ifdef _WIN64
    return ld_cmpb(x, y);
#else
    short sw;
    bool res;
    __asm
    {
        fld     tbyte ptr y
        fld     tbyte ptr x             // ST = x, ST1 = y
        fucomip ST(0),ST(1)
        setb    AL
        setnp   AH
        and     AL,AH
        mov     res,AL
        fstp    ST(0)
    }
    return res;
#endif
}
bool operator<=(longdouble x, longdouble y)
{
#ifdef _WIN64
    return ld_cmpbe(x, y);
#else
    short sw;
    bool res;
    __asm
    {
        fld     tbyte ptr y
        fld     tbyte ptr x             // ST = x, ST1 = y
        fucomip ST(0),ST(1)
        setbe   AL
        setnp   AH
        and     AL,AH
        mov     res,AL
        fstp    ST(0)
    }
    return res;
#endif
}
bool operator> (longdouble x, longdouble y)
{
#ifdef _WIN64
    return ld_cmpa(x, y);
#else
    short sw;
    bool res;
    __asm
    {
        fld     tbyte ptr y
        fld     tbyte ptr x             // ST = x, ST1 = y
        fucomip ST(0),ST(1)
        seta    AL
        setnp   AH
        and     AL,AH
        mov     res,AL
        fstp    ST(0)
    }
    return res;
#endif
}
bool operator>=(longdouble x, longdouble y)
{
#ifdef _WIN64
    return ld_cmpae(x, y);
#else
    short sw;
    bool res;
    __asm
    {
        fld     tbyte ptr y
        fld     tbyte ptr x             // ST = x, ST1 = y
        fucomip ST(0),ST(1)
        setae   AL
        setnp   AH
        and     AL,AH
        mov     res,AL
        fstp    ST(0)
    }
    return res;
#endif
}
bool operator==(longdouble x, longdouble y)
{
#ifdef _WIN64
    return ld_cmpe(x, y);
#else
    short sw;
    bool res;
    __asm
    {
        fld     tbyte ptr y
        fld     tbyte ptr x             // ST = x, ST1 = y
        fucomip ST(0),ST(1)
        sete    AL
        setnp   AH
        and     AL,AH
        mov     res,AL
        fstp    ST(0)
    }
    return res;
#endif
}
bool operator!=(longdouble x, longdouble y)
{
#ifdef _WIN64
    return ld_cmpne(x, y);
#else
    short sw;
    bool res;
    __asm
    {
        fld     tbyte ptr y
        fld     tbyte ptr x             // ST = x, ST1 = y
        fucomip ST(0),ST(1)
        setne   AL
        setp    AH
        or      AL,AH
        mov     res,AL
        fstp    ST(0)
    }
    return res;
#endif
}


int _isnan(longdouble ld)
{
    return (ld.exponent == 0x7fff && ld.mantissa != 0);
}

longdouble fabsl(longdouble ld)
{
    ld.sign = 0;
    return ld;
}

longdouble sqrtl(longdouble ld)
{
#ifdef _WIN64
    return ld_sqrt(ld);
#else
    longdouble res;
    __asm
    {
        fld tbyte ptr ld;
        fsqrt;
        fstp tbyte ptr res;
    }
    return res;
#endif
}

longdouble sinl (longdouble ld)
{
#ifdef _WIN64
    return ld_sin(ld);
#else
    longdouble res;
    __asm
    {
        fld tbyte ptr ld;
        fsin; // exact for |x|<=PI/4
        fstp tbyte ptr res
    }
    return res;
#endif
}
longdouble cosl (longdouble ld)
{
#ifdef _WIN64
    return ld_cos(ld);
#else
    longdouble res;
    __asm
    {
        fld tbyte ptr ld;
        fcos; // exact for |x|<=PI/4
        fstp tbyte ptr res;
    }
    return res;
#endif
}
longdouble tanl (longdouble ld)
{
#ifdef _WIN64
    return ld_tan(ld);
#else
    longdouble res;
    __asm
    {
        fld tbyte ptr ld;
        fptan;
        fstp ST(0); // always 1
        fstp tbyte ptr res;
    }
    return res;
#endif
}

longdouble fmodl(longdouble x, longdouble y)
{
#ifdef _WIN64
    return ld_mod(x, y);
#else
    short sw;
    longdouble res;
    __asm
    {
        fld     tbyte ptr y
        fld     tbyte ptr x             // ST = x, ST1 = y
FM1:    // We don't use fprem1 because for some inexplicable
        // reason we get -5 when we do _modulo(15, 10)
        fprem                           // ST = ST % ST1
        fstsw   word ptr sw
        fwait
        mov     AH,byte ptr sw+1        // get msb of status word in AH
        sahf                            // transfer to flags
        jp      FM1                     // continue till ST < ST1
        fstp    ST(1)                   // leave remainder on stack
        fstp    tbyte ptr res;
    }
    return res;
#endif
}

//////////////////////////////////////////////////////////////

longdouble ld_qnan = { 0x8000000000000000ULL, 0x7fff, 0 };
longdouble ld_snan = { 0x0000000000000001ULL, 0x7fff, 0 };
longdouble ld_inf  = { 0x0000000000000000ULL, 0x7fff, 0 };

longdouble ld_zero  = { 0, 0, 0 };
longdouble ld_one   = { 0x8000000000000000ULL, 0x3fff, 0 };
longdouble ld_pi    = { 0xc90fdaa22168c235ULL, 0x4000, 0 };
longdouble ld_log2t = { 0xd49a784bcd1b8afeULL, 0x4000, 0 };
longdouble ld_log2e = { 0xb8aa3b295c17f0bcULL, 0x3fff, 0 };
longdouble ld_log2  = { 0x9a209a84fbcff799ULL, 0x3ffd, 0 };
longdouble ld_ln2   = { 0xb17217f7d1cf79acULL, 0x3ffe, 0 };

longdouble ld_pi2     = ld_pi*2;
longdouble ld_piOver2 = ld_pi*0.5;
longdouble ld_piOver4 = ld_pi*0.25;

//////////////////////////////////////////////////////////////

#define LD_TYPE_OTHER    0
#define LD_TYPE_ZERO     1
#define LD_TYPE_INFINITE 2
#define LD_TYPE_SNAN     3
#define LD_TYPE_QNAN     4

int ld_type(longdouble x)
{
    if(x.exponent == 0)
        return x.mantissa == 0 ? LD_TYPE_ZERO : LD_TYPE_OTHER; // dnormal if not zero
    if(x.exponent != 0x7fff)
        return LD_TYPE_OTHER;
    if(x.mantissa == 0)
        return LD_TYPE_INFINITE;
    if(x.mantissa & (1LL << 63))
        return LD_TYPE_QNAN;
    return LD_TYPE_SNAN;
}

int ld_sprint(char* str, int fmt, longdouble x)
{
    // fmt is 'a','A','f' or 'g'
    if(fmt != 'a' && fmt != 'A')
    {
        char format[] = { '%', fmt, 0 };
        return sprintf(str, format, ld_read(&x));
    }

    unsigned short exp = x.exponent;
    unsigned long long mantissa = x.mantissa;

    switch(ld_type(x))
    {
    case LD_TYPE_ZERO:
        return sprintf(str, "0x0.0L");
    case LD_TYPE_QNAN:
    case LD_TYPE_SNAN:
        return sprintf(str, "NAN");
    case LD_TYPE_INFINITE:
        return sprintf(str, x.sign ? "-INF" : "INF");
    }

    int len = 0;
    if(x.sign)
        str[len++] = '-';
    len += sprintf(str + len, mantissa & (1LL << 63) ? "0x1." : "0x0.");
    mantissa = mantissa << 1;
    while(mantissa)
    {
        int dig = (mantissa >> 60) & 0xf;
        dig += dig < 10 ? '0' : fmt - 10;
        str[len++] = dig;
        mantissa = mantissa << 4;
    }
    str[len++] = 'p';
    if(exp < 0x3fff)
    {
        str[len++] = '-';
        exp = 0x3fff - exp;
    }
    else
    {
        str[len++] = '+';
        exp = exp - 0x3fff;
    }
    int exppos = len;
    for(int i = 12; i >= 0; i -= 4)
    {
        int dig = (exp >> i) & 0xf;
        if(dig != 0 || len > exppos || i == 0)
            str[len++] = dig + (dig < 10 ? '0' : fmt - 10);
    }
    str[len] = 0;
    return len;
}

//////////////////////////////////////////////////////////////

#if UNITTEST
static bool unittest()
{
    char buffer[32];
    ld_sprint(buffer, 'a', ld_pi);
    assert(strcmp(buffer, "0x1.921fb54442d1846ap+1") == 0);

    longdouble ldb = ldouble(0.4);
    long long b = ldb;
    assert(b == 0);

    b = ldouble(0.9);
    assert(b == 0);

    long long x = 0x12345678abcdef78LL;
    longdouble ldx = ldouble(x);
    assert(ldx > 0);
    long long y = ldx;
    assert(x == y);

    x = -0x12345678abcdef78LL;
    ldx = ldouble(x);
    assert(ldx < 0);
    y = ldx;
    assert(x == y);

    unsigned long long u = 0x12345678abcdef78LL;
    longdouble ldu = ldouble(u);
    assert(ldu > 0);
    unsigned long long v = ldu;
    assert(u == v);

    u = 0xf234567812345678ULL;
    ldu = ldouble(u);
    assert(ldu > 0);
    v = ldu;
    assert(u == v);

    u = 0xf2345678;
    ldu = ldouble(u);
    ldu = ldu * ldu;
    ldu = sqrt(ldu);
    v = ldu;
    assert(u == v);

    u = 0x123456789A;
    ldu = ldouble(u);
    ldu = ldu * (1LL << 23);
    v = ldu;
    u = u * (1LL << 23);
    assert(u == v);

    return true;
}

static bool runUnittest = unittest();

#endif // UNITTEST

#endif // _MSC_VER