| .file "libm_sincosl.s" |
| |
| |
| // Copyright (c) 2000 - 2004, Intel Corporation |
| // All rights reserved. |
| // |
| // 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 |
| // 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: |
| // 05/13/02 Initial version of sincosl (based on libm's sinl and cosl) |
| // 02/10/03 Reordered header: .section, .global, .proc, .align; |
| // used data8 for long double table values |
| // 10/13/03 Corrected .file name |
| // 02/11/04 cisl is moved to the separate file. |
| // 10/26/04 Avoided using r14-31 as scratch so not clobbered by dynamic loader |
| // |
| //********************************************************************* |
| // |
| // Function: Combined sincosl routine with 3 different API's |
| // |
| // API's |
| //============================================================== |
| // 1) void sincosl(long double, long double*s, long double*c) |
| // 2) __libm_sincosl - internal LIBM function, that accepts |
| // argument in f8 and returns cosine through f8, sine through f9 |
| // |
| // |
| //********************************************************************* |
| // |
| // Resources Used: |
| // |
| // Floating-Point Registers: f8 (Input x and cosl return value), |
| // f9 (sinl returned) |
| // f32-f121 |
| // |
| // General Purpose Registers: |
| // r32-r61 |
| // |
| // Predicate Registers: p6-p15 |
| // |
| //********************************************************************* |
| // |
| // IEEE Special Conditions: |
| // |
| // Denormal fault raised on denormal inputs |
| // Overflow exceptions do not occur |
| // Underflow exceptions raised when appropriate for sincosl |
| // (No specialized error handling for this routine) |
| // Inexact raised when appropriate by algorithm |
| // |
| // sincosl(SNaN) = QNaN, QNaN |
| // sincosl(QNaN) = QNaN, QNaN |
| // sincosl(inf) = QNaN, QNaN |
| // sincosl(+/-0) = +/-0, 1 |
| // |
| //********************************************************************* |
| // |
| // Mathematical Description |
| // ======================== |
| // |
| // The computation of FSIN and FCOS performed in parallel. |
| // |
| // Arg = N pi/2 + alpha, |alpha| <= pi/4. |
| // |
| // cosl( Arg ) = sinl( (N+1) pi/2 + alpha ), |
| // |
| // therefore, the code for computing sine will produce cosine as long |
| // as 1 is added to N immediately after the argument reduction |
| // process. |
| // |
| // Let M = N if sine |
| // N+1 if cosine. |
| // |
| // Now, given |
| // |
| // Arg = M pi/2 + alpha, |alpha| <= pi/4, |
| // |
| // let I = M mod 4, or I be the two lsb of M when M is represented |
| // as 2's complement. I = [i_0 i_1]. Then |
| // |
| // sinl( Arg ) = (-1)^i_0 sinl( alpha ) if i_1 = 0, |
| // = (-1)^i_0 cosl( alpha ) if i_1 = 1. |
| // |
| // For example: |
| // if M = -1, I = 11 |
| // sin ((-pi/2 + alpha) = (-1) cos (alpha) |
| // if M = 0, I = 00 |
| // sin (alpha) = sin (alpha) |
| // if M = 1, I = 01 |
| // sin (pi/2 + alpha) = cos (alpha) |
| // if M = 2, I = 10 |
| // sin (pi + alpha) = (-1) sin (alpha) |
| // if M = 3, I = 11 |
| // sin ((3/2)pi + alpha) = (-1) cos (alpha) |
| // |
| // The value of alpha is obtained by argument reduction and |
| // represented by two working precision numbers r and c where |
| // |
| // alpha = r + c accurately. |
| // |
| // The reduction method is described in a previous write up. |
| // The argument reduction scheme identifies 4 cases. For Cases 2 |
| // and 4, because |alpha| is small, sinl(r+c) and cosl(r+c) can be |
| // computed very easily by 2 or 3 terms of the Taylor series |
| // expansion as follows: |
| // |
| // Case 2: |
| // ------- |
| // |
| // sinl(r + c) = r + c - r^3/6 accurately |
| // cosl(r + c) = 1 - 2^(-67) accurately |
| // |
| // Case 4: |
| // ------- |
| // |
| // sinl(r + c) = r + c - r^3/6 + r^5/120 accurately |
| // cosl(r + c) = 1 - r^2/2 + r^4/24 accurately |
| // |
| // The only cases left are Cases 1 and 3 of the argument reduction |
| // procedure. These two cases will be merged since after the |
| // argument is reduced in either cases, we have the reduced argument |
| // represented as r + c and that the magnitude |r + c| is not small |
| // enough to allow the usage of a very short approximation. |
| // |
| // The required calculation is either |
| // |
| // sinl(r + c) = sinl(r) + correction, or |
| // cosl(r + c) = cosl(r) + correction. |
| // |
| // Specifically, |
| // |
| // sinl(r + c) = sinl(r) + c sin'(r) + O(c^2) |
| // = sinl(r) + c cos (r) + O(c^2) |
| // = sinl(r) + c(1 - r^2/2) accurately. |
| // Similarly, |
| // |
| // cosl(r + c) = cosl(r) - c sinl(r) + O(c^2) |
| // = cosl(r) - c(r - r^3/6) accurately. |
| // |
| // We therefore concentrate on accurately calculating sinl(r) and |
| // cosl(r) for a working-precision number r, |r| <= pi/4 to within |
| // 0.1% or so. |
| // |
| // The greatest challenge of this task is that the second terms of |
| // the Taylor series |
| // |
| // r - r^3/3! + r^r/5! - ... |
| // |
| // and |
| // |
| // 1 - r^2/2! + r^4/4! - ... |
| // |
| // are not very small when |r| is close to pi/4 and the rounding |
| // errors will be a concern if simple polynomial accumulation is |
| // used. When |r| < 2^-3, however, the second terms will be small |
| // enough (6 bits or so of right shift) that a normal Horner |
| // recurrence suffices. Hence there are two cases that we consider |
| // in the accurate computation of sinl(r) and cosl(r), |r| <= pi/4. |
| // |
| // Case small_r: |r| < 2^(-3) |
| // -------------------------- |
| // |
| // Since Arg = M pi/4 + r + c accurately, and M mod 4 is [i_0 i_1], |
| // we have |
| // |
| // sinl(Arg) = (-1)^i_0 * sinl(r + c) if i_1 = 0 |
| // = (-1)^i_0 * cosl(r + c) if i_1 = 1 |
| // |
| // can be accurately approximated by |
| // |
| // sinl(Arg) = (-1)^i_0 * [sinl(r) + c] if i_1 = 0 |
| // = (-1)^i_0 * [cosl(r) - c*r] if i_1 = 1 |
| // |
| // because |r| is small and thus the second terms in the correction |
| // are unneccessary. |
| // |
| // Finally, sinl(r) and cosl(r) are approximated by polynomials of |
| // moderate lengths. |
| // |
| // sinl(r) = r + S_1 r^3 + S_2 r^5 + ... + S_5 r^11 |
| // cosl(r) = 1 + C_1 r^2 + C_2 r^4 + ... + C_5 r^10 |
| // |
| // We can make use of predicates to selectively calculate |
| // sinl(r) or cosl(r) based on i_1. |
| // |
| // Case normal_r: 2^(-3) <= |r| <= pi/4 |
| // ------------------------------------ |
| // |
| // This case is more likely than the previous one if one considers |
| // r to be uniformly distributed in [-pi/4 pi/4]. Again, |
| // |
| // sinl(Arg) = (-1)^i_0 * sinl(r + c) if i_1 = 0 |
| // = (-1)^i_0 * cosl(r + c) if i_1 = 1. |
| // |
| // Because |r| is now larger, we need one extra term in the |
| // correction. sinl(Arg) can be accurately approximated by |
| // |
| // sinl(Arg) = (-1)^i_0 * [sinl(r) + c(1-r^2/2)] if i_1 = 0 |
| // = (-1)^i_0 * [cosl(r) - c*r*(1 - r^2/6)] i_1 = 1. |
| // |
| // Finally, sinl(r) and cosl(r) are approximated by polynomials of |
| // moderate lengths. |
| // |
| // sinl(r) = r + PP_1_hi r^3 + PP_1_lo r^3 + |
| // PP_2 r^5 + ... + PP_8 r^17 |
| // |
| // cosl(r) = 1 + QQ_1 r^2 + QQ_2 r^4 + ... + QQ_8 r^16 |
| // |
| // where PP_1_hi is only about 16 bits long and QQ_1 is -1/2. |
| // The crux in accurate computation is to calculate |
| // |
| // r + PP_1_hi r^3 or 1 + QQ_1 r^2 |
| // |
| // accurately as two pieces: U_hi and U_lo. The way to achieve this |
| // is to obtain r_hi as a 10 sig. bit number that approximates r to |
| // roughly 8 bits or so of accuracy. (One convenient way is |
| // |
| // r_hi := frcpa( frcpa( r ) ).) |
| // |
| // This way, |
| // |
| // r + PP_1_hi r^3 = r + PP_1_hi r_hi^3 + |
| // PP_1_hi (r^3 - r_hi^3) |
| // = [r + PP_1_hi r_hi^3] + |
| // [PP_1_hi (r - r_hi) |
| // (r^2 + r_hi r + r_hi^2) ] |
| // = U_hi + U_lo |
| // |
| // Since r_hi is only 10 bit long and PP_1_hi is only 16 bit long, |
| // PP_1_hi * r_hi^3 is only at most 46 bit long and thus computed |
| // exactly. Furthermore, r and PP_1_hi r_hi^3 are of opposite sign |
| // and that there is no more than 8 bit shift off between r and |
| // PP_1_hi * r_hi^3. Hence the sum, U_hi, is representable and thus |
| // calculated without any error. Finally, the fact that |
| // |
| // |U_lo| <= 2^(-8) |U_hi| |
| // |
| // says that U_hi + U_lo is approximating r + PP_1_hi r^3 to roughly |
| // 8 extra bits of accuracy. |
| // |
| // Similarly, |
| // |
| // 1 + QQ_1 r^2 = [1 + QQ_1 r_hi^2] + |
| // [QQ_1 (r - r_hi)(r + r_hi)] |
| // = U_hi + U_lo. |
| // |
| // Summarizing, we calculate r_hi = frcpa( frcpa( r ) ). |
| // |
| // If i_1 = 0, then |
| // |
| // U_hi := r + PP_1_hi * r_hi^3 |
| // U_lo := PP_1_hi * (r - r_hi) * (r^2 + r*r_hi + r_hi^2) |
| // poly := PP_1_lo r^3 + PP_2 r^5 + ... + PP_8 r^17 |
| // correction := c * ( 1 + C_1 r^2 ) |
| // |
| // Else ...i_1 = 1 |
| // |
| // U_hi := 1 + QQ_1 * r_hi * r_hi |
| // U_lo := QQ_1 * (r - r_hi) * (r + r_hi) |
| // poly := QQ_2 * r^4 + QQ_3 * r^6 + ... + QQ_8 r^16 |
| // correction := -c * r * (1 + S_1 * r^2) |
| // |
| // End |
| // |
| // Finally, |
| // |
| // V := poly + ( U_lo + correction ) |
| // |
| // / U_hi + V if i_0 = 0 |
| // result := | |
| // \ (-U_hi) - V if i_0 = 1 |
| // |
| // It is important that in the last step, negation of U_hi is |
| // performed prior to the subtraction which is to be performed in |
| // the user-set rounding mode. |
| // |
| // |
| // Algorithmic Description |
| // ======================= |
| // |
| // The argument reduction algorithm shares the same code between FSIN and FCOS. |
| // The argument reduction description given |
| // previously is repeated below. |
| // |
| // |
| // Step 0. Initialization. |
| // |
| // 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 64 |
| |
| LOCAL_OBJECT_START(FSINCOSL_CONSTANTS) |
| |
| sincosl_table_p: |
| //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 0xDBD171A1, 0x8D848E89, 0x0000BFBF,0x00000000 // d_1 |
| //data4 0x18A66F8E, 0xD5394C36, 0x0000BF7C,0x00000000 // d_2 |
| 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) |
| //data4 0x2168C234, 0xC90FDAA2, 0x00003FFE,0x00000000 // pi_by_4 |
| //data4 0x6EC6B45A, 0xA397E504, 0x00003FE7,0x00000000 // Inv_P_0 |
| 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) |
| //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 |
| 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) |
| //data4 0xC2B0FE52, 0xD56232EF, 0x00003FD2 // QQ_8 |
| //data4 0x2B48DCA6, 0xC9C99ABA, 0x0000BFDA // QQ_7 |
| //data4 0x9C716658, 0x8F76C650, 0x00003FE2 // QQ_6 |
| //data4 0xFDA8D0FC, 0x93F27DBA, 0x0000BFE9 // QQ_5 |
| //data4 0xAAAAAAAA, 0xAAAAAAAA, 0x0000BFFC // 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 |
| 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) |
| //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 |
| 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) |
| //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 |
| 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_ResultS = f9 |
| FR_ResultC = 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_r6 = f32 |
| FR_r7 = f33 |
| 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_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_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_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_d_2 = f93 |
| FR_P_0 = f95 |
| FR_C_lo = f96 |
| FR_PP_1 = f97 |
| FR_PP_1_lo = f98 |
| FR_ArgPrime = f99 |
| FR_inexact = f100 |
| |
| FR_Neg_Two_to_M3 = f109 |
| FR_Two_to_M3 = f110 |
| |
| FR_poly_hiS = f66 |
| FR_poly_hiC = f112 |
| |
| FR_poly_loS = f67 |
| FR_poly_loC = f113 |
| |
| FR_polyS = f92 |
| FR_polyC = f114 |
| |
| FR_cS = FR_c |
| FR_cC = f115 |
| |
| FR_corrS = f91 |
| FR_corrC = f116 |
| |
| FR_U_hiC = f117 |
| FR_U_loC = f118 |
| |
| FR_VS = f75 |
| FR_VC = f119 |
| |
| FR_FirstS = f120 |
| FR_FirstC = f121 |
| |
| FR_U_hiS = FR_U_hi |
| FR_U_loS = FR_U_lo |
| |
| FR_Tmp = f94 |
| |
| |
| |
| |
| sincos_pResSin = r34 |
| sincos_pResCos = r35 |
| |
| GR_exp_m2_to_m3= r36 |
| GR_N_Inc = r37 |
| GR_Cis = 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_N_SignS = r45 |
| GR_N_SignC = r46 |
| GR_N_SinCos = r47 |
| |
| GR_sig_inv_pi = r48 |
| GR_rshf_2to64 = r49 |
| GR_exp_2tom64 = r50 |
| GR_rshf = r51 |
| GR_ad_p = r52 |
| GR_ad_d = r53 |
| GR_ad_pp = r54 |
| GR_ad_qq = r55 |
| GR_ad_c = r56 |
| GR_ad_s = r57 |
| GR_ad_ce = r58 |
| GR_ad_se = r59 |
| GR_ad_m14 = r60 |
| GR_ad_s1 = r61 |
| |
| // For unwind support |
| GR_SAVE_B0 = r39 |
| GR_SAVE_GP = r40 |
| GR_SAVE_PFS = r41 |
| |
| |
| .section .text |
| |
| GLOBAL_IEEE754_ENTRY(sincosl) |
| { .mlx ///////////////////////////// 1 ///////////////// |
| alloc r32 = ar.pfs,3,27,2,0 |
| movl GR_sig_inv_pi = 0xa2f9836e4e44152a // significand of 1/pi |
| } |
| { .mlx |
| mov GR_N_Inc = 0x0 |
| movl GR_rshf_2to64 = 0x47e8000000000000 // 1.1000 2^(63+64) |
| };; |
| |
| { .mfi ///////////////////////////// 2 ///////////////// |
| 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 |
| mov GR_Cis = 0x0 |
| fnorm.s1 FR_norm_x = FR_Input_X // Normalize x |
| br.cond.sptk _COMMON_SINCOSL |
| };; |
| GLOBAL_IEEE754_END(sincosl) |
| |
| GLOBAL_LIBM_ENTRY(__libm_sincosl) |
| { .mlx ///////////////////////////// 1 ///////////////// |
| alloc r32 = ar.pfs,3,27,2,0 |
| movl GR_sig_inv_pi = 0xa2f9836e4e44152a // significand of 1/pi |
| } |
| { .mlx |
| mov GR_N_Inc = 0x0 |
| movl GR_rshf_2to64 = 0x47e8000000000000 // 1.1000 2^(63+64) |
| };; |
| |
| { .mfi ///////////////////////////// 2 ///////////////// |
| 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 |
| mov GR_Cis = 0x1 |
| fnorm.s1 FR_norm_x = FR_Input_X // Normalize x |
| nop.b 0 |
| };; |
| |
| _COMMON_SINCOSL: |
| { .mfi ///////////////////////////// 3 ///////////////// |
| setf.sig FR_inv_pi_2to63 = GR_sig_inv_pi // Form 1/pi * 2^63 |
| nop.f 0 |
| 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 ///////////////////////////// 4 ///////////////// |
| 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 0 |
| };; |
| |
| { .mfi ///////////////////////////// 5 ///////////////// |
| 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 0 |
| } |
| { .mib |
| mov GR_exp_mask = 0x1ffff // Exponent mask |
| nop.i 0 |
| (p6) br.cond.spnt SINCOSL_SPECIAL // Branch if x natval, nan, inf |
| };; |
| |
| { .mfi ///////////////////////////// 6 ///////////////// |
| setf.exp FR_2tom64 = GR_exp_2tom64 // Form 2^-64 for scaling N_float |
| nop.f 0 |
| add GR_ad_d = 0x70, GR_ad_p // Point to constant table d |
| } |
| { .mib |
| setf.d FR_rshf = GR_rshf // Form right shift const 1.1000 * 2^63 |
| mov GR_exp_m2_to_m3 = 0x2fffc // Form -(2^-3) |
| (p7) br.cond.spnt SINCOSL_DENORMAL // Branch if x denormal |
| };; |
| |
| SINCOSL_COMMON2: |
| { .mfi ///////////////////////////// 7 ///////////////// |
| and GR_exp_x = GR_exp_mask, GR_signexp_x // Get exponent of x |
| fclass.nm p8, p0 = FR_Input_X, 0x1FF // Test x unsupported type |
| mov GR_exp_2_to_63 = 0xffff + 63 // Exponent of 2^63 |
| } |
| { .mib |
| add GR_ad_pp = 0x40, GR_ad_d // Point to constant table pp |
| mov GR_exp_2_to_24 = 0xffff + 24 // Exponent of 2^24 |
| (p10) br.cond.spnt SINCOSL_ZERO // Branch if x zero |
| };; |
| |
| { .mfi ///////////////////////////// 8 ///////////////// |
| ldfe FR_Inv_pi_by_2 = [GR_ad_p], 16 // Load 2/pi |
| fcmp.eq.s0 p15, p0 = FR_Input_X, f0 // Dummy to set denormal |
| add GR_ad_qq = 0xa0, GR_ad_pp // Point to constant table qq |
| } |
| { .mfi |
| ldfe FR_Pi_by_4 = [GR_ad_d], 16 // Load pi/4 for range test |
| nop.f 0 |
| cmp.ge p10,p0 = GR_exp_x, GR_exp_2_to_63 // Is |x| >= 2^63 |
| };; |
| |
| { .mfi ///////////////////////////// 9 ///////////////// |
| ldfe FR_P_0 = [GR_ad_p], 16 // Load P_0 for pi/4 <= |x| < 2^63 |
| fmerge.s FR_abs_x = f1, FR_norm_x // |x| |
| add GR_ad_c = 0x90, GR_ad_qq // Point to constant table c |
| } |
| { .mfi |
| ldfe FR_Inv_P_0 = [GR_ad_d], 16 // Load 1/P_0 for pi/4 <= |x| < 2^63 |
| nop.f 0 |
| cmp.ge p7,p0 = GR_exp_x, GR_exp_2_to_24 // Is |x| >= 2^24 |
| };; |
| |
| { .mfi ///////////////////////////// 10 ///////////////// |
| ldfe FR_P_1 = [GR_ad_p], 16 // Load P_1 for pi/4 <= |x| < 2^63 |
| nop.f 0 |
| add GR_ad_s = 0x50, GR_ad_c // Point to constant table s |
| } |
| { .mfi |
| ldfe FR_PP_8 = [GR_ad_pp], 16 // Load PP_8 for 2^-3 < |r| < pi/4 |
| nop.f 0 |
| nop.i 0 |
| };; |
| |
| { .mfi ///////////////////////////// 11 ///////////////// |
| ldfe FR_P_2 = [GR_ad_p], 16 // Load P_2 for pi/4 <= |x| < 2^63 |
| nop.f 0 |
| add GR_ad_ce = 0x40, GR_ad_c // Point to end of constant table c |
| } |
| { .mfi |
| ldfe FR_QQ_8 = [GR_ad_qq], 16 // Load QQ_8 for 2^-3 < |r| < pi/4 |
| nop.f 0 |
| nop.i 0 |
| };; |
| |
| { .mfi ///////////////////////////// 12 ///////////////// |
| ldfe FR_QQ_7 = [GR_ad_qq], 16 // Load QQ_7 for 2^-3 < |r| < pi/4 |
| fma.s1 FR_N_float_signif = FR_Input_X, FR_inv_pi_2to63, FR_rshf_2to64 |
| add GR_ad_se = 0x40, GR_ad_s // Point to end of constant table s |
| } |
| { .mib |
| ldfe FR_PP_7 = [GR_ad_pp], 16 // Load PP_7 for 2^-3 < |r| < pi/4 |
| mov GR_ad_s1 = GR_ad_s // Save pointer to S_1 |
| (p10) br.cond.spnt SINCOSL_ARG_TOO_LARGE // Branch if |x| >= 2^63 |
| // Use Payne-Hanek Reduction |
| };; |
| |
| { .mfi ///////////////////////////// 13 ///////////////// |
| ldfe FR_P_3 = [GR_ad_p], 16 // Load P_3 for pi/4 <= |x| < 2^63 |
| fmerge.se FR_r = FR_norm_x, FR_norm_x // r = x, in case |x| < pi/4 |
| add GR_ad_m14 = 0x50, GR_ad_s // Point to constant table m14 |
| } |
| { .mfb |
| ldfps FR_Two_to_M3, FR_Neg_Two_to_M3 = [GR_ad_d], 8 |
| fma.s1 FR_rsq = FR_norm_x, FR_norm_x, f0 // rsq = x*x, in case |x| < pi/4 |
| (p7) br.cond.spnt SINCOSL_LARGER_ARG // Branch if 2^24 <= |x| < 2^63 |
| // Use pre-reduction |
| };; |
| |
| { .mmf ///////////////////////////// 14 ///////////////// |
| 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 ///////////////////////////// 15 ///////////////// |
| 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 0 |
| };; |
| |
| // Here if 0 < |x| < 2^24 |
| { .mfi ///////////////////////////// 17 ///////////////// |
| 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 0 |
| } |
| { .mfi |
| ldfe FR_C_5 = [GR_ad_ce], -16 // Load C_5 if i_1=1 |
| fms.s1 FR_N_float = FR_N_float_signif, FR_2tom64, FR_rshf |
| nop.i 0 |
| };; |
| |
| { .mmi ///////////////////////////// 18 ///////////////// |
| 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 0 |
| };; |
| |
| // |
| // N = Arg * 2/pi |
| // Check if Arg < pi/4 |
| // |
| // |
| // Case 2: Convert integer N_fix back to normalized floating-point value. |
| // Case 1: p8 is only affected when p6 is set |
| // |
| // |
| // Grab the integer part of N and call it N_fix |
| // |
| { .mfi ///////////////////////////// 19 ///////////////// |
| (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 = 0x0 // N_IncS if |x| < pi/4 |
| };; |
| |
| // If |x| < pi/4, r = x and c = 0 |
| // lf |x| < pi/4, is x < 2**(-3). |
| // r = Arg |
| // c = 0 |
| { .mmi ///////////////////////////// 20 ///////////////// |
| (p7) getf.sig GR_N_Inc = FR_N_float_signif |
| nop.m 0 |
| (p6) cmp.lt.unc p8,p0 = GR_exp_x, GR_exp_2_to_m3 // Is |x| < 2^-3 |
| };; |
| |
| // |
| // lf |x| < pi/4, is -2**(-3)< x < 2**(-3) - set p8. |
| // If |x| >= pi/4, |
| // Create the right N for |x| < pi/4 and otherwise |
| // Case 2: Place integer part of N in GP register |
| // |
| |
| { .mbb ///////////////////////////// 21 ///////////////// |
| nop.m 0 |
| (p8) br.cond.spnt SINCOSL_SMALL_R_0 // Branch if 0 < |x| < 2^-3 |
| (p6) br.cond.spnt SINCOSL_NORMAL_R_0 // Branch if 2^-3 <= |x| < pi/4 |
| };; |
| |
| // Here if pi/4 <= |x| < 2^24 |
| { .mfi |
| ldfs FR_Neg_Two_to_M67 = [GR_ad_d], 8 // Load -2^-67 |
| fnma.s1 FR_s = FR_N_float, FR_P_1, FR_Input_X // s = -N * P_1 + Arg |
| nop.i 0 |
| } |
| { .mfi |
| nop.m 0 |
| fma.s1 FR_w = FR_N_float, FR_P_2, f0 // w = N * P_2 |
| nop.i 0 |
| };; |
| |
| { .mfi |
| nop.m 0 |
| fms.s1 FR_r = FR_s, f1, FR_w // r = s - w, assume |s| >= 2^-33 |
| nop.i 0 |
| };; |
| |
| { .mfi |
| nop.m 0 |
| fcmp.lt.s1 p7, p6 = FR_s, FR_Two_to_M33 |
| nop.i 0 |
| };; |
| |
| { .mfi |
| nop.m 0 |
| (p7) fcmp.gt.s1 p7, p6 = FR_s, FR_Neg_Two_to_M33 // p6 if |s| >= 2^-33, else p7 |
| nop.i 0 |
| };; |
| |
| { .mfi |
| nop.m 0 |
| fms.s1 FR_c = FR_s, f1, FR_r // c = s - r, for |s| >= 2^-33 |
| nop.i 0 |
| } |
| { .mfi |
| nop.m 0 |
| fma.s1 FR_rsq = FR_r, FR_r, f0 // rsq = r * r, for |s| >= 2^-33 |
| nop.i 0 |
| };; |
| |
| { .mfi |
| nop.m 0 |
| (p7) fma.s1 FR_w = FR_N_float, FR_P_3, f0 |
| nop.i 0 |
| };; |
| |
| { .mmf |
| ldfe FR_C_1 = [GR_ad_pp], 16 // Load C_1 if i_1=0 |
| 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 0 |
| (p6) fcmp.lt.unc.s1 p8, p13 = FR_r, FR_Two_to_M3 // If big s, test r with 2^-3 |
| nop.i 0 |
| };; |
| |
| { .mfi |
| nop.m 0 |
| (p7) fma.s1 FR_U_1 = FR_N_float, FR_P_2, FR_w |
| nop.i 0 |
| };; |
| |
| // |
| // For big s: r = s - w: No futher reduction is necessary |
| // For small s: w = N * P_3 (change sign) More reduction |
| // |
| { .mfi |
| nop.m 0 |
| (p8) fcmp.gt.s1 p8, p13 = FR_r, FR_Neg_Two_to_M3 // If big s, p8 if |r| < 2^-3 |
| nop.i 0 |
| };; |
| |
| { .mfi |
| nop.m 0 |
| fma.s1 FR_polyS = FR_rsq, FR_PP_8, FR_PP_7 // poly = rsq*PP_8+PP_7 |
| nop.i 0 |
| } |
| { .mfi |
| nop.m 0 |
| fma.s1 FR_polyC = FR_rsq, FR_QQ_8, FR_QQ_7 // poly = rsq*QQ_8+QQ_7 |
| nop.i 0 |
| };; |
| |
| { .mfi |
| nop.m 0 |
| (p7) fms.s1 FR_r = FR_s, f1, FR_U_1 |
| nop.i 0 |
| };; |
| |
| { .mfi |
| nop.m 0 |
| (p6) fma.s1 FR_r_cubed = FR_r, FR_rsq, f0 // rcubed = r * rsq |
| nop.i 0 |
| };; |
| |
| { .mfi |
| // |
| // For big s: Is |r| < 2**(-3)? |
| // For big s: c = S - r |
| // For small s: U_1 = N * P_2 + w |
| // |
| // If p8 is set, prepare to branch to Small_R. |
| // If p9 is set, prepare to branch to Normal_R. |
| // For big s, r is complete here. |
| // |
| // |
| // For big s: c = c + w (w has not been negated.) |
| // For small s: r = S - U_1 |
| // |
| nop.m 0 |
| (p6) fms.s1 FR_c = FR_c, f1, FR_w |
| nop.i 0 |
| } |
| { .mbb |
| nop.m 0 |
| (p8) br.cond.spnt SINCOSL_SMALL_R_1 // Branch if |s|>=2^-33, |r| < 2^-3, |
| // and pi/4 <= |x| < 2^24 |
| (p13) br.cond.sptk SINCOSL_NORMAL_R_1 // Branch if |s|>=2^-33, |r| >= 2^-3, |
| // and pi/4 <= |x| < 2^24 |
| };; |
| |
| SINCOSL_S_TINY: |
| // |
| // Here if |s| < 2^-33, and pi/4 <= |x| < 2^24 |
| // |
| { .mfi |
| and GR_N_SinCos = 0x1, GR_N_Inc |
| fms.s1 FR_U_2 = FR_N_float, FR_P_2, FR_U_1 |
| tbit.z p8,p12 = GR_N_Inc, 0 |
| };; |
| |
| |
| // |
| // For small s: U_2 = N * P_2 - U_1 |
| // S_1 stored constant - grab the one stored with the |
| // coefficients. |
| // |
| { .mfi |
| ldfe FR_S_1 = [GR_ad_s1], 16 |
| fma.s1 FR_polyC = f0, f1, FR_Neg_Two_to_M67 |
| sub GR_N_SignS = GR_N_Inc, GR_N_SinCos |
| } |
| { .mfi |
| add GR_N_SignC = GR_N_Inc, GR_N_SinCos |
| nop.f 0 |
| nop.i 0 |
| };; |
| |
| { .mfi |
| nop.m 0 |
| fms.s1 FR_s = FR_s, f1, FR_r |
| (p8) tbit.z.unc p10,p11 = GR_N_SignC, 1 |
| } |
| { .mfi |
| nop.m 0 |
| fma.s1 FR_rsq = FR_r, FR_r, f0 |
| nop.i 0 |
| };; |
| |
| { .mfi |
| nop.m 0 |
| fma.s1 FR_U_2 = FR_U_2, f1, FR_w |
| (p8) tbit.z.unc p8,p9 = GR_N_SignS, 1 |
| };; |
| |
| { .mfi |
| nop.m 0 |
| fmerge.se FR_FirstS = FR_r, FR_r |
| (p12) tbit.z.unc p14,p15 = GR_N_SignC, 1 |
| } |
| { .mfi |
| nop.m 0 |
| fma.s1 FR_FirstC = f0, f1, f1 |
| nop.i 0 |
| };; |
| |
| { .mfi |
| nop.m 0 |
| fms.s1 FR_c = FR_s, f1, FR_U_1 |
| (p12) tbit.z.unc p12,p13 = GR_N_SignS, 1 |
| };; |
| |
| { .mfi |
| nop.m 0 |
| fma.s1 FR_r = FR_S_1, FR_r, f0 |
| nop.i 0 |
| };; |
| |
| { .mfi |
| nop.m 0 |
| fma.s0 FR_S_1 = FR_S_1, FR_S_1, f0 |
| nop.i 0 |
| };; |
| |
| { .mfi |
| nop.m 0 |
| fms.s1 FR_c = FR_c, f1, FR_U_2 |
| nop.i 0 |
| };; |
| |
| .pred.rel "mutex",p9,p15 |
| { .mfi |
| nop.m 0 |
| (p9) fms.s0 FR_FirstS = f1, f0, FR_FirstS |
| nop.i 0 |
| } |
| { .mfi |
| nop.m 0 |
| (p15) fms.s0 FR_FirstS = f1, f0, FR_FirstS |
| nop.i 0 |
| };; |
| |
| .pred.rel "mutex",p11,p13 |
| { .mfi |
| nop.m 0 |
| (p11) fms.s0 FR_FirstC = f1, f0, FR_FirstC |
| nop.i 0 |
| } |
| { .mfi |
| nop.m 0 |
| (p13) fms.s0 FR_FirstC = f1, f0, FR_FirstC |
| nop.i 0 |
| };; |
| |
| { .mfi |
| nop.m 0 |
| fma.s1 FR_polyS = FR_r, FR_rsq, FR_c |
| nop.i 0 |
| };; |
| |
| |
| .pred.rel "mutex",p8,p9 |
| { .mfi |
| nop.m 0 |
| (p8) fma.s0 FR_ResultS = FR_FirstS, f1, FR_polyS |
| nop.i 0 |
| } |
| { .mfi |
| nop.m 0 |
| (p9) fms.s0 FR_ResultS = FR_FirstS, f1, FR_polyS |
| nop.i 0 |
| };; |
| |
| .pred.rel "mutex",p10,p11 |
| { .mfi |
| nop.m 0 |
| (p10) fma.s0 FR_ResultC = FR_FirstC, f1, FR_polyC |
| nop.i 0 |
| } |
| { .mfi |
| nop.m 0 |
| (p11) fms.s0 FR_ResultC = FR_FirstC, f1, FR_polyC |
| nop.i 0 |
| };; |
| |
| |
| |
| .pred.rel "mutex",p12,p13 |
| { .mfi |
| nop.m 0 |
| (p12) fma.s0 FR_ResultS = FR_FirstC, f1, FR_polyC |
| nop.i 0 |
| } |
| { .mfi |
| nop.m 0 |
| (p13) fms.s0 FR_ResultS = FR_FirstC, f1, FR_polyC |
| nop.i 0 |
| };; |
| |
| .pred.rel "mutex",p14,p15 |
| { .mfi |
| nop.m 0 |
| (p14) fma.s0 FR_ResultC = FR_FirstS, f1, FR_polyS |
| nop.i 0 |
| } |
| { .mfb |
| cmp.eq p10, p0 = 0x1, GR_Cis |
| (p15) fms.s0 FR_ResultC = FR_FirstS, f1, FR_polyS |
| (p10) br.ret.sptk b0 |
| };; |
| |
| { .mmb // exit for sincosl |
| stfe [sincos_pResSin] = FR_ResultS |
| stfe [sincos_pResCos] = FR_ResultC |
| br.ret.sptk b0 |
| };; |
| |
| |
| |
| |
| |
| |
| SINCOSL_LARGER_ARG: |
| // |
| // 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 // N_0 = Arg * Inv_P_0 |
| nop.i 0 |
| };; |
| |
| { .mmi |
| ldfps FR_Two_to_M14, FR_Neg_Two_to_M14 = [GR_ad_m14] |
| nop.m 0 |
| nop.i 0 |
| };; |
| |
| { .mfi |
| ldfe FR_d_2 = [GR_ad_p], 16 // Load d_2 for |x| >= 2^24 path |
| nop.f 0 |
| nop.i 0 |
| };; |
| |
| { .mfi |
| nop.m 0 |
| fcvt.fx.s1 FR_N_0_fix = FR_N_0 // N_0_fix = integer part of N_0 |
| nop.i 0 |
| };; |
| |
| { .mfi |
| nop.m 0 |
| fcvt.xf FR_N_0 = FR_N_0_fix // Make N_0 the integer part |
| nop.i 0 |
| };; |
| |
| { .mfi |
| nop.m 0 |
| fnma.s1 FR_ArgPrime = FR_N_0, FR_P_0, FR_Input_X // Arg'=-N_0*P_0+Arg |
| nop.i 0 |
| } |
| { .mfi |
| nop.m 0 |
| fma.s1 FR_w = FR_N_0, FR_d_1, f0 // w = N_0 * d_1 |
| nop.i 0 |
| };; |
| |
| |
| { .mfi |
| nop.m 0 |
| fma.s1 FR_N_float = FR_ArgPrime, FR_Inv_pi_by_2, f0 // N = A' * 2/pi |
| nop.i 0 |
| };; |
| |
| { .mfi |
| nop.m 0 |
| fcvt.fx.s1 FR_N_fix = FR_N_float // N_fix is the integer part |
| nop.i 0 |
| };; |
| |
| { .mfi |
| nop.m 0 |
| fcvt.xf FR_N_float = FR_N_fix |
| nop.i 0 |
| };; |
| |
| { .mfi |
| getf.sig GR_N_Inc = FR_N_fix // N is the integer part of |
| // the reduced-reduced argument |
| nop.f 0 |
| nop.i 0 |
| };; |
| |
| { .mfi |
| nop.m 0 |
| fnma.s1 FR_s = FR_N_float, FR_P_1, FR_ArgPrime // s = -N*P_1 + Arg' |
| nop.i 0 |
| } |
| { .mfi |
| nop.m 0 |
| fnma.s1 FR_w = FR_N_float, FR_P_2, FR_w // w = -N*P_2 + w |
| nop.i 0 |
| };; |
| |
| // |
| // For |s| > 2**(-14) r = S + w (r complete) |
| // Else U_hi = N_0 * d_1 |
| // |
| { .mfi |
| nop.m 0 |
| fcmp.lt.unc.s1 p9, p8 = FR_s, FR_Two_to_M14 |
| nop.i 0 |
| };; |
| |
| { .mfi |
| nop.m 0 |
| (p9) fcmp.gt.s1 p9, p8 = FR_s, FR_Neg_Two_to_M14 // p9 if |s| < 2^-14 |
| nop.i 0 |
| };; |
| |
| // |
| // Either S <= -2**(-14) or S >= 2**(-14) |
| // or -2**(-14) < s < 2**(-14) |
| // |
| { .mfi |
| nop.m 0 |
| (p9) fma.s1 FR_V_hi = FR_N_float, FR_P_2, f0 |
| nop.i 0 |
| } |
| { .mfi |
| nop.m 0 |
| (p9) fma.s1 FR_U_hi = FR_N_0, FR_d_1, f0 |
| nop.i 0 |
| };; |
| |
| { .mfi |
| nop.m 0 |
| (p8) fma.s1 FR_r = FR_s, f1, FR_w |
| nop.i 0 |
| } |
| { .mfi |
| nop.m 0 |
| (p9) fma.s1 FR_w = FR_N_float, FR_P_3, f0 |
| nop.i 0 |
| };; |
| |
| // |
| // We need abs of both U_hi and V_hi - don't |
| // worry about switched sign of V_hi. |
| // |
| // 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. |
| // |
| { .mfi |
| nop.m 0 |
| (p9) fms.s1 FR_A = FR_U_hi, f1, FR_V_hi |
| nop.i 0 |
| } |
| { .mfi |
| nop.m 0 |
| (p9) fnma.s1 FR_V_lo = FR_N_float, FR_P_2, FR_V_hi |
| nop.i 0 |
| };; |
| |
| { .mfi |
| nop.m 0 |
| (p9) fmerge.s FR_V_hiabs = f0, FR_V_hi |
| nop.i 0 |
| } |
| { .mfi |
| nop.m 0 |
| (p9) fms.s1 FR_U_lo = FR_N_0, FR_d_1, FR_U_hi // For small s: U_lo=N_0*d_1-U_hi |
| nop.i 0 |
| };; |
| |
| // |
| // 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. |
| // |
| { .mfi |
| nop.m 0 |
| (p9) fmerge.s FR_U_hiabs = f0, FR_U_hi |
| nop.i 0 |
| } |
| { .mfi |
| nop.m 0 |
| (p8) fms.s1 FR_c = FR_s, f1, FR_r // For big s: c = S - r |
| nop.i 0 |
| };; |
| |
| // |
| // 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. |
| // |
| { .mfi |
| nop.m 0 |
| (p8) fcmp.lt.unc.s1 p12, p13 = FR_r, FR_Two_to_M3 |
| nop.i 0 |
| };; |
| |
| { .mfi |
| nop.m 0 |
| (p12) fcmp.gt.s1 p12, p13 = FR_r, FR_Neg_Two_to_M3 |
| nop.i 0 |
| };; |
| |
| { .mfi |
| nop.m 0 |
| (p8) fma.s1 FR_c = FR_c, f1, FR_w |
| nop.i 0 |
| } |
| { .mfb |
| nop.m 0 |
| (p9) fms.s1 FR_w = FR_N_0, FR_d_2, FR_w |
| (p12) br.cond.spnt SINCOSL_SMALL_R // Branch if |r| < 2^-3 |
| // and 2^24 <= |x| < 2^63 |
| };; |
| |
| { .mib |
| nop.m 0 |
| nop.i 0 |
| (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 |
| // |
| // 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. |
| // |
| { .mfi |
| and GR_N_SinCos = 0x1, GR_N_Inc |
| fcmp.ge.unc.s1 p6, p7 = FR_U_hiabs, FR_V_hiabs |
| tbit.z p8,p12 = GR_N_Inc, 0 |
| } |
| { .mfi |
| nop.m 0 |
| fma.s1 FR_t = FR_U_lo, f1, FR_V_lo // C_hi = S + A |
| nop.i 0 |
| };; |
| |
| { .mfi |
| sub GR_N_SignS = GR_N_Inc, GR_N_SinCos |
| (p6) fms.s1 FR_a = FR_U_hi, f1, FR_A |
| add GR_N_SignC = GR_N_Inc, GR_N_SinCos |
| } |
| { .mfi |
| nop.m 0 |
| (p7) fma.s1 FR_a = FR_V_hi, f1, FR_A |
| nop.i 0 |
| };; |
| |
| { .mmf |
| ldfe FR_C_1 = [GR_ad_c], 16 |
| ldfe FR_S_1 = [GR_ad_s], 16 |
| fma.s1 FR_C_hi = FR_s, f1, FR_A |
| };; |
| |
| { .mmi |
| ldfe FR_C_2 = [GR_ad_c], 64 |
| ldfe FR_S_2 = [GR_ad_s], 64 |
| (p8) tbit.z.unc p10,p11 = GR_N_SignC, 1 |
| };; |
| |
| // |
| // 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. |
| // |
| // For larger u than v: a = U_hi - A |
| // Else a = V_hi - A (do an add to account for missing (-) on V_hi |
| // |
| { .mfi |
| nop.m 0 |
| fma.s1 FR_t = FR_t, f1, FR_w // t = t + w |
| (p8) tbit.z.unc p8,p9 = GR_N_SignS, 1 |
| } |
| { .mfi |
| nop.m 0 |
| (p6) fms.s1 FR_a = FR_a, f1, FR_V_hi |
| nop.i 0 |
| };; |
| |
| // |
| // 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. |
| // |
| { .mfi |
| nop.m 0 |
| fms.s1 FR_C_lo = FR_s, f1, FR_C_hi |
| (p12) tbit.z.unc p14,p15 = GR_N_SignC, 1 |
| } |
| { .mfi |
| nop.m 0 |
| (p7) fms.s1 FR_a = FR_U_hi, f1, FR_a |
| nop.i 0 |
| };; |
| |
| { .mfi |
| nop.m 0 |
| fma.s1 FR_C_lo = FR_C_lo, f1, FR_A // C_lo = (S - C_hi) + A |
| (p12) tbit.z.unc p12,p13 = GR_N_SignS, 1 |
| } |
| { .mfi |
| nop.m 0 |
| fma.s1 FR_t = FR_t, f1, FR_a // t = t + a |
| nop.i 0 |
| };; |
| |
| { .mfi |
| nop.m 0 |
| fma.s1 FR_r = FR_C_hi, f1, FR_C_lo |
| nop.i 0 |
| };; |
| |
| { .mfi |
| nop.m 0 |
| fma.s1 FR_C_lo = FR_C_lo, f1, FR_t // C_lo = C_lo + t |
| nop.i 0 |
| };; |
| |
| |
| { .mfi |
| nop.m 0 |
| fma.s1 FR_rsq = FR_r, FR_r, f0 |
| nop.i 0 |
| } |
| { .mfi |
| nop.m 0 |
| fms.s1 FR_c = FR_C_hi, f1, FR_r |
| nop.i 0 |
| };; |
| |
| { .mfi |
| nop.m 0 |
| fma.s1 FR_FirstS = f0, f1, FR_r |
| nop.i 0 |
| } |
| { .mfi |
| nop.m 0 |
| fma.s1 FR_FirstC = f0, f1, f1 |
| nop.i 0 |
| };; |
| |
| { .mfi |
| nop.m 0 |
| fma.s1 FR_polyS = FR_rsq, FR_S_2, FR_S_1 |
| nop.i 0 |
| } |
| { .mfi |
| nop.m 0 |
| fma.s1 FR_polyC = FR_rsq, FR_C_2, FR_C_1 |
| nop.i 0 |
| };; |
| |
| { .mfi |
| nop.m 0 |
| fma.s1 FR_r_cubed = FR_rsq, FR_r, f0 |
| nop.i 0 |
| } |
| { .mfi |
| nop.m 0 |
| fma.s1 FR_c = FR_c, f1, FR_C_lo |
| nop.i 0 |
| };; |
| |
| .pred.rel "mutex",p9,p15 |
| { .mfi |
| nop.m 0 |
| (p9) fms.s0 FR_FirstS = f1, f0, FR_FirstS |
| nop.i 0 |
| } |
| { .mfi |
| nop.m 0 |
| (p15) fms.s0 FR_FirstS = f1, f0, FR_FirstS |
| nop.i 0 |
| };; |
| |
| .pred.rel "mutex",p11,p13 |
| { .mfi |
| nop.m 0 |
| (p11) fms.s0 FR_FirstC = f1, f0, FR_FirstC |
| nop.i 0 |
| } |
| { .mfi |
| nop.m 0 |
| (p13) fms.s0 FR_FirstC = f1, f0, FR_FirstC |
| nop.i 0 |
| };; |
| |
| |
| { .mfi |
| nop.m 0 |
| fma.s1 FR_polyS = FR_r_cubed, FR_polyS, FR_c |
| nop.i 0 |
| } |
| { .mfi |
| nop.m 0 |
| fma.s1 FR_polyC = FR_rsq, FR_polyC, f0 |
| nop.i 0 |
| };; |
| |
| |
| |
| .pred.rel "mutex",p8,p9 |
| { .mfi |
| nop.m 0 |
| (p8) fma.s0 FR_ResultS = FR_FirstS, f1, FR_polyS |
| nop.i 0 |
| } |
| { .mfi |
| nop.m 0 |
| (p9) fms.s0 FR_ResultS = FR_FirstS, f1, FR_polyS |
| nop.i 0 |
| };; |
| |
| .pred.rel "mutex",p10,p11 |
| { .mfi |
| nop.m 0 |
| (p10) fma.s0 FR_ResultC = FR_FirstC, f1, FR_polyC |
| nop.i 0 |
| } |
| { .mfi |
| nop.m 0 |
| (p11) fms.s0 FR_ResultC = FR_FirstC, f1, FR_polyC |
| nop.i 0 |
| };; |
| |
| |
| |
| .pred.rel "mutex",p12,p13 |
| { .mfi |
| nop.m 0 |
| (p12) fma.s0 FR_ResultS = FR_FirstC, f1, FR_polyC |
| nop.i 0 |
| } |
| { .mfi |
| nop.m 0 |
| (p13) fms.s0 FR_ResultS = FR_FirstC, f1, FR_polyC |
| nop.i 0 |
| };; |
| |
| .pred.rel "mutex",p14,p15 |
| { .mfi |
| nop.m 0 |
| (p14) fma.s0 FR_ResultC = FR_FirstS, f1, FR_polyS |
| nop.i 0 |
| } |
| { .mfb |
| cmp.eq p10, p0 = 0x1, GR_Cis |
| (p15) fms.s0 FR_ResultC = FR_FirstS, f1, FR_polyS |
| (p10) br.ret.sptk b0 |
| };; |
| |
| |
| { .mmb // exit for sincosl |
| stfe [sincos_pResSin] = FR_ResultS |
| stfe [sincos_pResCos] = FR_ResultC |
| br.ret.sptk b0 |
| };; |
| |
| |
| |
| SINCOSL_SMALL_R: |
| // |
| // Here if |r| < 2^-3 |
| // |
| // Enter with r, c, and N_Inc computed |
| // |
| { .mfi |
| nop.m 0 |
| fma.s1 FR_rsq = FR_r, FR_r, f0 // rsq = r * r |
| nop.i 0 |
| };; |
| |
| { .mmi |
| ldfe FR_S_5 = [GR_ad_se], -16 // Load S_5 |
| ldfe FR_C_5 = [GR_ad_ce], -16 // Load C_5 |
| nop.i 0 |
| };; |
| |
| { .mmi |
| ldfe FR_S_4 = [GR_ad_se], -16 // Load S_4 |
| ldfe FR_C_4 = [GR_ad_ce], -16 // Load C_4 |
| nop.i 0 |
| };; |
| |
| SINCOSL_SMALL_R_0: |
| // Entry point for 2^-3 < |x| < pi/4 |
| SINCOSL_SMALL_R_1: |
| // Entry point for pi/4 < |x| < 2^24 and |r| < 2^-3 |
| { .mfi |
| ldfe FR_S_3 = [GR_ad_se], -16 // Load S_3 |
| fma.s1 FR_r6 = FR_rsq, FR_rsq, f0 // Z = rsq * rsq |
| tbit.z p7,p11 = GR_N_Inc, 0 |
| } |
| { .mfi |
| ldfe FR_C_3 = [GR_ad_ce], -16 // Load C_3 |
| nop.f 0 |
| and GR_N_SinCos = 0x1, GR_N_Inc |
| };; |
| |
| { .mfi |
| ldfe FR_S_2 = [GR_ad_se], -16 // Load S_2 |
| fnma.s1 FR_cC = FR_c, FR_r, f0 // c = -c * r |
| sub GR_N_SignS = GR_N_Inc, GR_N_SinCos |
| } |
| { .mfi |
| ldfe FR_C_2 = [GR_ad_ce], -16 // Load C_2 |
| nop.f 0 |
| add GR_N_SignC = GR_N_Inc, GR_N_SinCos |
| };; |
| |
| { .mmi |
| ldfe FR_S_1 = [GR_ad_se], -16 // Load S_1 |
| ldfe FR_C_1 = [GR_ad_ce], -16 // Load C_1 |
| (p7) tbit.z.unc p9,p10 = GR_N_SignC, 1 |
| };; |
| |
| { .mfi |
| nop.m 0 |
| fma.s1 FR_r7 = FR_r6, FR_r, f0 // Z = Z * r |
| (p7) tbit.z.unc p7,p8 = GR_N_SignS, 1 |
| };; |
| |
| { .mfi |
| nop.m 0 |
| fma.s1 FR_poly_loS = FR_rsq, FR_S_5, FR_S_4 // poly_lo=rsq*S_5+S_4 |
| (p11) tbit.z.unc p13,p14 = GR_N_SignC, 1 |
| } |
| { .mfi |
| nop.m 0 |
| fma.s1 FR_poly_loC = FR_rsq, FR_C_5, FR_C_4 // poly_lo=rsq*C_5+C_4 |
| nop.i 0 |
| };; |
| |
| { .mfi |
| nop.m 0 |
| fma.s1 FR_poly_hiS = FR_rsq, FR_S_2, FR_S_1 // poly_hi=rsq*S_2+S_1 |
| (p11) tbit.z.unc p11,p12 = GR_N_SignS, 1 |
| } |
| { .mfi |
| nop.m 0 |
| fma.s1 FR_poly_hiC = FR_rsq, FR_C_2, FR_C_1 // poly_hi=rsq*C_2+C_1 |
| nop.i 0 |
| };; |
| |
| { .mfi |
| nop.m 0 |
| fma.s0 FR_FirstS = FR_r, f1, f0 |
| nop.i 0 |
| } |
| { .mfi |
| nop.m 0 |
| fma.s0 FR_FirstC = f1, f1, f0 |
| nop.i 0 |
| };; |
| |
| |
| { .mfi |
| nop.m 0 |
| fma.s1 FR_r6 = FR_r6, FR_rsq, f0 |
| nop.i 0 |
| } |
| { .mfi |
| nop.m 0 |
| fma.s1 FR_r7 = FR_r7, FR_rsq, f0 |
| nop.i 0 |
| };; |
| |
| { .mfi |
| nop.m 0 |
| fma.s1 FR_poly_loS = FR_rsq, FR_poly_loS, FR_S_3 // p_lo=p_lo*rsq+S_3 |
| nop.i 0 |
| } |
| { .mfi |
| nop.m 0 |
| fma.s1 FR_poly_loC = FR_rsq, FR_poly_loC, FR_C_3 // p_lo=p_lo*rsq+C_3 |
| nop.i 0 |
| };; |
| |
| { .mfi |
| nop.m 0 |
| fma.s0 FR_inexact = FR_S_4, FR_S_4, f0 // Dummy op to set inexact |
| nop.i 0 |
| };; |
| |
| { .mfi |
| nop.m 0 |
| fma.s1 FR_poly_hiS = FR_poly_hiS, FR_rsq, f0 // p_hi=p_hi*rsq |
| nop.i 0 |
| } |
| { .mfi |
| nop.m 0 |
| fma.s1 FR_poly_hiC = FR_poly_hiC, FR_rsq, f0 // p_hi=p_hi*rsq |
| nop.i 0 |
| };; |
| |
| .pred.rel "mutex",p8,p14 |
| { .mfi |
| nop.m 0 |
| (p8) fms.s0 FR_FirstS = f1, f0, FR_FirstS |
| nop.i 0 |
| } |
| { .mfi |
| nop.m 0 |
| (p14) fms.s0 FR_FirstS = f1, f0, FR_FirstS |
| nop.i 0 |
| };; |
| |
| .pred.rel "mutex",p10,p12 |
| { .mfi |
| nop.m 0 |
| (p10) fms.s0 FR_FirstC = f1, f0, FR_FirstC |
| nop.i 0 |
| } |
| { .mfi |
| nop.m 0 |
| (p12) fms.s0 FR_FirstC = f1, f0, FR_FirstC |
| nop.i 0 |
| };; |
| |
| { .mfi |
| nop.m 0 |
| fma.s1 FR_polyS = FR_r7, FR_poly_loS, FR_cS // poly=Z*poly_lo+c |
| nop.i 0 |
| } |
| { .mfi |
| nop.m 0 |
| fma.s1 FR_polyC = FR_r6, FR_poly_loC, FR_cC // poly=Z*poly_lo+c |
| nop.i 0 |
| };; |
| |
| { .mfi |
| nop.m 0 |
| fma.s1 FR_poly_hiS = FR_r, FR_poly_hiS, f0 // p_hi=r*p_hi |
| nop.i 0 |
| };; |
| |
| |
| { .mfi |
| nop.m 0 |
| fma.s1 FR_polyS = FR_polyS, f1, FR_poly_hiS |
| nop.i 0 |
| } |
| { .mfi |
| nop.m 0 |
| fma.s1 FR_polyC = FR_polyC, f1, FR_poly_hiC |
| nop.i 0 |
| };; |
| |
| .pred.rel "mutex",p7,p8 |
| { .mfi |
| nop.m 0 |
| (p7) fma.s0 FR_ResultS = FR_FirstS, f1, FR_polyS |
| nop.i 0 |
| } |
| { .mfi |
| nop.m 0 |
| (p8) fms.s0 FR_ResultS = FR_FirstS, f1, FR_polyS |
| nop.i 0 |
| };; |
| |
| .pred.rel "mutex",p9,p10 |
| { .mfi |
| nop.m 0 |
| (p9) fma.s0 FR_ResultC = FR_FirstC, f1, FR_polyC |
| nop.i 0 |
| } |
| { .mfi |
| nop.m 0 |
| (p10) fms.s0 FR_ResultC = FR_FirstC, f1, FR_polyC |
| nop.i 0 |
| };; |
| |
| .pred.rel "mutex",p11,p12 |
| { .mfi |
| nop.m 0 |
| (p11) fma.s0 FR_ResultS = FR_FirstC, f1, FR_polyC |
| nop.i 0 |
| } |
| { .mfi |
| nop.m 0 |
| (p12) fms.s0 FR_ResultS = FR_FirstC, f1, FR_polyC |
| nop.i 0 |
| };; |
| |
| .pred.rel "mutex",p13,p14 |
| { .mfi |
| nop.m 0 |
| (p13) fma.s0 FR_ResultC = FR_FirstS, f1, FR_polyS |
| nop.i 0 |
| } |
| { .mfb |
| cmp.eq p15, p0 = 0x1, GR_Cis |
| (p14) fms.s0 FR_ResultC = FR_FirstS, f1, FR_polyS |
| (p15) br.ret.sptk b0 |
| };; |
| |
| |
| { .mmb // exit for sincosl |
| stfe [sincos_pResSin] = FR_ResultS |
| stfe [sincos_pResCos] = FR_ResultC |
| br.ret.sptk b0 |
| };; |
| |
| |
| |
| |
| |
| |
| SINCOSL_NORMAL_R: |
| // |
| // Here if 2^-3 <= |r| < pi/4 |
| // THIS IS THE MAIN PATH |
| // |
| // 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 |
| nop.i 0 |
| } |
| { .mfi |
| ldfe FR_QQ_6 = [GR_ad_qq], 16 // Load QQ_6 |
| nop.f 0 |
| nop.i 0 |
| };; |
| |
| { .mmi |
| ldfe FR_PP_5 = [GR_ad_pp], 16 // Load PP_5 |
| ldfe FR_QQ_5 = [GR_ad_qq], 16 // Load QQ_5 |
| nop.i 0 |
| };; |
| |
| |
| |
| SINCOSL_NORMAL_R_0: |
| // Entry for 2^-3 < |x| < pi/4 |
| .pred.rel "mutex",p9,p10 |
| { .mmf |
| ldfe FR_C_1 = [GR_ad_pp], 16 // Load C_1 |
| ldfe FR_S_1 = [GR_ad_qq], 16 // Load S_1 |
| frcpa.s1 FR_r_hi, p6 = f1, FR_r // r_hi = frcpa(r) |
| };; |
| |
| { .mfi |
| nop.m 0 |
| fma.s1 FR_polyS = FR_rsq, FR_PP_8, FR_PP_7 // poly = rsq*PP_8+PP_7 |
| nop.i 0 |
| } |
| { .mfi |
| nop.m 0 |
| fma.s1 FR_polyC = FR_rsq, FR_QQ_8, FR_QQ_7 // poly = rsq*QQ_8+QQ_7 |
| nop.i 0 |
| };; |
| |
| { .mfi |
| nop.m 0 |
| fma.s1 FR_r_cubed = FR_r, FR_rsq, f0 // rcubed = r * rsq |
| nop.i 0 |
| };; |
| |
| |
| SINCOSL_NORMAL_R_1: |
| // Entry for pi/4 <= |x| < 2^24 |
| .pred.rel "mutex",p9,p10 |
| { .mmf |
| ldfe FR_PP_1 = [GR_ad_pp], 16 // Load PP_1_hi |
| ldfe FR_QQ_1 = [GR_ad_qq], 16 // Load QQ_1 |
| frcpa.s1 FR_r_hi, p6 = f1, FR_r_hi // r_hi = frpca(frcpa(r)) |
| };; |
| |
| { .mfi |
| ldfe FR_PP_4 = [GR_ad_pp], 16 // Load PP_4 |
| fma.s1 FR_polyS = FR_rsq, FR_polyS, FR_PP_6 // poly = rsq*poly+PP_6 |
| and GR_N_SinCos = 0x1, GR_N_Inc |
| } |
| { .mfi |
| ldfe FR_QQ_4 = [GR_ad_qq], 16 // Load QQ_4 |
| fma.s1 FR_polyC = FR_rsq, FR_polyC, FR_QQ_6 // poly = rsq*poly+QQ_6 |
| nop.i 0 |
| };; |
| |
| { .mfi |
| nop.m 0 |
| fma.s1 FR_corrS = FR_C_1, FR_rsq, f0 // corr = C_1 * rsq |
| sub GR_N_SignS = GR_N_Inc, GR_N_SinCos |
| } |
| { .mfi |
| nop.m 0 |
| fma.s1 FR_corrC = FR_S_1, FR_r_cubed, FR_r // corr = S_1 * r^3 + r |
| add GR_N_SignC = GR_N_Inc, GR_N_SinCos |
| };; |
| |
| { .mfi |
| ldfe FR_PP_3 = [GR_ad_pp], 16 // Load PP_3 |
| fma.s1 FR_r_hi_sq = FR_r_hi, FR_r_hi, f0 // r_hi_sq = r_hi * r_hi |
| tbit.z p7,p11 = GR_N_Inc, 0 |
| } |
| { .mfi |
| ldfe FR_QQ_3 = [GR_ad_qq], 16 // Load QQ_3 |
| fms.s1 FR_r_lo = FR_r, f1, FR_r_hi // r_lo = r - r_hi |
| nop.i 0 |
| };; |
| |
| { .mfi |
| ldfe FR_PP_2 = [GR_ad_pp], 16 // Load PP_2 |
| fma.s1 FR_polyS = FR_rsq, FR_polyS, FR_PP_5 // poly = rsq*poly+PP_5 |
| (p7) tbit.z.unc p9,p10 = GR_N_SignC, 1 |
| } |
| { .mfi |
| ldfe FR_QQ_2 = [GR_ad_qq], 16 // Load QQ_2 |
| fma.s1 FR_polyC = FR_rsq, FR_polyC, FR_QQ_5 // poly = rsq*poly+QQ_5 |
| nop.i 0 |
| };; |
| |
| { .mfi |
| ldfe FR_PP_1_lo = [GR_ad_pp], 16 // Load PP_1_lo |
| fma.s1 FR_corrS = FR_corrS, FR_c, FR_c // corr = corr * c + c |
| (p7) tbit.z.unc p7,p8 = GR_N_SignS, 1 |
| } |
| { .mfi |
| nop.m 0 |
| fnma.s1 FR_corrC = FR_corrC, FR_c, f0 // corr = -corr * c |
| nop.i 0 |
| };; |
| |
| { .mfi |
| nop.m 0 |
| fma.s1 FR_U_loS = FR_r, FR_r_hi, FR_r_hi_sq // U_lo = r*r_hi+r_hi_sq |
| (p11) tbit.z.unc p13,p14 = GR_N_SignC, 1 |
| } |
| { .mfi |
| nop.m 0 |
| fma.s1 FR_U_loC = FR_r_hi, f1, FR_r // U_lo = r_hi + r |
| nop.i 0 |
| };; |
| |
| { .mfi |
| nop.m 0 |
| fma.s1 FR_U_hiS = FR_r_hi, FR_r_hi_sq, f0 // U_hi = r_hi*r_hi_sq |
| (p11) tbit.z.unc p11,p12 = GR_N_SignS, 1 |
| } |
| { .mfi |
| nop.m 0 |
| fma.s1 FR_U_hiC = FR_QQ_1, FR_r_hi_sq, f1 // U_hi = QQ_1*r_hi_sq+1 |
| nop.i 0 |
| };; |
| |
| { .mfi |
| nop.m 0 |
| fma.s1 FR_polyS = FR_rsq, FR_polyS, FR_PP_4 // poly = poly*rsq+PP_4 |
| nop.i 0 |
| } |
| { .mfi |
| nop.m 0 |
| fma.s1 FR_polyC = FR_rsq, FR_polyC, FR_QQ_4 // poly = poly*rsq+QQ_4 |
| nop.i 0 |
| };; |
| |
| { .mfi |
| nop.m 0 |
| fma.s1 FR_U_loS = FR_r, FR_r, FR_U_loS // U_lo = r * r + U_lo |
| nop.i 0 |
| } |
| { .mfi |
| nop.m 0 |
| fma.s1 FR_U_loC = FR_r_lo, FR_U_loC, f0 // U_lo = r_lo * U_lo |
| nop.i 0 |
| };; |
| |
| { .mfi |
| nop.m 0 |
| fma.s1 FR_U_hiS = FR_PP_1, FR_U_hiS, f0 // U_hi = PP_1 * U_hi |
| nop.i 0 |
| };; |
| |
| { .mfi |
| nop.m 0 |
| fma.s1 FR_polyS = FR_rsq, FR_polyS, FR_PP_3 // poly = poly*rsq+PP_3 |
| nop.i 0 |
| } |
| { .mfi |
| nop.m 0 |
| fma.s1 FR_polyC = FR_rsq, FR_polyC, FR_QQ_3 // poly = poly*rsq+QQ_3 |
| nop.i 0 |
| };; |
| |
| { .mfi |
| nop.m 0 |
| fma.s1 FR_U_loS = FR_r_lo, FR_U_loS, f0 // U_lo = r_lo * U_lo |
| nop.i 0 |
| } |
| { .mfi |
| nop.m 0 |
| fma.s1 FR_U_loC = FR_QQ_1,FR_U_loC, f0 // U_lo = QQ_1 * U_lo |
| nop.i 0 |
| };; |
| |
| { .mfi |
| nop.m 0 |
| fma.s1 FR_U_hiS = FR_r, f1, FR_U_hiS // U_hi = r + U_hi |
| nop.i 0 |
| };; |
| |
| { .mfi |
| nop.m 0 |
| fma.s1 FR_polyS = FR_rsq, FR_polyS, FR_PP_2 // poly = poly*rsq+PP_2 |
| nop.i 0 |
| } |
| { .mfi |
| nop.m 0 |
| fma.s1 FR_polyC = FR_rsq, FR_polyC, FR_QQ_2 // poly = poly*rsq+QQ_2 |
| nop.i 0 |
| };; |
| |
| { .mfi |
| nop.m 0 |
| fma.s1 FR_U_loS = FR_PP_1, FR_U_loS, f0 // U_lo = PP_1 * U_lo |
| nop.i 0 |
| };; |
| |
| { .mfi |
| nop.m 0 |
| fma.s1 FR_polyS = FR_rsq, FR_polyS, FR_PP_1_lo // poly =poly*rsq+PP1lo |
| nop.i 0 |
| } |
| { .mfi |
| nop.m 0 |
| fma.s1 FR_polyC = FR_rsq, FR_polyC, f0 // poly = poly*rsq |
| nop.i 0 |
| };; |
| |
| |
| .pred.rel "mutex",p8,p14 |
| { .mfi |
| nop.m 0 |
| (p8) fms.s0 FR_U_hiS = f1, f0, FR_U_hiS |
| nop.i 0 |
| } |
| { .mfi |
| nop.m 0 |
| (p14) fms.s0 FR_U_hiS = f1, f0, FR_U_hiS |
| nop.i 0 |
| };; |
| |
| .pred.rel "mutex",p10,p12 |
| { .mfi |
| nop.m 0 |
| (p10) fms.s0 FR_U_hiC = f1, f0, FR_U_hiC |
| nop.i 0 |
| } |
| { .mfi |
| nop.m 0 |
| (p12) fms.s0 FR_U_hiC = f1, f0, FR_U_hiC |
| nop.i 0 |
| };; |
| |
| |
| { .mfi |
| nop.m 0 |
| fma.s1 FR_VS = FR_U_loS, f1, FR_corrS // V = U_lo + corr |
| nop.i 0 |
| } |
| { .mfi |
| nop.m 0 |
| fma.s1 FR_VC = FR_U_loC, f1, FR_corrC // V = U_lo + corr |
| nop.i 0 |
| };; |
| |
| { .mfi |
| nop.m 0 |
| fma.s0 FR_inexact = FR_PP_5, FR_PP_4, f0 // Dummy op to set inexact |
| nop.i 0 |
| };; |
| |
| |
| { .mfi |
| nop.m 0 |
| fma.s1 FR_polyS = FR_r_cubed, FR_polyS, f0 // poly = poly*r^3 |
| nop.i 0 |
| } |
| { .mfi |
| nop.m 0 |
| fma.s1 FR_polyC = FR_rsq, FR_polyC, f0 // poly = poly*rsq |
| nop.i 0 |
| };; |
| |
| |
| { .mfi |
| nop.m 0 |
| fma.s1 FR_VS = FR_polyS, f1, FR_VS // V = poly + V |
| nop.i 0 |
| } |
| { .mfi |
| nop.m 0 |
| fma.s1 FR_VC = FR_polyC, f1, FR_VC // V = poly + V |
| nop.i 0 |
| };; |
| |
| |
| |
| .pred.rel "mutex",p7,p8 |
| { .mfi |
| nop.m 0 |
| (p7) fma.s0 FR_ResultS = FR_U_hiS, f1, FR_VS |
| nop.i 0 |
| } |
| { .mfi |
| nop.m 0 |
| (p8) fms.s0 FR_ResultS = FR_U_hiS, f1, FR_VS |
| nop.i 0 |
| };; |
| |
| .pred.rel "mutex",p9,p10 |
| { .mfi |
| nop.m 0 |
| (p9) fma.s0 FR_ResultC = FR_U_hiC, f1, FR_VC |
| nop.i 0 |
| } |
| { .mfi |
| nop.m 0 |
| (p10) fms.s0 FR_ResultC = FR_U_hiC, f1, FR_VC |
| nop.i 0 |
| };; |
| |
| |
| |
| .pred.rel "mutex",p11,p12 |
| { .mfi |
| nop.m 0 |
| (p11) fma.s0 FR_ResultS = FR_U_hiC, f1, FR_VC |
| nop.i 0 |
| } |
| { .mfi |
| nop.m 0 |
| (p12) fms.s0 FR_ResultS = FR_U_hiC, f1, FR_VC |
| nop.i 0 |
| };; |
| |
| .pred.rel "mutex",p13,p14 |
| { .mfi |
| nop.m 0 |
| (p13) fma.s0 FR_ResultC = FR_U_hiS, f1, FR_VS |
| nop.i 0 |
| } |
| { .mfb |
| cmp.eq p15, p0 = 0x1, GR_Cis |
| (p14) fms.s0 FR_ResultC = FR_U_hiS, f1, FR_VS |
| (p15) br.ret.sptk b0 |
| };; |
| |
| { .mmb // exit for sincosl |
| stfe [sincos_pResSin] = FR_ResultS |
| stfe [sincos_pResCos] = FR_ResultC |
| br.ret.sptk b0 |
| };; |
| |
| |
| |
| |
| |
| SINCOSL_ZERO: |
| |
| { .mfi |
| nop.m 0 |
| fmerge.s FR_ResultS = FR_Input_X, FR_Input_X // If sin, result = input |
| nop.i 0 |
| } |
| { .mfb |
| cmp.eq p15, p0 = 0x1, GR_Cis |
| fma.s0 FR_ResultC = f1, f1, f0 // If cos, result=1.0 |
| (p15) br.ret.sptk b0 |
| };; |
| |
| { .mmb // exit for sincosl |
| stfe [sincos_pResSin] = FR_ResultS |
| stfe [sincos_pResCos] = FR_ResultC |
| br.ret.sptk b0 |
| };; |
| |
| |
| SINCOSL_DENORMAL: |
| { .mmb |
| getf.exp GR_signexp_x = FR_norm_x // Get sign and exponent of x |
| nop.m 999 |
| br.cond.sptk SINCOSL_COMMON2 // Return to common code |
| } |
| ;; |
| |
| |
| SINCOSL_SPECIAL: |
| // |
| // Path for Arg = +/- QNaN, SNaN, Inf |
| // Invalid can be raised. SNaNs |
| // become QNaNs |
| // |
| { .mfi |
| cmp.eq p15, p0 = 0x1, GR_Cis |
| fmpy.s0 FR_ResultS = FR_Input_X, f0 |
| nop.i 0 |
| } |
| { .mfb |
| nop.m 0 |
| fmpy.s0 FR_ResultC = FR_Input_X, f0 |
| (p15) br.ret.sptk b0 |
| };; |
| |
| { .mmb // exit for sincosl |
| stfe [sincos_pResSin] = FR_ResultS |
| stfe [sincos_pResCos] = FR_ResultC |
| br.ret.sptk b0 |
| };; |
| |
| GLOBAL_LIBM_END(__libm_sincosl) |
| |
| |
| // ******************************************************************* |
| // ******************************************************************* |
| // ******************************************************************* |
| // |
| // 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 r62-63. 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 |
| nop.f 0 |
| .save ar.pfs,GR_SAVE_PFS |
| mov GR_SAVE_PFS=ar.pfs // Save ar.pfs |
| };; |
| |
| { .mmi |
| 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 |
| setf.exp FR_Neg_Two_to_M3 = GR_exp_m2_to_m3 // Form -(2^-3) |
| nop.i 0 |
| br.call.sptk b0=__libm_pi_by_2_reduce# |
| };; |
| |
| { .mfi |
| mov GR_N_Inc = 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 0 |
| (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# |
| |
| |
| |
