e_powf.s

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

data8 0x3fe4bea2a5bdbe84 // log(1/frcpa(1+233/256))=  +6.48271e-001
data8 0x3fe4cdf28f10ac44 // log(1/frcpa(1+234/256))=  +6.50140e-001
data8 0x3fe4dd49cf994058 // log(1/frcpa(1+235/256))=  +6.52013e-001
data8 0x3fe4eca86e64a680 // log(1/frcpa(1+236/256))=  +6.53889e-001
data8 0x3fe503c43cd8eb68 // log(1/frcpa(1+237/256))=  +6.56710e-001
data8 0x3fe513356667fc54 // log(1/frcpa(1+238/256))=  +6.58595e-001
data8 0x3fe522ae0738a3d4 // log(1/frcpa(1+239/256))=  +6.60483e-001
data8 0x3fe5322e26867854 // log(1/frcpa(1+240/256))=  +6.62376e-001
data8 0x3fe541b5cb979808 // log(1/frcpa(1+241/256))=  +6.64271e-001
data8 0x3fe55144fdbcbd60 // log(1/frcpa(1+242/256))=  +6.66171e-001
data8 0x3fe560dbc45153c4 // log(1/frcpa(1+243/256))=  +6.68074e-001
data8 0x3fe5707a26bb8c64 // log(1/frcpa(1+244/256))=  +6.69980e-001
data8 0x3fe587f60ed5b8fc // log(1/frcpa(1+245/256))=  +6.72847e-001
data8 0x3fe597a7977c8f30 // log(1/frcpa(1+246/256))=  +6.74763e-001
data8 0x3fe5a760d634bb88 // log(1/frcpa(1+247/256))=  +6.76682e-001
data8 0x3fe5b721d295f10c // log(1/frcpa(1+248/256))=  +6.78605e-001
data8 0x3fe5c6ea94431ef8 // log(1/frcpa(1+249/256))=  +6.80532e-001
data8 0x3fe5d6bb22ea86f4 // log(1/frcpa(1+250/256))=  +6.82462e-001
data8 0x3fe5e6938645d38c // log(1/frcpa(1+251/256))=  +6.84397e-001
data8 0x3fe5f673c61a2ed0 // log(1/frcpa(1+252/256))=  +6.86335e-001
data8 0x3fe6065bea385924 // log(1/frcpa(1+253/256))=  +6.88276e-001
data8 0x3fe6164bfa7cc068 // log(1/frcpa(1+254/256))=  +6.90222e-001
data8 0x3fe62643fecf9740 // log(1/frcpa(1+255/256))=  +6.92171e-001
LOCAL_OBJECT_END(pow_Tt)


// Table 1 is 2^(index_1/128) where
// index_1 goes from 0 to 15
LOCAL_OBJECT_START(pow_tbl1)
data8 0x8000000000000000 , 0x00003FFF
data8 0x80B1ED4FD999AB6C , 0x00003FFF
data8 0x8164D1F3BC030773 , 0x00003FFF
data8 0x8218AF4373FC25EC , 0x00003FFF
data8 0x82CD8698AC2BA1D7 , 0x00003FFF
data8 0x8383594EEFB6EE37 , 0x00003FFF
data8 0x843A28C3ACDE4046 , 0x00003FFF
data8 0x84F1F656379C1A29 , 0x00003FFF
data8 0x85AAC367CC487B15 , 0x00003FFF
data8 0x8664915B923FBA04 , 0x00003FFF
data8 0x871F61969E8D1010 , 0x00003FFF
data8 0x87DB357FF698D792 , 0x00003FFF
data8 0x88980E8092DA8527 , 0x00003FFF
data8 0x8955EE03618E5FDD , 0x00003FFF
data8 0x8A14D575496EFD9A , 0x00003FFF
data8 0x8AD4C6452C728924 , 0x00003FFF
LOCAL_OBJECT_END(pow_tbl1)


