s_cosl.s

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

data8 0xAAAA000000000000, 0x0000BFFC // PP_1_hi
data8 0xB8EF1D2ABAF69EEA, 0x00003FEC // PP_4
data8 0xD00D00D00D03BB69, 0x0000BFF2 // PP_3
data8 0x8888888888888962, 0x00003FF8 // PP_2
data8 0xAAAAAAAAAAAB0000, 0x0000BFEC // PP_1_lo
LOCAL_OBJECT_END(sincosl_table_pp)

LOCAL_OBJECT_START(sincosl_table_qq)
data8 0xD56232EFC2B0FE52, 0x00003FD2 // QQ_8
data8 0xC9C99ABA2B48DCA6, 0x0000BFDA // QQ_7
data8 0x8F76C6509C716658, 0x00003FE2 // QQ_6
data8 0x93F27DBAFDA8D0FC, 0x0000BFE9 // QQ_5
data8 0xAAAAAAAAAAAAAAAA, 0x0000BFFC // S_1
data8 0x8000000000000000, 0x0000BFFE // QQ_1
data8 0xD00D00D00C6E5041, 0x00003FEF // QQ_4
data8 0xB60B60B60B607F60, 0x0000BFF5 // QQ_3
data8 0xAAAAAAAAAAAAAA9B, 0x00003FFA // QQ_2
LOCAL_OBJECT_END(sincosl_table_qq)

LOCAL_OBJECT_START(sincosl_table_c)
data8 0xFFFFFFFFFFFFFFFE, 0x0000BFFD // C_1
data8 0xAAAAAAAAAAAA719F, 0x00003FFA // C_2
data8 0xB60B60B60356F994, 0x0000BFF5 // C_3
data8 0xD00CFFD5B2385EA9, 0x00003FEF // C_4
data8 0x93E4BD18292A14CD, 0x0000BFE9 // C_5
LOCAL_OBJECT_END(sincosl_table_c)

LOCAL_OBJECT_START(sincosl_table_s)
data8 0xAAAAAAAAAAAAAAAA, 0x0000BFFC // S_1
data8 0x88888888888868DB, 0x00003FF8 // S_2
data8 0xD00D00D0055EFD4B, 0x0000BFF2 // S_3
data8 0xB8EF1C5D839730B9, 0x00003FEC // S_4
data8 0xD71EA3A4E5B3F492, 0x0000BFE5 // S_5
data4 0x38800000, 0xB8800000                        // two**-14 and -two**-14
LOCAL_OBJECT_END(sincosl_table_s)

FR_Input_X        = f8
FR_Result         = f8

FR_r              = f8
FR_c              = f9

FR_norm_x         = f9
FR_inv_pi_2to63   = f10
FR_rshf_2to64     = f11
FR_2tom64         = f12
FR_rshf           = f13
FR_N_float_signif = f14
FR_abs_x          = f15
FR_Pi_by_4        = f34
FR_Two_to_M14     = f35
FR_Neg_Two_to_M14 = f36
FR_Two_to_M33     = f37
FR_Neg_Two_to_M33 = f38
FR_Neg_Two_to_M67 = f39
FR_Inv_pi_by_2    = f40
FR_N_float        = f41
FR_N_fix          = f42
FR_P_1            = f43
FR_P_2            = f44
FR_P_3            = f45
FR_s              = f46
FR_w              = f47
FR_d_2            = f48
FR_tmp_result     = f49
FR_Z              = f50
FR_A              = f51
FR_a              = f52
FR_t              = f53
FR_U_1            = f54
FR_U_2            = f55
FR_C_1            = f56
FR_C_2            = f57
FR_C_3            = f58
FR_C_4            = f59
FR_C_5            = f60
FR_S_1            = f61
FR_S_2            = f62
FR_S_3            = f63
FR_S_4            = f64
FR_S_5            = f65
FR_poly_hi        = f66
FR_poly_lo        = f67
FR_r_hi           = f68
FR_r_lo           = f69
FR_rsq            = f70
FR_r_cubed        = f71
FR_C_hi           = f72
FR_N_0            = f73
FR_d_1            = f74
FR_V              = f75
FR_V_hi           = f75
FR_V_lo           = f76
FR_U_hi           = f77
FR_U_lo           = f78
FR_U_hiabs        = f79
FR_V_hiabs        = f80
FR_PP_8           = f81
FR_QQ_8           = f101
FR_PP_7           = f82
FR_QQ_7           = f102
FR_PP_6           = f83
FR_QQ_6           = f103
FR_PP_5           = f84
FR_QQ_5           = f104
FR_PP_4           = f85
FR_QQ_4           = f105
FR_PP_3           = f86
FR_QQ_3           = f106
FR_PP_2           = f87
FR_QQ_2           = f107
FR_QQ_1           = f108
FR_r_hi_sq        = f88
FR_N_0_fix        = f89
FR_Inv_P_0        = f90
FR_corr           = f91
FR_poly           = f92
FR_Neg_Two_to_M3  = f93
FR_Two_to_M3      = f94
FR_P_0            = f95
FR_C_lo           = f96
FR_PP_1           = f97
FR_PP_1_lo        = f98
FR_ArgPrime       = f99
FR_inexact        = f100

