libm_support.h

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

        (r)->ex   =   (y)->ex;                              \
        (r)->ldhi = -((y)->ldhi);                           \
        (r)->ldlo = -((y)->ldlo);                           \
    } else {                                                \
        /* |y|<<|x| */                                      \
        *(r) = *(x);                                        \
    }
#elif defined(SIZE_INT_64)
#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)->significand = 0x8000000000000000;     \
            (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| */                                      \
        (r)->ex   =   (y)->ex;                              \
        (r)->ldhi = -((y)->ldhi);                           \
        (r)->ldlo = -((y)->ldlo);                           \
    } else {                                                \
        /* |y|<<|x| */                                      \
        *(r) = *(x);                                        \
    }
#endif

/* 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 for any finite x and y                           */
#define __LIBM_SUBL_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_SUBL_K80(r,x,y, t1)

/* Multiplication: x*y                                      */
/* The result is sum rhi+rlo                                */
/* Here t32 is the constant 2^32+1                          */
/* Temporary variables: t1,t2,t3,t4,t5,t6                   */
/* All variables are in long double precision               */
/* Correct if no over/underflow (algorithm by T.J.Dekker)   */
#define __LIBM_MULL1_K80(rhi,rlo,x,y,                       \
                                     t32,t1,t2,t3,t4,t5,t6) \
    t1=(x)*(t32); t3=x-t1; t3=t3+t1; t4=x-t3;               \
    t1=(y)*(t32); t5=y-t1; t5=t5+t1; t6=y-t5;               \
    t1=(t3)*(t5);                                           \
    t2=(t3)*(t6)+(t4)*(t5);                                 \
    rhi=t1+t2;                                              \
    rlo=t1-rhi; rlo=rlo+t2; rlo=rlo+(t4*t6);

/* Multiplication: (xhi+xlo)*(yhi+ylo)                      */
/* The result is sum rhi+rlo                                */
/* Here t32 is the constant 2^32+1                          */
/* Temporary variables: t1,t2,t3,t4,t5,t6,t7,t8             */
/* All variables are in long double precision               */
/* Correct if no over/underflow (algorithm by T.J.Dekker)   */
#define __LIBM_MULL2_K80(rhi,rlo,xhi,xlo,yhi,ylo,           \
                               t32,t1,t2,t3,t4,t5,t6,t7,t8) \
    __LIBM_MULL1_K80(t7,t8,xhi,yhi, t32,t1,t2,t3,t4,t5,t6)  \
    t1=(xhi)*(ylo)+(xlo)*(yhi); t1=t1+t8;                   \
    rhi=t7+t1;                                              \
    rlo=t7-rhi; rlo=rlo+t1;

/* Multiplication: r=x*y                                    */
/* Variables r,x,y are pointers to struct ker80,            */
/* all other variables are in long double precision         */
/* Here t32 is the constant 2^32+1                          */
/* Temporary variables: t1,t2,t3,t4,t5,t6,t7,t8             */
/* Correct if x and y belong to interval [2^-8000;2^8000]   */
#define __LIBM_MULL_K80(r,x,y, t32,t1,t2,t3,t4,t5,t6,t7,t8) \
    (r)->ex = (x)->ex + (y)->ex;                            \
    __LIBM_MULL2_K80((r)->ldhi,(r)->ldlo,                   \
        (x)->ldhi,(x)->ldlo,(y)->ldhi,(y)->ldlo,            \
        t32,t1,t2,t3,t4,t5,t6,t7,t8)

/* Multiplication: r=x*y                                    */
/* Variables r,x,y are pointers to struct ker80,            */
/* all other variables are in long double precision         */
/* Here t32 is the constant 2^32+1                          */
/* Temporary variables: t1,t2,t3,t4,t5,t6,t7,t8             */
/* Correct for any finite x and y                           */
#define __LIBM_MULL_NORM_K80(r,x,y,                         \
                               t32,t1,t2,t3,t4,t5,t6,t7,t8) \
    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_MULL_K80(r,x,y, t32,t1,t2,t3,t4,t5,t6,t7,t8)

/* Division: (xhi+xlo)/(yhi+ylo)                            */
/* The result is sum rhi+rlo                                */
/* Here t32 is the constant 2^32+1                          */
/* Temporary variables: t1,t2,t3,t4,t5,t6,t7,t8,t9          */
/* All variables are in long double precision               */
/* Correct if no over/underflow (algorithm by T.J.Dekker)   */
#define __LIBM_DIVL2_K80(rhi,rlo,xhi,xlo,yhi,ylo,           \
                            t32,t1,t2,t3,t4,t5,t6,t7,t8,t9) \
    t7=(xhi)/(yhi);                                         \
    __LIBM_MULL1_K80(t8,t9,t7,yhi, t32,t1,t2,t3,t4,t5,t6)   \
    t1=xhi-t8; t1=t1-t9; t1=t1+xlo; t1=t1-(t7)*(ylo);       \
    t1=(t1)/(yhi);                                          \
    rhi=t7+t1;                                              \
    rlo=t7-rhi; rlo=rlo+t1;

/* Division: r=x/y                                          */
/* Variables r,x,y are pointers to struct ker80,            */
/* all other variables are in long double precision         */
/* Here t32 is the constant 2^32+1                          */
/* Temporary variables: t1,t2,t3,t4,t5,t6,t7,t8,t9          */
/* Correct if x and y belong to interval [2^-8000;2^8000]   */
#define __LIBM_DIVL_K80(r,x,y,                              \
                            t32,t1,t2,t3,t4,t5,t6,t7,t8,t9) \
    (r)->ex = (x)->ex - (y)->ex;                            \
    __LIBM_DIVL2_K80( (r)->ldhi,(r)->ldlo,                  \
        (x)->ldhi,(x)->ldlo,(y)->ldhi,(y)->ldlo,            \
        t32,t1,t2,t3,t4,t5,t6,t7,t8,t9)

