/***
* ieeemisc.c - IEEE miscellaneous recommended functions
*
*       Copyright (c) 1992-2001, Microsoft Corporation. All rights reserved.
*
*Purpose:
*
*Revision History:
*        5-04-92  GDP   written
*        8-13-96  JWM   Order of tests in _fpclass() rearranged, since
*                       "if (x==0.0)" now uses FP hardware.
*       11-25-00  PML   IA64 _logb and _isnan supplied by libm .s code.
*
*******************************************************************************/

#include <trans.h>
#include <math.h>
#include <float.h>


/***
* _copysign - copy sign
*
*Purpose:
*   copysign(x,y) returns x with the sign of y. Hence, abs(x) := copysign
*   even if x is NaN [IEEE std 854-1987 Appendix]
*
*
*Entry:
*
*Exit:
*
*Exceptions:
*   No exceptions, even if one of the arguments is NaN.
*
*   (Currently the i386 compiler returns doubles on the fp stack
*   so the fld instruction at the end will cause an invalid operation
*   if x is NaN. However this compiler calling convention will change
*   soon)
*
*******************************************************************************/

double _copysign (double x, double y)
{
    double retval;
    *D_LO(retval) = *D_LO(x);
    *D_HI(retval) = *D_HI(x) & ~(1<<31) |
                    *D_HI(y) &  (1<<31) ;

    return retval;
}



/***
* _chgsign - change sign
*
*Purpose:
*  x is copied with its sign reversed, not 0-x; the distinction is germane
*  when x is +0, -0, or NaN
*
*Entry:
*
*Exit:
*
*Exceptions:
*   No exceptions, even if x is NaN.
*
*   (Currently the i386 compiler returns doubles on the fp stack
*   so the fld instruction at the end will cause an invalid operation
*   if x is NaN. However this compiler calling convention will change
*   soon)
*
*******************************************************************************/

double _chgsign (double x)
{
    double retval;

    *D_LO(retval) = *D_LO(x);
    *D_HI(retval) = *D_HI(x) & ~(1 << 31)  |
                    ~*D_HI(x) & (1<<31);

    return retval;
}


/***
* _scalb - scale by power of 2
*
*Purpose:
*   _scalb(x,n) returns x * 2^n for integral values of n without
*   computing 2^n
*   Special case:
*      If x is infinity or zero, _scaleb returns x
*
*
*Entry:
*   double x
*   int   n
*
*Exit:
*
*Exceptions:
*   Invalid operation, Overflow, Underflow
*
*******************************************************************************/

double _scalb(double x, long n)
{
    //
    // It turns out that our implementation of ldexp matces the IEEE
    // description of _scalb. The only problem with calling ldexp
    // is that if an exception occurs, the operation code reported
    // to the handler will be the one that corresponds to ldexp
    // (i.e., we do not define a new operation code for _scalb
    //

    return ldexp(x,n);
}


#if !defined(_M_IA64)

/***
* _logb - extract exponent
*
*Purpose:
*   _logb(x) returns the unbiased exponent of x, a signed integer in the
*   format of x, except that logb(NaN) is a NaN, logb(+INF) is +INF,and
*   logb(0) is is -INF and signals the division by zero exception.
*   For x positive and finite, 1<= abs(scalb(x, -logb(x))) < 2
*
*
*Entry:
*   double x
*   int   n
*
*Exit:
*
*Exceptions:
*   Invalid operation, Division by zero
*
*******************************************************************************/
double _logb(double x)
{
    uintptr_t savedcw;
    int exp;
    double retval;

    /* save user fp control word */
    savedcw = _maskfp();

    /* check for infinity or NAN */
    if (IS_D_SPECIAL(x)){
        switch (_sptype(x)) {
        case T_PINF:
        case T_NINF:
            RETURN(savedcw, x);
        case T_QNAN:
            return _handle_qnan1(OP_LOGB, x, savedcw);
        default: //T_SNAN
            return _except1(FP_I, OP_LOGB, x, _s2qnan(x), savedcw);
        }
    }

    if (x == 0) {
         return _except1(FP_Z, OP_LOGB, x, -D_INF, savedcw);
    }

    (void) _decomp(x, &exp);

    //
    // x == man * 2^exp, where .5 <= man < 1. According to the spec
    // of this function, we should compute the exponent so that
    // 1<=man<2, i.e., we should decrement the computed exp by one
    //

    retval = (double) (exp - 1);

    RETURN(savedcw, retval);

}

#endif  // !defined(_M_IA64)



/***
* _nextafter - next representable neighbor
*
*Purpose:
*  _nextafter(x,y) returns the next representable neighbor of x in
*  the direction toward y. The following special cases arise: if
*  x=y, then the result is x without any exception being signaled;
*  otherwise, if either x or y is a quiet NaN, then the result is
*  one or the other of the input NaNs. Overflow is sibnaled when x
*  is finite but nextafter(x,y) is infinite; underflow is signaled
*  when nextafter(x,y) lies strictly between -2^Emin, 2^Emin; in
*  both cases, inexact is signaled.
*
*
*Entry:
*
*Exit:
*
*Exceptions:
*   O, U, I, P
*
*******************************************************************************/

