summaryrefslogtreecommitdiff
path: root/sysdeps/ia64/fpu/libm_sincos_large.S
diff options
context:
space:
mode:
Diffstat (limited to 'sysdeps/ia64/fpu/libm_sincos_large.S')
-rw-r--r--sysdeps/ia64/fpu/libm_sincos_large.S2757
1 files changed, 2757 insertions, 0 deletions
diff --git a/sysdeps/ia64/fpu/libm_sincos_large.S b/sysdeps/ia64/fpu/libm_sincos_large.S
new file mode 100644
index 0000000000..b09d3693a6
--- /dev/null
+++ b/sysdeps/ia64/fpu/libm_sincos_large.S
@@ -0,0 +1,2757 @@
+.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 16
+
+LOCAL_OBJECT_START(FSINCOS_CONSTANTS)
+
+data4 0x4B800000 // two**24
+data4 0xCB800000 // -two**24
+data4 0x00000000 // pad
+data4 0x00000000 // pad
+data8 0xA2F9836E4E44152A, 0x00003FFE // Inv_pi_by_2
+data8 0xC84D32B0CE81B9F1, 0x00004016 // P_0
+data8 0xC90FDAA22168C235, 0x00003FFF // P_1
+data8 0xECE675D1FC8F8CBB, 0x0000BFBD // P_2
+data8 0xB7ED8FBBACC19C60, 0x0000BF7C // P_3
+data4 0x5F000000 // two**63
+data4 0xDF000000 // -two**63
+data4 0x00000000 // pad
+data4 0x00000000 // pad
+data8 0xA397E5046EC6B45A, 0x00003FE7 // Inv_P_0
+data8 0x8D848E89DBD171A1, 0x0000BFBF // d_1
+data8 0xD5394C3618A66F8E, 0x0000BF7C // d_2
+data8 0xC90FDAA22168C234, 0x00003FFE // pi_by_4
+data8 0xC90FDAA22168C234, 0x0000BFFE // neg_pi_by_4
+data4 0x3E000000 // two**-3
+data4 0xBE000000 // -two**-3
+data4 0x00000000 // pad
+data4 0x00000000 // pad
+data4 0x2F000000 // two**-33
+data4 0xAF000000 // -two**-33
+data4 0x9E000000 // -two**-67
+data4 0x00000000 // pad
+data8 0xCC8ABEBCA21C0BC9, 0x00003FCE // PP_8
+data8 0xD7468A05720221DA, 0x0000BFD6 // PP_7
+data8 0xB092382F640AD517, 0x00003FDE // PP_6
+data8 0xD7322B47D1EB75A4, 0x0000BFE5 // PP_5
+data8 0xFFFFFFFFFFFFFFFE, 0x0000BFFD // C_1
+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
+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
+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
+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 // two**-14
+data4 0xB8800000 // -two**-14
+LOCAL_OBJECT_END(FSINCOS_CONSTANTS)
+
+// sin and cos registers
+
+// FR
+FR_Input_X = f8
+
+FR_r = f8
+FR_c = f9
+
+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_d_2 = f48
+FR_prelim = 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_Neg_Two_to_M3 = 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
+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
+
+GR_SAVE_B0 = r39
+GR_SAVE_GP = r40
+GR_SAVE_PFS = r41
+
+// sincos combined routine registers
+
+// GR
+GR_SINCOS_SAVE_PFS = r32
+GR_SINCOS_SAVE_B0 = r33
+GR_SINCOS_SAVE_GP = r34
+
+// FR
+FR_SINCOS_ARG = f100
+FR_SINCOS_RES_SIN = f101
+
+
+.section .text
+
+
+GLOBAL_LIBM_ENTRY(__libm_sincos_large)
+
+{ .mfi
+ alloc GR_SINCOS_SAVE_PFS = ar.pfs,0,3,0,0
+ fma.s1 FR_SINCOS_ARG = f8, f1, f0 // Save argument for sin and cos
+ mov GR_SINCOS_SAVE_B0 = b0
+};;
+
+{ .mfb
+ mov GR_SINCOS_SAVE_GP = gp
+ nop.f 0
+ br.call.sptk b0 = __libm_sin_large // Call sin
+};;
+
+{ .mfi
+ nop.m 0
+ fma.s1 FR_SINCOS_RES_SIN = f8, f1, f0 // Save sin result
+ nop.i 0
+};;
+
+{ .mfb
+ nop.m 0
+ fma.s1 f8 = FR_SINCOS_ARG, f1, f0 // Arg for cos
+ br.call.sptk b0 = __libm_cos_large // Call cos
+};;
+
+{ .mfi
+ mov gp = GR_SINCOS_SAVE_GP
+ fma.s1 f9 = FR_SINCOS_RES_SIN, f1, f0 // Out sin result
+ mov b0 = GR_SINCOS_SAVE_B0
+};;
+
+{ .mib
+ nop.m 0
+ mov ar.pfs = GR_SINCOS_SAVE_PFS
+ br.ret.sptk b0 // sincos_large exit
+};;
+
+GLOBAL_LIBM_END(__libm_sincos_large)
+
+
+
+
+GLOBAL_LIBM_ENTRY(__libm_sin_large)
+
+{ .mlx
+alloc GR_Table_Base = ar.pfs,0,12,2,0
+ movl GR_Sin_or_Cos = 0x0 ;;
+}
+
+{ .mmi
+ nop.m 999
+ addl GR_Table_Base = @ltoff(FSINCOS_CONSTANTS#), gp
+ nop.i 999
+}
+;;
+
+{ .mmi
+ ld8 GR_Table_Base = [GR_Table_Base]
+ nop.m 999
+ nop.i 999
+}
+;;
+
+
+{ .mib
+ nop.m 999
+ nop.i 999
+ br.cond.sptk SINCOS_CONTINUE ;;
+}
+
+GLOBAL_LIBM_END(__libm_sin_large)
+
+GLOBAL_LIBM_ENTRY(__libm_cos_large)
+
+{ .mlx
+alloc GR_Table_Base= ar.pfs,0,12,2,0
+ movl GR_Sin_or_Cos = 0x1 ;;
+}
+
+{ .mmi
+ nop.m 999
+ addl GR_Table_Base = @ltoff(FSINCOS_CONSTANTS#), gp
+ nop.i 999
+}
+;;
+
+{ .mmi
+ ld8 GR_Table_Base = [GR_Table_Base]
+ nop.m 999
+ nop.i 999
+}
+;;
+
+//
+// Load Table Address
+//
+SINCOS_CONTINUE:
+
+{ .mmi
+ add GR_Table_Base1 = 96, GR_Table_Base
+ ldfs FR_Two_to_24 = [GR_Table_Base], 4
+ nop.i 999
+}
+;;
+
+{ .mmi
+ nop.m 999
+//
+// Load 2**24, load 2**63.
+//
+ ldfs FR_Neg_Two_to_24 = [GR_Table_Base], 12
+ mov r41 = ar.pfs ;;
+}
+
+{ .mfi
+ 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
+//
+ fclass.m.unc p6, p8 = FR_Input_X, 0x1E3
+ mov r40 = gp ;;
+}
+
+{ .mfi
+ nop.m 999
+ fclass.nm.unc p8, p0 = FR_Input_X, 0x1FF
+// GR_Sin_or_Cos denotes
+ mov r39 = b0
+}
+
+{ .mfb
+ ldfs FR_Neg_Two_to_63 = [GR_Table_Base1], 12
+ fclass.m.unc p10, p0 = FR_Input_X, 0x007
+(p6) br.cond.spnt SINCOS_SPECIAL ;;
+}
+
+{ .mib
+ nop.m 999
+ nop.i 999
+(p8) br.cond.spnt SINCOS_SPECIAL ;;
+}
+
+{ .mib
+ nop.m 999
+ nop.i 999
+//
+// Branch if +/- NaN, Inf.
+// Load -2**24, load -2**63.
+//
+(p10) br.cond.spnt SINCOS_ZERO ;;
+}
+
+{ .mmb
+ ldfe FR_Inv_pi_by_2 = [GR_Table_Base], 16
+ ldfe FR_Inv_P_0 = [GR_Table_Base1], 16
+ nop.b 999 ;;
+}
+
+{ .mmb
+ nop.m 999
+ ldfe FR_d_1 = [GR_Table_Base1], 16
+ nop.b 999 ;;
+}
+//
+// Raise possible denormal operand flag with useful fcmp
+// Is x <= -2**63
+// Load Inv_P_0 for pre-reduction
+// Load Inv_pi_by_2
+//
+
+{ .mmb
+ ldfe FR_P_0 = [GR_Table_Base], 16
+ ldfe FR_d_2 = [GR_Table_Base1], 16
+ nop.b 999 ;;
+}
+//
+// Load P_0
+// Load d_1
+// Is x >= 2**63
+// Is x <= -2**24?
+//
+
+{ .mmi
+ ldfe FR_P_1 = [GR_Table_Base], 16 ;;
+//
+// Load P_1
+// Load d_2
+// Is x >= 2**24?
+//
+ ldfe FR_P_2 = [GR_Table_Base], 16
+ nop.i 999 ;;
+}
+
+{ .mmf
+ nop.m 999
+ ldfe FR_P_3 = [GR_Table_Base], 16
+ fcmp.le.unc.s1 p7, p8 = FR_Input_X, FR_Neg_Two_to_24
+}
+
+{ .mfi
+ nop.m 999
+//
+// 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
+//
+ fcmp.le.unc.s1 p10, p11 = FR_Input_X, FR_Neg_Two_to_63
+ nop.i 999 ;;
+}
+
+{ .mfi
+ nop.m 999
+(p8) fcmp.ge.s1 p7, p0 = FR_Input_X, FR_Two_to_24
+ nop.i 999
+}
+
+{ .mfi
+ ldfe FR_Pi_by_4 = [GR_Table_Base1], 16
+(p11) fcmp.ge.s1 p10, p0 = FR_Input_X, FR_Two_to_63
+ nop.i 999 ;;
+}
+
+{ .mmi
+ ldfe FR_Neg_Pi_by_4 = [GR_Table_Base1], 16 ;;
+ ldfs FR_Two_to_M3 = [GR_Table_Base1], 4
+ nop.i 999 ;;
+}
+
+{ .mib
+ 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 SINCOS_ARG_TOO_LARGE ;;
+}
+
+{ .mib
+ nop.m 999
+ nop.i 999
+//
+// Branch out if x >= 2**63. Use Payne-Hanek Reduction
+//
+(p7) br.cond.spnt SINCOS_LARGER_ARG ;;
+}
+
+{ .mfi
+ nop.m 999
+//
+// Branch if Arg <= -2**24 or Arg >= 2**24 and use pre-reduction.
+//
+ fma.s1 FR_N_float = FR_Input_X, FR_Inv_pi_by_2, f0
+ nop.i 999 ;;
+}
+
+{ .mfi
+ nop.m 999
+ 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
+//
+ 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 SINCOS_SMALL_R ;;
+}
+
+{ .mib
+ nop.m 999
+ nop.i 999
+(p6) br.cond.sptk SINCOS_NORMAL_R ;;
+}
+//
+// if |x| < pi/4, branch based on |x| < 2**(-3) or otherwise.
+//
+//
+// In this branch, |x| >= pi/4.
+//
+
+{ .mfi
+ ldfs FR_Neg_Two_to_M67 = [GR_Table_Base1], 8
+//
+// Load -2**(-67)
+//
+ fnma.s1 FR_s = FR_N_float, FR_P_1, FR_Input_X
+//
+// w = N * P_2
+// s = -N * P_1 + Arg
+//
+ add GR_N_Inc = GR_N_Inc, GR_Sin_or_Cos
+}
+
+{ .mfi
+ nop.m 999
+ 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
+//
+ 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 SINCOS_SMALL_R ;;
+}
+
+{ .mib
+ nop.m 999
+ nop.i 999
+(p9) br.cond.sptk SINCOS_NORMAL_R ;;
+}
+
+{ .mfi
+(p7) add GR_Table_Base1 = 224, GR_Table_Base1
+//
+// Branch to SINCOS_SMALL_R or SINCOS_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
+(p7) fmerge.se FR_prelim = FR_r, FR_r
+ nop.i 999 ;;
+}
+
+{ .mfi
+ nop.m 999
+//(p10) fma.s1 FR_Input_X = f0, f1, f1
+(p10) fma.s1 FR_prelim = 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
+(p12) fnma.s1 FR_prelim = FR_prelim, 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.d.s1 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.s1 FR_Input_X = FR_prelim, f1, FR_poly
+ nop.i 999 ;;
+}
+
+{ .mfb
+ nop.m 999
+(p12) fms.s1 FR_Input_X = FR_prelim, f1, FR_poly
+//
+// if (i_0 == 0), Result = Result + poly
+// else Result = Result - poly
+//
+ br.ret.sptk b0 ;;
+}
+SINCOS_LARGER_ARG:
+
+{ .mfi
+ nop.m 999
+ 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
+ addl GR_Table_Base = @ltoff(FSINCOS_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
+ add GR_Table_Base = 688, GR_Table_Base ;;
+ ldfs FR_Two_to_M14 = [GR_Table_Base], 4
+ nop.i 999 ;;
+}
+
+{ .mfi
+ 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)
+//
+ 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
+//
+ fcvt.xf FR_N_0 = FR_N_0_fix
+ nop.i 999 ;;
+}
+
+{ .mfi
+ nop.m 999
+//
+// Make N_0 the integer part
+//
+ fnma.s1 FR_ArgPrime = FR_N_0, FR_P_0, FR_Input_X
+ nop.i 999
+}
+
+{ .mfi
+ nop.m 999
+ 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
+//
+ fma.s1 FR_N_float = FR_ArgPrime, FR_Inv_pi_by_2, f0
+ nop.i 999 ;;
+}
+
+{ .mfi
+ nop.m 999
+//
+// N = A' * 2/pi
+//
+ fcvt.fx.s1 FR_N_fix = FR_N_float
+ nop.i 999 ;;
+}
+
+{ .mfi
+ nop.m 999
+//
+// N_fix is the integer part
+//
+ fcvt.xf FR_N_float = FR_N_fix
+ nop.i 999 ;;
+}
+
+{ .mfi
+ getf.sig GR_N_Inc = FR_N_fix
+ nop.f 999
+ nop.i 999 ;;
+}
+
+{ .mii
+ nop.m 999
+ nop.i 999 ;;
+ 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
+//
+ fnma.s1 FR_s = FR_N_float, FR_P_1, FR_ArgPrime
+ nop.i 999
+}
+
+{ .mfi
+ nop.m 999
+ 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
+//
+ 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 ;;
+}
+
+{ .mfi
+ nop.m 999
+(p9) fms.s1 FR_U_lo = FR_N_0, FR_d_1, FR_U_hi
+ nop.i 999 ;;
+}
+
+{ .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 SINCOS_SMALL_R ;;
+}
+
+{ .mib
+ nop.m 999
+ nop.i 999
+(p13) br.cond.sptk SINCOS_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
+ 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
+ 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
+ addl GR_Table_Base = @ltoff(FSINCOS_CONSTANTS#), gp
+ nop.i 999
+}
+;;
+
+{ .mmi
+ ld8 GR_Table_Base = [GR_Table_Base]
+ nop.m 999
+ nop.i 999
+}
+;;
+
+
+{ .mfi
+ 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
+ ldfe FR_C_1 = [GR_Table_Base], 16 ;;
+ ldfe FR_C_2 = [GR_Table_Base], 64
+ nop.i 999 ;;
+}
+
+{ .mmf
+ nop.m 999
+//
+// c = c + C_lo finished.
+// Load C_2
+//
+ ldfe FR_S_1 = [GR_Table_Base], 16
+//
+// C_lo = S - C_hi
+//
+ fma.s1 FR_t = FR_t, f1, FR_w ;;
+}
+//
+// 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
+ ldfe FR_S_2 = [GR_Table_Base], 64
+//
+// t = t + w
+//
+(p10) fms.s1 FR_a = FR_a, f1, FR_V_hi
+ 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
+//
+ 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
+ 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.
+//
+ fma.s1 FR_C_lo = FR_C_lo, f1, FR_A
+ nop.i 999 ;;
+}
+
+{ .mfi
+ nop.m 999
+//
+// C_lo = (S - C_hi) + A
+//
+ fma.s1 FR_t = FR_t, f1, FR_a
+ nop.i 999 ;;
+}
+
+{ .mfi
+ nop.m 999
+//
+// t = t + a
+//
+ 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
+//
+ fma.s1 FR_r = FR_C_hi, f1, FR_C_lo
+ nop.i 999 ;;
+}
+
+{ .mfi
+ nop.m 999
+//
+// Load S_2
+//
+ 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
+//
+ 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
+(p11) fma.s1 FR_prelim = f0, f1, FR_r
+ nop.i 999 ;;
+}
+
+{ .mfi
+ nop.m 999
+//(p12) fma.s1 FR_Input_X = f0, f1, f1
+(p12) fma.s1 FR_prelim = 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 ?
+//
+ fma.s1 FR_r_cubed = FR_rsq, FR_r, f0
+ nop.i 999 ;;
+}
+
+{ .mfi
+ nop.m 999
+//
+// c = C_hi - r
+// Load C_1
+//
+ 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
+(p10) fms.s1 FR_prelim = f0, f1, FR_prelim
+ 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.s1 FR_Input_X = FR_prelim, f1, FR_poly
+ nop.i 999 ;;
+}
+
+{ .mfb
+ nop.m 999
+(p10) fms.s1 FR_Input_X = FR_prelim, f1, FR_poly
+//
+// if i_0 == 0: Result = Result + poly
+// else Result = Result - poly
+//
+ br.ret.sptk b0 ;;
+}
+SINCOS_SMALL_R:
+
+{ .mii
+ nop.m 999
+ 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.
+//
+ cmp.eq.unc p9, p10 = 0x0, GR_i_1 ;;
+}
+
+{ .mfi
+ nop.m 999
+ fma.s1 FR_rsq = FR_r, FR_r, f0
+ 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
+ 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
+ addl GR_Table_Base = @ltoff(FSINCOS_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
+//
+ 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.d.s1 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
+//
+ 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.d.s1 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
+//
+ 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
+//
+ fma.s1 FR_poly = FR_poly, f1, FR_poly_hi
+ nop.i 999 ;;
+}
+
+{ .mfi
+ nop.m 999
+(p12) fms.s1 FR_Input_X = FR_r, f1, FR_poly
+ nop.i 999
+}
+
+{ .mfb
+ nop.m 999
+//
+// poly = poly + poly_hi
+//
+(p11) fma.s1 FR_Input_X = FR_r, f1, FR_poly
+//
+// if (i_0 == 0) Result = r + poly
+// if (i_0 != 0) Result = r - poly
+//
+ br.ret.sptk b0 ;;
+}
+SINCOS_NORMAL_R:
+
+{ .mii
+ nop.m 999
+ extr.u GR_i_1 = GR_N_Inc, 0, 1 ;;
+//
+// Set table_ptr1 and table_ptr2 to base address of
+// constant table.
+ cmp.eq.unc p9, p10 = 0x0, GR_i_1 ;;
+}
+
+{ .mfi
+ nop.m 999
+ fma.s1 FR_rsq = FR_r, FR_r, f0
+ extr.u GR_i_0 = GR_N_Inc, 1, 1 ;;
+}
+
+{ .mfi
+ nop.m 999
+ frcpa.s1 FR_r_hi, p6 = f1, FR_r
+ 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
+ addl GR_Table_Base = @ltoff(FSINCOS_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
+(p12) fms.s1 FR_prelim = f0, f1, f1
+(p9) add GR_Table_Base = 224, GR_Table_Base ;;
+}
+
+{ .mmf
+ nop.m 999
+(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 ;;
+(p11) fma.s1 FR_prelim = f0, f1, f1 ;;
+}
+
+{ .mmf
+(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
+ fma.s1 FR_r_cubed = FR_r, FR_rsq, f0 ;;
+}
+
+{ .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
+//
+ 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 ;;
+}
+
+{ .mmi
+(p10) ldfe FR_QQ_4 = [GR_Table_Base], 16 ;;
+(p9) ldfe FR_PP_1 = [GR_Table_Base], 16
+ nop.i 999 ;;
+}
+
+{ .mmf
+(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
+(p10) fma.s1 FR_poly = FR_rsq, FR_QQ_8, FR_QQ_7 ;;
+}
+//
+// if (i_1=0) poly = rsq * poly + PP_5
+// else poly = rsq * poly + QQ_5
+// Load PP_4 or QQ_4
+//
+
+{ .mmf
+(p9) ldfe FR_PP_3 = [GR_Table_Base], 16
+(p10) ldfe FR_QQ_2 = [GR_Table_Base], 16
+//
+// r_hi = frcpa(frcpa(r)).
+// r_cube = r * FR_rsq.
+//
+(p9) fma.s1 FR_poly = FR_rsq, FR_PP_8, FR_PP_7 ;;
+}
+//
+// 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 ;;
+}
+
+{ .mmf
+ nop.m 999
+(p9) ldfe FR_PP_1_lo = [GR_Table_Base], 16
+(p10) fma.s1 FR_corr = FR_S_1, FR_r_cubed, FR_r
+}
+
+{ .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
+//
+ 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
+//
+ 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.d.s1 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.d.s1 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
+//
+ 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
+//
+ fma.s1 FR_V = FR_poly, f1, FR_V
+ nop.i 999 ;;
+}
+
+{ .mfi
+ nop.m 999
+//(p12) fms.s1 FR_Input_X = FR_Input_X, FR_U_hi, FR_V
+(p12) fms.s1 FR_Input_X = FR_prelim, FR_U_hi, FR_V
+ nop.i 999
+}
+
+{ .mfb
+ nop.m 999
+//
+// V = V + poly
+//
+//(p11) fma.s1 FR_Input_X = FR_Input_X, FR_U_hi, FR_V
+(p11) fma.s1 FR_Input_X = FR_prelim, FR_U_hi, FR_V
+//
+// if (i_0==0) Result = Result * U_hi + V
+// else Result = Result * U_hi - V
+//
+ 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
+//
+SINCOS_ZERO:
+
+{ .mmb
+ cmp.eq.unc p6, p7 = 0x1, GR_Sin_or_Cos
+ nop.m 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
+ br.ret.sptk b0 ;;
+}
+
+SINCOS_SPECIAL:
+
+//
+// Path for Arg = +/- QNaN, SNaN, Inf
+// Invalid can be raised. SNaNs
+// become QNaNs
+//
+
+{ .mfb
+ nop.m 999
+ fmpy.s1 FR_Input_X = FR_Input_X, f0
+ br.ret.sptk b0 ;;
+}
+GLOBAL_LIBM_END(__libm_cos_large)
+
+
+// *******************************************************************
+// *******************************************************************
+// *******************************************************************
+//
+// Special Code to handle very large argument case.
+// Call int __libm_pi_by_2_reduce(x,r,c) for |arguments| >= 2**63
+// The interface is custom:
+// On input:
+// (Arg or x) is in f8
+// On output:
+// r is in f8
+// c is in f9
+// N is in r8
+// Be sure to allocate at least 2 GP registers as output registers for
+// __libm_pi_by_2_reduce. This routine uses r49-50. These are used as
+// scratch registers within the __libm_pi_by_2_reduce routine (for speed).
+//
+// We know also that __libm_pi_by_2_reduce preserves f10-15, f71-127. We
+// use this to eliminate save/restore of key fp registers in this calling
+// function.
+//
+// *******************************************************************
+// *******************************************************************
+// *******************************************************************
+
+LOCAL_LIBM_ENTRY(__libm_callout_2)
+SINCOS_ARG_TOO_LARGE:
+
+.prologue
+// Readjust Table ptr
+{ .mfi
+ adds GR_Table_Base1 = -16, GR_Table_Base1
+ nop.f 999
+.save ar.pfs,GR_SAVE_PFS
+ mov GR_SAVE_PFS=ar.pfs // Save ar.pfs
+};;
+
+{ .mmi
+ ldfs FR_Two_to_M3 = [GR_Table_Base1],4
+ mov GR_SAVE_GP=gp // Save gp
+.save b0, GR_SAVE_B0
+ mov GR_SAVE_B0=b0 // Save b0
+};;
+
+.body
+//
+// Call argument reduction with x in f8
+// Returns with N in r8, r in f8, c in f9
+// Assumes f71-127 are preserved across the call
+//
+{ .mib
+ ldfs FR_Neg_Two_to_M3 = [GR_Table_Base1],0
+ nop.i 0
+ br.call.sptk b0=__libm_pi_by_2_reduce#
+};;
+
+{ .mfi
+ add GR_N_Inc = GR_Sin_or_Cos,r8
+ fcmp.lt.unc.s1 p6, p0 = FR_r, FR_Two_to_M3
+ mov b0 = GR_SAVE_B0 // Restore return address
+};;
+
+{ .mfi
+ mov gp = GR_SAVE_GP // Restore gp
+(p6) fcmp.gt.unc.s1 p6, p0 = FR_r, FR_Neg_Two_to_M3
+ mov ar.pfs = GR_SAVE_PFS // Restore ar.pfs
+};;
+
+{ .mbb
+ nop.m 999
+(p6) br.cond.spnt SINCOS_SMALL_R // Branch if |r| < 1/4
+ br.cond.sptk SINCOS_NORMAL_R ;; // Branch if 1/4 <= |r| < pi/4
+}
+
+LOCAL_LIBM_END(__libm_callout_2)
+
+.type __libm_pi_by_2_reduce#,@function
+.global __libm_pi_by_2_reduce#
+