/* Division: r=x/y                                          */
/* Variables r,x,y are pointers to struct ker80,            */
/* all other variables are in long double precision         */
/* Here t32 is the constant 2^32+1                          */
/* Temporary variables: t1,t2,t3,t4,t5,t6,t7,t8             */
/* Correct for any finite x and y                           */
#define __LIBM_DIVL_NORM_K80(r,x,y,                         \
                            t32,t1,t2,t3,t4,t5,t6,t7,t8,t9) \
    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_DIVL_K80(r,x,y, t32,t1,t2,t3,t4,t5,t6,t7,t8,t9)

/* Square root: sqrt(xhi+xlo)                               */
/* The result is sum rhi+rlo                                */
/* Here t32 is the constant 2^32+1                          */
/*      half is the constant 0.5                            */
/* Temporary variables: t1,t2,t3,t4,t5,t6,t7,t8,t9          */
/* All variables are in long double precision               */
/* Correct for positive xhi+xlo (algorithm by T.J.Dekker)   */
#define __LIBM_SQRTL2_NORM_K80(rhi,rlo,xhi,xlo,             \
                       t32,half,t1,t2,t3,t4,t5,t6,t7,t8,t9) \
    t7=sqrtl(xhi);                                          \
    __LIBM_MULL1_K80(t8,t9,t7,t7, t32,t1,t2,t3,t4,t5,t6)    \
    t1=xhi-t8; t1=t1-t9; t1=t1+xlo; t1=(t1)*(half);         \
    t1=(t1)/(t7);                                           \
    rhi=t7+t1;                                              \
    rlo=t7-rhi; rlo=rlo+t1;

/* Square root: r=sqrt(x)                                   */
/* Variables r,x,y are pointers to struct ker80,            */
/* all other variables are in long double precision         */
/* Here t32 is the constant 2^32+1                          */
/*      half is the constant 0.5                            */
/* Temporary variables: t1,t2,t3,t4,t5,t6,t7,t8,t9          */
/* Correct if x belongs to interval [2^-16000;2^16000]      */
#define __LIBM_SQRTL_K80(r,x,                               \
                       t32,half,t1,t2,t3,t4,t5,t6,t7,t8,t9) \
    if ( ((x)->ex & 1) == 1 ) {                             \
        (x)->ex    = (x)->ex + 1;                           \
        (x)->ldhi *= half;                                  \
        (x)->ldlo *= half;                                  \
    }                                                       \
    (r)->ex = (x)->ex >> 1;                                 \
    __LIBM_SQRTL2_NORM_K80( (r)->ldhi,(r)->ldlo,            \
        (x)->ldhi,(x)->ldlo,                                \
        t32,half,t1,t2,t3,t4,t5,t6,t7,t8,t9)

/* Square root: r=sqrt(x)                                   */
/* Variables r,x,y are pointers to struct ker80,            */
/* all other variables are in long double precision         */
/* Here t32 is the constant 2^32+1                          */
/*      half is the constant 0.5                            */
/* Temporary variables: t1,t2,t3,t4,t5,t6,t7,t8,t9          */
/* Correct for any positive x                               */
#define __LIBM_SQRTL_NORM_K80(r,x,                          \
                       t32,half,t1,t2,t3,t4,t5,t6,t7,t8,t9) \
    if ( ((x)->fphi.exponent-BIAS_80<-16000) ||             \
         ((x)->fphi.exponent-BIAS_80>+16000) )              \
    {                                                       \
        __libm_normalizel_k80(x);                           \
    }                                                       \
    __LIBM_SQRTL_K80(r,x, t32,half,t1,t2,t3,t4,t5,t6,t7,t8,t9)


#ifdef __INTEL_COMPILER
#define ALIGN(n) __declspec(align(n))
#else /* __INTEL_COMPILER */
#define ALIGN(n)
#endif /* __INTEL_COMPILER */

/* macros to form a long double value in hex representation (unsigned short type) */

#if (defined(__unix__) && defined(__i386__))
# define LDOUBLE_ALIGN 12	/* IA32 Linux: 12-byte alignment */
#else	/*__linux__ & IA32*/
# define LDOUBLE_ALIGN 16	/* EFI2/IA32 Win or IPF Win/Linux: 16-byte alignment */
#endif	/*__linux__ & IA32*/

#if (LDOUBLE_ALIGN == 16)
#define _XPD_ ,0x0000,0x0000,0x0000
#else /*12*/
#define _XPD_ ,0x0000
#endif

#define LDOUBLE_HEX(w4,w3,w2,w1,w0) 0x##w0,0x##w1,0x##w2,0x##w3,0x##w4 _XPD_ /*LITTLE_ENDIAN*/

/* macros to sign-expand low 'num' bits of 'val' to native integer */

#if defined(SIZE_INT_32)
# define SIGN_EXPAND(val,num)  ((int)(val) << (32-(num))) >> (32-(num)) /* sign expand of 'num' LSBs */
#elif defined(SIZE_INT_64)
# define SIGN_EXPAND(val,num)  ((int)(val) << (64-(num))) >> (64-(num)) /* sign expand of 'num' LSBs */
#endif

/* macros to form pointers to FP number on-the-fly */

#define FP32(f)  ((struct fp32 *)&f)
#define FP64(d)  ((struct fp64 *)&d)