mirror of
https://sourceware.org/git/glibc.git
synced 2024-12-12 22:30:12 +00:00
30891f35fa
We stopped adding "Contributed by" or similar lines in sources in 2012 in favour of git logs and keeping the Contributors section of the glibc manual up to date. Removing these lines makes the license header a bit more consistent across files and also removes the possibility of error in attribution when license blocks or files are copied across since the contributed-by lines don't actually reflect reality in those cases. Move all "Contributed by" and similar lines (Written by, Test by, etc.) into a new file CONTRIBUTED-BY to retain record of these contributions. These contributors are also mentioned in manual/contrib.texi, so we just maintain this additional record as a courtesy to the earlier developers. The following scripts were used to filter a list of files to edit in place and to clean up the CONTRIBUTED-BY file respectively. These were not added to the glibc sources because they're not expected to be of any use in future given that this is a one time task: https://gist.github.com/siddhesh/b5ecac94eabfd72ed2916d6d8157e7dc https://gist.github.com/siddhesh/15ea1f5e435ace9774f485030695ee02 Reviewed-by: Carlos O'Donell <carlos@redhat.com>
2526 lines
62 KiB
ArmAsm
2526 lines
62 KiB
ArmAsm
.file "libm_sincosl.s"
|
|
|
|
|
|
// Copyright (c) 2000 - 2004, Intel Corporation
|
|
// All rights reserved.
|
|
//
|
|
//
|
|
// 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 unnecessary.
|
|
//
|
|
// 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)
|
|
libm_alias_ldouble_other (__sincos, sincos)
|
|
|
|
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#
|