import later fedora-branch tweaks

[glibc.git] / sysdeps / ia64 / fpu / s_cosl.S

.file "sincosl.s"

// Copyright (C) 2000, 2001, Intel Corporation
// All rights reserved.
// 
// Contributed 2/2/2000 by John Harrison, Ted Kubaska, Bob Norin, Shane Story,
// and Ping Tak Peter Tang of the Computational Software Lab, 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://developer.intel.com/opensource.
//
// *********************************************************************
//
// History: 
// 2/02/2000 (hand-optimized)
// 4/04/00  Unwind support added
//
// *********************************************************************
//
// Function:   Combined sinl(x) and cosl(x), where
//
//             sinl(x) = sine(x), for double-extended precision x values
//             cosl(x) = cosine(x), for double-extended precision x values
//
// *********************************************************************
//
// Resources Used:
//
//    Floating-Point Registers: f8 (Input and Return Value) 
//                              f32-f99
//
//    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
//
//    sinl(SNaN) = QNaN
//    sinl(QNaN) = QNaN
//    sinl(inf) = QNaN 
//    sinl(+/-0) = +/-0
//    cosl(inf) = QNaN 
//    cosl(SNaN) = QNaN
//    cosl(QNaN) = QNaN
//    cosl(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
//
//  cosl( Arg ) = sinl( (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
//
//  sinl( Arg ) = (-1)^i_0  sinl( alpha )       if i_1 = 0,
//             = (-1)^i_0  cosl( 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, sinl(r+c) and cosl(r+c) can be 
//  computed very easily by 2 or 3 terms of the Taylor series 
//  expansion as follows:
//
//  Case 2:
//  -------
//
//  sinl(r + c) = r + c - r^3/6 accurately
//  cosl(r + c) = 1 - 2^(-67)   accurately
//
//  Case 4:
//  -------
//
//  sinl(r + c) = r + c - r^3/6 + r^5/120       accurately
//  cosl(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
//
//  sinl(r + c)  =  sinl(r)  +  correction,  or
//  cosl(r + c)  =  cosl(r)  +  correction.
//
//  Specifically,
//
//      sinl(r + c) = sinl(r) + c sin'(r) + O(c^2)
//                 = sinl(r) + c cos (r) + O(c^2)
//                 = sinl(r) + c(1 - r^2/2)  accurately.
//  Similarly,
//
//      cosl(r + c) = cosl(r) - c sinl(r) + O(c^2)
//                 = cosl(r) - c(r - r^3/6)  accurately.
//
//  We therefore concentrate on accurately calculating sinl(r) and 
//  cosl(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 sinl(r) and cosl(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
//
//      sinl(Arg) = (-1)^i_0 * sinl(r + c)      if i_1 = 0
//               = (-1)^i_0 * cosl(r + c)       if i_1 = 1
//
//  can be accurately approximated by
//
//  sinl(Arg) = (-1)^i_0 * [sinl(r) + c]        if i_1 = 0
//           = (-1)^i_0 * [cosl(r) - c*r] if i_1 = 1
//
//  because |r| is small and thus the second terms in the correction 
//  are unneccessary.
//
//  Finally, sinl(r) and cosl(r) are approximated by polynomials of 
//  moderate lengths.
//
//  sinl(r) =  r + S_1 r^3 + S_2 r^5 + ... + S_5 r^11
//  cosl(r) =  1 + C_1 r^2 + C_2 r^4 + ... + C_5 r^10
//
//  We can make use of predicates to selectively calculate 
//  sinl(r) or cosl(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,
// 
//  sinl(Arg) = (-1)^i_0 * sinl(r + c)  if i_1 = 0
//           = (-1)^i_0 * cosl(r + c)   if i_1 = 1.
//
//  Because |r| is now larger, we need one extra term in the 
//  correction. sinl(Arg) can be accurately approximated by
//
//  sinl(Arg) = (-1)^i_0 * [sinl(r) + c(1-r^2/2)]      if i_1 = 0
//           = (-1)^i_0 * [cosl(r) - c*r*(1 - r^2/6)]    i_1 = 1.
//
//  Finally, sinl(r) and cosl(r) are approximated by polynomials of 
//  moderate lengths.
//
//      sinl(r) =  r + PP_1_hi r^3 + PP_1_lo r^3 + 
//                    PP_2 r^5 + ... + PP_8 r^17
//
//      cosl(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.
//

#include "libm_support.h" 

#ifdef _LIBC
.rodata
#else
.data
#endif
.align 64 

FSINCOSL_CONSTANTS:
ASM_TYPE_DIRECTIVE(FSINCOSL_CONSTANTS,@object)
data4 0x4B800000, 0xCB800000, 0x00000000,0x00000000 // two**24, -two**24
data4 0x4E44152A, 0xA2F9836E, 0x00003FFE,0x00000000 // Inv_pi_by_2
data4 0xCE81B9F1, 0xC84D32B0, 0x00004016,0x00000000 // P_0 
data4 0x2168C235, 0xC90FDAA2, 0x00003FFF,0x00000000 // P_1 
data4 0xFC8F8CBB, 0xECE675D1, 0x0000BFBD,0x00000000 // P_2 
data4 0xACC19C60, 0xB7ED8FBB, 0x0000BF7C,0x00000000 // P_3 
data4 0x5F000000, 0xDF000000, 0x00000000,0x00000000 // two_to_63, -two_to_63
data4 0x6EC6B45A, 0xA397E504, 0x00003FE7,0x00000000 // Inv_P_0 
data4 0xDBD171A1, 0x8D848E89, 0x0000BFBF,0x00000000 // d_1 
data4 0x18A66F8E, 0xD5394C36, 0x0000BF7C,0x00000000 // d_2 
data4 0x2168C234, 0xC90FDAA2, 0x00003FFE,0x00000000 // pi_by_4 
data4 0x2168C234, 0xC90FDAA2, 0x0000BFFE,0x00000000 // neg_pi_by_4 
data4 0x3E000000, 0xBE000000, 0x00000000,0x00000000 // two**-3, -two**-3
data4 0x2F000000, 0xAF000000, 0x9E000000,0x00000000 // two**-33, -two**-33, -two**-67
data4 0xA21C0BC9, 0xCC8ABEBC, 0x00003FCE,0x00000000 // PP_8 
data4 0x720221DA, 0xD7468A05, 0x0000BFD6,0x00000000 // PP_7 
data4 0x640AD517, 0xB092382F, 0x00003FDE,0x00000000 // PP_6 
data4 0xD1EB75A4, 0xD7322B47, 0x0000BFE5,0x00000000 // PP_5 
data4 0xFFFFFFFE, 0xFFFFFFFF, 0x0000BFFD,0x00000000 // C_1 
data4 0x00000000, 0xAAAA0000, 0x0000BFFC,0x00000000 // PP_1_hi 
data4 0xBAF69EEA, 0xB8EF1D2A, 0x00003FEC,0x00000000 // PP_4 
data4 0x0D03BB69, 0xD00D00D0, 0x0000BFF2,0x00000000 // PP_3 
data4 0x88888962, 0x88888888, 0x00003FF8,0x00000000 // PP_2
data4 0xAAAB0000, 0xAAAAAAAA, 0x0000BFEC,0x00000000 // PP_1_lo 
data4 0xC2B0FE52, 0xD56232EF, 0x00003FD2,0x00000000 // QQ_8
data4 0x2B48DCA6, 0xC9C99ABA, 0x0000BFDA,0x00000000 // QQ_7
data4 0x9C716658, 0x8F76C650, 0x00003FE2,0x00000000 // QQ_6
data4 0xFDA8D0FC, 0x93F27DBA, 0x0000BFE9,0x00000000 // QQ_5
data4 0xAAAAAAAA, 0xAAAAAAAA, 0x0000BFFC,0x00000000 // S_1 
data4 0x00000000, 0x80000000, 0x0000BFFE,0x00000000 // QQ_1 
data4 0x0C6E5041, 0xD00D00D0, 0x00003FEF,0x00000000 // QQ_4 
data4 0x0B607F60, 0xB60B60B6, 0x0000BFF5,0x00000000 // QQ_3 
data4 0xAAAAAA9B, 0xAAAAAAAA, 0x00003FFA,0x00000000 // QQ_2 
data4 0xFFFFFFFE, 0xFFFFFFFF, 0x0000BFFD,0x00000000 // C_1 
data4 0xAAAA719F, 0xAAAAAAAA, 0x00003FFA,0x00000000 // C_2 
data4 0x0356F994, 0xB60B60B6, 0x0000BFF5,0x00000000 // C_3
data4 0xB2385EA9, 0xD00CFFD5, 0x00003FEF,0x00000000 // C_4 
data4 0x292A14CD, 0x93E4BD18, 0x0000BFE9,0x00000000 // C_5
data4 0xAAAAAAAA, 0xAAAAAAAA, 0x0000BFFC,0x00000000 // S_1 
data4 0x888868DB, 0x88888888, 0x00003FF8,0x00000000 // S_2 
data4 0x055EFD4B, 0xD00D00D0, 0x0000BFF2,0x00000000 // S_3 
data4 0x839730B9, 0xB8EF1C5D, 0x00003FEC,0x00000000 // S_4
data4 0xE5B3F492, 0xD71EA3A4, 0x0000BFE5,0x00000000 // S_5
data4 0x38800000, 0xB8800000, 0x00000000            // two**-14, -two**-14
ASM_SIZE_DIRECTIVE(FSINCOSL_CONSTANTS)

