libm_support.h

glibc 2.9,最新版的C语言库函数

#if !(defined(opensource))
typedef          __int32  INT32;
typedef   signed __int32 SINT32;
typedef unsigned __int32 UINT32;

typedef          __int64  INT64;
typedef   signed __int64 SINT64;
typedef unsigned __int64 UINT64;
#else
typedef          int  INT32;
typedef   signed int SINT32;
typedef unsigned int UINT32;

typedef          long long  INT64;
typedef   signed long long SINT64;
typedef unsigned long long UINT64;
#endif

#if (defined(_WIN32) || defined(_WIN64))        /* Windows */
# define I64CONST(bits) 0x##bits##i64
# define U64CONST(bits) 0x##bits##ui64
#elif (defined(__linux__) && defined(_M_IA64))  /* Linux,64 */
# define I64CONST(bits) 0x##bits##L
# define U64CONST(bits) 0x##bits##uL
#else                                           /* Linux,32 */
# define I64CONST(bits) 0x##bits##LL
# define U64CONST(bits) 0x##bits##uLL
#endif

struct ker80 {
    union {
        long double ldhi;
        struct fp80 fphi;
    };
    union {
        long double ldlo;
        struct fp80 fplo;
    };
    int ex;
};

/* Addition: x+y                                            */
/* The result is sum rhi+rlo                                */
/* Temporary variables: t1                                  */
/* All variables are in long double precision               */
/* Correct if no overflow (algorithm by D.Knuth)           */
#define __LIBM_ADDL1_K80( rhi,rlo,x,y, t1 )                 \
    rhi = x   + y;                                          \
    rlo = rhi - x;                                          \
    t1  = rhi - rlo;                                        \
    rlo = y   - rlo;                                        \
    t1  = x   - t1;                                         \
    rlo = rlo + t1;

/* Addition: (xhi+xlo) + (yhi+ylo)                          */
/* The result is sum rhi+rlo                                */
/* Temporary variables: t1                                  */
/* All variables are in long double precision               */
/* Correct if no overflow (algorithm by T.J.Dekker)         */
#define __LIBM_ADDL2_K80( rhi,rlo,xhi,xlo,yhi,ylo, t1 )     \
    rlo = xhi+yhi;                                          \
    if ( VALUE_GT_80(FP80(xhi),FP80(yhi)) ) {               \
        t1=xhi-rlo;t1=t1+yhi;t1=t1+ylo;t1=t1+xlo;           \
    } else {                                                \
        t1=yhi-rlo;t1=t1+xhi;t1=t1+xlo;t1=t1+ylo;           \
    }                                                       \
    rhi=rlo+t1;                                             \
    rlo=rlo-rhi;rlo=rlo+t1;

/* Addition: r=x+y                                          */
/* Variables r,x,y are pointers to struct ker80,            */
/* all other variables are in long double precision         */
/* Temporary variables: t1                                  */
/* Correct if x and y belong to interval [2^-8000;2^8000],  */
/* or when one or both of them are zero                     */
#if   defined(SIZE_INT_32)
#define __LIBM_ADDL_K80(r,x,y, t1)                          \
    if ( ((y)->ex+(y)->fphi.exponent-134 <                  \
          (x)->ex+(x)->fphi.exponent)       &&              \
         ((x)->ex+(x)->fphi.exponent <                      \
          (y)->ex+(y)->fphi.exponent+134)   &&              \
         !SIGNIFICAND_ZERO_80(&((x)->fphi)) &&              \
         !SIGNIFICAND_ZERO_80(&((y)->fphi)) )               \
    {                                                       \
        /* y/2^134 < x < y*2^134,               */          \
        /* and x,y are nonzero finite numbers   */          \
        if ( (x)->ex != (y)->ex ) {                         \
            /* adjust x->ex to y->ex */                     \
            /* t1 = 2^(x->ex - y->ex) */                    \
            FP80(t1)->sign = 0;                             \
            FP80(t1)->exponent = BIAS_80 + (x)->ex-(y)->ex; \
            /*  exponent is correct because             */  \
            /*  |x->ex - y->ex| =                       */  \
            /*  = |  (x->ex + x->fphi.exponent) -       */  \
            /*      -(y->ex + y->fphi.exponent) +       */  \
            /*              + y->fphi.exponent  -       */  \
            /*              - x->fphi.exponent     | <  */  \
            /*  < |  (x->ex+x->fphi.exponent) -         */  \
            /*      -(y->ex+y->fphi.exponent)      | +  */  \
            /*   +|  y->fphi.exponent -                 */  \
            /*      -x->fphi.exponent              | <  */  \
            /*  < 134 + 16000                           */  \
            FP80(t1)->hi_significand = 0x80000000;          \
            FP80(t1)->lo_significand = 0x00000000;          \
            (x)->ex = (y)->ex;                              \
            (x)->ldhi *= t1;                                \
            (x)->ldlo *= t1;                                \
        }                                                   \
        /* r==x+y */                                        \
        (r)->ex = (y)->ex;                                  \
        __LIBM_ADDL2_K80( (r)->ldhi,(r)->ldlo,              \
            (x)->ldhi,(x)->ldlo, (y)->ldhi,(y)->ldlo, t1 ); \
    } else if ( SIGNIFICAND_ZERO_80(&((x)->fphi)) ||        \
             ((y)->ex+(y)->fphi.exponent-BIAS_80 - 134 >=   \
              (x)->ex+(x)->fphi.exponent-BIAS_80) )         \
    {                                                       \
        /* |x|<<|y| */                                      \
        *(r) = *(y);                                        \
    } else {                                                \
        /* |y|<<|x| */                                      \
        *(r) = *(x);                                        \
    }