// Table 2 is 2^(index_1/8) where
// index_2 goes from 0 to 7
LOCAL_OBJECT_START(pow_tbl2)
data8 0x8000000000000000 , 0x00003FFF
data8 0x8B95C1E3EA8BD6E7 , 0x00003FFF
data8 0x9837F0518DB8A96F , 0x00003FFF
data8 0xA5FED6A9B15138EA , 0x00003FFF
data8 0xB504F333F9DE6484 , 0x00003FFF
data8 0xC5672A115506DADD , 0x00003FFF
data8 0xD744FCCAD69D6AF4 , 0x00003FFF
data8 0xEAC0C6E7DD24392F , 0x00003FFF
LOCAL_OBJECT_END(pow_tbl2)

.section .text
GLOBAL_LIBM_ENTRY(powf)

// Get exponent of x.  Will be used to calculate K.
{ .mfi
          getf.exp     pow_GR_signexp_X = f8
          fms.s1 POW_Xm1 = f8,f1,f1     // Will be used for r1 if x>0
          mov           pow_GR_17ones   = 0x1FFFF
}
{ .mfi
          addl          pow_AD_P        = @ltoff(pow_table_P), gp
          fma.s1 POW_Xp1 = f8,f1,f1     // Will be used for r1 if x<0
          nop.i 999
}
;;

// Get significand of x.  Will be used to get index to fetch T, Tt.
{ .mfi
          getf.sig      pow_GR_sig_X    = f8
          frcpa.s1      POW_B, p6       = f1,f8
          mov           pow_GR_exp_half = 0xFFFE   // Exponent for 0.5
}
{ .mfi
          ld8 pow_AD_P = [pow_AD_P]
          fma.s1        POW_NORM_X      = f8,f1,f0
          mov          pow_GR_exp_2tom8 = 0xFFF7
}
;;

// DOUBLE 0x10033  exponent limit at which y is an integer
{ .mfi
          nop.m 999
          fcmp.lt.s1 p8,p9 = f8, f0     // Test for x<0
          addl pow_GR_10033             = 0x10033, r0
}
{ .mfi
          mov           pow_GR_16ones   = 0xFFFF
          fma.s1        POW_NORM_Y      = f9,f1,f0
          nop.i 999
}
;;

// p13 = TRUE ==> X is unorm
{ .mfi
          setf.exp      POW_Q0_half     = pow_GR_exp_half  // Form 0.5
          fclass.m  p13,p0              = f8, 0x0b  // Test for x unorm
          adds          pow_AD_Tt       = pow_Tt - pow_table_P,  pow_AD_P
}
{ .mfi
          adds          pow_AD_Q        = pow_table_Q - pow_table_P,  pow_AD_P
          nop.f 999
          nop.i 999
}
;;

// p14 = TRUE ==> X is ZERO
{ .mfi
          ldfe          POW_P2          = [pow_AD_Q], 16
          fclass.m  p14,p0              = f8, 0x07
          nop.i 999
}
// Note POW_Xm1 and POW_r1 are used interchangably
{ .mfb
          nop.m 999
(p8)      fnma.s1        POW_Xm1        = POW_Xp1,f1,f0
(p13)     br.cond.spnt POW_X_DENORM
}
;;

// Continue normal and denormal paths here
POW_COMMON:
// p11 = TRUE ==> Y is a NAN
{ .mfi
          and           pow_GR_exp_X    = pow_GR_signexp_X, pow_GR_17ones
          fclass.m  p11,p0              = f9, 0xc3
          nop.i 999
}
{ .mfi
          nop.m 999
          fms.s1        POW_r           = POW_B, POW_NORM_X,f1
          mov pow_GR_y_zero = 0
}
;;

// Get exponent of |x|-1 to use in comparison to 2^-8
{ .mmi
          getf.exp  pow_GR_signexp_Xm1  = POW_Xm1
          sub       pow_GR_true_exp_X   = pow_GR_exp_X, pow_GR_16ones
          extr.u        pow_GR_offset   = pow_GR_sig_X, 55, 8
}
;;

{ .mfi
          alloc         r32=ar.pfs,2,19,4,0
          fcvt.fx.s1   POW_int_Y        = POW_NORM_Y
          shladd pow_AD_Tt = pow_GR_offset, 3, pow_AD_Tt
}
{ .mfi
          setf.sig POW_int_K            = pow_GR_true_exp_X
          nop.f 999
          nop.i 999
}
;;