GR_exp_m2_to_m3= r36
GR_N_Inc       = r37
GR_Sin_or_Cos  = r38
GR_signexp_x   = r40
GR_exp_x       = r40
GR_exp_mask    = r41
GR_exp_2_to_63 = r42
GR_exp_2_to_m3 = r43
GR_exp_2_to_24 = r44

GR_sig_inv_pi  = r45
GR_rshf_2to64  = r46
GR_exp_2tom64  = r47
GR_rshf        = r48
GR_ad_p        = r49
GR_ad_d        = r50
GR_ad_pp       = r51
GR_ad_qq       = r52
GR_ad_c        = r53
GR_ad_s        = r54
GR_ad_ce       = r55
GR_ad_se       = r56
GR_ad_m14      = r57
GR_ad_s1       = r58

// Added for unwind support

GR_SAVE_B0     = r39
GR_SAVE_GP     = r40
GR_SAVE_PFS    = r41


.section .text

GLOBAL_IEEE754_ENTRY(sinl)
{ .mlx
      alloc r32 = ar.pfs,0,27,2,0
      movl GR_sig_inv_pi = 0xa2f9836e4e44152a // significand of 1/pi
}
{ .mlx
      mov GR_Sin_or_Cos = 0x0
      movl GR_rshf_2to64 = 0x47e8000000000000 // 1.1000 2^(63+64)
}
;;