FR_Input_X        = f8 
FR_Neg_Two_to_M3  = f32 
FR_Two_to_63      = f32 
FR_Two_to_24      = f33 
FR_Pi_by_4        = f33 
FR_Two_to_M14     = f34 
FR_Two_to_M33     = f35 
FR_Neg_Two_to_24  = f36 
FR_Neg_Pi_by_4    = f36 
FR_Neg_Two_to_M14 = 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_c              = f48 
FR_r              = 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           = f81 
FR_PP_7           = f82 
FR_QQ_7           = f82 
FR_PP_6           = f83 
FR_QQ_6           = f83 
FR_PP_5           = f84 
FR_QQ_5           = f84 
FR_PP_4           = f85 
FR_QQ_4           = f85 
FR_PP_3           = f86 
FR_QQ_3           = f86 
FR_PP_2           = f87 
FR_QQ_2           = f87 
FR_QQ_1           = f88 
FR_N_0_fix        = f89 
FR_Inv_P_0        = f90 
FR_corr           = f91 
FR_poly           = f92 
FR_d_2            = f93 
FR_Two_to_M3      = f94 
FR_Neg_Two_to_63  = f94 
FR_P_0            = f95 
FR_C_lo           = f96 
FR_PP_1           = f97 
FR_PP_1_lo        = f98 
FR_ArgPrime       = f99 

GR_Table_Base  = r32 
GR_Table_Base1 = r33 
GR_i_0         = r34
GR_i_1         = r35
GR_N_Inc       = r36 
GR_Sin_or_Cos  = r37 

// Added for unwind support

GR_SAVE_B0     = r39
GR_SAVE_GP     = r40
GR_SAVE_PFS    = r41


.global sinl#
.global cosl#
#ifdef _LIBC
.global __sinl#
.global __cosl#
#endif

.section .text
.proc sinl#
#ifdef _LIBC
.proc __sinl#
#endif
.align 64 
sinl:
#ifdef _LIBC
__sinl:
#endif
{ .mlx
alloc GR_Table_Base = ar.pfs,0,12,2,0
(p0)   movl GR_Sin_or_Cos = 0x0 ;;
}