// p12 = TRUE if Y is ZERO
// Compute xsq to decide later if |x|=1
{ .mfi
          ldfe          POW_P1          = [pow_AD_P], 16
          fclass.m      p12,p0          = f9, 0x07
          nop.i 999
}
{ .mfb
          ldfe          POW_P0          = [pow_AD_Q], 16
          fma.s1        POW_xsq = POW_NORM_X, POW_NORM_X, f0
(p11)     br.cond.spnt  POW_Y_NAN       // Branch if y=nan
}
;;

{ .mmf
          getf.exp  pow_GR_signexp_Y    = POW_NORM_Y
          ldfd  POW_T                   = [pow_AD_Tt]
          fma.s1        POW_rsq         = POW_r, POW_r,f0
}
;;

// p11 = TRUE ==> X is a NAN
{ .mfi
          ldfpd         POW_log2_hi, POW_log2_lo  = [pow_AD_Q], 16
          fclass.m      p11,p0          = POW_NORM_X, 0xc3
          nop.i 999
}
{ .mfi
          ldfe          POW_inv_log2_by_128 = [pow_AD_P], 16
          fma.s1 POW_delta              = f0,f0,f0 // delta=0 in case |x| near 1
(p12)     mov pow_GR_y_zero = 1
}
;;

{ .mfi
          ldfd   POW_Q2                 = [pow_AD_P], 16
          fnma.s1 POW_twoV              = POW_r, POW_Q0_half,f1
          and       pow_GR_exp_Xm1      = pow_GR_signexp_Xm1, pow_GR_17ones
}
{ .mfi
          nop.m 999
          fma.s1 POW_U                  = POW_NORM_Y,POW_r,f0
          nop.i 999
}
;;

// Determine if we will use the |x| near 1 path (p6) or normal path (p7)
{ .mfi
          nop.m 999
          fcvt.xf POW_K                 = POW_int_K
          cmp.lt p6,p7                  = pow_GR_exp_Xm1, pow_GR_exp_2tom8
}
{ .mfb
          nop.m 999
          fma.s1 POW_G                  = f0,f0,f0  // G=0 in case |x| near 1
(p11)     br.cond.spnt  POW_X_NAN       // Branch if x=nan and y not nan
}
;;

// If on the x near 1 path, assign r1 to r
{ .mfi
          ldfpd  POW_Q1, POW_RSHF       = [pow_AD_P], 16
(p6)      fma.s1    POW_r               = POW_r1, f1, f0
          nop.i 999
}
{ .mfb
          nop.m 999
(p6)      fma.s1    POW_rsq             = POW_r1, POW_r1, f0
(p14)     br.cond.spnt POW_X_0          // Branch if x zero and y not nan
}
;;

{ .mfi
          getf.sig pow_GR_sig_int_Y     = POW_int_Y
(p6)      fnma.s1 POW_twoV              = POW_r1, POW_Q0_half,f1
          and pow_GR_exp_Y              = pow_GR_signexp_Y, pow_GR_17ones
}
{ .mfb
          andcm pow_GR_sign_Y           = pow_GR_signexp_Y, pow_GR_17ones
(p6)      fma.s1 POW_U                  = POW_NORM_Y,POW_r1,f0
(p12)     br.cond.spnt POW_Y_0   // Branch if y=zero, x not zero or nan
}
;;

{ .mfi
          ldfe      POW_log2_by_128_lo  = [pow_AD_P], 16
(p7)      fma.s1 POW_Z2                 = POW_twoV, POW_U, f0
          nop.i 999
}
{ .mfi
          ldfe      POW_log2_by_128_hi  = [pow_AD_Q], 16
          nop.f 999
          nop.i 999
}
;;

{ .mfi
          nop.m 999
          fcvt.xf   POW_float_int_Y     = POW_int_Y
          nop.i 999
}
{ .mfi
          nop.m 999
(p7)      fma.s1 POW_G                  = POW_K, POW_log2_hi, POW_T
          adds          pow_AD_tbl1     = pow_tbl1 - pow_Tt,  pow_AD_Q
}
;;

// p11 = TRUE ==> X is NEGATIVE but not inf
{ .mfi
          nop.m 999
          fclass.m  p11,p0              = POW_NORM_X, 0x1a
          nop.i 999
}
{ .mfi
          nop.m 999
(p7)      fma.s1 POW_delta              = POW_K, POW_log2_lo, f0
          adds pow_AD_tbl2              = pow_tbl2 - pow_tbl1,  pow_AD_tbl1
}
;;

