diff options
Diffstat (limited to 'sysdeps/ia64/fpu/s_cosl.S')
| -rw-r--r-- | sysdeps/ia64/fpu/s_cosl.S | 2760 |
1 files changed, 1303 insertions, 1457 deletions
diff --git a/sysdeps/ia64/fpu/s_cosl.S b/sysdeps/ia64/fpu/s_cosl.S index 2755580..8d71e50 100644 --- a/sysdeps/ia64/fpu/s_cosl.S +++ b/sysdeps/ia64/fpu/s_cosl.S @@ -1,10 +1,10 @@ .file "sincosl.s" -// Copyright (C) 2000, 2001, Intel Corporation + +// Copyright (c) 2000 - 2004, 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. +// +// Contributed 2000 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 @@ -20,76 +20,82 @@ // * 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 + +// 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 +// 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 +// 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. -// +// 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. +// problem reports or change requests be submitted to it directly at +// http://www.intel.com/software/products/opensource/libraries/num.htm. // -// ********************************************************************* +//********************************************************************* // -// History: -// 2/02/2000 (hand-optimized) -// 4/04/00 Unwind support added +// History: +// 02/02/00 (hand-optimized) +// 04/04/00 Unwind support added +// 07/30/01 Improved speed on all paths +// 08/20/01 Fixed bundling typo +// 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 +// 10/13/03 Corrected final .endp name to match .proc +// 10/26/04 Avoided using r14-31 as scratch so not clobbered by dynamic loader // -// ********************************************************************* +//********************************************************************* // // 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) +// 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) +// r32-r58 // // 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 +// 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(inf) = QNaN // sinl(+/-0) = +/-0 -// cosl(inf) = QNaN +// 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 +// 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. @@ -98,62 +104,62 @@ // // 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 +// 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. +// 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 +// 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, +// 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 +// if M = -1, I = 11 // sin ((-pi/2 + alpha) = (-1) cos (alpha) -// if M = 0, I = 00 +// if M = 0, I = 00 // sin (alpha) = sin (alpha) -// if M = 1, I = 01 +// if M = 1, I = 01 // sin (pi/2 + alpha) = cos (alpha) -// if M = 2, I = 10 +// if M = 2, I = 10 // sin (pi + alpha) = (-1) sin (alpha) -// if M = 3, I = 11 +// if M = 3, I = 11 // sin ((3/2)pi + alpha) = (-1) cos (alpha) // -// The value of alpha is obtained by argument reduction and +// 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 +// 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 +// 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 +// 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 +// 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 @@ -163,32 +169,32 @@ // // 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. +// 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. +// 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 +// 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 greatest challenge of this task is that the second terms of // the Taylor series -// -// r - r^3/3! + r^r/5! - ... +// +// r - r^3/3! + r^r/5! - ... // // and // -// 1 - r^2/2! + r^4/4! - ... +// 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 +// 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) @@ -197,88 +203,88 @@ // 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 +// 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 +// 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 +// 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 +// 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. +// 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 +// 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 +// 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 +// 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 +// 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 +// 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 +// 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 +// 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 +// 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| +// |U_lo| <= 2^(-8) |U_hi| // -// says that U_hi + U_lo is approximating r + PP_1_hi r^3 to roughly +// 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 ) ). +// 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 // @@ -297,35 +303,35 @@ // End // // Finally, -// -// V := poly + ( U_lo + correction ) +// +// V := poly + ( U_lo + correction ) // // / U_hi + V if i_0 = 0 -// result := | +// 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. +// 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 +// 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. +// 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 +// * 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. @@ -335,18 +341,18 @@ // // Step 2. Reduction of moderate arguments. // -// If |Arg| < pi/4 ...quick branch -// N_fix := N_inc (integer) +// 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) +// 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) +// 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 @@ -358,8 +364,8 @@ // Step 3. Case_1_reduce. // // r := s + w -// c := (s - r) + w ...observe order -// +// c := (s - r) + w ...observe order +// // Step 4. Case_1_complete // // ...At this point, the reduced argument alpha is @@ -375,17 +381,17 @@ // // 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_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 +// 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_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 +// correction := -c*(r + S_1*FR_rsq*r) ...any order // Endif // -// V := poly + (U_lo + correction) ...observe order +// V := poly + (U_lo + correction) ...observe order // // result := (i_0 == 0? 1.0 : -1.0) // @@ -397,7 +403,7 @@ // 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. // @@ -412,7 +418,7 @@ // 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) +// poly_hi := FR_rsq*(C_1 + FR_rsq*C_2) // correction := -c*r // result := 1 // Endif @@ -429,15 +435,15 @@ // // Step 7. Case_2_reduce. // -// ...Refer to the write up for argument reduction for +// ...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 := 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 // @@ -446,29 +452,29 @@ // ...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 +// ...|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 +// 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. @@ -476,13 +482,13 @@ // 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 +// N_fix := N_fix + N_inc // // s := Arg' - N * P_1 // w := w - N * P_2 @@ -494,15 +500,15 @@ // 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 +// 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 +// 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 @@ -517,12 +523,12 @@ // 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 +// ...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_hi := s + A t := t + a +// C_lo := (s - C_hi) + A // C_lo := C_lo + t // // ...Final steps for reduction @@ -548,156 +554,192 @@ // 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 + +RODATA +.align 16 + +LOCAL_OBJECT_START(FSINCOSL_CONSTANTS) + +sincosl_table_p: +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 +data8 0x8D848E89DBD171A1, 0x0000BFBF // d_1 +data8 0xD5394C3618A66F8E, 0x0000BF7C // d_2 +LOCAL_OBJECT_END(FSINCOSL_CONSTANTS) + +LOCAL_OBJECT_START(sincosl_table_d) +data8 0xC90FDAA22168C234, 0x00003FFE // pi_by_4 +data8 0xA397E5046EC6B45A, 0x00003FE7 // Inv_P_0 +data4 0x3E000000, 0xBE000000 // 2^-3 and -2^-3 +data4 0x2F000000, 0xAF000000 // 2^-33 and -2^-33 +data4 0x9E000000, 0x00000000 // -2^-67 +data4 0x00000000, 0x00000000 // pad +LOCAL_OBJECT_END(sincosl_table_d) + +LOCAL_OBJECT_START(sincosl_table_pp) +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 +LOCAL_OBJECT_END(sincosl_table_pp) + +LOCAL_OBJECT_START(sincosl_table_qq) +data8 0xD56232EFC2B0FE52, 0x00003FD2 // QQ_8 +data8 0xC9C99ABA2B48DCA6, 0x0000BFDA // QQ_7 +data8 0x8F76C6509C716658, 0x00003FE2 // QQ_6 +data8 0x93F27DBAFDA8D0FC, 0x0000BFE9 // QQ_5 +data8 0xAAAAAAAAAAAAAAAA, 0x0000BFFC // S_1 +data8 0x8000000000000000, 0x0000BFFE // QQ_1 +data8 0xD00D00D00C6E5041, 0x00003FEF // QQ_4 +data8 0xB60B60B60B607F60, 0x0000BFF5 // QQ_3 +data8 0xAAAAAAAAAAAAAA9B, 0x00003FFA // QQ_2 +LOCAL_OBJECT_END(sincosl_table_qq) + +LOCAL_OBJECT_START(sincosl_table_c) +data8 0xFFFFFFFFFFFFFFFE, 0x0000BFFD // C_1 +data8 0xAAAAAAAAAAAA719F, 0x00003FFA // C_2 +data8 0xB60B60B60356F994, 0x0000BFF5 // C_3 +data8 0xD00CFFD5B2385EA9, 0x00003FEF // C_4 +data8 0x93E4BD18292A14CD, 0x0000BFE9 // C_5 +LOCAL_OBJECT_END(sincosl_table_c) + +LOCAL_OBJECT_START(sincosl_table_s) +data8 0xAAAAAAAAAAAAAAAA, 0x0000BFFC // S_1 +data8 0x88888888888868DB, 0x00003FF8 // S_2 +data8 0xD00D00D0055EFD4B, 0x0000BFF2 // S_3 +data8 0xB8EF1C5D839730B9, 0x00003FEC // S_4 +data8 0xD71EA3A4E5B3F492, 0x0000BFE5 // S_5 +data4 0x38800000, 0xB8800000 // two**-14 and -two**-14 +LOCAL_OBJECT_END(sincosl_table_s) + +FR_Input_X = f8 +FR_Result = f8 + +FR_r = f8 +FR_c = f9 + +FR_norm_x = f9 +FR_inv_pi_2to63 = f10 +FR_rshf_2to64 = f11 +FR_2tom64 = f12 +FR_rshf = f13 +FR_N_float_signif = f14 +FR_abs_x = f15 +FR_Pi_by_4 = f34 +FR_Two_to_M14 = f35 +FR_Neg_Two_to_M14 = f36 +FR_Two_to_M33 = f37 +FR_Neg_Two_to_M33 = f38 +FR_Neg_Two_to_M67 = f39 +FR_Inv_pi_by_2 = f40 +FR_N_float = f41 +FR_N_fix = f42 +FR_P_1 = f43 +FR_P_2 = f44 +FR_P_3 = f45 +FR_s = f46 +FR_w = f47 +FR_d_2 = f48 +FR_tmp_result = f49 +FR_Z = f50 +FR_A = f51 +FR_a = f52 +FR_t = f53 +FR_U_1 = f54 +FR_U_2 = f55 +FR_C_1 = f56 +FR_C_2 = f57 +FR_C_3 = f58 +FR_C_4 = f59 +FR_C_5 = f60 +FR_S_1 = f61 +FR_S_2 = f62 +FR_S_3 = f63 +FR_S_4 = f64 +FR_S_5 = f65 +FR_poly_hi = f66 +FR_poly_lo = f67 +FR_r_hi = f68 +FR_r_lo = f69 +FR_rsq = f70 +FR_r_cubed = f71 +FR_C_hi = f72 +FR_N_0 = f73 +FR_d_1 = f74 +FR_V = f75 +FR_V_hi = f75 +FR_V_lo = f76 +FR_U_hi = f77 +FR_U_lo = f78 +FR_U_hiabs = f79 +FR_V_hiabs = f80 +FR_PP_8 = f81 +FR_QQ_8 = f101 +FR_PP_7 = f82 +FR_QQ_7 = f102 +FR_PP_6 = f83 +FR_QQ_6 = f103 +FR_PP_5 = f84 +FR_QQ_5 = f104 +FR_PP_4 = f85 +FR_QQ_4 = f105 +FR_PP_3 = f86 +FR_QQ_3 = f106 +FR_PP_2 = f87 +FR_QQ_2 = f107 +FR_QQ_1 = f108 +FR_r_hi_sq = f88 +FR_N_0_fix = f89 +FR_Inv_P_0 = f90 +FR_corr = f91 +FR_poly = f92 +FR_Neg_Two_to_M3 = f93 +FR_Two_to_M3 = f94 +FR_P_0 = f95 +FR_C_lo = f96 +FR_PP_1 = f97 +FR_PP_1_lo = f98 +FR_ArgPrime = f99 +FR_inexact = f100 + +GR_exp_m2_to_m3= r36 +GR_N_Inc = r37 +GR_Sin_or_Cos = r38 +GR_signexp_x = r40 +GR_exp_x = r40 +GR_exp_mask = r41 +GR_exp_2_to_63 = r42 +GR_exp_2_to_m3 = r43 +GR_exp_2_to_24 = r44 + +GR_sig_inv_pi = r45 +GR_rshf_2to64 = r46 +GR_exp_2tom64 = r47 +GR_rshf = r48 +GR_ad_p = r49 +GR_ad_d = r50 +GR_ad_pp = r51 +GR_ad_qq = r52 +GR_ad_c = r53 +GR_ad_s = r54 +GR_ad_ce = r55 +GR_ad_se = r56 +GR_ad_m14 = r57 +GR_ad_s1 = r58 // Added for unwind support @@ -706,386 +748,377 @@ 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 + +GLOBAL_IEEE754_ENTRY(sinl) { .mlx -alloc GR_Table_Base = ar.pfs,0,12,2,0 -(p0) movl GR_Sin_or_Cos = 0x0 ;; + alloc r32 = ar.pfs,0,27,2,0 + movl GR_sig_inv_pi = 0xa2f9836e4e44152a // significand of 1/pi } - -{ .mmi - nop.m 999 -(p0) addl GR_Table_Base = @ltoff(FSINCOSL_CONSTANTS#), gp - nop.i 999 +{ .mlx + mov GR_Sin_or_Cos = 0x0 + movl GR_rshf_2to64 = 0x47e8000000000000 // 1.1000 2^(63+64) } ;; -{ .mmb - ld8 GR_Table_Base = [GR_Table_Base] +{ .mfi + addl GR_ad_p = @ltoff(FSINCOSL_CONSTANTS#), gp + fclass.m p6, p0 = FR_Input_X, 0x1E3 // Test x natval, nan, inf + mov GR_exp_2_to_m3 = 0xffff - 3 // Exponent of 2^-3 +} +{ .mfb nop.m 999 -(p0) br.cond.sptk L(SINCOSL_CONTINUE) ;; + fnorm.s1 FR_norm_x = FR_Input_X // Normalize x + br.cond.sptk SINCOSL_CONTINUE } ;; +GLOBAL_IEEE754_END(sinl) -.endp sinl# -ASM_SIZE_DIRECTIVE(sinl#) - -.section .text -.proc cosl# -cosl: -#ifdef _LIBC -.proc __cosl# -__cosl: -#endif +GLOBAL_IEEE754_ENTRY(cosl) +{ .mlx + alloc r32 = ar.pfs,0,27,2,0 + movl GR_sig_inv_pi = 0xa2f9836e4e44152a // significand of 1/pi +} { .mlx -alloc GR_Table_Base= ar.pfs,0,12,2,0 -(p0) movl GR_Sin_or_Cos = 0x1 ;; + mov GR_Sin_or_Cos = 0x1 + movl GR_rshf_2to64 = 0x47e8000000000000 // 1.1000 2^(63+64) } ;; -{ .mmi +{ .mfi + addl GR_ad_p = @ltoff(FSINCOSL_CONSTANTS#), gp + fclass.m p6, p0 = FR_Input_X, 0x1E3 // Test x natval, nan, inf + mov GR_exp_2_to_m3 = 0xffff - 3 // Exponent of 2^-3 +} +{ .mfi nop.m 999 -(p0) addl GR_Table_Base = @ltoff(FSINCOSL_CONSTANTS#), gp + fnorm.s1 FR_norm_x = FR_Input_X // Normalize x nop.i 999 } ;; -{ .mmb - ld8 GR_Table_Base = [GR_Table_Base] - nop.m 999 - nop.b 999 +SINCOSL_CONTINUE: +{ .mfi + setf.sig FR_inv_pi_2to63 = GR_sig_inv_pi // Form 1/pi * 2^63 + nop.f 999 + mov GR_exp_2tom64 = 0xffff - 64 // Scaling constant to compute N +} +{ .mlx + setf.d FR_rshf_2to64 = GR_rshf_2to64 // Form const 1.1000 * 2^(63+64) + movl GR_rshf = 0x43e8000000000000 // Form const 1.1000 * 2^63 } ;; +{ .mfi + ld8 GR_ad_p = [GR_ad_p] // Point to Inv_pi_by_2 + fclass.m p7, p0 = FR_Input_X, 0x0b // Test x denormal + nop.i 999 +} +;; - -// -// 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 ;; +{ .mfi + getf.exp GR_signexp_x = FR_Input_X // Get sign and exponent of x + fclass.m p10, p0 = FR_Input_X, 0x007 // Test x zero + nop.i 999 } -{ .mmi - nop.m 0 -// -// Load 2**24, load 2**63. -// -(p0) ldfs FR_Neg_Two_to_24 = [GR_Table_Base], 12 - nop.i 0 +{ .mib + mov GR_exp_mask = 0x1ffff // Exponent mask + nop.i 999 +(p6) br.cond.spnt SINCOSL_SPECIAL // Branch if x natval, nan, inf } +;; + { .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) ;; + setf.exp FR_2tom64 = GR_exp_2tom64 // Form 2^-64 for scaling N_float + nop.f 0 + add GR_ad_d = 0x70, GR_ad_p // Point to constant table d } { .mib - nop.m 999 - nop.i 999 -(p8) br.cond.spnt L(SINCOSL_SPECIAL) ;; + setf.d FR_rshf = GR_rshf // Form right shift const 1.1000 * 2^63 + mov GR_exp_m2_to_m3 = 0x2fffc // Form -(2^-3) +(p7) br.cond.spnt SINCOSL_DENORMAL // Branch if x denormal } -{ .mib - nop.m 999 - nop.i 999 -// -// Branch if +/- NaN, Inf. -// Load -2**24, load -2**63. -// -(p10) br.cond.spnt L(SINCOSL_ZERO) ;; +;; + +SINCOSL_COMMON: +{ .mfi + and GR_exp_x = GR_exp_mask, GR_signexp_x // Get exponent of x + fclass.nm p8, p0 = FR_Input_X, 0x1FF // Test x unsupported type + mov GR_exp_2_to_63 = 0xffff + 63 // Exponent of 2^63 } -{ .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 ;; +{ .mib + add GR_ad_pp = 0x40, GR_ad_d // Point to constant table pp + mov GR_exp_2_to_24 = 0xffff + 24 // Exponent of 2^24 +(p10) br.cond.spnt SINCOSL_ZERO // Branch if x zero } -{ .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 ;; +;; + +{ .mfi + ldfe FR_Inv_pi_by_2 = [GR_ad_p], 16 // Load 2/pi + fcmp.eq.s0 p15, p0 = FR_Input_X, f0 // Dummy to set denormal + add GR_ad_qq = 0xa0, GR_ad_pp // Point to constant table qq } -{ .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 ;; +{ .mfi + ldfe FR_Pi_by_4 = [GR_ad_d], 16 // Load pi/4 for range test + nop.f 999 + cmp.ge p10,p0 = GR_exp_x, GR_exp_2_to_63 // Is |x| >= 2^63 } -// -// 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 ;; + ldfe FR_P_0 = [GR_ad_p], 16 // Load P_0 for pi/4 <= |x| < 2^63 + fmerge.s FR_abs_x = f1, FR_norm_x // |x| + add GR_ad_c = 0x90, GR_ad_qq // Point to constant table c } -{ .mbb -(p0) ldfe FR_P_3 = [GR_Table_Base], 16 - nop.b 999 - nop.b 999 ;; +{ .mfi + ldfe FR_Inv_P_0 = [GR_ad_d], 16 // Load 1/P_0 for pi/4 <= |x| < 2^63 + nop.f 999 + cmp.ge p7,p0 = GR_exp_x, GR_exp_2_to_24 // Is |x| >= 2^24 } +;; + { .mfi - nop.m 999 -(p8) fcmp.ge.s1 p7, p0 = FR_Input_X, FR_Two_to_24 - nop.i 999 + ldfe FR_P_1 = [GR_ad_p], 16 // Load P_1 for pi/4 <= |x| < 2^63 + nop.f 999 + add GR_ad_s = 0x50, GR_ad_c // Point to constant table s } { .mfi -(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 ;; + ldfe FR_PP_8 = [GR_ad_pp], 16 // Load PP_8 for 2^-3 < |r| < pi/4 + nop.f 999 + 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 + ldfe FR_P_2 = [GR_ad_p], 16 // Load P_2 for pi/4 <= |x| < 2^63 + nop.f 999 + add GR_ad_ce = 0x40, GR_ad_c // Point to end of constant table c } { .mfi - nop.m 999 -(p11) fcmp.ge.s1 p10, p0 = FR_Input_X, FR_Two_to_63 - nop.i 999 ;; + ldfe FR_QQ_8 = [GR_ad_qq], 16 // Load QQ_8 for 2^-3 < |r| < pi/4 + nop.f 999 + 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) ;; +;; + +{ .mfi + ldfe FR_QQ_7 = [GR_ad_qq], 16 // Load QQ_7 for 2^-3 < |r| < pi/4 + fma.s1 FR_N_float_signif = FR_Input_X, FR_inv_pi_2to63, FR_rshf_2to64 + add GR_ad_se = 0x40, GR_ad_s // Point to end of constant table s } { .mib - nop.m 999 - nop.i 999 -// -// Branch out if x >= 2**63. Use Payne-Hanek Reduction -// -(p7) br.cond.spnt L(SINCOSL_LARGER_ARG) ;; + ldfe FR_PP_7 = [GR_ad_pp], 16 // Load PP_7 for 2^-3 < |r| < pi/4 + mov GR_ad_s1 = GR_ad_s // Save pointer to S_1 +(p10) br.cond.spnt SINCOSL_ARG_TOO_LARGE // Branch if |x| >= 2^63 + // Use Payne-Hanek Reduction } +;; + { .mfi - 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 ;; + ldfe FR_P_3 = [GR_ad_p], 16 // Load P_3 for pi/4 <= |x| < 2^63 + fmerge.se FR_r = FR_norm_x, FR_norm_x // r = x, in case |x| < pi/4 + add GR_ad_m14 = 0x50, GR_ad_s // Point to constant table m14 } -{ .mfi - nop.m 999 -(p0) fcmp.lt.unc.s1 p6, p7 = FR_Input_X, FR_Pi_by_4 - nop.i 999 ;; +{ .mfb + ldfps FR_Two_to_M3, FR_Neg_Two_to_M3 = [GR_ad_d], 8 + fma.s1 FR_rsq = FR_norm_x, FR_norm_x, f0 // rsq = x*x, in case |x| < pi/4 +(p7) br.cond.spnt SINCOSL_LARGER_ARG // Branch if 2^24 <= |x| < 2^63 + // Use pre-reduction +} +;; + +{ .mmf + ldfe FR_PP_6 = [GR_ad_pp], 16 // Load PP_6 for normal path + ldfe FR_QQ_6 = [GR_ad_qq], 16 // Load QQ_6 for normal path + fmerge.se FR_c = f0, f0 // c = 0 in case |x| < pi/4 } +;; + +{ .mmf + ldfe FR_PP_5 = [GR_ad_pp], 16 // Load PP_5 for normal path + ldfe FR_QQ_5 = [GR_ad_qq], 16 // Load QQ_5 for normal path + nop.f 999 +} +;; + +// Here if 0 < |x| < 2^24 { .mfi - 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 ;; + ldfe FR_S_5 = [GR_ad_se], -16 // Load S_5 if i_1=0 + fcmp.lt.s1 p6, p7 = FR_abs_x, FR_Pi_by_4 // Test |x| < pi/4 + nop.i 999 } { .mfi - nop.m 999 + ldfe FR_C_5 = [GR_ad_ce], -16 // Load C_5 if i_1=1 + fms.s1 FR_N_float = FR_N_float_signif, FR_2tom64, FR_rshf + nop.i 999 +} +;; + +{ .mmi + ldfe FR_S_4 = [GR_ad_se], -16 // Load S_4 if i_1=0 + ldfe FR_C_4 = [GR_ad_ce], -16 // Load C_4 if i_1=1 + nop.i 999 +} +;; + // // N = Arg * 2/pi // Check if Arg < pi/4 // -(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 +{ .mfi +(p7) ldfps FR_Two_to_M33, FR_Neg_Two_to_M33 = [GR_ad_d], 8 +(p6) fma.s1 FR_r_cubed = FR_r, FR_rsq, f0 // r^3 if |x| < pi/4 +(p6) mov GR_N_Inc = GR_Sin_or_Cos // N_Inc if |x| < pi/4 +} +;; + +// If |x| < pi/4, r = x and c = 0 // lf |x| < pi/4, is x < 2**(-3). -// r = Arg +// 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 ;; +{ .mmi +(p7) getf.sig GR_N_Inc = FR_N_float_signif +(p6) cmp.lt.unc p8,p0 = GR_exp_x, GR_exp_2_to_m3 // Is |x| < 2^-3 +(p6) tbit.z p9,p10 = GR_N_Inc, 0 // p9 if i_1=0, N mod 4 = 0,1 + // p10 if i_1=1, N mod 4 = 2,3 } -{ .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 +// 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) ;; + + +{ .mbb + nop.m 999 +(p8) br.cond.spnt SINCOSL_SMALL_R_0 // Branch if 0 < |x| < 2^-3 +(p6) br.cond.spnt SINCOSL_NORMAL_R_0 // Branch if 2^-3 <= |x| < pi/4 } -{ .mib - nop.m 999 - nop.i 999 -(p6) br.cond.sptk L(SINCOSL_NORMAL_R) ;; +;; + +// Here if pi/4 <= |x| < 2^24 +{ .mfi + ldfs FR_Neg_Two_to_M67 = [GR_ad_d], 8 // Load -2^-67 + fnma.s1 FR_s = FR_N_float, FR_P_1, FR_Input_X // s = -N * P_1 + Arg + add GR_N_Inc = GR_N_Inc, GR_Sin_or_Cos // Adjust N_Inc for sin/cos } -// -// 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 + nop.m 999 + fma.s1 FR_w = FR_N_float, FR_P_2, f0 // w = N * P_2 + nop.i 999 } +;; + { .mfi - nop.m 999 -(p0) fma.s1 FR_w = FR_N_float, FR_P_2, f0 - nop.i 999 ;; + nop.m 999 + fms.s1 FR_r = FR_s, f1, FR_w // r = s - w, assume |s| >= 2^-33 + tbit.z p9,p10 = GR_N_Inc, 0 // p9 if i_1=0, N mod 4 = 0,1 + // p10 if i_1=1, N mod 4 = 2,3 } +;; + { .mfi - nop.m 999 -// -// 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 ;; + nop.m 999 + fcmp.lt.s1 p7, p6 = FR_s, FR_Two_to_M33 + nop.i 999 } +;; + { .mfi - nop.m 999 -(p7) fcmp.gt.s1 p7, p6 = FR_s, FR_Neg_Two_to_M33 - nop.i 999 ;; + nop.m 999 +(p7) fcmp.gt.s1 p7, p6 = FR_s, FR_Neg_Two_to_M33 // p6 if |s| >= 2^-33, else p7 + nop.i 999 } +;; + { .mfi - nop.m 999 -// 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 + nop.m 999 + fms.s1 FR_c = FR_s, f1, FR_r // c = s - r, for |s| >= 2^-33 + nop.i 999 } { .mfi - nop.m 999 -(p7) fma.s1 FR_w = FR_N_float, FR_P_3, f0 - nop.i 999 ;; + nop.m 999 + fma.s1 FR_rsq = FR_r, FR_r, f0 // rsq = r * r, for |s| >= 2^-33 + nop.i 999 } +;; + { .mfi - nop.m 999 -(p7) fma.s1 FR_U_1 = FR_N_float, FR_P_2, FR_w - nop.i 999 + nop.m 999 +(p7) fma.s1 FR_w = FR_N_float, FR_P_3, f0 + nop.i 999 } +;; + +{ .mmf +(p9) ldfe FR_C_1 = [GR_ad_pp], 16 // Load C_1 if i_1=0 +(p10) ldfe FR_S_1 = [GR_ad_qq], 16 // Load S_1 if i_1=1 + frcpa.s1 FR_r_hi, p15 = f1, FR_r // r_hi = frcpa(r) +} +;; + { .mfi - nop.m 999 -(p6) fms.s1 FR_c = FR_s, f1, FR_r - nop.i 999 ;; + nop.m 999 +(p6) fcmp.lt.unc.s1 p8, p13 = FR_r, FR_Two_to_M3 // If big s, test r with 2^-3 + nop.i 999 } +;; + { .mfi - nop.m 999 -// -// For big s: r = s - w: No futher reduction is necessary + nop.m 999 +(p7) fma.s1 FR_U_1 = FR_N_float, FR_P_2, FR_w + nop.i 999 +} +;; + +// +// For big s: r = s - w: No futher reduction is necessary // For small s: w = N * P_3 (change sign) More reduction // -(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, p13 = FR_r, FR_Neg_Two_to_M3 // If big s, p8 if |r| < 2^-3 + nop.i 999 ;; } + { .mfi - nop.m 999 -(p8) fcmp.gt.s1 p8, p9 = FR_r, FR_Neg_Two_to_M3 - nop.i 999 ;; + nop.m 999 +(p9) fma.s1 FR_poly = FR_rsq, FR_PP_8, FR_PP_7 // poly = rsq*PP_8+PP_7 if i_1=0 + nop.i 999 } { .mfi - nop.m 999 + nop.m 999 +(p10) fma.s1 FR_poly = FR_rsq, FR_QQ_8, FR_QQ_7 // poly = rsq*QQ_8+QQ_7 if i_1=1 + nop.i 999 +} +;; + +{ .mfi + nop.m 999 (p7) fms.s1 FR_r = FR_s, f1, FR_U_1 - nop.i 999 + nop.i 999 } -{ .mfb - nop.m 999 +;; + +{ .mfi + nop.m 999 +(p6) fma.s1 FR_r_cubed = FR_r, FR_rsq, f0 // rcubed = r * rsq + nop.i 999 +} +;; + +{ .mfi // // For big s: Is |r| < 2**(-3)? // For big s: c = S - r @@ -1095,355 +1128,356 @@ L(SINCOSL_CONTINUE): // 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) ;; + nop.m 999 +(p6) fms.s1 FR_c = FR_c, f1, FR_w + nop.i 999 } -{ .mib - nop.m 999 - nop.i 999 -(p9) br.cond.sptk L(SINCOSL_NORMAL_R) ;; +{ .mbb + nop.m 999 +(p8) br.cond.spnt SINCOSL_SMALL_R_1 // Branch if |s|>=2^-33, |r| < 2^-3, + // and pi/4 <= |x| < 2^24 +(p13) br.cond.sptk SINCOSL_NORMAL_R_1 // Branch if |s|>=2^-33, |r| >= 2^-3, + // and pi/4 <= |x| < 2^24 } -{ .mfi -(p7) add GR_Table_Base1 = 224, GR_Table_Base1 +;; + +SINCOSL_S_TINY: // -// Branch to SINCOSL_SMALL_R or SINCOSL_NORMAL_R +// Here if |s| < 2^-33, and pi/4 <= |x| < 2^24 +// +{ .mfi + fms.s1 FR_U_2 = FR_N_float, FR_P_2, FR_U_1 // -(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 + 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 ;; + tbit.z p9,p10 = GR_N_Inc, 0 // p9 if i_1=0, N mod 4 = 0,1 + // p10 if i_1=1, N mod 4 = 2,3 } -// +;; + +// // For small s: U_2 = N * P_2 - U_1 // S_1 stored constant - grab the one stored with the // coefficients. -// +// { .mfi -(p7) ldfe FR_S_1 = [GR_Table_Base1], 16 + ldfe FR_S_1 = [GR_ad_s1], 16 // // Check if i_1 and i_0 != 0 // -(p10) fma.s1 FR_poly = f0, f1, FR_Neg_Two_to_M67 -(p7) cmp.eq.unc p11, p12 = 0x0, GR_i_0 ;; +(p10) fma.s1 FR_poly = f0, f1, FR_Neg_Two_to_M67 + tbit.z p11,p12 = GR_N_Inc, 1 // p11 if i_0=0, N mod 4 = 0,2 + // p12 if i_0=1, N mod 4 = 1,3 } +;; + { .mfi - nop.m 999 -(p7) fms.s1 FR_s = FR_s, f1, FR_r - nop.i 999 + nop.m 999 + fms.s1 FR_s = FR_s, f1, FR_r + nop.i 999 } { .mfi - nop.m 999 -// + 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 ;; + 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 + nop.m 999 + 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 ;; + nop.m 999 + fmerge.se FR_tmp_result = FR_r, FR_r + nop.i 999 ;; } { .mfi - nop.m 999 -(p10) fma.s1 FR_Input_X = f0, f1, f1 - nop.i 999 ;; + nop.m 999 +(p10) fma.s1 FR_tmp_result = f0, f1, f1 + nop.i 999 ;; } { .mfi - nop.m 999 -// + 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 ;; + fms.s1 FR_c = FR_s, f1, FR_U_1 + nop.i 999 ;; } { .mfi - nop.m 999 -// + 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 +(p12) fnma.s1 FR_tmp_result = FR_tmp_result, f1, f0 + nop.i 999 } { .mfi - nop.m 999 -(p7) fma.s1 FR_r = FR_S_1, FR_r, f0 - nop.i 999 ;; + nop.m 999 + 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 ;; + nop.m 999 + fma.s0 FR_S_1 = FR_S_1, FR_S_1, f0 + nop.i 999 ;; } { .mfi - nop.m 999 + nop.m 999 // // If i_1 != 0, poly = 2**(-67) // -(p7) fms.s1 FR_c = FR_c, f1, FR_U_2 - nop.i 999 ;; + fms.s1 FR_c = FR_c, f1, FR_U_2 + nop.i 999 ;; } { .mfi - nop.m 999 -// + nop.m 999 +// // c = c - U_2 -// +// (p9) fma.s1 FR_poly = FR_r, FR_rsq, FR_c - nop.i 999 ;; + nop.i 999 ;; } { .mfi - nop.m 999 + nop.m 999 // // i_0 != 0, so Result = -Result // -(p11) fma.s0 FR_Input_X = FR_Input_X, f1, FR_poly - nop.i 999 ;; +(p11) fma.s0 FR_Result = FR_tmp_result, f1, FR_poly + nop.i 999 ;; } { .mfb - nop.m 999 -(p12) fms.s0 FR_Input_X = FR_Input_X, f1, FR_poly + nop.m 999 +(p12) fms.s0 FR_Result = FR_tmp_result, 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 + br.ret.sptk b0 // Exit if |s| < 2^-33, and pi/4 <= |x| < 2^24 } ;; -// This path for argument > 2*24 -// Adjust table_ptr1 to beginning of table. +SINCOSL_LARGER_ARG: // - -{ .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 +// Here if 2^24 <= |x| < 2^63 +// +{ .mfi + ldfe FR_d_1 = [GR_ad_p], 16 // Load d_1 for |x| >= 2^24 path + fma.s1 FR_N_0 = FR_Input_X, FR_Inv_P_0, f0 + nop.i 999 } ;; - -// -// Point to 2*-14 +// // N_0 = Arg * Inv_P_0 // +// Load values 2**(-14) and -2**(-14) { .mmi -(p0) add GR_Table_Base = 688, GR_Table_Base ;; -(p0) ldfs FR_Two_to_M14 = [GR_Table_Base], 4 - nop.i 999 ;; + ldfps FR_Two_to_M14, FR_Neg_Two_to_M14 = [GR_ad_m14] + nop.i 999 ;; } { .mfi -(p0) ldfs FR_Neg_Two_to_M14 = [GR_Table_Base], 0 - nop.f 999 - nop.i 999 ;; + ldfe FR_d_2 = [GR_ad_p], 16 // Load d_2 for |x| >= 2^24 path + nop.f 999 + nop.i 999 ;; } { .mfi - nop.m 999 + nop.m 999 // -// Load values 2**(-14) and -2**(-14) // -(p0) fcvt.fx.s1 FR_N_0_fix = FR_N_0 - nop.i 999 ;; + fcvt.fx.s1 FR_N_0_fix = FR_N_0 + nop.i 999 ;; } { .mfi - nop.m 999 + nop.m 999 // // N_0_fix = integer part of N_0 // -(p0) fcvt.xf FR_N_0 = FR_N_0_fix - nop.i 999 ;; + fcvt.xf FR_N_0 = FR_N_0_fix + nop.i 999 ;; } { .mfi - nop.m 999 + 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 + 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 ;; + nop.m 999 + fma.s1 FR_w = FR_N_0, FR_d_1, f0 + nop.i 999 ;; } { .mfi - nop.m 999 + 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 ;; + fma.s1 FR_N_float = FR_ArgPrime, FR_Inv_pi_by_2, f0 + nop.i 999 ;; } { .mfi - nop.m 999 + nop.m 999 // -// N = A' * 2/pi +// N = A' * 2/pi // -(p0) fcvt.fx.s1 FR_N_fix = FR_N_float - nop.i 999 ;; + fcvt.fx.s1 FR_N_fix = FR_N_float + nop.i 999 ;; } { .mfi - nop.m 999 + nop.m 999 // -// N_fix is the integer part +// N_fix is the integer part // -(p0) fcvt.xf FR_N_float = FR_N_fix - nop.i 999 ;; + 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 ;; + 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 ;; + nop.m 999 + nop.i 999 ;; + add GR_N_Inc = GR_N_Inc, GR_Sin_or_Cos ;; } { .mfi - nop.m 999 + 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 + 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 ;; + nop.m 999 + fnma.s1 FR_w = FR_N_float, FR_P_2, FR_w + nop.i 999 ;; } { .mfi - nop.m 999 + 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 ;; + 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 ;; + nop.m 999 +(p9) fcmp.gt.s1 p9, p8 = FR_s, FR_Neg_Two_to_M14 // p9 if |s| < 2^-14 + nop.i 999 ;; } + { .mfi - nop.m 999 + 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 + nop.i 999 } { .mfi - nop.m 999 + nop.m 999 (p9) fma.s1 FR_U_hi = FR_N_0, FR_d_1, f0 - nop.i 999 ;; + nop.i 999 ;; } { .mfi - nop.m 999 + 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 + nop.i 999 } { .mfi - nop.m 999 + nop.m 999 (p9) fma.s1 FR_w = FR_N_float, FR_P_3, f0 - nop.i 999 ;; + nop.i 999 ;; } { .mfi - nop.m 999 + 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 + nop.i 999 } { .mfi - nop.m 999 + nop.m 999 // -// Big s: finish up c = (S - r) + w (c complete) +// 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 ;; + nop.i 999 ;; } { .mmf - nop.m 999 - nop.m 999 + 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 + nop.m 999 (p9) fmerge.s FR_V_hiabs = f0, FR_V_hi - nop.i 999 ;; + nop.i 999 ;; } +//{ .mfb +//(p9) fmerge.s f8= FR_V_lo,FR_V_lo +//(p9) br.ret.sptk b0 +//} +//;; { .mfi - nop.m 999 + 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 + nop.i 999 } { .mfi - nop.m 999 + nop.m 999 // -// For big s: Is |r| < 2**(-3) +// 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 ;; +(p8) fms.s1 FR_c = FR_s, f1, FR_r + nop.i 999 ;; } { .mfi - nop.m 999 + nop.m 999 // // For small S: V_hi = N * P_2 // w = N * P_3 @@ -1451,104 +1485,99 @@ L(SINCOSL_LARGER_ARG): // so (-) missing for V_hi and w. // (p8) fcmp.lt.unc.s1 p12, p13 = FR_r, FR_Two_to_M3 - nop.i 999 ;; + nop.i 999 ;; } { .mfi - nop.m 999 + nop.m 999 (p12) fcmp.gt.s1 p12, p13 = FR_r, FR_Neg_Two_to_M3 - nop.i 999 ;; + nop.i 999 ;; } { .mfi - nop.m 999 + nop.m 999 (p8) fma.s1 FR_c = FR_c, f1, FR_w - nop.i 999 + nop.i 999 } { .mfb - nop.m 999 + nop.m 999 (p9) fms.s1 FR_w = FR_N_0, FR_d_2, FR_w -(p12) br.cond.spnt L(SINCOSL_SMALL_R) ;; +(p12) br.cond.spnt SINCOSL_SMALL_R // Branch if |r| < 2^-3 + // and 2^24 <= |x| < 2^63 } +;; + { .mib - nop.m 999 - nop.i 999 -(p13) br.cond.sptk L(SINCOSL_NORMAL_R) ;; + nop.m 999 + nop.i 999 +(p13) br.cond.sptk SINCOSL_NORMAL_R // Branch if |r| >= 2^-3 + // and 2^24 <= |x| < 2^63 } +;; + +SINCOSL_LARGER_S_TINY: +// +// Here if |s| < 2^-14, and 2^24 <= |x| < 2^63 +// { .mfi - nop.m 999 -// -// Big s: Vector off when |r| < 2**(-3). Recall that p8 will be true. + 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 + fcmp.ge.unc.s1 p7, p8 = FR_U_hiabs, FR_V_hiabs } { .mfi - nop.m 999 + 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 ;; + fma.s1 FR_t = FR_U_lo, f1, FR_V_lo } +;; + { .mfi - nop.m 999 + nop.m 999 // -// t = U_lo + V_lo +// t = U_lo + V_lo // // -(p10) fms.s1 FR_a = FR_U_hi, f1, FR_A - nop.i 999 ;; +(p7) 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 + nop.m 999 +(p8) fma.s1 FR_a = FR_V_hi, f1, FR_A + 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 ;; + nop.m 999 + 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 ;; + ldfe FR_C_1 = [GR_ad_c], 16 ;; + ldfe FR_C_2 = [GR_ad_c], 64 + nop.i 999 ;; } // // c = c + C_lo finished. // Load C_2 // { .mfi -(p0) ldfe FR_S_1 = [GR_Table_Base], 16 + ldfe FR_S_1 = [GR_ad_s], 16 // -// C_lo = S - C_hi +// C_lo = S - C_hi // -(p0) fma.s1 FR_t = FR_t, f1, FR_w - nop.i 999 ;; + fma.s1 FR_t = FR_t, f1, FR_w + nop.i 999 ;; } // // r and c have been computed. @@ -1558,855 +1587,695 @@ L(SINCOSL_LARGER_ARG): // Load S_1 // { .mfi -(p0) ldfe FR_S_2 = [GR_Table_Base], 64 + ldfe FR_S_2 = [GR_ad_s], 64 // -// t = t + w +// t = t + w // -(p10) fms.s1 FR_a = FR_a, f1, FR_V_hi -(p0) cmp.eq.unc p9, p10 = 0x0, GR_i_0 ;; +(p7) fms.s1 FR_a = FR_a, f1, FR_V_hi + tbit.z p9,p10 = GR_N_Inc, 0 // p9 if i_1=0, N mod 4 = 0,1 + // p10 if i_1=1, N mod 4 = 2,3 } +;; { .mfi - nop.m 999 + 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 ;; + 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 ;; + nop.m 999 +(p8) fms.s1 FR_a = FR_U_hi, f1, FR_a + tbit.z p11,p12 = GR_N_Inc, 1 // p11 if i_0=0, N mod 4 = 0,2 + // p12 if i_0=1, N mod 4 = 1,3 } +;; + { .mfi - nop.m 999 + 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 ;; + fma.s1 FR_C_lo = FR_C_lo, f1, FR_A + nop.i 999 ;; } { .mfi - nop.m 999 + nop.m 999 // -// C_lo = (S - C_hi) + A +// C_lo = (S - C_hi) + A // -(p0) fma.s1 FR_t = FR_t, f1, FR_a - nop.i 999 ;; + fma.s1 FR_t = FR_t, f1, FR_a + nop.i 999 ;; } { .mfi - nop.m 999 + nop.m 999 // -// t = t + a +// t = t + a // -(p0) fma.s1 FR_C_lo = FR_C_lo, f1, FR_t - nop.i 999 ;; + fma.s1 FR_C_lo = FR_C_lo, f1, FR_t + nop.i 999 ;; } { .mfi - nop.m 999 + 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 ;; + fma.s1 FR_r = FR_C_hi, f1, FR_C_lo + nop.i 999 ;; } { .mfi - nop.m 999 + nop.m 999 // // Load S_2 // -(p0) fma.s1 FR_rsq = FR_r, FR_r, f0 - nop.i 999 + fma.s1 FR_rsq = FR_r, FR_r, f0 + nop.i 999 } { .mfi - nop.m 999 + 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 ;; + fms.s1 FR_c = FR_C_hi, f1, FR_r + nop.i 999 ;; } { .mfi - nop.m 999 + 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 ;; +(p9) fma.s1 FR_tmp_result = f0, f1, FR_r + nop.i 999 ;; } { .mfi - nop.m 999 -(p12) fma.s1 FR_Input_X = f0, f1, f1 - nop.i 999 ;; + nop.m 999 +(p10) fma.s1 FR_tmp_result = f0, f1, f1 + nop.i 999 ;; } { .mfi - nop.m 999 + nop.m 999 // -// Compute r_cube = FR_rsq * r +// Compute r_cube = FR_rsq * r // -(p11) fma.s1 FR_poly = FR_rsq, FR_S_2, FR_S_1 - nop.i 999 ;; +(p9) 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 + nop.m 999 +(p10) fma.s1 FR_poly = FR_rsq, FR_C_2, FR_C_1 + nop.i 999 } { .mfi - nop.m 999 + 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 ;; + fma.s1 FR_r_cubed = FR_rsq, FR_r, f0 + nop.i 999 ;; } { .mfi - nop.m 999 + nop.m 999 // // c = C_hi - r // Load C_1 // -(p0) fma.s1 FR_c = FR_c, f1, FR_C_lo - nop.i 999 + fma.s1 FR_c = FR_c, f1, FR_C_lo + nop.i 999 } { .mfi - nop.m 999 + 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 ;; +(p12) fms.s1 FR_tmp_result = f0, f1, FR_tmp_result + nop.i 999 ;; } { .mfi - nop.m 999 + 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 ;; +(p9) 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 ;; + nop.m 999 +(p10) fma.s1 FR_poly = FR_rsq, FR_poly, f0 + nop.i 999 ;; } { .mfi - nop.m 999 + nop.m 999 // -// if i_0 !=0: Result = -Result +// if i_0 !=0: Result = -Result // -(p9) fma.s0 FR_Input_X = FR_Input_X, f1, FR_poly - nop.i 999 ;; +(p11) fma.s0 FR_Result = FR_tmp_result, f1, FR_poly + nop.i 999 ;; } { .mfb - nop.m 999 -(p10) fms.s0 FR_Input_X = FR_Input_X, f1, FR_poly + nop.m 999 +(p12) fms.s0 FR_Result = FR_tmp_result, f1, FR_poly // // if i_0 == 0: Result = Result + poly // else Result = Result - poly // -(p0) br.ret.sptk b0 ;; + br.ret.sptk b0 // Exit for |s| < 2^-14, and 2^24 <= |x| < 2^63 } -L(SINCOSL_SMALL_R): -{ .mii - nop.m 999 -(p0) extr.u GR_i_1 = GR_N_Inc, 0, 1 ;; +;; + + +SINCOSL_SMALL_R: +// +// Here if |r| < 2^-3 // +// Enter with r, c, and N_Inc computed // // 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 + nop.m 999 + fma.s1 FR_rsq = FR_r, FR_r, f0 // rsq = r * r + tbit.z p9,p10 = GR_N_Inc, 0 // p9 if i_1=0, N mod 4 = 0,1 + // p10 if i_1=1, N mod 4 = 2,3 } ;; -// ****************************************************************** -// ****************************************************************** -// ****************************************************************** -// 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 +(p9) ldfe FR_S_5 = [GR_ad_se], -16 // Load S_5 if i_1=0 +(p10) ldfe FR_C_5 = [GR_ad_ce], -16 // Load C_5 if i_1=1 nop.i 999 } ;; { .mmi - ld8 GR_Table_Base = [GR_Table_Base] - nop.m 999 +(p9) ldfe FR_S_4 = [GR_ad_se], -16 // Load S_4 if i_1=0 +(p10) ldfe FR_C_4 = [GR_ad_ce], -16 // Load C_4 if i_1=1 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 ;; +SINCOSL_SMALL_R_0: +// Entry point for 2^-3 < |x| < pi/4 +.pred.rel "mutex",p9,p10 +SINCOSL_SMALL_R_1: +// Entry point for pi/4 < |x| < 2^24 and |r| < 2^-3 +.pred.rel "mutex",p9,p10 +{ .mfi +(p9) ldfe FR_S_3 = [GR_ad_se], -16 // Load S_3 if i_1=0 + fma.s1 FR_Z = FR_rsq, FR_rsq, f0 // Z = rsq * rsq + nop.i 999 } -// -// 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 ;; +{ .mfi +(p10) ldfe FR_C_3 = [GR_ad_ce], -16 // Load C_3 if i_1=1 +(p10) fnma.s1 FR_c = FR_c, FR_r, f0 // c = -c * r if i_1=0 + 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 ;; +(p9) ldfe FR_S_2 = [GR_ad_se], -16 // Load S_2 if i_1=0 +(p10) ldfe FR_C_2 = [GR_ad_ce], -16 // Load C_2 if i_1=1 +(p10) fmerge.s FR_r = f1, f1 } +;; + { .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 ;; +(p9) ldfe FR_S_1 = [GR_ad_se], -16 // Load S_1 if i_1=0 +(p10) ldfe FR_C_1 = [GR_ad_ce], -16 // Load C_1 if i_1=1 + 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 + nop.m 999 +(p9) fma.s1 FR_Z = FR_Z, FR_r, f0 // Z = Z * r if i_1=0 + 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 ;; + nop.m 999 +(p9) fma.s1 FR_poly_lo = FR_rsq, FR_S_5, FR_S_4 // poly_lo=rsq*S_5+S_4 if i_1=0 + 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 ;; + nop.m 999 +(p10) fma.s1 FR_poly_lo = FR_rsq, FR_C_5, FR_C_4 // poly_lo=rsq*C_5+C_4 if i_1=1 + nop.i 999 } +;; + { .mfi - nop.m 999 -(p9) fma.s1 FR_poly_hi = FR_rsq, FR_S_2, FR_S_1 - nop.i 999 + nop.m 999 +(p9) fma.s1 FR_poly_hi = FR_rsq, FR_S_2, FR_S_1 // poly_hi=rsq*S_2+S_1 if i_1=0 + 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 ;; + nop.m 999 +(p10) fma.s1 FR_poly_hi = FR_rsq, FR_C_2, FR_C_1 // poly_hi=rsq*C_2+C_1 if i_1=1 + 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 ;; + nop.m 999 + fma.s1 FR_Z = FR_Z, FR_rsq, f0 // Z = Z * rsq + nop.i 999 } +;; + { .mfi - nop.m 999 -(p10) fma.s1 FR_poly_hi = FR_rsq, FR_C_2, FR_C_1 - nop.i 999 + nop.m 999 +(p9) fma.s1 FR_poly_lo = FR_rsq, FR_poly_lo, FR_S_3 // p_lo=p_lo*rsq+S_3, i_1=0 + nop.i 999 } { .mfi - nop.m 999 -(p10) fma.s0 FR_C_1 = FR_C_1, FR_C_1, f0 - nop.i 999 ;; + nop.m 999 +(p10) fma.s1 FR_poly_lo = FR_rsq, FR_poly_lo, FR_C_3 // p_lo=p_lo*rsq+C_3, i_1=1 + nop.i 999 } +;; + { .mfi - nop.m 999 -(p9) fma.s1 FR_poly_hi = FR_poly_hi, FR_rsq, f0 - nop.i 999 + nop.m 999 +(p9) fma.s0 FR_inexact = FR_S_4, FR_S_4, f0 // Dummy op to set inexact + tbit.z p11,p12 = GR_N_Inc, 1 // p11 if i_0=0, N mod 4 = 0,2 + // p12 if i_0=1, N mod 4 = 1,3 } { .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 ;; + nop.m 999 +(p10) fma.s0 FR_inexact = FR_C_1, FR_C_1, f0 // Dummy op to set inexact + nop.i 999 } +;; + { .mfi - nop.m 999 -(p10) fma.s1 FR_poly_hi = FR_poly_hi, FR_rsq, f0 - nop.i 999 ;; + nop.m 999 +(p9) fma.s1 FR_poly_hi = FR_poly_hi, FR_rsq, f0 // p_hi=p_hi*rsq if i_1=0 + 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 + nop.m 999 +(p10) fma.s1 FR_poly_hi = FR_poly_hi, FR_rsq, f0 // p_hi=p_hi*rsq if i_1=1 + 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 ;; + nop.m 999 + fma.s1 FR_poly = FR_Z, FR_poly_lo, FR_c // poly=Z*poly_lo+c + nop.i 999 } +;; + { .mfi - nop.m 999 -(p12) fms.s1 FR_r = f0, f1, FR_r - nop.i 999 ;; + nop.m 999 +(p9) fma.s1 FR_poly_hi = FR_r, FR_poly_hi, f0 // p_hi=r*p_hi if i_1=0 + 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 ;; + nop.m 999 +(p12) fms.s1 FR_r = f0, f1, FR_r // r = -r if i_0=1 + nop.i 999 } +;; + { .mfi - nop.m 999 -(p12) fms.s0 FR_Input_X = FR_r, f1, FR_poly - nop.i 999 + nop.m 999 + fma.s1 FR_poly = FR_poly, f1, FR_poly_hi // poly=poly+poly_hi + 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 ;; + nop.m 999 +(p11) fma.s0 FR_Result = FR_r, f1, FR_poly + nop.i 999 } -{ .mfi - nop.m 999 -(p0) frcpa.s1 FR_r_hi, p6 = f1, FR_r -(p0) cmp.eq.unc p11, p12 = 0x0, GR_i_0 +{ .mfb + nop.m 999 +(p12) fms.s0 FR_Result = FR_r, f1, FR_poly + br.ret.sptk b0 // Exit for |r| < 2^-3 } ;; -// ****************************************************************** -// ****************************************************************** -// ****************************************************************** + +SINCOSL_NORMAL_R: // -// 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. +// Here if 2^-3 <= |r| < pi/4 +// THIS IS THE MAIN PATH // - -{ .mmi - nop.m 999 -(p0) addl GR_Table_Base = @ltoff(FSINCOSL_CONSTANTS#), gp +// Enter with r, c, and N_Inc having been computed +// +{ .mfi + ldfe FR_PP_6 = [GR_ad_pp], 16 // Load PP_6 + fma.s1 FR_rsq = FR_r, FR_r, f0 // rsq = r * r + tbit.z p9,p10 = GR_N_Inc, 0 // p9 if i_1=0, N mod 4 = 0,1 + // p10 if i_1=1, N mod 4 = 2,3 +} +{ .mfi + ldfe FR_QQ_6 = [GR_ad_qq], 16 // Load QQ_6 + nop.f 999 nop.i 999 } ;; { .mmi - ld8 GR_Table_Base = [GR_Table_Base] - nop.m 999 +(p9) ldfe FR_PP_5 = [GR_ad_pp], 16 // Load PP_5 if i_1=0 +(p10) ldfe FR_QQ_5 = [GR_ad_qq], 16 // Load QQ_5 if i_1=1 nop.i 999 } ;; +SINCOSL_NORMAL_R_0: +// Entry for 2^-3 < |x| < pi/4 +.pred.rel "mutex",p9,p10 +{ .mmf +(p9) ldfe FR_C_1 = [GR_ad_pp], 16 // Load C_1 if i_1=0 +(p10) ldfe FR_S_1 = [GR_ad_qq], 16 // Load S_1 if i_1=1 + frcpa.s1 FR_r_hi, p6 = f1, FR_r // r_hi = frcpa(r) +} +;; { .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 ;; + nop.m 999 +(p9) fma.s1 FR_poly = FR_rsq, FR_PP_8, FR_PP_7 // poly = rsq*PP_8+PP_7 if i_1=0 + nop.i 999 } { .mfi -(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 ;; + nop.m 999 +(p10) fma.s1 FR_poly = FR_rsq, FR_QQ_8, FR_QQ_7 // poly = rsq*QQ_8+QQ_7 if i_1=1 + nop.i 999 } +;; + { .mfi - nop.m 999 -(p0) fma.s1 FR_r_cubed = FR_r, FR_rsq, f0 - nop.i 999 ;; + nop.m 999 + fma.s1 FR_r_cubed = FR_r, FR_rsq, f0 // rcubed = r * rsq + nop.i 999 } +;; + + +SINCOSL_NORMAL_R_1: +// Entry for pi/4 <= |x| < 2^24 +.pred.rel "mutex",p9,p10 { .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 +(p9) ldfe FR_PP_1 = [GR_ad_pp], 16 // Load PP_1_hi if i_1=0 +(p10) ldfe FR_QQ_1 = [GR_ad_qq], 16 // Load QQ_1 if i_1=1 + frcpa.s1 FR_r_hi, p6 = f1, FR_r_hi // r_hi = frpca(frcpa(r)) } +;; + { .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 ;; +(p9) ldfe FR_PP_4 = [GR_ad_pp], 16 // Load PP_4 if i_1=0 +(p9) fma.s1 FR_poly = FR_rsq, FR_poly, FR_PP_6 // poly = rsq*poly+PP_6 if i_1=0 + 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 ;; +(p10) ldfe FR_QQ_4 = [GR_ad_qq], 16 // Load QQ_4 if i_1=1 +(p10) fma.s1 FR_poly = FR_rsq, FR_poly, FR_QQ_6 // poly = rsq*poly+QQ_6 if i_1=1 + nop.i 999 } +;; + { .mfi - nop.m 999 -(p10) fma.s1 FR_corr = FR_S_1, FR_r_cubed, FR_r - nop.i 999 + nop.m 999 +(p9) fma.s1 FR_corr = FR_C_1, FR_rsq, f0 // corr = C_1 * rsq if i_1=0 + nop.i 999 } { .mfi - nop.m 999 -(p10) fma.s1 FR_poly = FR_rsq, FR_poly, FR_QQ_6 - nop.i 999 ;; + nop.m 999 +(p10) fma.s1 FR_corr = FR_S_1, FR_r_cubed, FR_r // corr = S_1 * r^3 + r if i_1=1 + 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 ;; +(p9) ldfe FR_PP_3 = [GR_ad_pp], 16 // Load PP_3 if i_1=0 + fma.s1 FR_r_hi_sq = FR_r_hi, FR_r_hi, f0 // r_hi_sq = r_hi * r_hi + 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 ;; +(p10) ldfe FR_QQ_3 = [GR_ad_qq], 16 // Load QQ_3 if i_1=1 + fms.s1 FR_r_lo = FR_r, f1, FR_r_hi // r_lo = r - r_hi + nop.i 999 } +;; + { .mfi - nop.m 999 -(p10) fma.s1 FR_U_lo = FR_r_hi, f1, FR_r - nop.i 999 +(p9) ldfe FR_PP_2 = [GR_ad_pp], 16 // Load PP_2 if i_1=0 +(p9) fma.s1 FR_poly = FR_rsq, FR_poly, FR_PP_5 // poly = rsq*poly+PP_5 if i_1=0 + 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 ;; +(p10) ldfe FR_QQ_2 = [GR_ad_qq], 16 // Load QQ_2 if i_1=1 +(p10) fma.s1 FR_poly = FR_rsq, FR_poly, FR_QQ_5 // poly = rsq*poly+QQ_5 if i_1=1 + 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 +(p9) ldfe FR_PP_1_lo = [GR_ad_pp], 16 // Load PP_1_lo if i_1=0 +(p9) fma.s1 FR_corr = FR_corr, FR_c, FR_c // corr = corr * c + c if i_1=0 + 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 ;; + nop.m 999 +(p10) fnma.s1 FR_corr = FR_corr, FR_c, f0 // corr = -corr * c if i_1=1 + nop.i 999 } +;; + { .mfi - nop.m 999 -(p10) fma.s1 FR_poly = FR_rsq, FR_poly, FR_QQ_5 - nop.i 999 + nop.m 999 +(p9) fma.s1 FR_U_lo = FR_r, FR_r_hi, FR_r_hi_sq // U_lo = r*r_hi+r_hi_sq, i_1=0 + nop.i 999 } { .mfi - nop.m 999 -(p10) fnma.s1 FR_corr = FR_corr, FR_c, f0 - nop.i 999 ;; + nop.m 999 +(p10) fma.s1 FR_U_lo = FR_r_hi, f1, FR_r // U_lo = r_hi + r if i_1=1 + 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 ;; + nop.m 999 +(p9) fma.s1 FR_U_hi = FR_r_hi, FR_r_hi_sq, f0 // U_hi = r_hi*r_hi_sq if i_1=0 + nop.i 999 } { .mfi - nop.m 999 -(p9) fma.s1 FR_poly = FR_rsq, FR_poly, FR_PP_5 - nop.i 999 + nop.m 999 +(p10) fma.s1 FR_U_hi = FR_QQ_1, FR_r_hi_sq, f1 // U_hi = QQ_1*r_hi_sq+1, i_1=1 + 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 ;; + nop.m 999 +(p9) fma.s1 FR_poly = FR_rsq, FR_poly, FR_PP_4 // poly = poly*rsq+PP_4 if i_1=0 + nop.i 999 } { .mfi - nop.m 999 -(p9) fma.s1 FR_U_lo = FR_r, FR_r, FR_U_lo - nop.i 999 ;; + nop.m 999 +(p10) fma.s1 FR_poly = FR_rsq, FR_poly, FR_QQ_4 // poly = poly*rsq+QQ_4 if i_1=1 + nop.i 999 } +;; + { .mfi - nop.m 999 -(p10) fma.s1 FR_U_lo = FR_r_lo, FR_U_lo, f0 - nop.i 999 ;; + nop.m 999 +(p9) fma.s1 FR_U_lo = FR_r, FR_r, FR_U_lo // U_lo = r * r + U_lo if i_1=0 + 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 + nop.m 999 +(p10) fma.s1 FR_U_lo = FR_r_lo, FR_U_lo, f0 // U_lo = r_lo * U_lo if i_1=1 + nop.i 999 } +;; + { .mfi - nop.m 999 -(p10) fma.s1 FR_poly = FR_rsq, FR_poly, FR_QQ_4 - nop.i 999 ;; + nop.m 999 +(p9) fma.s1 FR_U_hi = FR_PP_1, FR_U_hi, f0 // U_hi = PP_1 * U_hi if i_1=0 + nop.i 999 } +;; + { .mfi - nop.m 999 -(p10) fma.s1 FR_U_hi = FR_QQ_1, FR_U_hi, f1 - nop.i 999 + nop.m 999 +(p9) fma.s1 FR_poly = FR_rsq, FR_poly, FR_PP_3 // poly = poly*rsq+PP_3 if i_1=0 + nop.i 999 } { .mfi - nop.m 999 -(p10) fma.s0 FR_QQ_5 = FR_QQ_5, FR_QQ_5, f0 - nop.i 999 ;; + nop.m 999 +(p10) fma.s1 FR_poly = FR_rsq, FR_poly, FR_QQ_3 // poly = poly*rsq+QQ_3 if i_1=1 + 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 ;; + nop.m 999 +(p9) fma.s1 FR_U_lo = FR_r_lo, FR_U_lo, f0 // U_lo = r_lo * U_lo if i_1=0 + nop.i 999 } { .mfi - nop.m 999 -(p9) fma.s1 FR_U_lo = FR_r_lo, FR_U_lo, f0 - nop.i 999 ;; + nop.m 999 +(p10) fma.s1 FR_U_lo = FR_QQ_1,FR_U_lo, f0 // U_lo = QQ_1 * U_lo if i_1=1 + nop.i 999 } +;; + { .mfi - nop.m 999 -(p10) fma.s1 FR_U_lo = FR_QQ_1,FR_U_lo, f0 - nop.i 999 ;; + nop.m 999 +(p9) fma.s1 FR_U_hi = FR_r, f1, FR_U_hi // U_hi = r + U_hi if i_1=0 + nop.i 999 } +;; + { .mfi - nop.m 999 -(p9) fma.s1 FR_U_hi = FR_PP_1, FR_U_hi, f0 - nop.i 999 + nop.m 999 +(p9) fma.s1 FR_poly = FR_rsq, FR_poly, FR_PP_2 // poly = poly*rsq+PP_2 if i_1=0 + nop.i 999 } { .mfi - nop.m 999 -(p10) fma.s1 FR_poly = FR_rsq, FR_poly, FR_QQ_3 - nop.i 999 ;; + nop.m 999 +(p10) fma.s1 FR_poly = FR_rsq, FR_poly, FR_QQ_2 // poly = poly*rsq+QQ_2 if i_1=1 + 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 ;; + nop.m 999 +(p9) fma.s1 FR_U_lo = FR_PP_1, FR_U_lo, f0 // U_lo = PP_1 * U_lo if i_1=0 + 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 ;; + nop.m 999 +(p9) fma.s1 FR_poly = FR_rsq, FR_poly, FR_PP_1_lo // poly =poly*rsq+PP1lo i_1=0 + 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 + nop.m 999 +(p10) fma.s1 FR_poly = FR_rsq, FR_poly, f0 // poly = poly*rsq if i_1=1 + nop.i 999 } +;; + { .mfi - nop.m 999 -(p10) fma.s1 FR_poly = FR_rsq, FR_poly, FR_QQ_2 - nop.i 999 ;; + nop.m 999 + fma.s1 FR_V = FR_U_lo, f1, FR_corr // V = U_lo + corr + tbit.z p11,p12 = GR_N_Inc, 1 // p11 if i_0=0, N mod 4 = 0,2 + // p12 if i_0=1, N mod 4 = 1,3 } +;; + { .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 ;; + nop.m 999 +(p9) fma.s0 FR_inexact = FR_PP_5, FR_PP_4, f0 // Dummy op to set inexact + 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 ;; + nop.m 999 +(p10) fma.s0 FR_inexact = FR_QQ_5, FR_QQ_5, f0 // Dummy op to set inexact + nop.i 999 } +;; + { .mfi - nop.m 999 -(p10) fma.s1 FR_poly = FR_rsq, FR_poly, f0 - nop.i 999 ;; + nop.m 999 +(p9) fma.s1 FR_poly = FR_r_cubed, FR_poly, f0 // poly = poly*r^3 if i_1=0 + 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 ;; + nop.m 999 +(p10) fma.s1 FR_poly = FR_rsq, FR_poly, f0 // poly = poly*rsq if i_1=1 + nop.i 999 } +;; + { .mfi - nop.m 999 -(p10) fma.s1 FR_poly = FR_rsq, FR_poly, f0 - nop.i 999 ;; + nop.m 999 +(p11) fma.s1 FR_tmp_result = f0, f1, f1// tmp_result=+1.0 if i_0=0 + 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 ;; + nop.m 999 +(p12) fms.s1 FR_tmp_result = f0, f1, f1// tmp_result=-1.0 if i_0=1 + 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 ;; + nop.m 999 + fma.s1 FR_V = FR_poly, f1, FR_V // V = poly + V + nop.i 999 } +;; + +// If i_0 = 0 Result = U_hi + V +// If i_0 = 1 Result = -U_hi - V { .mfi - nop.m 999 -(p12) fms.s0 FR_Input_X = FR_Input_X, FR_U_hi, FR_V - nop.i 999 + nop.m 999 +(p11) fma.s0 FR_Result = FR_tmp_result, 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 -// + nop.m 999 +(p12) fms.s0 FR_Result = FR_tmp_result, FR_U_hi, FR_V + br.ret.sptk b0 // Exit for 2^-3 <= |r| < pi/4 +} +;; -L(SINCOSL_ZERO): -{ .mbb -(p0) cmp.eq.unc p6, p7 = 0x1, GR_Sin_or_Cos - nop.b 999 - nop.b 999 ;; +SINCOSL_ZERO: +// Here if x = 0 +{ .mfi + cmp.eq.unc p6, p7 = 0x1, GR_Sin_or_Cos + nop.f 999 + nop.i 999 } +;; + { .mfi - nop.m 999 -(p7) fmerge.s FR_Input_X = FR_Input_X, FR_Input_X - nop.i 999 + nop.m 999 +(p7) fmerge.s FR_Result = FR_Input_X, FR_Input_X // If sin, result = input + nop.i 999 } { .mfb - nop.m 999 -(p6) fmerge.s FR_Input_X = f1, f1 -(p0) br.ret.sptk b0 ;; + nop.m 999 +(p6) fma.s0 FR_Result = f1, f1, f0 // If cos, result=1.0 + br.ret.sptk b0 // Exit for x=0 +} +;; + + +SINCOSL_DENORMAL: +{ .mmb + getf.exp GR_signexp_x = FR_norm_x // Get sign and exponent of x + nop.m 999 + br.cond.sptk SINCOSL_COMMON // Return to common code } -L(SINCOSL_SPECIAL): +;; + +SINCOSL_SPECIAL: { .mfb nop.m 999 // @@ -2414,106 +2283,83 @@ L(SINCOSL_SPECIAL): // Invalid can be raised. SNaNs // become QNaNs // -(p0) fmpy.s0 FR_Input_X = FR_Input_X, f0 -(p0) br.ret.sptk b0 ;; + fmpy.s0 FR_Result = FR_Input_X, f0 + 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 -// +GLOBAL_IEEE754_END(cosl) -.proc __libm_callout -__libm_callout: -L(SINCOSL_ARG_TOO_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 r59-60. 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) +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 + setf.exp FR_Two_to_M3 = GR_exp_2_to_m3 // Form 2^-3 + 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 - 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 + setf.exp FR_Neg_Two_to_M3 = GR_exp_m2_to_m3 // Form -(2^-3) nop.i 0 -(p0) br.call.sptk b0=__libm_pi_by_2_reduce# ;; + 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 + 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 };; -{ .mib + +{ .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 - 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) + +{ .mbb + nop.m 999 +(p6) br.cond.spnt SINCOSL_SMALL_R // Branch if |r|< 2^-3 for |x| >= 2^63 + br.cond.sptk SINCOSL_NORMAL_R // Branch if |r|>=2^-3 for |x| >= 2^63 +};; + +LOCAL_LIBM_END(__libm_callout) .type __libm_pi_by_2_reduce#,@function .global __libm_pi_by_2_reduce# |