{ .mfi
      addl           GR_ad_p   = @ltoff(FSINCOSL_CONSTANTS#), gp
      fclass.m p6, p0 =  FR_Input_X, 0x1E3 // Test x natval, nan, inf
      mov GR_exp_2_to_m3 = 0xffff - 3      // Exponent of 2^-3
}
{ .mfb
      nop.m 999
      fnorm.s1 FR_norm_x = FR_Input_X      // Normalize x
      br.cond.sptk SINCOSL_CONTINUE
}
;;

GLOBAL_IEEE754_END(sinl)

GLOBAL_IEEE754_ENTRY(cosl)
{ .mlx
      alloc r32 = ar.pfs,0,27,2,0
      movl GR_sig_inv_pi = 0xa2f9836e4e44152a // significand of 1/pi
}
{ .mlx
      mov GR_Sin_or_Cos = 0x1
      movl GR_rshf_2to64 = 0x47e8000000000000 // 1.1000 2^(63+64)
}
;;

{ .mfi
      addl           GR_ad_p   = @ltoff(FSINCOSL_CONSTANTS#), gp
      fclass.m p6, p0 =  FR_Input_X, 0x1E3 // Test x natval, nan, inf
      mov GR_exp_2_to_m3 = 0xffff - 3      // Exponent of 2^-3
}
{ .mfi
      nop.m 999
      fnorm.s1 FR_norm_x = FR_Input_X      // Normalize x
      nop.i 999
}
;;

SINCOSL_CONTINUE:
{ .mfi
      setf.sig FR_inv_pi_2to63 = GR_sig_inv_pi // Form 1/pi * 2^63
      nop.f 999
      mov GR_exp_2tom64 = 0xffff - 64      // Scaling constant to compute N
}
{ .mlx
      setf.d FR_rshf_2to64 = GR_rshf_2to64    // Form const 1.1000 * 2^(63+64)
      movl GR_rshf = 0x43e8000000000000       // Form const 1.1000 * 2^63
}
;;

{ .mfi
      ld8 GR_ad_p = [GR_ad_p]              // Point to Inv_pi_by_2
      fclass.m p7, p0 = FR_Input_X, 0x0b   // Test x denormal
      nop.i 999
}
;;

{ .mfi
      getf.exp GR_signexp_x = FR_Input_X   // Get sign and exponent of x
      fclass.m p10, p0 = FR_Input_X, 0x007 // Test x zero
      nop.i 999
}
{ .mib
      mov GR_exp_mask = 0x1ffff            // Exponent mask
      nop.i 999
(p6)  br.cond.spnt SINCOSL_SPECIAL         // Branch if x natval, nan, inf
}
;;

{ .mfi
      setf.exp FR_2tom64 = GR_exp_2tom64   // Form 2^-64 for scaling N_float
      nop.f 0
      add GR_ad_d = 0x70, GR_ad_p          // Point to constant table d
}
{ .mib
      setf.d FR_rshf = GR_rshf         // Form right shift const 1.1000 * 2^63
      mov  GR_exp_m2_to_m3 = 0x2fffc       // Form -(2^-3)
(p7)  br.cond.spnt SINCOSL_DENORMAL        // Branch if x denormal
}
;;

SINCOSL_COMMON:
{ .mfi
      and GR_exp_x = GR_exp_mask, GR_signexp_x // Get exponent of x
      fclass.nm p8, p0 = FR_Input_X, 0x1FF // Test x unsupported type
      mov GR_exp_2_to_63 = 0xffff + 63     // Exponent of 2^63
}
{ .mib
      add GR_ad_pp = 0x40, GR_ad_d         // Point to constant table pp
      mov GR_exp_2_to_24 = 0xffff + 24     // Exponent of 2^24
(p10) br.cond.spnt SINCOSL_ZERO            // Branch if x zero
}
;;

{ .mfi
      ldfe FR_Inv_pi_by_2 = [GR_ad_p], 16  // Load 2/pi
      fcmp.eq.s0 p15, p0 = FR_Input_X, f0  // Dummy to set denormal
      add GR_ad_qq = 0xa0, GR_ad_pp        // Point to constant table qq
}
{ .mfi
      ldfe FR_Pi_by_4 = [GR_ad_d], 16      // Load pi/4 for range test
      nop.f 999
      cmp.ge p10,p0 = GR_exp_x, GR_exp_2_to_63   // Is |x| >= 2^63
}
;;

{ .mfi
      ldfe FR_P_0 = [GR_ad_p], 16          // Load P_0 for pi/4 <= |x| < 2^63
      fmerge.s FR_abs_x = f1, FR_norm_x    // |x|
      add GR_ad_c = 0x90, GR_ad_qq         // Point to constant table c
}
{ .mfi
      ldfe FR_Inv_P_0 = [GR_ad_d], 16      // Load 1/P_0 for pi/4 <= |x| < 2^63
      nop.f 999
      cmp.ge p7,p0 = GR_exp_x, GR_exp_2_to_24   // Is |x| >= 2^24
}
;;

{ .mfi
      ldfe FR_P_1 = [GR_ad_p], 16          // Load P_1 for pi/4 <= |x| < 2^63
      nop.f 999
      add GR_ad_s = 0x50, GR_ad_c          // Point to constant table s
}
{ .mfi
      ldfe FR_PP_8 = [GR_ad_pp], 16        // Load PP_8 for 2^-3 < |r| < pi/4
      nop.f 999
      nop.i 999
}
;;

{ .mfi
      ldfe FR_P_2 = [GR_ad_p], 16          // Load P_2 for pi/4 <= |x| < 2^63
      nop.f 999
      add GR_ad_ce = 0x40, GR_ad_c         // Point to end of constant table c
}
{ .mfi
      ldfe FR_QQ_8 = [GR_ad_qq], 16        // Load QQ_8 for 2^-3 < |r| < pi/4
      nop.f 999
      nop.i 999
}
;;

{ .mfi
      ldfe FR_QQ_7 = [GR_ad_qq], 16        // Load QQ_7 for 2^-3 < |r| < pi/4
      fma.s1        FR_N_float_signif = FR_Input_X, FR_inv_pi_2to63, FR_rshf_2to64
      add GR_ad_se = 0x40, GR_ad_s         // Point to end of constant table s
}
{ .mib
      ldfe FR_PP_7 = [GR_ad_pp], 16        // Load PP_7 for 2^-3 < |r| < pi/4
      mov GR_ad_s1 = GR_ad_s               // Save pointer to S_1
(p10) br.cond.spnt SINCOSL_ARG_TOO_LARGE   // Branch if |x| >= 2^63
                                           // Use Payne-Hanek Reduction
}
;;

{ .mfi
      ldfe FR_P_3 = [GR_ad_p], 16          // Load P_3 for pi/4 <= |x| < 2^63
      fmerge.se FR_r = FR_norm_x, FR_norm_x // r = x, in case |x| < pi/4
      add GR_ad_m14 = 0x50, GR_ad_s        // Point to constant table m14
}
{ .mfb
      ldfps FR_Two_to_M3, FR_Neg_Two_to_M3 = [GR_ad_d], 8
      fma.s1 FR_rsq = FR_norm_x, FR_norm_x, f0 // rsq = x*x, in case |x| < pi/4
(p7)  br.cond.spnt SINCOSL_LARGER_ARG      // Branch if 2^24 <= |x| < 2^63
                                           // Use pre-reduction
}
;;

{ .mmf
      ldfe FR_PP_6 = [GR_ad_pp], 16       // Load PP_6 for normal path
      ldfe FR_QQ_6 = [GR_ad_qq], 16       // Load QQ_6 for normal path
      fmerge.se FR_c = f0, f0             // c = 0 in case |x| < pi/4
}
;;

{ .mmf
      ldfe FR_PP_5 = [GR_ad_pp], 16       // Load PP_5 for normal path
      ldfe FR_QQ_5 = [GR_ad_qq], 16       // Load QQ_5 for normal path
      nop.f 999
}
;;

// Here if 0 < |x| < 2^24
{ .mfi
      ldfe FR_S_5 = [GR_ad_se], -16       // Load S_5 if i_1=0
      fcmp.lt.s1  p6, p7 = FR_abs_x, FR_Pi_by_4  // Test |x| < pi/4
      nop.i 999
}
{ .mfi
      ldfe FR_C_5 = [GR_ad_ce], -16       // Load C_5 if i_1=1
      fms.s1 FR_N_float = FR_N_float_signif, FR_2tom64, FR_rshf
      nop.i 999
}
;;

{ .mmi
      ldfe FR_S_4 = [GR_ad_se], -16       // Load S_4 if i_1=0
      ldfe FR_C_4 = [GR_ad_ce], -16       // Load C_4 if i_1=1
      nop.i 999
}
;;

//
//     N  = Arg * 2/pi
//     Check if Arg < pi/4
//
//
//     Case 2: Convert integer N_fix back to normalized floating-point value.
//     Case 1: p8 is only affected  when p6 is set
//
//
//     Grab the integer part of N and call it N_fix
//
{ .mfi
(p7)  ldfps FR_Two_to_M33, FR_Neg_Two_to_M33 = [GR_ad_d], 8
(p6)  fma.s1 FR_r_cubed = FR_r, FR_rsq, f0        // r^3 if |x| < pi/4
(p6)  mov GR_N_Inc = GR_Sin_or_Cos                // N_Inc if |x| < pi/4
}
;;

//     If |x| < pi/4, r = x and c = 0
//     lf |x| < pi/4, is x < 2**(-3).
//     r = Arg
//     c = 0
{ .mmi
(p7)  getf.sig        GR_N_Inc = FR_N_float_signif
(p6)  cmp.lt.unc p8,p0 = GR_exp_x, GR_exp_2_to_m3   // Is |x| < 2^-3
(p6)  tbit.z p9,p10 = GR_N_Inc, 0         // p9  if i_1=0, N mod 4 = 0,1
                                          // p10 if i_1=1, N mod 4 = 2,3
}
;;

//
//     lf |x| < pi/4, is -2**(-3)< x < 2**(-3) - set p8.
//     If |x| >= pi/4,
//     Create the right N for |x| < pi/4 and otherwise
//     Case 2: Place integer part of N in GP register
//


{ .mbb
      nop.m 999
(p8)  br.cond.spnt SINCOSL_SMALL_R_0    // Branch if 0 < |x| < 2^-3
(p6)  br.cond.spnt SINCOSL_NORMAL_R_0   // Branch if 2^-3 <= |x| < pi/4
}
;;

// Here if pi/4 <= |x| < 2^24
{ .mfi
      ldfs FR_Neg_Two_to_M67 = [GR_ad_d], 8     // Load -2^-67
      fnma.s1 FR_s = FR_N_float, FR_P_1, FR_Input_X // s = -N * P_1  + Arg
      add GR_N_Inc = GR_N_Inc, GR_Sin_or_Cos    // Adjust N_Inc for sin/cos
}
{ .mfi
      nop.m 999
      fma.s1 FR_w = FR_N_float, FR_P_2, f0      // w = N * P_2
      nop.i 999
}
;;

{ .mfi
      nop.m 999
      fms.s1 FR_r = FR_s, f1, FR_w        // r = s - w, assume |s| >= 2^-33
      tbit.z p9,p10 = GR_N_Inc, 0         // p9  if i_1=0, N mod 4 = 0,1
                                          // p10 if i_1=1, N mod 4 = 2,3
}
;;

{ .mfi
      nop.m 999
      fcmp.lt.s1 p7, p6 = FR_s, FR_Two_to_M33
      nop.i 999
}
;;

{ .mfi
      nop.m 999
(p7)  fcmp.gt.s1 p7, p6 = FR_s, FR_Neg_Two_to_M33 // p6 if |s| >= 2^-33, else p7
      nop.i 999
}
;;

{ .mfi
      nop.m 999
      fms.s1 FR_c = FR_s, f1, FR_r             // c = s - r, for |s| >= 2^-33
      nop.i 999
}
{ .mfi
      nop.m 999
      fma.s1 FR_rsq = FR_r, FR_r, f0           // rsq = r * r, for |s| >= 2^-33
      nop.i 999
}
;;

{ .mfi
      nop.m 999
(p7)  fma.s1 FR_w = FR_N_float, FR_P_3, f0
      nop.i 999
}
;;

{ .mmf
(p9)  ldfe FR_C_1 = [GR_ad_pp], 16     // Load C_1 if i_1=0
(p10) ldfe FR_S_1 = [GR_ad_qq], 16     // Load S_1 if i_1=1
      frcpa.s1 FR_r_hi, p15 = f1, FR_r  // r_hi = frcpa(r)
}
;;

{ .mfi
      nop.m 999
(p6)  fcmp.lt.unc.s1 p8, p13 = FR_r, FR_Two_to_M3 // If big s, test r with 2^-3
      nop.i 999
}
;;

{ .mfi
      nop.m 999
(p7)  fma.s1 FR_U_1 = FR_N_float, FR_P_2, FR_w
      nop.i 999
}
;;

//
//     For big s: r = s - w: No futher reduction is necessary
//     For small s: w = N * P_3 (change sign) More reduction
//
{ .mfi
        nop.m 999
(p8)   fcmp.gt.s1 p8, p13 = FR_r, FR_Neg_Two_to_M3 // If big s, p8 if |r| < 2^-3
        nop.i 999 ;;
}

{ .mfi
      nop.m 999
(p9)  fma.s1 FR_poly = FR_rsq, FR_PP_8, FR_PP_7 // poly = rsq*PP_8+PP_7 if i_1=0
      nop.i 999
}
{ .mfi
      nop.m 999
(p10) fma.s1 FR_poly = FR_rsq, FR_QQ_8, FR_QQ_7 // poly = rsq*QQ_8+QQ_7 if i_1=1
      nop.i 999
}
;;

{ .mfi
        nop.m 999
(p7)   fms.s1 FR_r = FR_s, f1, FR_U_1
        nop.i 999
}
;;

{ .mfi
      nop.m 999
(p6)  fma.s1 FR_r_cubed = FR_r, FR_rsq, f0  // rcubed = r * rsq
      nop.i 999
}
;;

{ .mfi
//
//     For big s: Is |r| < 2**(-3)?
//     For big s: c = S - r
//     For small s: U_1 = N * P_2 + w
//
//     If p8 is set, prepare to branch to Small_R.
//     If p9 is set, prepare to branch to Normal_R.
//     For big s,  r is complete here.
//
//
//     For big s: c = c + w (w has not been negated.)
//     For small s: r = S - U_1
//
      nop.m 999
(p6)  fms.s1 FR_c = FR_c, f1, FR_w
      nop.i 999
}
{ .mbb
      nop.m 999
(p8)  br.cond.spnt    SINCOSL_SMALL_R_1  // Branch if |s|>=2^-33, |r| < 2^-3,
                                         // and pi/4 <= |x| < 2^24
(p13) br.cond.sptk    SINCOSL_NORMAL_R_1 // Branch if |s|>=2^-33, |r| >= 2^-3,
                                         // and pi/4 <= |x| < 2^24
}
;;

SINCOSL_S_TINY:
//
// Here if |s| < 2^-33, and pi/4 <= |x| < 2^24
//
{ .mfi
       fms.s1 FR_U_2 = FR_N_float, FR_P_2, FR_U_1
//
//     c = S - U_1
//     r = S_1 * r
//
//
}
;;

{ .mmi
        nop.m 999
//
//     Get [i_0,i_1] - two lsb of N_fix_gr.
//     Do dummy fmpy so inexact is always set.
//
      tbit.z p9,p10 = GR_N_Inc, 0      // p9  if i_1=0, N mod 4 = 0,1
                                       // p10 if i_1=1, N mod 4 = 2,3
}
;;

//
//     For small s: U_2 = N * P_2 - U_1
//     S_1 stored constant - grab the one stored with the
//     coefficients.
//
{ .mfi
       ldfe FR_S_1 = [GR_ad_s1], 16
//
//     Check if i_1 and i_0  != 0
//
(p10)  fma.s1        FR_poly = f0, f1, FR_Neg_Two_to_M67
      tbit.z p11,p12 = GR_N_Inc, 1     // p11 if i_0=0, N mod 4 = 0,2