double _nextafter(double x, double y)
{
    uintptr_t savedcw;
    double result;

    /* save user fp control word */
    savedcw = _maskfp();

    /* check for infinity or NAN */
    if (IS_D_SPECIAL(x) || IS_D_SPECIAL(y)){
        if (IS_D_SNAN(x) || IS_D_SNAN(y)){
            return _except2(FP_I,OP_NEXTAFTER,x,y,_d_snan2(x,y),savedcw);
        }
        if (IS_D_QNAN(x) || IS_D_QNAN(y)){
            return _handle_qnan2(OP_NEXTAFTER,x,y,savedcw);
        }

        //
        // infinite arguments are not treated as special cases
        //
    }

    if (y == x) {

        //
        // no exceptions are raised in this case
        //

        RETURN(savedcw, x);
    }

    if (x == 0) {

        *D_LO(result) = 1;

        if (y > x) {
            *D_HI(result) = 0;
        }

        else {

            //
            // result should be negative
            //

            *D_HI(result) = (unsigned long)(1<<31);
        }

    }


    //
    // At this point x!=y, and x!=0. x can be treated as a 64bit
    // integer in sign/magnitude representation. To get the next
    // representable neighbor we add or subtract one from this
    // integer. (Note that for boundary cases like x==INF, need to
    // add one will never occur --this would mean that y should
    // be greater than INF, which is impossible)
    //

    if (x > 0 && y < x ||
        x < 0 && y > x) {

        //
        // decrease value by one
        //

        *D_LO(result) = *D_LO(x) - 1;
        *D_HI(result) = *D_HI(x);

        if (*D_LO(x) == 0) {

            //
            // a borrow should propagate to the high order dword
            //

            (*D_HI(result)) --;
        }
    }

    else if (x > 0 && y > x ||
             x < 0 && y < x) {

        //
        // increase value by one
        //

        *D_LO(result) = *D_LO(x) + 1;
        *D_HI(result) = *D_HI(x);

        if (*D_LO(result) == 0) {

            //
            // a carry should propagate to the high order dword
            //

            (*D_HI(result)) ++;
        }
    }


    //
    // check if an exception should be raised
    //


    if ( IS_D_DENORM(result) ) {

        //
        // should signal underflow and inexact
        // and provide a properly scaled value
        //

        double mant;
        int exp;

        mant = _decomp(result, &exp);
        result = _set_exp(mant, exp+IEEE_ADJUST);

        return _except2(FP_U|FP_P,OP_NEXTAFTER,x,y,result,savedcw);
    }



    if ( IS_D_INF(result) || IS_D_MINF(result) ) {

        //
        // should signal overflow and inexact
        // and provide a properly scaled value
        //

        double mant;
        int exp;

        mant = _decomp(result, &exp);
        result = _set_exp(mant, exp-IEEE_ADJUST);

        return _except2(FP_O|FP_P,OP_NEXTAFTER,x,y,result,savedcw);
    }


    RETURN(savedcw, result);
}




/***
* _finite -
*
*Purpose:
*   finite(x) returns the value TRUE if -INF < x < +INF and returns
*   false otherwise [IEEE std]
*
*Entry:
*
*Exit:
*
*Exceptions:
*
*   This routine is treated as a nonarithmetic operation, therefore
*   it does not signal any floating point exceptions
*
*******************************************************************************/

int _finite(double x)
{
    if (IS_D_SPECIAL(x)) {

        //
        // x is INF or NaN
        //

        return 0;
    }
    return 1;
}



#if !defined(_M_IA64)

/***
* _isnan -
*
*Purpose:
*   isnan(x) returns the value TRUE if x is a NaN, and returns FALSE
*   otherwise.
*
*
*Entry:
*
*Exit:
*
*Exceptions:
*
*   This routine is treated as a nonarithmetic operation, therefore
*   it does not signal any floating point exceptions
*
*******************************************************************************/

int _isnan(double x)
{
    if (IS_D_SNAN(x) || IS_D_QNAN(x)) {
        return 1;
    }
    return 0;
}

#endif  // !defined(_M_IA64)


/***
*double _fpclass(double x) - floating point class
*
*Purpose:
*   Compute the floating point class of a number, according
*   to the recommendations of the IEEE std. 754
*
*Entry:
*
*Exit:
*
*Exceptions:
*   This function is never exceptional, even when the argument is SNAN
*
*******************************************************************************/

int _fpclass(double x)
{
    int sign;

    if (IS_D_SPECIAL(x)){
        switch (_sptype(x)) {
        case T_PINF:
            return _FPCLASS_PINF;
        case T_NINF:
            return _FPCLASS_NINF;
        case T_QNAN:
            return _FPCLASS_QNAN;
        default: //T_SNAN
            return _FPCLASS_SNAN;
        }
    }
    sign = (*D_EXP(x)) & 0x8000;

    if (IS_D_DENORM(x))
        return sign? _FPCLASS_ND : _FPCLASS_PD;

    if (x == 0.0)
        return sign? _FPCLASS_NZ : _FPCLASS_PZ;

    return sign? _FPCLASS_NN : _FPCLASS_PN;
}