{ .mfi
          nop.m 999
(p6)      fma.s1 POW_Z                  = POW_twoV, POW_U, f0
          nop.i 999
}
{ .mfi
          nop.m 999
          fma.s1 POW_v2                 = POW_P1, POW_r,  POW_P0
          nop.i 999
}
;;

// p11 = TRUE ==> X is NEGATIVE but not inf
//    p12 = TRUE ==> X is NEGATIVE  AND  Y  already even int
//    p13 = TRUE ==> X is NEGATIVE  AND  Y possible int
{ .mfi
          nop.m 999
(p7)      fma.s1 POW_Z                  = POW_NORM_Y, POW_G, POW_Z2
(p11)     cmp.gt.unc  p12,p13           = pow_GR_exp_Y, pow_GR_10033
}
{ .mfi
          nop.m 999
          fma.s1 POW_Gpr                = POW_G, f1, POW_r
          nop.i 999
}
;;

{ .mfi
          nop.m 999
          fma.s1 POW_Yrcub              = POW_rsq, POW_U, f0
          nop.i 999
}
{ .mfi
          nop.m 999
          fma.s1 POW_p                  = POW_rsq, POW_P2, POW_v2
          nop.i 999
}
;;

// Test if x inf
{ .mfi
          nop.m 999
          fclass.m p15,p0 = POW_NORM_X,  0x23
          nop.i 999
}
// By adding RSHF (1.1000...*2^63) we put integer part in rightmost significand
{ .mfi
          nop.m 999
          fma.s1 POW_W1  = POW_Z, POW_inv_log2_by_128, POW_RSHF
          nop.i 999
}
;;

// p13 = TRUE ==> X is NEGATIVE  AND  Y possible int
//     p10 = TRUE ==> X is NEG and Y is an int
//     p12 = TRUE ==> X is NEG and Y is not an int
{ .mfi
          nop.m 999
(p13)     fcmp.eq.unc.s1 p10,p12        = POW_float_int_Y,  POW_NORM_Y
          mov pow_GR_xneg_yodd = 0
}
{ .mfi
          nop.m 999
          fma.s1 POW_Y_Gpr              = POW_NORM_Y, POW_Gpr, f0
          nop.i 999
}
;;

// p11 = TRUE ==> X is +1.0
{ .mfi
          nop.m 999
          fcmp.eq.s1 p11,p0 = POW_NORM_X, f1
          nop.i 999
}
;;

// Extract rounded integer from rightmost significand of POW_W1
// By subtracting RSHF we get rounded integer POW_Nfloat
{ .mfi
          getf.sig pow_GR_int_N        = POW_W1
          fms.s1 POW_Nfloat  = POW_W1, f1, POW_RSHF
          nop.i 999
}
{ .mfb
          nop.m 999
          fma.s1 POW_Z3                 = POW_p, POW_Yrcub, f0
(p12)     br.cond.spnt POW_X_NEG_Y_NONINT  // Branch if x neg, y not integer
}
;;

// p7  = TRUE ==> Y is +1.0
// p12 = TRUE ==> X is NEGATIVE  AND Y is an odd integer
{ .mfi
          getf.exp pow_GR_signexp_Y_Gpr = POW_Y_Gpr
          fcmp.eq.s1 p7,p0 = POW_NORM_Y, f1  // Test for y=1.0
(p10)     tbit.nz.unc  p12,p0           = pow_GR_sig_int_Y,0
}
{ .mfb
          nop.m 999
(p11)     fma.s.s0 f8 = f1,f1,f0    // If x=1, result is +1
(p15)     br.cond.spnt POW_X_INF
}
;;

// Test x and y and flag denormal
{ .mfi
          nop.m 999
          fcmp.eq.s0 p15,p0 = f8,f9
          nop.i 999
}
{ .mfb
          nop.m 999
          fma.s1 POW_e3                 = POW_NORM_Y, POW_delta, f0
(p11)     br.ret.spnt b0            // Early exit if x=1.0, result is +1
}
;;

{ .mfi