#elif defined(SIZE_INT_64)
#define __LIBM_ADDL_K80(r,x,y, t1)                          \
    if ( ((y)->ex+(y)->fphi.exponent-134 <                  \
          (x)->ex+(x)->fphi.exponent)       &&              \
         ((x)->ex+(x)->fphi.exponent <                      \
          (y)->ex+(y)->fphi.exponent+134)   &&              \
         !SIGNIFICAND_ZERO_80(&((x)->fphi)) &&              \
         !SIGNIFICAND_ZERO_80(&((y)->fphi)) )               \
    {                                                       \
        /* y/2^134 < x < y*2^134,               */          \
        /* and x,y are nonzero finite numbers   */          \
        if ( (x)->ex != (y)->ex ) {                         \
            /* adjust x->ex to y->ex */                     \
            /* t1 = 2^(x->ex - y->ex) */                    \
            FP80(t1)->sign = 0;                             \
            FP80(t1)->exponent = BIAS_80 + (x)->ex-(y)->ex; \
            /*  exponent is correct because             */  \
            /*  |x->ex - y->ex| =                       */  \
            /*  = |  (x->ex + x->fphi.exponent) -       */  \
            /*      -(y->ex + y->fphi.exponent) +       */  \
            /*              + y->fphi.exponent  -       */  \
            /*              - x->fphi.exponent     | <  */  \
            /*  < |  (x->ex+x->fphi.exponent) -         */  \
            /*      -(y->ex+y->fphi.exponent)      | +  */  \
            /*   +|  y->fphi.exponent -                 */  \
            /*      -x->fphi.exponent              | <  */  \
            /*  < 134 + 16000                           */  \
            FP80(t1)->significand = 0x8000000000000000;     \
            (x)->ex = (y)->ex;                              \
            (x)->ldhi *= t1;                                \
            (x)->ldlo *= t1;                                \
        }                                                   \
        /* r==x+y */                                        \
        (r)->ex = (y)->ex;                                  \
        __LIBM_ADDL2_K80( (r)->ldhi,(r)->ldlo,              \
            (x)->ldhi,(x)->ldlo, (y)->ldhi,(y)->ldlo, t1 ); \
    } else if ( SIGNIFICAND_ZERO_80(&((x)->fphi)) ||        \
             ((y)->ex+(y)->fphi.exponent-BIAS_80 - 134 >=   \
              (x)->ex+(x)->fphi.exponent-BIAS_80) )         \
    {                                                       \
        /* |x|<<|y| */                                      \
        *(r) = *(y);                                        \
    } else {                                                \
        /* |y|<<|x| */                                      \
        *(r) = *(x);                                        \
    }
#endif

/* Addition: r=x+y                                          */
/* Variables r,x,y are pointers to struct ker80,            */
/* all other variables are in long double precision         */
/* Temporary variables: t1                                  */
/* Correct for any finite x and y                           */
#define __LIBM_ADDL_NORM_K80(r,x,y, t1)                     \
    if ( ((x)->fphi.exponent-BIAS_80<-8000) ||              \
         ((x)->fphi.exponent-BIAS_80>+8000) ||              \
         ((y)->fphi.exponent-BIAS_80<-8000) ||              \
         ((y)->fphi.exponent-BIAS_80>+8000) )               \
    {                                                       \
        __libm_normalizel_k80(x);                           \
        __libm_normalizel_k80(y);                           \
    }                                                       \
    __LIBM_ADDL_K80(r,x,y, t1)