{ .mmi
      nop.m 999
(p0)  addl           GR_Table_Base   = @ltoff(FSINCOSL_CONSTANTS#), gp
      nop.i 999
}
;;

{ .mmb
      ld8 GR_Table_Base = [GR_Table_Base]
      nop.m 999
(p0)   br.cond.sptk L(SINCOSL_CONTINUE) ;;
}
;;


.endp sinl#
ASM_SIZE_DIRECTIVE(sinl#)

.section .text
.proc cosl#
cosl:
#ifdef _LIBC
.proc __cosl#
__cosl:
#endif
{ .mlx
alloc GR_Table_Base= ar.pfs,0,12,2,0
(p0)   movl GR_Sin_or_Cos = 0x1 ;;
}
;;

{ .mmi
      nop.m 999
(p0)  addl           GR_Table_Base   = @ltoff(FSINCOSL_CONSTANTS#), gp
      nop.i 999
}
;;

{ .mmb
      ld8 GR_Table_Base = [GR_Table_Base]
      nop.m 999
      nop.b 999
}
;;



//
//     Load Table Address
//

L(SINCOSL_CONTINUE): 
{ .mmi
(p0)   add GR_Table_Base1 = 96, GR_Table_Base
(p0)   ldfs     FR_Two_to_24 = [GR_Table_Base], 4
// GR_Sin_or_Cos denotes 
(p0)   mov   r39 = b0 ;;
}
{ .mmi
       nop.m 0
//
//     Load 2**24, load 2**63.
//
(p0)   ldfs     FR_Neg_Two_to_24 = [GR_Table_Base], 12
       nop.i 0
}
{ .mfi
(p0)   ldfs     FR_Two_to_63 = [GR_Table_Base1], 4
//
//     Check for unnormals - unsupported operands. We do not want
//     to generate denormal exception
//     Check for NatVals, QNaNs, SNaNs, +/-Infs
//     Check for EM unsupporteds
//     Check for Zero 
//
(p0)   fclass.m.unc  p6, p0 =  FR_Input_X, 0x1E3
       nop.i 0
};;
{ .mmf
        nop.m 999
(p0)   ldfs     FR_Neg_Two_to_63 = [GR_Table_Base1], 12
(p0)   fclass.nm.unc p8, p0 =  FR_Input_X, 0x1FF
}
{ .mfb
        nop.m 999
(p0)   fclass.m.unc p10, p0 = FR_Input_X, 0x007
(p6)   br.cond.spnt L(SINCOSL_SPECIAL) ;;
}
{ .mib
        nop.m 999
        nop.i 999
(p8)   br.cond.spnt L(SINCOSL_SPECIAL) ;;
}
{ .mib
        nop.m 999
        nop.i 999
//
//     Branch if +/- NaN, Inf.
//     Load -2**24, load -2**63.
//
(p10)  br.cond.spnt L(SINCOSL_ZERO) ;;
}
{ .mmb
(p0)   ldfe     FR_Inv_pi_by_2 = [GR_Table_Base], 16
(p0)   ldfe     FR_Inv_P_0 = [GR_Table_Base1], 16
        nop.b 999 ;;
}
{ .mmb
(p0)   ldfe             FR_d_1 = [GR_Table_Base1], 16
//
//     Raise possible denormal operand flag with useful fcmp
//     Is x <= -2**63
//     Load Inv_P_0 for pre-reduction
//     Load Inv_pi_by_2
//
(p0)   ldfe             FR_P_0 = [GR_Table_Base], 16
        nop.b 999 ;;
}
{ .mmb
(p0)   ldfe     FR_d_2 = [GR_Table_Base1], 16
//
//     Load P_0
//     Load d_1
//     Is x >= 2**63
//     Is x <= -2**24?
//
(p0)   ldfe     FR_P_1 = [GR_Table_Base], 16
        nop.b 999 ;;
}
//
//     Load P_1
//     Load d_2
//     Is x >= 2**24?
//
{ .mfi
(p0)   ldfe     FR_P_2 = [GR_Table_Base], 16
(p0)   fcmp.le.unc.s1   p7, p8 = FR_Input_X, FR_Neg_Two_to_24
        nop.i 999 ;;
}
{ .mbb
(p0)   ldfe     FR_P_3 = [GR_Table_Base], 16
        nop.b 999
        nop.b 999 ;;
}
{ .mfi
        nop.m 999
(p8)   fcmp.ge.s1 p7, p0 = FR_Input_X, FR_Two_to_24
        nop.i 999
}
{ .mfi
(p0)   ldfe     FR_Pi_by_4 = [GR_Table_Base1], 16
//
//     Branch if +/- zero.
//     Decide about the paths to take:
//     If -2**24 < FR_Input_X < 2**24 - CASE 1 OR 2 
//     OTHERWISE - CASE 3 OR 4 
//
(p0)   fcmp.le.unc.s0   p10, p11 = FR_Input_X, FR_Neg_Two_to_63
        nop.i 999 ;;
}
{ .mmi
(p0)   ldfe     FR_Neg_Pi_by_4 = [GR_Table_Base1], 16 ;;
(p0)   ldfs     FR_Two_to_M3 = [GR_Table_Base1], 4
        nop.i 999
}
{ .mfi
        nop.m 999
(p11)  fcmp.ge.s1       p10, p0 = FR_Input_X, FR_Two_to_63
        nop.i 999 ;;
}
{ .mib
(p0)   ldfs     FR_Neg_Two_to_M3 = [GR_Table_Base1], 12
        nop.i 999
//
//     Load P_2
//     Load P_3
//     Load pi_by_4
//     Load neg_pi_by_4
//     Load 2**(-3)
//     Load -2**(-3).
//
(p10)  br.cond.spnt L(SINCOSL_ARG_TOO_LARGE) ;;
}
{ .mib
        nop.m 999
        nop.i 999
//
//     Branch out if x >= 2**63. Use Payne-Hanek Reduction
//
(p7)   br.cond.spnt L(SINCOSL_LARGER_ARG) ;;
}
{ .mfi
        nop.m 999
// 
//     Branch if Arg <= -2**24 or Arg >= 2**24 and use pre-reduction.
//
(p0)   fma.s1   FR_N_float = FR_Input_X, FR_Inv_pi_by_2, f0
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
(p0)   fcmp.lt.unc.s1   p6, p7 = FR_Input_X, FR_Pi_by_4
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
// 
//     Select the case when |Arg| < pi/4 
//     Else Select the case when |Arg| >= pi/4 
//
(p0)   fcvt.fx.s1 FR_N_fix = FR_N_float
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
//
//     N  = Arg * 2/pi
//     Check if Arg < pi/4
//
(p6)   fcmp.gt.s1 p6, p7 = FR_Input_X, FR_Neg_Pi_by_4
        nop.i 999 ;;
}
//
//     Case 2: Convert integer N_fix back to normalized floating-point value.
//     Case 1: p8 is only affected  when p6 is set
//
{ .mfi
(p7)   ldfs FR_Two_to_M33 = [GR_Table_Base1], 4
//
//     Grab the integer part of N and call it N_fix
//
(p6)   fmerge.se FR_r = FR_Input_X, FR_Input_X
//     If |x| < pi/4, r = x and c = 0 
//     lf |x| < pi/4, is x < 2**(-3).
//     r = Arg 
//     c = 0
(p6)   mov GR_N_Inc = GR_Sin_or_Cos ;;
}
{ .mmf
        nop.m 999
(p7)   ldfs FR_Neg_Two_to_M33 = [GR_Table_Base1], 4
(p6)   fmerge.se FR_c = f0, f0
}
{ .mfi
        nop.m 999
(p6)   fcmp.lt.unc.s1   p8, p9 = FR_Input_X, FR_Two_to_M3
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
//
//     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
//
(p7)   fcvt.xf FR_N_float = FR_N_fix
        nop.i 999 ;;
}
{ .mmf
        nop.m 999
(p7)   getf.sig GR_N_Inc = FR_N_fix
(p8)   fcmp.gt.s1 p8, p0 = FR_Input_X, FR_Neg_Two_to_M3 ;;
}
{ .mib
        nop.m 999
        nop.i 999
//
//     Load 2**(-33), -2**(-33)
//
(p8)   br.cond.spnt L(SINCOSL_SMALL_R) ;;
}
{ .mib
        nop.m 999
        nop.i 999
(p6)   br.cond.sptk L(SINCOSL_NORMAL_R) ;;
}
//
//     if |x| < pi/4, branch based on |x| < 2**(-3) or otherwise.
//
//
//     In this branch, |x| >= pi/4.
// 
{ .mfi
(p0)   ldfs FR_Neg_Two_to_M67 = [GR_Table_Base1], 8
//
//     Load -2**(-67)
// 
(p0)   fnma.s1  FR_s = FR_N_float, FR_P_1, FR_Input_X
//
//     w = N * P_2
//     s = -N * P_1  + Arg
//
(p0)   add GR_N_Inc = GR_N_Inc, GR_Sin_or_Cos
}
{ .mfi
        nop.m 999
(p0)   fma.s1   FR_w = FR_N_float, FR_P_2, f0
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
// 
//     Adjust N_fix by N_inc to determine whether sine or
//     cosine is being calculated
//
(p0)   fcmp.lt.unc.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
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
//     Remember x >= pi/4.
//     Is s <= -2**(-33) or s >= 2**(-33) (p6)
//     or -2**(-33) < s < 2**(-33) (p7)
(p6)   fms.s1 FR_r = FR_s, f1, FR_w
        nop.i 999
}
{ .mfi
        nop.m 999
(p7)   fma.s1 FR_w = FR_N_float, FR_P_3, f0
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
(p7)   fma.s1 FR_U_1 = FR_N_float, FR_P_2, FR_w
        nop.i 999
}
{ .mfi
        nop.m 999
(p6)   fms.s1 FR_c = FR_s, f1, FR_r
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
// 
//     For big s: r = s - w: No futher reduction is necessary 
//     For small s: w = N * P_3 (change sign) More reduction
//
(p6)   fcmp.lt.unc.s1 p8, p9 = FR_r, FR_Two_to_M3
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
(p8)   fcmp.gt.s1 p8, p9 = FR_r, FR_Neg_Two_to_M3
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
(p7)   fms.s1 FR_r = FR_s, f1, FR_U_1
        nop.i 999
}
{ .mfb
        nop.m 999
//
//     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.
//
(p6)   fms.s1 FR_c = FR_c, f1, FR_w
// 
//     For big s: c = c + w (w has not been negated.)
//     For small s: r = S - U_1
//
(p8)   br.cond.spnt     L(SINCOSL_SMALL_R) ;;
}
{ .mib
        nop.m 999
        nop.i 999
(p9)   br.cond.sptk     L(SINCOSL_NORMAL_R) ;;
}
{ .mfi
(p7)   add GR_Table_Base1 = 224, GR_Table_Base1
//
//     Branch to SINCOSL_SMALL_R or SINCOSL_NORMAL_R
//
(p7)   fms.s1 FR_U_2 = FR_N_float, FR_P_2, FR_U_1
// 
//     c = S - U_1
//     r = S_1 * r
//
//
(p7)   extr.u   GR_i_1 = GR_N_Inc, 0, 1 ;;
}
{ .mmi
        nop.m 999
//
//     Get [i_0,i_1] - two lsb of N_fix_gr.
//     Do dummy fmpy so inexact is always set.
//
(p7)   cmp.eq.unc p9, p10 = 0x0, GR_i_1
(p7)   extr.u   GR_i_0 = GR_N_Inc, 1, 1 ;;
}
// 
//     For small s: U_2 = N * P_2 - U_1
//     S_1 stored constant - grab the one stored with the
//     coefficients.
// 
{ .mfi
(p7)   ldfe FR_S_1 = [GR_Table_Base1], 16
//
//     Check if i_1 and i_0  != 0
//
(p10)  fma.s1   FR_poly = f0, f1, FR_Neg_Two_to_M67
(p7)   cmp.eq.unc p11, p12 = 0x0, GR_i_0 ;;
}
{ .mfi
        nop.m 999
(p7)   fms.s1   FR_s = FR_s, f1, FR_r
        nop.i 999
}
{ .mfi
        nop.m 999
// 
//     S = S - r
//     U_2 = U_2 + w
//     load S_1
//
(p7)   fma.s1   FR_rsq = FR_r, FR_r, f0
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
(p7)   fma.s1   FR_U_2 = FR_U_2, f1, FR_w
        nop.i 999
}
{ .mfi
        nop.m 999
(p7)   fmerge.se FR_Input_X = FR_r, FR_r
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
(p10)  fma.s1 FR_Input_X = f0, f1, f1
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
// 
//     FR_rsq = r * r
//     Save r as the result.
//
(p7)   fms.s1   FR_c = FR_s, f1, FR_U_1
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
// 
//     if ( i_1 ==0) poly = c + S_1*r*r*r
//     else Result = 1
//
(p12)  fnma.s1 FR_Input_X = FR_Input_X, f1, f0
        nop.i 999
}
{ .mfi
        nop.m 999
(p7)   fma.s1   FR_r = FR_S_1, FR_r, f0
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
(p7)   fma.s0   FR_S_1 = FR_S_1, FR_S_1, f0
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
//
//     If i_1 != 0, poly = 2**(-67)
//
(p7)   fms.s1 FR_c = FR_c, f1, FR_U_2
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
// 
//     c = c - U_2
// 
(p9)   fma.s1 FR_poly = FR_r, FR_rsq, FR_c
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
//
//     i_0 != 0, so Result = -Result
//
(p11)  fma.s0 FR_Input_X = FR_Input_X, f1, FR_poly
        nop.i 999 ;;
}
{ .mfb
        nop.m 999
(p12)  fms.s0 FR_Input_X = FR_Input_X, f1, FR_poly
//
//     if (i_0 == 0),  Result = Result + poly
//     else            Result = Result - poly
//
(p0)    br.ret.sptk   b0 ;;
}
L(SINCOSL_LARGER_ARG): 
{ .mfi
        nop.m 999
(p0)   fma.s1 FR_N_0 = FR_Input_X, FR_Inv_P_0, f0
        nop.i 999
}
;;

//     This path for argument > 2*24 
//     Adjust table_ptr1 to beginning of table.
//

{ .mmi
      nop.m 999
(p0)  addl           GR_Table_Base   = @ltoff(FSINCOSL_CONSTANTS#), gp
      nop.i 999
}
;;

{ .mmi
      ld8 GR_Table_Base = [GR_Table_Base]
      nop.m 999
      nop.i 999
}
;;


// 
//     Point to  2*-14 
//     N_0 = Arg * Inv_P_0
//
{ .mmi
(p0)   add GR_Table_Base = 688, GR_Table_Base ;;
(p0)   ldfs FR_Two_to_M14 = [GR_Table_Base], 4
        nop.i 999 ;;
}
{ .mfi
(p0)   ldfs FR_Neg_Two_to_M14 = [GR_Table_Base], 0
        nop.f 999
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
//
//     Load values 2**(-14) and -2**(-14)
//
(p0)   fcvt.fx.s1 FR_N_0_fix = FR_N_0
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
//
//     N_0_fix  = integer part of N_0
//
(p0)   fcvt.xf FR_N_0 = FR_N_0_fix 
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
//
//     Make N_0 the integer part
//
(p0)   fnma.s1 FR_ArgPrime = FR_N_0, FR_P_0, FR_Input_X
        nop.i 999
}
{ .mfi
        nop.m 999
(p0)   fma.s1 FR_w = FR_N_0, FR_d_1, f0
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
//
//     Arg' = -N_0 * P_0 + Arg
//     w  = N_0 * d_1
//
(p0)   fma.s1 FR_N_float = FR_ArgPrime, FR_Inv_pi_by_2, f0
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
//
//     N = A' * 2/pi    
//
(p0)   fcvt.fx.s1 FR_N_fix = FR_N_float
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
//
//     N_fix is the integer part        
//
(p0)   fcvt.xf FR_N_float = FR_N_fix 
        nop.i 999 ;;
}
{ .mfi
(p0)   getf.sig GR_N_Inc = FR_N_fix
        nop.f 999
        nop.i 999 ;;
}
{ .mii
        nop.m 999
        nop.i 999 ;;
(p0)   add GR_N_Inc = GR_N_Inc, GR_Sin_or_Cos ;;
}
{ .mfi
        nop.m 999
//
//     N is the integer part of the reduced-reduced argument.
//     Put the integer in a GP register
//
(p0)   fnma.s1 FR_s = FR_N_float, FR_P_1, FR_ArgPrime
        nop.i 999
}
{ .mfi
        nop.m 999
(p0)   fnma.s1 FR_w = FR_N_float, FR_P_2, FR_w
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
//
//     s = -N*P_1 + Arg'
//     w = -N*P_2 + w
//     N_fix_gr = N_fix_gr + N_inc
//
(p0)   fcmp.lt.unc.s1 p9, p8 = FR_s, FR_Two_to_M14
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
(p9)   fcmp.gt.s1 p9, p8 = FR_s, FR_Neg_Two_to_M14
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
//
//     For |s|  > 2**(-14) r = S + w (r complete)
//     Else       U_hi = N_0 * d_1
//
(p9)   fma.s1 FR_V_hi = FR_N_float, FR_P_2, f0
        nop.i 999
}
{ .mfi
        nop.m 999
(p9)   fma.s1 FR_U_hi = FR_N_0, FR_d_1, f0
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
//
//     Either S <= -2**(-14) or S >= 2**(-14)
//     or -2**(-14) < s < 2**(-14)
//
(p8)   fma.s1 FR_r = FR_s, f1, FR_w
        nop.i 999
}
{ .mfi
        nop.m 999
(p9)   fma.s1 FR_w = FR_N_float, FR_P_3, f0
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
//
//     We need abs of both U_hi and V_hi - don't
//     worry about switched sign of V_hi.
//
(p9)   fms.s1 FR_A = FR_U_hi, f1, FR_V_hi
        nop.i 999
}
{ .mfi
        nop.m 999
//
//     Big s: finish up c = (S - r) + w (c complete)    
//     Case 4: A =  U_hi + V_hi
//     Note: Worry about switched sign of V_hi, so subtract instead of add.
//
(p9)   fnma.s1 FR_V_lo = FR_N_float, FR_P_2, FR_V_hi
        nop.i 999 ;;
}
{ .mmf
        nop.m 999
        nop.m 999
(p9)   fms.s1 FR_U_lo = FR_N_0, FR_d_1, FR_U_hi
}
{ .mfi
        nop.m 999
(p9)   fmerge.s FR_V_hiabs = f0, FR_V_hi
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
//     For big s: c = S - r
//     For small s do more work: U_lo = N_0 * d_1 - U_hi
//
(p9)   fmerge.s FR_U_hiabs = f0, FR_U_hi
        nop.i 999
}
{ .mfi
        nop.m 999
//
//     For big s: Is |r| < 2**(-3)      
//     For big s: if p12 set, prepare to branch to Small_R.
//     For big s: If p13 set, prepare to branch to Normal_R.
//
(p8)   fms.s1 FR_c = FR_s, f1, FR_r 
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
//
//     For small S: V_hi = N * P_2
//                  w = N * P_3
//     Note the product does not include the (-) as in the writeup
//     so (-) missing for V_hi and w.
//
(p8)   fcmp.lt.unc.s1 p12, p13 = FR_r, FR_Two_to_M3
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
(p12)  fcmp.gt.s1 p12, p13 = FR_r, FR_Neg_Two_to_M3
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
(p8)   fma.s1 FR_c = FR_c, f1, FR_w
        nop.i 999
}
{ .mfb
        nop.m 999
(p9)   fms.s1 FR_w = FR_N_0, FR_d_2, FR_w
(p12)  br.cond.spnt L(SINCOSL_SMALL_R) ;;
}
{ .mib
        nop.m 999
        nop.i 999
(p13)  br.cond.sptk L(SINCOSL_NORMAL_R) ;;
}
{ .mfi
        nop.m 999
// 
//     Big s: Vector off when |r| < 2**(-3).  Recall that p8 will be true. 
//     The remaining stuff is for Case 4.
//     Small s: V_lo = N * P_2 + U_hi (U_hi is in place of V_hi in writeup)
//     Note: the (-) is still missing for V_lo.
//     Small s: w = w + N_0 * d_2
//     Note: the (-) is now incorporated in w.
//
(p9)   fcmp.ge.unc.s1 p10, p11 = FR_U_hiabs, FR_V_hiabs
(p0)   extr.u   GR_i_1 = GR_N_Inc, 0, 1
}
{ .mfi
        nop.m 999
//
//     C_hi = S + A
//
(p9)   fma.s1 FR_t = FR_U_lo, f1, FR_V_lo
(p0)   extr.u   GR_i_0 = GR_N_Inc, 1, 1 ;;
}
{ .mfi
        nop.m 999
//
//     t = U_lo + V_lo 
//
//
(p10)  fms.s1 FR_a = FR_U_hi, f1, FR_A
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
(p11)  fma.s1 FR_a = FR_V_hi, f1, FR_A
        nop.i 999
}
;;

{ .mmi
      nop.m 999
(p0)  addl           GR_Table_Base   = @ltoff(FSINCOSL_CONSTANTS#), gp
      nop.i 999
}
;;

{ .mmi
      ld8 GR_Table_Base = [GR_Table_Base]
      nop.m 999
      nop.i 999
}
;;


{ .mfi
(p0)   add GR_Table_Base = 528, GR_Table_Base
//
//     Is U_hiabs >= V_hiabs?
//
(p9)   fma.s1 FR_C_hi = FR_s, f1, FR_A
        nop.i 999 ;;
}
{ .mmi
(p0)   ldfe FR_C_1 = [GR_Table_Base], 16 ;;
(p0)   ldfe FR_C_2 = [GR_Table_Base], 64
        nop.i 999 ;;
}
//
//     c = c + C_lo  finished.
//     Load  C_2
//
{ .mfi
(p0)   ldfe     FR_S_1 = [GR_Table_Base], 16
//
//     C_lo = S - C_hi 
//
(p0)   fma.s1 FR_t = FR_t, f1, FR_w
        nop.i 999 ;;
}
//
//     r and c have been computed.
//     Make sure ftz mode is set - should be automatic when using wre
//     |r| < 2**(-3)
//     Get [i_0,i_1] - two lsb of N_fix.
//     Load S_1
//
{ .mfi
(p0)   ldfe FR_S_2 = [GR_Table_Base], 64
//
//     t = t + w        
//
(p10)  fms.s1 FR_a = FR_a, f1, FR_V_hi
(p0)   cmp.eq.unc p9, p10 = 0x0, GR_i_0 ;;
}
{ .mfi
        nop.m 999
//
//     For larger u than v: a = U_hi - A
//     Else a = V_hi - A (do an add to account for missing (-) on V_hi
//
(p0)   fms.s1 FR_C_lo = FR_s, f1, FR_C_hi
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
(p11)  fms.s1 FR_a = FR_U_hi, f1, FR_a
(p0)   cmp.eq.unc p11, p12 = 0x0, GR_i_1 ;;
}
{ .mfi
        nop.m 999
//
//     If u > v: a = (U_hi - A)  + V_hi
//     Else      a = (V_hi - A)  + U_hi
//     In each case account for negative missing from V_hi.
//
(p0)   fma.s1 FR_C_lo = FR_C_lo, f1, FR_A
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
//
//     C_lo = (S - C_hi) + A    
//
(p0)   fma.s1 FR_t = FR_t, f1, FR_a
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
//
//     t = t + a 
//
(p0)   fma.s1 FR_C_lo = FR_C_lo, f1, FR_t
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
//
//     C_lo = C_lo + t
//     Adjust Table_Base to beginning of table
//
(p0)   fma.s1 FR_r = FR_C_hi, f1, FR_C_lo
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
//
//     Load S_2
//
(p0)   fma.s1 FR_rsq = FR_r, FR_r, f0
        nop.i 999
}
{ .mfi
        nop.m 999
//
//     Table_Base points to C_1
//     r = C_hi + C_lo
//
(p0)   fms.s1 FR_c = FR_C_hi, f1, FR_r
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
//
//     if i_1 ==0: poly = S_2 * FR_rsq + S_1
//     else        poly = C_2 * FR_rsq + C_1
//
(p11)  fma.s1 FR_Input_X = f0, f1, FR_r
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
(p12)  fma.s1 FR_Input_X = f0, f1, f1
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
//
//     Compute r_cube = FR_rsq * r      
//
(p11)  fma.s1 FR_poly = FR_rsq, FR_S_2, FR_S_1
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
(p12)  fma.s1 FR_poly = FR_rsq, FR_C_2, FR_C_1
        nop.i 999
}
{ .mfi
        nop.m 999
//
//     Compute FR_rsq = r * r
//     Is i_1 == 0 ?
//
(p0)   fma.s1 FR_r_cubed = FR_rsq, FR_r, f0
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
//
//     c = C_hi - r
//     Load  C_1
//
(p0)   fma.s1 FR_c = FR_c, f1, FR_C_lo
        nop.i 999
}
{ .mfi
        nop.m 999
//
//     if i_1 ==0: poly = r_cube * poly + c
//     else        poly = FR_rsq * poly
//
(p10)  fms.s1 FR_Input_X = f0, f1, FR_Input_X
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
//
//     if i_1 ==0: Result = r
//     else        Result = 1.0
//
(p11)  fma.s1 FR_poly = FR_r_cubed, FR_poly, FR_c
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
(p12)  fma.s1 FR_poly = FR_rsq, FR_poly, f0
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
//
//     if i_0 !=0: Result = -Result 
//
(p9)   fma.s0 FR_Input_X = FR_Input_X, f1, FR_poly
        nop.i 999 ;;
}
{ .mfb
        nop.m 999
(p10)  fms.s0 FR_Input_X = FR_Input_X, f1, FR_poly
//
//     if i_0 == 0: Result = Result + poly
//     else         Result = Result - poly
//
(p0)    br.ret.sptk   b0 ;;
}
L(SINCOSL_SMALL_R): 
{ .mii
        nop.m 999
(p0)    extr.u  GR_i_1 = GR_N_Inc, 0, 1 ;;
//
//
//      Compare both i_1 and i_0 with 0.
//      if i_1 == 0, set p9.
//      if i_0 == 0, set p11.
//
(p0)    cmp.eq.unc p9, p10 = 0x0, GR_i_1 ;;
}
{ .mfi
        nop.m 999
(p0)    fma.s1 FR_rsq = FR_r, FR_r, f0
(p0)    extr.u  GR_i_0 = GR_N_Inc, 1, 1 ;;
}
{ .mfi
        nop.m 999
//
//      Z = Z * FR_rsq 
//
(p10)   fnma.s1 FR_c = FR_c, FR_r, f0
(p0)    cmp.eq.unc p11, p12 = 0x0, GR_i_0
}
;;

// ******************************************************************
// ******************************************************************
// ******************************************************************
//      r and c have been computed.
//      We know whether this is the sine or cosine routine.
//      Make sure ftz mode is set - should be automatic when using wre
//      |r| < 2**(-3)
//
//      Set table_ptr1 to beginning of constant table.
//      Get [i_0,i_1] - two lsb of N_fix_gr.
//

{ .mmi
      nop.m 999
(p0)  addl           GR_Table_Base   = @ltoff(FSINCOSL_CONSTANTS#), gp
      nop.i 999
}
;;

{ .mmi
      ld8 GR_Table_Base = [GR_Table_Base]
      nop.m 999
      nop.i 999
}
;;


// 
//      Set table_ptr1 to point to S_5.
//      Set table_ptr1 to point to C_5.
//      Compute FR_rsq = r * r
//
{ .mfi
(p9)    add GR_Table_Base = 672, GR_Table_Base
(p10)   fmerge.s FR_r = f1, f1
(p10)   add GR_Table_Base = 592, GR_Table_Base ;;
}
// 
//      Set table_ptr1 to point to S_5.
//      Set table_ptr1 to point to C_5.
//
{ .mmi
(p9)    ldfe FR_S_5 = [GR_Table_Base], -16 ;;
//
//      if (i_1 == 0) load S_5
//      if (i_1 != 0) load C_5
//
(p9)    ldfe FR_S_4 = [GR_Table_Base], -16
        nop.i 999 ;;
}
{ .mmf
(p10)   ldfe FR_C_5 = [GR_Table_Base], -16
// 
//      Z = FR_rsq * FR_rsq
//
(p9)    ldfe FR_S_3 = [GR_Table_Base], -16
//
//      Compute FR_rsq = r * r
//      if (i_1 == 0) load S_4
//      if (i_1 != 0) load C_4
//
(p0)    fma.s1 FR_Z = FR_rsq, FR_rsq, f0 ;;
}
//
//      if (i_1 == 0) load S_3
//      if (i_1 != 0) load C_3
//
{ .mmi
(p9)    ldfe FR_S_2 = [GR_Table_Base], -16 ;;
//
//      if (i_1 == 0) load S_2
//      if (i_1 != 0) load C_2
//
(p9)    ldfe FR_S_1 = [GR_Table_Base], -16
        nop.i 999
}
{ .mmi
(p10)   ldfe FR_C_4 = [GR_Table_Base], -16 ;;
(p10)   ldfe FR_C_3 = [GR_Table_Base], -16
        nop.i 999 ;;
}
{ .mmi
(p10)   ldfe FR_C_2 = [GR_Table_Base], -16 ;;
(p10)   ldfe FR_C_1 = [GR_Table_Base], -16
        nop.i 999
}
{ .mfi
        nop.m 999
//
//      if (i_1 != 0):
//      poly_lo = FR_rsq * C_5 + C_4
//      poly_hi = FR_rsq * C_2 + C_1
//
(p9)    fma.s1 FR_Z = FR_Z, FR_r, f0
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
//
//      if (i_1 == 0) load S_1
//      if (i_1 != 0) load C_1
//
(p9)    fma.s1 FR_poly_lo = FR_rsq, FR_S_5, FR_S_4
        nop.i 999
}
{ .mfi
        nop.m 999
//
//      c = -c * r
//      dummy fmpy's to flag inexact.
//
(p9)    fma.s0 FR_S_4 = FR_S_4, FR_S_4, f0
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
//
//      poly_lo = FR_rsq * poly_lo + C_3
//      poly_hi = FR_rsq * poly_hi
//
(p0)    fma.s1  FR_Z = FR_Z, FR_rsq, f0
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
(p9)    fma.s1 FR_poly_hi = FR_rsq, FR_S_2, FR_S_1
        nop.i 999
}
{ .mfi
        nop.m 999
//
//      if (i_1 == 0):
//      poly_lo = FR_rsq * S_5 + S_4
//      poly_hi = FR_rsq * S_2 + S_1
//
(p10)   fma.s1 FR_poly_lo = FR_rsq, FR_C_5, FR_C_4
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
//
//      if (i_1 == 0):
//      Z = Z * r  for only one of the small r cases - not there
//      in original implementation notes.
// 
(p9)    fma.s1 FR_poly_lo = FR_rsq, FR_poly_lo, FR_S_3
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
(p10)   fma.s1 FR_poly_hi = FR_rsq, FR_C_2, FR_C_1
        nop.i 999
}
{ .mfi
        nop.m 999
(p10)   fma.s0 FR_C_1 = FR_C_1, FR_C_1, f0
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
(p9)    fma.s1 FR_poly_hi = FR_poly_hi, FR_rsq, f0
        nop.i 999
}
{ .mfi
        nop.m 999
//
//      poly_lo = FR_rsq * poly_lo + S_3
//      poly_hi = FR_rsq * poly_hi
//
(p10)   fma.s1 FR_poly_lo = FR_rsq, FR_poly_lo, FR_C_3
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
(p10)   fma.s1 FR_poly_hi = FR_poly_hi, FR_rsq, f0
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
//
//      if (i_1 == 0): dummy fmpy's to flag inexact
//      r = 1
//
(p9)    fma.s1 FR_poly_hi = FR_r, FR_poly_hi, f0
        nop.i 999
}
{ .mfi
        nop.m 999
//
//      poly_hi = r * poly_hi 
//
(p0)    fma.s1  FR_poly = FR_Z, FR_poly_lo, FR_c
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
(p12)   fms.s1  FR_r = f0, f1, FR_r
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
//
//      poly_hi = Z * poly_lo + c       
//      if i_0 == 1: r = -r     
//
(p0)    fma.s1  FR_poly = FR_poly, f1, FR_poly_hi
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
(p12)   fms.s0 FR_Input_X = FR_r, f1, FR_poly
        nop.i 999
}
{ .mfb
        nop.m 999
//
//      poly = poly + poly_hi   
//
(p11)   fma.s0 FR_Input_X = FR_r, f1, FR_poly
//
//      if (i_0 == 0) Result = r + poly
//      if (i_0 != 0) Result = r - poly
//
(p0)    br.ret.sptk   b0 ;;
}
L(SINCOSL_NORMAL_R): 
{ .mii
        nop.m 999
(p0)    extr.u  GR_i_1 = GR_N_Inc, 0, 1 ;;
//
//      Set table_ptr1 and table_ptr2 to base address of
//      constant table.
(p0)    cmp.eq.unc p9, p10 = 0x0, GR_i_1 ;;
}
{ .mfi
        nop.m 999
(p0)    fma.s1  FR_rsq = FR_r, FR_r, f0
(p0)    extr.u  GR_i_0 = GR_N_Inc, 1, 1 ;;
}
{ .mfi
        nop.m 999
(p0)    frcpa.s1 FR_r_hi, p6 = f1, FR_r
(p0)    cmp.eq.unc p11, p12 = 0x0, GR_i_0
}
;;

// ******************************************************************
// ******************************************************************
// ******************************************************************
//
//      r and c have been computed.
//      We known whether this is the sine or cosine routine.
//      Make sure ftz mode is set - should be automatic when using wre
//      Get [i_0,i_1] - two lsb of N_fix_gr alone.
//

{ .mmi
      nop.m 999
(p0)  addl           GR_Table_Base   = @ltoff(FSINCOSL_CONSTANTS#), gp
      nop.i 999
}
;;

{ .mmi
      ld8 GR_Table_Base = [GR_Table_Base]
      nop.m 999
      nop.i 999
}
;;


{ .mfi
(p10)   add GR_Table_Base = 384, GR_Table_Base
(p12)   fms.s1 FR_Input_X = f0, f1, f1
(p9)    add GR_Table_Base = 224, GR_Table_Base ;;
}
{ .mfi
(p10)   ldfe FR_QQ_8 = [GR_Table_Base], 16
//
//      if (i_1==0) poly = poly * FR_rsq + PP_1_lo
//      else        poly = FR_rsq * poly
//
(p11)   fma.s1 FR_Input_X = f0, f1, f1
        nop.i 999 ;;
}
{ .mmb
(p10)   ldfe FR_QQ_7 = [GR_Table_Base], 16
//
//      Adjust table pointers based on i_0 
//      Compute rsq = r * r
//
(p9)    ldfe FR_PP_8 = [GR_Table_Base], 16
        nop.b 999 ;;
}
{ .mfi
        nop.m 999
(p0)    fma.s1 FR_r_cubed = FR_r, FR_rsq, f0
        nop.i 999 ;;
}
{ .mmf
(p9)    ldfe FR_PP_7 = [GR_Table_Base], 16
(p10)   ldfe FR_QQ_6 = [GR_Table_Base], 16
//
//      Load PP_8 and QQ_8; PP_7 and QQ_7
//
(p0)    frcpa.s1 FR_r_hi, p6 = f1, FR_r_hi ;;
}
//
//      if (i_1==0) poly =   PP_7 + FR_rsq * PP_8.
//      else        poly =   QQ_7 + FR_rsq * QQ_8.
//
{ .mmb
(p9)    ldfe FR_PP_6 = [GR_Table_Base], 16
(p10)   ldfe FR_QQ_5 = [GR_Table_Base], 16
        nop.b 999 ;;
}
{ .mmb
(p9)    ldfe FR_PP_5 = [GR_Table_Base], 16
(p10)   ldfe FR_S_1 = [GR_Table_Base], 16
        nop.b 999 ;;
}
{ .mmb
(p10)   ldfe FR_QQ_1 = [GR_Table_Base], 16
(p9)    ldfe FR_C_1 = [GR_Table_Base], 16
        nop.b 999 ;;
}
{ .mmb
(p10)   ldfe FR_QQ_4 = [GR_Table_Base], 16
(p9)    ldfe FR_PP_1 = [GR_Table_Base], 16
        nop.b 999 ;;
}
{ .mmb
(p10)   ldfe FR_QQ_3 = [GR_Table_Base], 16
//
//      if (i_1=0) corr = corr + c*c
//      else       corr = corr * c 
//
(p9)    ldfe FR_PP_4 = [GR_Table_Base], 16
        nop.b 999 ;;
}
{ .mfi
        nop.m 999
(p10)   fma.s1 FR_poly = FR_rsq, FR_QQ_8, FR_QQ_7
        nop.i 999 ;;
}
//
//      if (i_1=0) poly = rsq * poly + PP_5 
//      else       poly = rsq * poly + QQ_5 
//      Load PP_4 or QQ_4
//
{ .mmi
(p9)    ldfe FR_PP_3 = [GR_Table_Base], 16 ;;
(p10)   ldfe FR_QQ_2 = [GR_Table_Base], 16
        nop.i 999
}
{ .mfi
        nop.m 999
//
//      r_hi =   frcpa(frcpa(r)).
//      r_cube = r * FR_rsq.
//
(p9)    fma.s1 FR_poly = FR_rsq, FR_PP_8, FR_PP_7
        nop.i 999 ;;
}
//
//      Do dummy multiplies so inexact is always set. 
//
{ .mfi
(p9)    ldfe FR_PP_2 = [GR_Table_Base], 16
//
//      r_lo = r - r_hi 
//
(p9)    fma.s1 FR_U_lo = FR_r_hi, FR_r_hi, f0
        nop.i 999 ;;
}
{ .mbb
(p9)    ldfe FR_PP_1_lo = [GR_Table_Base], 16
        nop.b 999
        nop.b 999 ;;
}
{ .mfi
        nop.m 999
(p10)   fma.s1 FR_corr = FR_S_1, FR_r_cubed, FR_r
        nop.i 999
}
{ .mfi
        nop.m 999
(p10)   fma.s1 FR_poly = FR_rsq, FR_poly, FR_QQ_6
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
//
//      if (i_1=0) U_lo = r_hi * r_hi
//      else       U_lo = r_hi + r
//
(p9)    fma.s1 FR_corr = FR_C_1, FR_rsq, f0
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
//
//      if (i_1=0) corr = C_1 * rsq
//      else       corr = S_1 * r_cubed + r
//
(p9)    fma.s1 FR_poly = FR_rsq, FR_poly, FR_PP_6
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
(p10)   fma.s1 FR_U_lo = FR_r_hi, f1, FR_r
        nop.i 999
}
{ .mfi
        nop.m 999
//
//      if (i_1=0) U_hi = r_hi + U_hi 
//      else       U_hi = QQ_1 * U_hi + 1
//
(p9)    fma.s1 FR_U_lo = FR_r, FR_r_hi, FR_U_lo
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
//
//      U_hi = r_hi * r_hi      
//
(p0)    fms.s1 FR_r_lo = FR_r, f1, FR_r_hi
        nop.i 999
}
{ .mfi
        nop.m 999
//
//      Load PP_1, PP_6, PP_5, and C_1
//      Load QQ_1, QQ_6, QQ_5, and S_1
//
(p0)    fma.s1 FR_U_hi = FR_r_hi, FR_r_hi, f0
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
(p10)   fma.s1 FR_poly = FR_rsq, FR_poly, FR_QQ_5
        nop.i 999
}
{ .mfi
        nop.m 999
(p10)   fnma.s1 FR_corr = FR_corr, FR_c, f0
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
//
//      if (i_1=0) U_lo = r * r_hi + U_lo 
//      else       U_lo = r_lo * U_lo
//
(p9)    fma.s1 FR_corr = FR_corr, FR_c, FR_c
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
(p9)    fma.s1 FR_poly = FR_rsq, FR_poly, FR_PP_5
        nop.i 999
}
{ .mfi
        nop.m 999
//
//      if (i_1 =0) U_hi = r + U_hi
//      if (i_1 =0) U_lo = r_lo * U_lo 
//      
//
(p9)    fma.s0 FR_PP_5 = FR_PP_5, FR_PP_4, f0
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
(p9)    fma.s1 FR_U_lo = FR_r, FR_r, FR_U_lo
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
(p10)   fma.s1 FR_U_lo = FR_r_lo, FR_U_lo, f0
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
//
//      if (i_1=0) poly = poly * rsq + PP_6
//      else       poly = poly * rsq + QQ_6 
//
(p9)    fma.s1 FR_U_hi = FR_r_hi, FR_U_hi, f0
        nop.i 999
}
{ .mfi
        nop.m 999
(p10)   fma.s1 FR_poly = FR_rsq, FR_poly, FR_QQ_4
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
(p10)   fma.s1 FR_U_hi = FR_QQ_1, FR_U_hi, f1
        nop.i 999
}
{ .mfi
        nop.m 999
(p10)   fma.s0 FR_QQ_5 = FR_QQ_5, FR_QQ_5, f0
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
//
//      if (i_1!=0) U_hi = PP_1 * U_hi  
//      if (i_1!=0) U_lo = r * r  + U_lo  
//      Load PP_3 or QQ_3
//
(p9)    fma.s1 FR_poly = FR_rsq, FR_poly, FR_PP_4
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
(p9)    fma.s1 FR_U_lo = FR_r_lo, FR_U_lo, f0
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
(p10)   fma.s1 FR_U_lo = FR_QQ_1,FR_U_lo, f0
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
(p9)    fma.s1 FR_U_hi = FR_PP_1, FR_U_hi, f0
        nop.i 999
}
{ .mfi
        nop.m 999
(p10)   fma.s1 FR_poly = FR_rsq, FR_poly, FR_QQ_3
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
//
//      Load PP_2, QQ_2
//
(p9)    fma.s1 FR_poly = FR_rsq, FR_poly, FR_PP_3
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
//
//      if (i_1==0) poly = FR_rsq * poly  + PP_3
//      else        poly = FR_rsq * poly  + QQ_3
//      Load PP_1_lo
//
(p9)    fma.s1 FR_U_lo = FR_PP_1, FR_U_lo, f0
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
//
//      if (i_1 =0) poly = poly * rsq + pp_r4
//      else        poly = poly * rsq + qq_r4
//
(p9)    fma.s1 FR_U_hi = FR_r, f1, FR_U_hi
        nop.i 999
}
{ .mfi
        nop.m 999
(p10)   fma.s1 FR_poly = FR_rsq, FR_poly, FR_QQ_2
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
//
//      if (i_1==0) U_lo =  PP_1_hi * U_lo
//      else        U_lo =  QQ_1 * U_lo
//
(p9)    fma.s1 FR_poly = FR_rsq, FR_poly, FR_PP_2
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
//
//      if (i_0==0)  Result = 1
//      else         Result = -1
//
(p0)    fma.s1 FR_V = FR_U_lo, f1, FR_corr
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
(p10)   fma.s1 FR_poly = FR_rsq, FR_poly, f0
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
//
//      if (i_1==0) poly =  FR_rsq * poly + PP_2
//      else poly =  FR_rsq * poly + QQ_2
// 
(p9)    fma.s1 FR_poly = FR_rsq, FR_poly, FR_PP_1_lo
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
(p10)   fma.s1 FR_poly = FR_rsq, FR_poly, f0
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
//
//      V = U_lo + corr
//
(p9)    fma.s1 FR_poly = FR_r_cubed, FR_poly, f0
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
//
//      if (i_1==0) poly = r_cube * poly
//      else        poly = FR_rsq * poly
//
(p0)    fma.s1  FR_V = FR_poly, f1, FR_V
        nop.i 999 ;;
}
{ .mfi
        nop.m 999
(p12)   fms.s0 FR_Input_X = FR_Input_X, FR_U_hi, FR_V
        nop.i 999
}
{ .mfb
        nop.m 999
//
//      V = V + poly    
//
(p11)   fma.s0 FR_Input_X = FR_Input_X, FR_U_hi, FR_V
//
//      if (i_0==0) Result = Result * U_hi + V
//      else        Result = Result * U_hi - V
//
(p0)    br.ret.sptk   b0 
};;

//
//      If cosine, FR_Input_X = 1
//      If sine, FR_Input_X = +/-Zero (Input FR_Input_X)
//      Results are exact, no exceptions
//

L(SINCOSL_ZERO):
{ .mbb
(p0)    cmp.eq.unc p6, p7 = 0x1, GR_Sin_or_Cos
        nop.b 999
        nop.b 999 ;;
}
{ .mfi
        nop.m 999
(p7)    fmerge.s FR_Input_X = FR_Input_X, FR_Input_X
        nop.i 999
}
{ .mfb
        nop.m 999
(p6)    fmerge.s FR_Input_X = f1, f1
(p0)    br.ret.sptk   b0 ;;
}
L(SINCOSL_SPECIAL):
{ .mfb
        nop.m 999
//
//      Path for Arg = +/- QNaN, SNaN, Inf
//      Invalid can be raised. SNaNs
//      become QNaNs
//
(p0)    fmpy.s0 FR_Input_X = FR_Input_X, f0
(p0)    br.ret.sptk   b0 ;;
}
.endp cosl#
ASM_SIZE_DIRECTIVE(cosl#)

//      Call int pi_by_2_reduce(double* x, double *y)
//      for |arguments| >= 2**63
//      Address to save r and c as double 
//
//             sp+32  -> f0
//      r45    sp+16  -> f0
//      r44 -> sp     -> InputX  
//      

.proc __libm_callout
__libm_callout:
L(SINCOSL_ARG_TOO_LARGE): 
.prologue
{ .mfi
        add   r45=-32,sp                        // Parameter: r address 
        nop.f 0
.save   ar.pfs,GR_SAVE_PFS
        mov  GR_SAVE_PFS=ar.pfs                 // Save ar.pfs
}
{ .mfi
.fframe 64
        add sp=-64,sp                           // Create new stack
        nop.f 0
        mov GR_SAVE_GP=gp                       // Save gp
};;
{ .mmi
        stfe [r45] = f0,16                      // Clear Parameter r on stack
        add  r44 = 16,sp                        // Parameter x address
.save   b0, GR_SAVE_B0
        mov GR_SAVE_B0=b0                       // Save b0
};;
.body
{ .mib
        stfe [r45] = f0,-16                     // Clear Parameter c on stack 
        nop.i 0
        nop.b 0
}
{ .mib
        stfe [r44] = FR_Input_X                 // Store Parameter x on stack
        nop.i 0
(p0)    br.call.sptk b0=__libm_pi_by_2_reduce# ;;
};;
{ .mii
(p0)    ldfe  FR_Input_X =[r44],16
//
//      Get r and c off stack
//
(p0)    adds  GR_Table_Base1 = -16, GR_Table_Base1
//
//      Get r and c off stack
//
(p0)    add   GR_N_Inc = GR_Sin_or_Cos,r8 ;;
}
{ .mmb
(p0)    ldfe  FR_r =[r45],16
//
//      Get X off the stack
//      Readjust Table ptr
//
(p0)    ldfs FR_Two_to_M3 = [GR_Table_Base1],4
        nop.b 999 ;;
}
{ .mmb
(p0)    ldfs FR_Neg_Two_to_M3 = [GR_Table_Base1],0
(p0)    ldfe  FR_c =[r45]
        nop.b 999 ;;
}
{ .mfi
.restore sp
        add   sp = 64,sp                       // Restore stack pointer
(p0)    fcmp.lt.unc.s1  p6, p0 = FR_r, FR_Two_to_M3
        mov   b0 = GR_SAVE_B0                  // Restore return address
};;
{ .mib
        mov   gp = GR_SAVE_GP                  // Restore gp
        mov   ar.pfs = GR_SAVE_PFS             // Restore ar.pfs
        nop.b 0
};;
{ .mfi
        nop.m 999
(p6)    fcmp.gt.unc.s1  p6, p0 = FR_r, FR_Neg_Two_to_M3
        nop.i 999 ;;
}
{ .mib
        nop.m 999
        nop.i 999
(p6)    br.cond.spnt L(SINCOSL_SMALL_R) ;;
}
{ .mib
        nop.m 999
        nop.i 999
(p0)    br.cond.sptk L(SINCOSL_NORMAL_R) ;;
}
.endp __libm_callout
ASM_SIZE_DIRECTIVE(__libm_callout)
.type   __libm_pi_by_2_reduce#,@function
.global __libm_pi_by_2_reduce#