.file "libm_sincos_large.s"// Copyright (c) 2002 - 2003, Intel Corporation// All rights reserved.//// Contributed 2002 by the Intel Numerics Group, Intel Corporation//// Redistribution and use in source and binary forms, with or without// modification, are permitted provided that the following conditions are// met://// * Redistributions of source code must retain the above copyright// notice, this list of conditions and the following disclaimer.//// * Redistributions in binary form must reproduce the above copyright// notice, this list of conditions and the following disclaimer in the// documentation and/or other materials provided with the distribution.//// * The name of Intel Corporation may not be used to endorse or promote// products derived from this software without specific prior written// permission.// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS// "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT// LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR// A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL INTEL OR ITS// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY// OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY OR TORT (INCLUDING// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.//// Intel Corporation is the author of this code, and requests that all// problem reports or change requests be submitted to it directly at// http://www.intel.com/software/products/opensource/libraries/num.htm.//// History//==============================================================// 02/15/02 Initial version// 05/13/02 Changed interface to __libm_pi_by_2_reduce// 02/10/03 Reordered header: .section, .global, .proc, .align;//          used data8 for long double table values// 05/15/03 Reformatted data tables////// Overview of operation//==============================================================//// These functions calculate the sin and cos for inputs// greater than 2^10//// __libm_sin_large#// __libm_cos_large#// They accept argument in f8// and return result in f8 without final rounding//// __libm_sincos_large#// It accepts argument in f8// and returns cos in f8 and sin in f9 without final rounding//////*********************************************************************//// Accuracy:       Within .7 ulps for 80-bit floating point values//                 Very accurate for double precision values////*********************************************************************//// Resources Used:////    Floating-Point Registers: f8 as Input Value, f8 and f9 as Return Values//                              f32-f103////    General Purpose Registers://      r32-r43//      r44-r45 (Used to pass arguments to pi_by_2 reduce routine)////    Predicate Registers:      p6-p13////*********************************************************************////  IEEE Special Conditions:////    Denormal  fault raised on denormal inputs//    Overflow exceptions do not occur//    Underflow exceptions raised when appropriate for sin//    (No specialized error handling for this routine)//    Inexact raised when appropriate by algorithm////    sin(SNaN) = QNaN//    sin(QNaN) = QNaN//    sin(inf) = QNaN//    sin(+/-0) = +/-0//    cos(inf) = QNaN//    cos(SNaN) = QNaN//    cos(QNaN) = QNaN//    cos(0) = 1////*********************************************************************////  Mathematical Description//  ========================////  The computation of FSIN and FCOS is best handled in one piece of//  code. The main reason is that given any argument Arg, computation//  of trigonometric functions first calculate N and an approximation//  to alpha where////  Arg = N pi/2 + alpha, |alpha| <= pi/4.////  Since////  cos( Arg ) = sin( (N+1) pi/2 + alpha ),////  therefore, the code for computing sine will produce cosine as long//  as 1 is added to N immediately after the argument reduction//  process.////  Let M = N if sine//      N+1 if cosine.////  Now, given////  Arg = M pi/2  + alpha, |alpha| <= pi/4,////  let I = M mod 4, or I be the two lsb of M when M is represented//  as 2's complement. I = [i_0 i_1]. Then////  sin( Arg ) = (-1)^i_0  sin( alpha ) if i_1 = 0,//             = (-1)^i_0  cos( alpha )     if i_1 = 1.////  For example://       if M = -1, I = 11//         sin ((-pi/2 + alpha) = (-1) cos (alpha)//       if M = 0, I = 00//         sin (alpha) = sin (alpha)//       if M = 1, I = 01//         sin (pi/2 + alpha) = cos (alpha)//       if M = 2, I = 10//         sin (pi + alpha) = (-1) sin (alpha)//       if M = 3, I = 11//         sin ((3/2)pi + alpha) = (-1) cos (alpha)////  The value of alpha is obtained by argument reduction and//  represented by two working precision numbers r and c where////  alpha =  r  +  c     accurately.////  The reduction method is described in a previous write up.//  The argument reduction scheme identifies 4 cases. For Cases 2//  and 4, because |alpha| is small, sin(r+c) and cos(r+c) can be//  computed very easily by 2 or 3 terms of the Taylor series//  expansion as follows:////  Case 2://  -------////  sin(r + c) = r + c - r^3/6  accurately//  cos(r + c) = 1 - 2^(-67)    accurately////  Case 4://  -------////  sin(r + c) = r + c - r^3/6 + r^5/120    accurately//  cos(r + c) = 1 - r^2/2 + r^4/24     accurately////  The only cases left are Cases 1 and 3 of the argument reduction//  procedure. These two cases will be merged since after the//  argument is reduced in either cases, we have the reduced argument//  represented as r + c and that the magnitude |r + c| is not small//  enough to allow the usage of a very short approximation.////  The required calculation is either////  sin(r + c)  =  sin(r)  +  correction,  or//  cos(r + c)  =  cos(r)  +  correction.////  Specifically,////  sin(r + c) = sin(r) + c sin'(r) + O(c^2)//         = sin(r) + c cos (r) + O(c^2)//         = sin(r) + c(1 - r^2/2)  accurately.//  Similarly,////  cos(r + c) = cos(r) - c sin(r) + O(c^2)//         = cos(r) - c(r - r^3/6)  accurately.////  We therefore concentrate on accurately calculating sin(r) and//  cos(r) for a working-precision number r, |r| <= pi/4 to within//  0.1% or so.////  The greatest challenge of this task is that the second terms of//  the Taylor series////  r - r^3/3! + r^r/5! - ...////  and////  1 - r^2/2! + r^4/4! - ...////  are not very small when |r| is close to pi/4 and the rounding//  errors will be a concern if simple polynomial accumulation is//  used. When |r| < 2^-3, however, the second terms will be small//  enough (6 bits or so of right shift) that a normal Horner//  recurrence suffices. Hence there are two cases that we consider//  in the accurate computation of sin(r) and cos(r), |r| <= pi/4.////  Case small_r: |r| < 2^(-3)//  --------------------------////  Since Arg = M pi/4 + r + c accurately, and M mod 4 is [i_0 i_1],//  we have////  sin(Arg) = (-1)^i_0 * sin(r + c)    if i_1 = 0//       = (-1)^i_0 * cos(r + c)    if i_1 = 1////  can be accurately approximated by////  sin(Arg) = (-1)^i_0 * [sin(r) + c]  if i_1 = 0//           = (-1)^i_0 * [cos(r) - c*r] if i_1 = 1////  because |r| is small and thus the second terms in the correction//  are unneccessary.////  Finally, sin(r) and cos(r) are approximated by polynomials of//  moderate lengths.////  sin(r) =  r + S_1 r^3 + S_2 r^5 + ... + S_5 r^11//  cos(r) =  1 + C_1 r^2 + C_2 r^4 + ... + C_5 r^10////  We can make use of predicates to selectively calculate//  sin(r) or cos(r) based on i_1.////  Case normal_r: 2^(-3) <= |r| <= pi/4//  ------------------------------------////  This case is more likely than the previous one if one considers//  r to be uniformly distributed in [-pi/4 pi/4]. Again,////  sin(Arg) = (-1)^i_0 * sin(r + c)    if i_1 = 0//           = (-1)^i_0 * cos(r + c)    if i_1 = 1.////  Because |r| is now larger, we need one extra term in the//  correction. sin(Arg) can be accurately approximated by////  sin(Arg) = (-1)^i_0 * [sin(r) + c(1-r^2/2)]      if i_1 = 0//           = (-1)^i_0 * [cos(r) - c*r*(1 - r^2/6)]    i_1 = 1.////  Finally, sin(r) and cos(r) are approximated by polynomials of//  moderate lengths.////  sin(r) =  r + PP_1_hi r^3 + PP_1_lo r^3 +//                PP_2 r^5 + ... + PP_8 r^17////  cos(r) =  1 + QQ_1 r^2 + QQ_2 r^4 + ... + QQ_8 r^16////  where PP_1_hi is only about 16 bits long and QQ_1 is -1/2.//  The crux in accurate computation is to calculate////  r + PP_1_hi r^3   or  1 + QQ_1 r^2////  accurately as two pieces: U_hi and U_lo. The way to achieve this//  is to obtain r_hi as a 10 sig. bit number that approximates r to//  roughly 8 bits or so of accuracy. (One convenient way is////  r_hi := frcpa( frcpa( r ) ).)////  This way,////  r + PP_1_hi r^3 =  r + PP_1_hi r_hi^3 +//                          PP_1_hi (r^3 - r_hi^3)//              =  [r + PP_1_hi r_hi^3]  +//             [PP_1_hi (r - r_hi)//                (r^2 + r_hi r + r_hi^2) ]//              =  U_hi  +  U_lo////  Since r_hi is only 10 bit long and PP_1_hi is only 16 bit long,//  PP_1_hi * r_hi^3 is only at most 46 bit long and thus computed//  exactly. Furthermore, r and PP_1_hi r_hi^3 are of opposite sign//  and that there is no more than 8 bit shift off between r and//  PP_1_hi * r_hi^3. Hence the sum, U_hi, is representable and thus//  calculated without any error. Finally, the fact that////  |U_lo| <= 2^(-8) |U_hi|////  says that U_hi + U_lo is approximating r + PP_1_hi r^3 to roughly//  8 extra bits of accuracy.////  Similarly,////  1 + QQ_1 r^2  =  [1 + QQ_1 r_hi^2]  +//                      [QQ_1 (r - r_hi)(r + r_hi)]//            =  U_hi  +  U_lo.////  Summarizing, we calculate r_hi = frcpa( frcpa( r ) ).////  If i_1 = 0, then////    U_hi := r + PP_1_hi * r_hi^3//    U_lo := PP_1_hi * (r - r_hi) * (r^2 + r*r_hi + r_hi^2)//    poly := PP_1_lo r^3 + PP_2 r^5 + ... + PP_8 r^17//    correction := c * ( 1 + C_1 r^2 )////  Else ...i_1 = 1////    U_hi := 1 + QQ_1 * r_hi * r_hi//    U_lo := QQ_1 * (r - r_hi) * (r + r_hi)//    poly := QQ_2 * r^4 + QQ_3 * r^6 + ... + QQ_8 r^16//    correction := -c * r * (1 + S_1 * r^2)////  End////  Finally,////  V := poly + ( U_lo + correction )////                 /    U_hi  +  V         if i_0 = 0//  result := |//                 \  (-U_hi) -  V         if i_0 = 1////  It is important that in the last step, negation of U_hi is//  performed prior to the subtraction which is to be performed in//  the user-set rounding mode.//////  Algorithmic Description//  =======================////  The argument reduction algorithm is tightly integrated into FSIN//  and FCOS which share the same code. The following is complete and//  self-contained. The argument reduction description given//  previously is repeated below.//////  Step 0. Initialization.////   If FSIN is invoked, set N_inc := 0; else if FCOS is invoked,//   set N_inc := 1.////  Step 1. Check for exceptional and special cases.////   * If Arg is +-0, +-inf, NaN, NaT, go to Step 10 for special//     handling.//   * If |Arg| < 2^24, go to Step 2 for reduction of moderate//     arguments. This is the most likely case.//   * If |Arg| < 2^63, go to Step 8 for pre-reduction of large//     arguments.//   * If |Arg| >= 2^63, go to Step 10 for special handling.////  Step 2. Reduction of moderate arguments.////  If |Arg| < pi/4     ...quick branch//     N_fix := N_inc   (integer)//     r     := Arg//     c     := 0.0//     Branch to Step 4, Case_1_complete//  Else        ...cf. argument reduction//     N     := Arg * two_by_PI (fp)//     N_fix := fcvt.fx( N )    (int)//     N     := fcvt.xf( N_fix )//     N_fix := N_fix + N_inc//     s     := Arg - N * P_1   (first piece of pi/2)//     w     := -N * P_2    (second piece of pi/2)////     If |s| >= 2^(-33)//        go to Step 3, Case_1_reduce//     Else//        go to Step 7, Case_2_reduce//     Endif//  Endif////  Step 3. Case_1_reduce.////  r := s + w//  c := (s - r) + w    ...observe order////  Step 4. Case_1_complete////  ...At this point, the reduced argument alpha is//  ...accurately represented as r + c.//  If |r| < 2^(-3), go to Step 6, small_r.////  Step 5. Normal_r.////  Let [i_0 i_1] by the 2 lsb of N_fix.//  FR_rsq  := r * r//  r_hi := frcpa( frcpa( r ) )//  r_lo := r - r_hi////  If i_1 = 0, then//    poly := r*FR_rsq*(PP_1_lo + FR_rsq*(PP_2 + ... FR_rsq*PP_8))//    U_hi := r + PP_1_hi*r_hi*r_hi*r_hi    ...any order//    U_lo := PP_1_hi*r_lo*(r*r + r*r_hi + r_hi*r_hi)//    correction := c + c*C_1*FR_rsq        ...any order//  Else//    poly := FR_rsq*FR_rsq*(QQ_2 + FR_rsq*(QQ_3 + ... + FR_rsq*QQ_8))//    U_hi := 1 + QQ_1 * r_hi * r_hi        ...any order//    U_lo := QQ_1 * r_lo * (r + r_hi)//    correction := -c*(r + S_1*FR_rsq*r)   ...any order//  Endif////  V := poly + (U_lo + correction) ...observe order////  result := (i_0 == 0?   1.0 : -1.0)////  Last instruction in user-set rounding mode////  result := (i_0 == 0?   result*U_hi + V ://                        result*U_hi - V)////  Return////  Step 6. Small_r.////  ...Use flush to zero mode without causing exception//    Let [i_0 i_1] be the two lsb of N_fix.////  FR_rsq := r * r////  If i_1 = 0 then//     z := FR_rsq*FR_rsq; z := FR_rsq*z *r//     poly_lo := S_3 + FR_rsq*(S_4 + FR_rsq*S_5)//     poly_hi := r*FR_rsq*(S_1 + FR_rsq*S_2)//     correction := c//     result := r//  Else//     z := FR_rsq*FR_rsq; z := FR_rsq*z//     poly_lo := C_3 + FR_rsq*(C_4 + FR_rsq*C_5)//     poly_hi := FR_rsq*(C_1 + FR_rsq*C_2)//     correction := -c*r//     result := 1//  Endif////  poly := poly_hi + (z * poly_lo + correction)////  If i_0 = 1, result := -result////  Last operation. Perform in user-set rounding mode////  result := (i_0 == 0?     result + poly ://                          result - poly )//  Return////  Step 7. Case_2_reduce.////  ...Refer to the write up for argument reduction for//  ...rationale. The reduction algorithm below is taken from//  ...argument reduction description and integrated this.////  w := N*P_3//  U_1 := N*P_2 + w        ...FMA//  U_2 := (N*P_2 - U_1) + w    ...2 FMA//  ...U_1 + U_2 is  N*(P_2+P_3) accurately////  r := s - U_1//  c := ( (s - r) - U_1 ) - U_2////  ...The mathematical sum r + c approximates the reduced//  ...argument accurately. Note that although compared to//  ...Case 1, this case requires much more work to reduce//  ...the argument, the subsequent calculation needed for//  ...any of the trigonometric function is very little because//  ...|alpha| < 1.01*2^(-33) and thus two terms of the//  ...Taylor series expansion suffices.////  If i_1 = 0 then//     poly := c + S_1 * r * r * r  ...any order//     result := r//  Else//     poly := -2^(-67)//     result := 1.0//  Endif////  If i_0 = 1, result := -result////  Last operation. Perform in user-set rounding mode////  result := (i_0 == 0?     result + poly ://                           result - poly )////  Return//////  Step 8. Pre-reduction of large arguments.////  ...Again, the following reduction procedure was described//  ...in the separate write up for argument reduction, which//  ...is tightly integrated here.//  N_0 := Arg * Inv_P_0//  N_0_fix := fcvt.fx( N_0 )//  N_0 := fcvt.xf( N_0_fix)//  Arg' := Arg - N_0 * P_0//  w := N_0 * d_1//  N := Arg' * two_by_PI//  N_fix := fcvt.fx( N )//  N := fcvt.xf( N_fix )//  N_fix := N_fix + N_inc////  s := Arg' - N * P_1//  w := w - N * P_2////  If |s| >= 2^(-14)//     go to Step 3//  Else//     go to Step 9//  Endif////  Step 9. Case_4_reduce.////    ...first obtain N_0*d_1 and -N*P_2 accurately//   U_hi := N_0 * d_1      V_hi := -N*P_2//   U_lo := N_0 * d_1 - U_hi   V_lo := -N*P_2 - U_hi   ...FMAs////   ...compute the contribution from N_0*d_1 and -N*P_3//   w := -N*P_3//   w := w + N_0*d_2//   t := U_lo + V_lo + w       ...any order////   ...at this point, the mathematical value//   ...s + U_hi + V_hi  + t approximates the true reduced argument//   ...accurately. Just need to compute this accurately.////   ...Calculate U_hi + V_hi accurately://   A := U_hi + V_hi//   if |U_hi| >= |V_hi| then//      a := (U_hi - A) + V_hi//   else//      a := (V_hi - A) + U_hi//   endif//   ...order in computing "a" must be observed. This branch is//   ...best implemented by predicates.//   ...A + a  is U_hi + V_hi accurately. Moreover, "a" is//   ...much smaller than A: |a| <= (1/2)ulp(A).////   ...Just need to calculate   s + A + a + t//   C_hi := s + A      t := t + a//   C_lo := (s - C_hi) + A//   C_lo := C_lo + t////   ...Final steps for reduction//   r := C_hi + C_lo//   c := (C_hi - r) + C_lo////   ...At this point, we have r and c//   ...And all we need is a couple of terms of the corresponding//   ...Taylor series.////   If i_1 = 0//      poly := c + r*FR_rsq*(S_1 + FR_rsq*S_2)//      result := r//   Else//      poly := FR_rsq*(C_1 + FR_rsq*C_2)//      result := 1//   Endif////   If i_0 = 1, result := -result////   Last operation. Perform in user-set rounding mode////   result := (i_0 == 0?     result + poly ://                            result - poly )//   Return////   Large Arguments: For arguments above 2**63, a Payne-Hanek//   style argument reduction is used and pi_by_2 reduce is called.//RODATA.align 16LOCAL_OBJECT_START(FSINCOS_CONSTANTS)data4 0x4B800000 // two**24data4 0xCB800000 // -two**24data4 0x00000000 // paddata4 0x00000000 // paddata8 0xA2F9836E4E44152A, 0x00003FFE // Inv_pi_by_2data8 0xC84D32B0CE81B9F1, 0x00004016 // P_0data8 0xC90FDAA22168C235, 0x00003FFF // P_1data8 0xECE675D1FC8F8CBB, 0x0000BFBD // P_2data8 0xB7ED8FBBACC19C60, 0x0000BF7C // P_3data4 0x5F000000 // two**63data4 0xDF000000 // -two**63data4 0x00000000 // paddata4 0x00000000 // paddata8 0xA397E5046EC6B45A, 0x00003FE7 // Inv_P_0data8 0x8D848E89DBD171A1, 0x0000BFBF // d_1data8 0xD5394C3618A66F8E, 0x0000BF7C // d_2data8 0xC90FDAA22168C234, 0x00003FFE // pi_by_4data8 0xC90FDAA22168C234, 0x0000BFFE // neg_pi_by_4data4 0x3E000000 // two**-3data4 0xBE000000 // -two**-3data4 0x00000000 // paddata4 0x00000000 // paddata4 0x2F000000 // two**-33data4 0xAF000000 // -two**-33data4 0x9E000000 // -two**-67data4 0x00000000 // paddata8 0xCC8ABEBCA21C0BC9, 0x00003FCE // PP_8data8 0xD7468A05720221DA, 0x0000BFD6 // PP_7data8 0xB092382F640AD517, 0x00003FDE // PP_6data8 0xD7322B47D1EB75A4, 0x0000BFE5 // PP_5data8 0xFFFFFFFFFFFFFFFE, 0x0000BFFD // C_1data8 0xAAAA000000000000, 0x0000BFFC // PP_1_hidata8 0xB8EF1D2ABAF69EEA, 0x00003FEC // PP_4data8 0xD00D00D00D03BB69, 0x0000BFF2 // PP_3data8 0x8888888888888962, 0x00003FF8 // PP_2data8 0xAAAAAAAAAAAB0000, 0x0000BFEC // PP_1_lodata8 0xD56232EFC2B0FE52, 0x00003FD2 // QQ_8data8 0xC9C99ABA2B48DCA6, 0x0000BFDA // QQ_7data8 0x8F76C6509C716658, 0x00003FE2 // QQ_6data8 0x93F27DBAFDA8D0FC, 0x0000BFE9 // QQ_5data8 0xAAAAAAAAAAAAAAAA, 0x0000BFFC // S_1data8 0x8000000000000000, 0x0000BFFE // QQ_1data8 0xD00D00D00C6E5041, 0x00003FEF // QQ_4data8 0xB60B60B60B607F60, 0x0000BFF5 // QQ_3data8 0xAAAAAAAAAAAAAA9B, 0x00003FFA // QQ_2data8 0xFFFFFFFFFFFFFFFE, 0x0000BFFD // C_1data8 0xAAAAAAAAAAAA719F, 0x00003FFA // C_2data8 0xB60B60B60356F994, 0x0000BFF5 // C_3data8 0xD00CFFD5B2385EA9, 0x00003FEF // C_4data8 0x93E4BD18292A14CD, 0x0000BFE9 // C_5data8 0xAAAAAAAAAAAAAAAA, 0x0000BFFC // S_1data8 0x88888888888868DB, 0x00003FF8 // S_2data8 0xD00D00D0055EFD4B, 0x0000BFF2 // S_3data8 0xB8EF1C5D839730B9, 0x00003FEC // S_4data8 0xD71EA3A4E5B3F492, 0x0000BFE5 // S_5data4 0x38800000 // two**-14data4 0xB8800000 // -two**-14LOCAL_OBJECT_END(FSINCOS_CONSTANTS)// sin and cos registers// FRFR_Input_X        = f8FR_r              = f8FR_c              = f9FR_Two_to_63      = f32FR_Two_to_24      = f33FR_Pi_by_4        = f33FR_Two_to_M14     = f34FR_Two_to_M33     = f35FR_Neg_Two_to_24  = f36FR_Neg_Pi_by_4    = f36FR_Neg_Two_to_M14 = f37FR_Neg_Two_to_M33 = f38FR_Neg_Two_to_M67 = f39FR_Inv_pi_by_2    = f40FR_N_float        = f41FR_N_fix          = f42FR_P_1            = f43FR_P_2            = f44FR_P_3            = f45FR_s              = f46FR_w              = f47FR_d_2            = f48FR_prelim         = f49FR_Z              = f50FR_A              = f51FR_a              = f52FR_t              = f53FR_U_1            = f54FR_U_2            = f55FR_C_1            = f56FR_C_2            = f57FR_C_3            = f58FR_C_4            = f59FR_C_5            = f60FR_S_1            = f61FR_S_2            = f62FR_S_3            = f63FR_S_4            = f64FR_S_5            = f65FR_poly_hi        = f66FR_poly_lo        = f67FR_r_hi           = f68FR_r_lo           = f69FR_rsq            = f70FR_r_cubed        = f71FR_C_hi           = f72FR_N_0            = f73FR_d_1            = f74