/* Subtraction: x-y                                         */
/* The result is sum rhi+rlo                                */
/* Temporary variables: t1                                  */
/* All variables are in long double precision               */
/* Correct if no overflow (algorithm by D.Knuth)           */
#define __LIBM_SUBL1_K80( rhi, rlo, x, y, t1 )              \
    rhi = x   - y;                                          \
    rlo = rhi - x;                                          \
    t1  = rhi - rlo;                                        \
    rlo = y   + rlo;                                        \
    t1  = x   - t1;                                         \
    rlo = t1  - rlo;

/* Subtraction: (xhi+xlo) - (yhi+ylo)                       */
/* The result is sum rhi+rlo                                */
/* Temporary variables: t1                                  */
/* All variables are in long double precision               */
/* Correct if no overflow (algorithm by T.J.Dekker)         */
#define __LIBM_SUBL2_K80( rhi,rlo,xhi,xlo,yhi,ylo, t1 )     \
    rlo = xhi-yhi;                                          \
    if ( VALUE_GT_80(FP80(xhi),FP80(yhi)) ) {               \
        t1=xhi-rlo;t1=t1-yhi;t1=t1-ylo;t1=t1+xlo;           \
    } else {                                                \
        t1=yhi+rlo;t1=xhi-t1;t1=t1+xlo;t1=t1-ylo;           \
    }                                                       \
    rhi=rlo+t1;                                             \
    rlo=rlo-rhi;rlo=rlo+t1;

/* Subtraction: r=x-y                                       */
/* Variables r,x,y are pointers to struct ker80,            */
/* all other variables are in long double precision         */
/* Temporary variables: t1                                  */
/* Correct if x and y belong to interval [2^-8000;2^8000],  */
/* or when one or both of them are zero                     */
#if   defined(SIZE_INT_32)
#define __LIBM_SUBL_K80(r,x,y, t1)                          \
    if ( ((y)->ex+(y)->fphi.exponent-134 <                  \
          (x)->ex+(x)->fphi.exponent)       &&              \
         ((x)->ex+(x)->fphi.exponent <                      \
          (y)->ex+(y)->fphi.exponent+134)   &&              \
         !SIGNIFICAND_ZERO_80(&((x)->fphi)) &&              \
         !SIGNIFICAND_ZERO_80(&((y)->fphi)) )               \
    {                                                       \
        /* y/2^134 < x < y*2^134,               */          \
        /* and x,y are nonzero finite numbers   */          \
        if ( (x)->ex != (y)->ex ) {                         \
            /* adjust x->ex to y->ex */                     \
            /* t1 = 2^(x->ex - y->ex) */                    \
            FP80(t1)->sign = 0;                             \
            FP80(t1)->exponent = BIAS_80 + (x)->ex-(y)->ex; \
            /*  exponent is correct because             */  \
            /*  |x->ex - y->ex| =                       */  \
            /*  = |  (x->ex + x->fphi.exponent) -       */  \
            /*      -(y->ex + y->fphi.exponent) +       */  \
            /*              + y->fphi.exponent  -       */  \
            /*              - x->fphi.exponent     | <  */  \
            /*  < |  (x->ex+x->fphi.exponent) -         */  \
            /*      -(y->ex+y->fphi.exponent)      | +  */  \
            /*   +|  y->fphi.exponent -                 */  \
            /*      -x->fphi.exponent              | <  */  \
            /*  < 134 + 16000                           */  \
            FP80(t1)->hi_significand = 0x80000000;          \
            FP80(t1)->lo_significand = 0x00000000;          \
            (x)->ex = (y)->ex;                              \
            (x)->ldhi *= t1;                                \
            (x)->ldlo *= t1;                                \
        }                                                   \
        /* r==x+y */                                        \
        (r)->ex = (y)->ex;                                  \
        __LIBM_SUBL2_K80( (r)->ldhi,(r)->ldlo,              \
            (x)->ldhi,(x)->ldlo, (y)->ldhi,(y)->ldlo, t1 ); \
    } else if ( SIGNIFICAND_ZERO_80(&((x)->fphi)) ||        \
             ((y)->ex+(y)->fphi.exponent-BIAS_80 - 134 >=   \
              (x)->ex+(x)->fphi.exponent-BIAS_80) )         \
    {                                                       \
        /* |x|<<|y| */                                      \