4 // Copyright (c) 2000 - 2004, Intel Corporation
5 // All rights reserved.
7 // Contributed 2000 by the Intel Numerics Group, Intel Corporation
9 // Redistribution and use in source and binary forms, with or without
10 // modification, are permitted provided that the following conditions are
13 // * Redistributions of source code must retain the above copyright
14 // notice, this list of conditions and the following disclaimer.
16 // * Redistributions in binary form must reproduce the above copyright
17 // notice, this list of conditions and the following disclaimer in the
18 // documentation and/or other materials provided with the distribution.
20 // * The name of Intel Corporation may not be used to endorse or promote
21 // products derived from this software without specific prior written
24 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
25 // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
26 // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
27 // A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL INTEL OR ITS
28 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
29 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
30 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
31 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY
32 // OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY OR TORT (INCLUDING
33 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
34 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 // Intel Corporation is the author of this code, and requests that all
37 // problem reports or change requests be submitted to it directly at
38 // http://www.intel.com/software/products/opensource/libraries/num.htm.
40 //*********************************************************************
43 // 02/02/00 (hand-optimized)
44 // 04/04/00 Unwind support added
45 // 07/30/01 Improved speed on all paths
46 // 08/20/01 Fixed bundling typo
47 // 05/13/02 Changed interface to __libm_pi_by_2_reduce
48 // 02/10/03 Reordered header: .section, .global, .proc, .align;
49 // used data8 for long double table values
50 // 10/13/03 Corrected final .endp name to match .proc
51 // 10/26/04 Avoided using r14-31 as scratch so not clobbered by dynamic loader
53 //*********************************************************************
55 // Function: Combined sinl(x) and cosl(x), where
57 // sinl(x) = sine(x), for double-extended precision x values
58 // cosl(x) = cosine(x), for double-extended precision x values
60 //*********************************************************************
64 // Floating-Point Registers: f8 (Input and Return Value)
67 // General Purpose Registers:
70 // Predicate Registers: p6-p13
72 //*********************************************************************
74 // IEEE Special Conditions:
76 // Denormal fault raised on denormal inputs
77 // Overflow exceptions do not occur
78 // Underflow exceptions raised when appropriate for sin
79 // (No specialized error handling for this routine)
80 // Inexact raised when appropriate by algorithm
91 //*********************************************************************
93 // Mathematical Description
94 // ========================
96 // The computation of FSIN and FCOS is best handled in one piece of
97 // code. The main reason is that given any argument Arg, computation
98 // of trigonometric functions first calculate N and an approximation
101 // Arg = N pi/2 + alpha, |alpha| <= pi/4.
105 // cosl( Arg ) = sinl( (N+1) pi/2 + alpha ),
107 // therefore, the code for computing sine will produce cosine as long
108 // as 1 is added to N immediately after the argument reduction
116 // Arg = M pi/2 + alpha, |alpha| <= pi/4,
118 // let I = M mod 4, or I be the two lsb of M when M is represented
119 // as 2's complement. I = [i_0 i_1]. Then
121 // sinl( Arg ) = (-1)^i_0 sinl( alpha ) if i_1 = 0,
122 // = (-1)^i_0 cosl( alpha ) if i_1 = 1.
126 // sin ((-pi/2 + alpha) = (-1) cos (alpha)
128 // sin (alpha) = sin (alpha)
130 // sin (pi/2 + alpha) = cos (alpha)
132 // sin (pi + alpha) = (-1) sin (alpha)
134 // sin ((3/2)pi + alpha) = (-1) cos (alpha)
136 // The value of alpha is obtained by argument reduction and
137 // represented by two working precision numbers r and c where
139 // alpha = r + c accurately.
141 // The reduction method is described in a previous write up.
142 // The argument reduction scheme identifies 4 cases. For Cases 2
143 // and 4, because |alpha| is small, sinl(r+c) and cosl(r+c) can be
144 // computed very easily by 2 or 3 terms of the Taylor series
145 // expansion as follows:
150 // sinl(r + c) = r + c - r^3/6 accurately
151 // cosl(r + c) = 1 - 2^(-67) accurately
156 // sinl(r + c) = r + c - r^3/6 + r^5/120 accurately
157 // cosl(r + c) = 1 - r^2/2 + r^4/24 accurately
159 // The only cases left are Cases 1 and 3 of the argument reduction
160 // procedure. These two cases will be merged since after the
161 // argument is reduced in either cases, we have the reduced argument
162 // represented as r + c and that the magnitude |r + c| is not small
163 // enough to allow the usage of a very short approximation.
165 // The required calculation is either
167 // sinl(r + c) = sinl(r) + correction, or
168 // cosl(r + c) = cosl(r) + correction.
172 // sinl(r + c) = sinl(r) + c sin'(r) + O(c^2)
173 // = sinl(r) + c cos (r) + O(c^2)
174 // = sinl(r) + c(1 - r^2/2) accurately.
177 // cosl(r + c) = cosl(r) - c sinl(r) + O(c^2)
178 // = cosl(r) - c(r - r^3/6) accurately.
180 // We therefore concentrate on accurately calculating sinl(r) and
181 // cosl(r) for a working-precision number r, |r| <= pi/4 to within
184 // The greatest challenge of this task is that the second terms of
187 // r - r^3/3! + r^r/5! - ...
191 // 1 - r^2/2! + r^4/4! - ...
193 // are not very small when |r| is close to pi/4 and the rounding
194 // errors will be a concern if simple polynomial accumulation is
195 // used. When |r| < 2^-3, however, the second terms will be small
196 // enough (6 bits or so of right shift) that a normal Horner
197 // recurrence suffices. Hence there are two cases that we consider
198 // in the accurate computation of sinl(r) and cosl(r), |r| <= pi/4.
200 // Case small_r: |r| < 2^(-3)
201 // --------------------------
203 // Since Arg = M pi/4 + r + c accurately, and M mod 4 is [i_0 i_1],
206 // sinl(Arg) = (-1)^i_0 * sinl(r + c) if i_1 = 0
207 // = (-1)^i_0 * cosl(r + c) if i_1 = 1
209 // can be accurately approximated by
211 // sinl(Arg) = (-1)^i_0 * [sinl(r) + c] if i_1 = 0
212 // = (-1)^i_0 * [cosl(r) - c*r] if i_1 = 1
214 // because |r| is small and thus the second terms in the correction
217 // Finally, sinl(r) and cosl(r) are approximated by polynomials of
220 // sinl(r) = r + S_1 r^3 + S_2 r^5 + ... + S_5 r^11
221 // cosl(r) = 1 + C_1 r^2 + C_2 r^4 + ... + C_5 r^10
223 // We can make use of predicates to selectively calculate
224 // sinl(r) or cosl(r) based on i_1.
226 // Case normal_r: 2^(-3) <= |r| <= pi/4
227 // ------------------------------------
229 // This case is more likely than the previous one if one considers
230 // r to be uniformly distributed in [-pi/4 pi/4]. Again,
232 // sinl(Arg) = (-1)^i_0 * sinl(r + c) if i_1 = 0
233 // = (-1)^i_0 * cosl(r + c) if i_1 = 1.
235 // Because |r| is now larger, we need one extra term in the
236 // correction. sinl(Arg) can be accurately approximated by
238 // sinl(Arg) = (-1)^i_0 * [sinl(r) + c(1-r^2/2)] if i_1 = 0
239 // = (-1)^i_0 * [cosl(r) - c*r*(1 - r^2/6)] i_1 = 1.
241 // Finally, sinl(r) and cosl(r) are approximated by polynomials of
244 // sinl(r) = r + PP_1_hi r^3 + PP_1_lo r^3 +
245 // PP_2 r^5 + ... + PP_8 r^17
247 // cosl(r) = 1 + QQ_1 r^2 + QQ_2 r^4 + ... + QQ_8 r^16
249 // where PP_1_hi is only about 16 bits long and QQ_1 is -1/2.
250 // The crux in accurate computation is to calculate
252 // r + PP_1_hi r^3 or 1 + QQ_1 r^2
254 // accurately as two pieces: U_hi and U_lo. The way to achieve this
255 // is to obtain r_hi as a 10 sig. bit number that approximates r to
256 // roughly 8 bits or so of accuracy. (One convenient way is
258 // r_hi := frcpa( frcpa( r ) ).)
262 // r + PP_1_hi r^3 = r + PP_1_hi r_hi^3 +
263 // PP_1_hi (r^3 - r_hi^3)
264 // = [r + PP_1_hi r_hi^3] +
265 // [PP_1_hi (r - r_hi)
266 // (r^2 + r_hi r + r_hi^2) ]
269 // Since r_hi is only 10 bit long and PP_1_hi is only 16 bit long,
270 // PP_1_hi * r_hi^3 is only at most 46 bit long and thus computed
271 // exactly. Furthermore, r and PP_1_hi r_hi^3 are of opposite sign
272 // and that there is no more than 8 bit shift off between r and
273 // PP_1_hi * r_hi^3. Hence the sum, U_hi, is representable and thus
274 // calculated without any error. Finally, the fact that
276 // |U_lo| <= 2^(-8) |U_hi|
278 // says that U_hi + U_lo is approximating r + PP_1_hi r^3 to roughly
279 // 8 extra bits of accuracy.
283 // 1 + QQ_1 r^2 = [1 + QQ_1 r_hi^2] +
284 // [QQ_1 (r - r_hi)(r + r_hi)]
287 // Summarizing, we calculate r_hi = frcpa( frcpa( r ) ).
291 // U_hi := r + PP_1_hi * r_hi^3
292 // U_lo := PP_1_hi * (r - r_hi) * (r^2 + r*r_hi + r_hi^2)
293 // poly := PP_1_lo r^3 + PP_2 r^5 + ... + PP_8 r^17
294 // correction := c * ( 1 + C_1 r^2 )
298 // U_hi := 1 + QQ_1 * r_hi * r_hi
299 // U_lo := QQ_1 * (r - r_hi) * (r + r_hi)
300 // poly := QQ_2 * r^4 + QQ_3 * r^6 + ... + QQ_8 r^16
301 // correction := -c * r * (1 + S_1 * r^2)
307 // V := poly + ( U_lo + correction )
309 // / U_hi + V if i_0 = 0
311 // \ (-U_hi) - V if i_0 = 1
313 // It is important that in the last step, negation of U_hi is
314 // performed prior to the subtraction which is to be performed in
315 // the user-set rounding mode.
318 // Algorithmic Description
319 // =======================
321 // The argument reduction algorithm is tightly integrated into FSIN
322 // and FCOS which share the same code. The following is complete and
323 // self-contained. The argument reduction description given
324 // previously is repeated below.
327 // Step 0. Initialization.
329 // If FSIN is invoked, set N_inc := 0; else if FCOS is invoked,
332 // Step 1. Check for exceptional and special cases.
334 // * If Arg is +-0, +-inf, NaN, NaT, go to Step 10 for special
336 // * If |Arg| < 2^24, go to Step 2 for reduction of moderate
337 // arguments. This is the most likely case.
338 // * If |Arg| < 2^63, go to Step 8 for pre-reduction of large
340 // * If |Arg| >= 2^63, go to Step 10 for special handling.
342 // Step 2. Reduction of moderate arguments.
344 // If |Arg| < pi/4 ...quick branch
345 // N_fix := N_inc (integer)
348 // Branch to Step 4, Case_1_complete
349 // Else ...cf. argument reduction
350 // N := Arg * two_by_PI (fp)
351 // N_fix := fcvt.fx( N ) (int)
352 // N := fcvt.xf( N_fix )
353 // N_fix := N_fix + N_inc
354 // s := Arg - N * P_1 (first piece of pi/2)
355 // w := -N * P_2 (second piece of pi/2)
358 // go to Step 3, Case_1_reduce
360 // go to Step 7, Case_2_reduce
364 // Step 3. Case_1_reduce.
367 // c := (s - r) + w ...observe order
369 // Step 4. Case_1_complete
371 // ...At this point, the reduced argument alpha is
372 // ...accurately represented as r + c.
373 // If |r| < 2^(-3), go to Step 6, small_r.
377 // Let [i_0 i_1] by the 2 lsb of N_fix.
379 // r_hi := frcpa( frcpa( r ) )
383 // poly := r*FR_rsq*(PP_1_lo + FR_rsq*(PP_2 + ... FR_rsq*PP_8))
384 // U_hi := r + PP_1_hi*r_hi*r_hi*r_hi ...any order
385 // U_lo := PP_1_hi*r_lo*(r*r + r*r_hi + r_hi*r_hi)
386 // correction := c + c*C_1*FR_rsq ...any order
388 // poly := FR_rsq*FR_rsq*(QQ_2 + FR_rsq*(QQ_3 + ... + FR_rsq*QQ_8))
389 // U_hi := 1 + QQ_1 * r_hi * r_hi ...any order
390 // U_lo := QQ_1 * r_lo * (r + r_hi)
391 // correction := -c*(r + S_1*FR_rsq*r) ...any order
394 // V := poly + (U_lo + correction) ...observe order
396 // result := (i_0 == 0? 1.0 : -1.0)
398 // Last instruction in user-set rounding mode
400 // result := (i_0 == 0? result*U_hi + V :
407 // ...Use flush to zero mode without causing exception
408 // Let [i_0 i_1] be the two lsb of N_fix.
413 // z := FR_rsq*FR_rsq; z := FR_rsq*z *r
414 // poly_lo := S_3 + FR_rsq*(S_4 + FR_rsq*S_5)
415 // poly_hi := r*FR_rsq*(S_1 + FR_rsq*S_2)
419 // z := FR_rsq*FR_rsq; z := FR_rsq*z
420 // poly_lo := C_3 + FR_rsq*(C_4 + FR_rsq*C_5)
421 // poly_hi := FR_rsq*(C_1 + FR_rsq*C_2)
422 // correction := -c*r
426 // poly := poly_hi + (z * poly_lo + correction)
428 // If i_0 = 1, result := -result
430 // Last operation. Perform in user-set rounding mode
432 // result := (i_0 == 0? result + poly :
436 // Step 7. Case_2_reduce.
438 // ...Refer to the write up for argument reduction for
439 // ...rationale. The reduction algorithm below is taken from
440 // ...argument reduction description and integrated this.
443 // U_1 := N*P_2 + w ...FMA
444 // U_2 := (N*P_2 - U_1) + w ...2 FMA
445 // ...U_1 + U_2 is N*(P_2+P_3) accurately
448 // c := ( (s - r) - U_1 ) - U_2
450 // ...The mathematical sum r + c approximates the reduced
451 // ...argument accurately. Note that although compared to
452 // ...Case 1, this case requires much more work to reduce
453 // ...the argument, the subsequent calculation needed for
454 // ...any of the trigonometric function is very little because
455 // ...|alpha| < 1.01*2^(-33) and thus two terms of the
456 // ...Taylor series expansion suffices.
459 // poly := c + S_1 * r * r * r ...any order
466 // If i_0 = 1, result := -result
468 // Last operation. Perform in user-set rounding mode
470 // result := (i_0 == 0? result + poly :
476 // Step 8. Pre-reduction of large arguments.
478 // ...Again, the following reduction procedure was described
479 // ...in the separate write up for argument reduction, which
480 // ...is tightly integrated here.
482 // N_0 := Arg * Inv_P_0
483 // N_0_fix := fcvt.fx( N_0 )
484 // N_0 := fcvt.xf( N_0_fix)
486 // Arg' := Arg - N_0 * P_0
488 // N := Arg' * two_by_PI
489 // N_fix := fcvt.fx( N )
490 // N := fcvt.xf( N_fix )
491 // N_fix := N_fix + N_inc
493 // s := Arg' - N * P_1
502 // Step 9. Case_4_reduce.
504 // ...first obtain N_0*d_1 and -N*P_2 accurately
505 // U_hi := N_0 * d_1 V_hi := -N*P_2
506 // U_lo := N_0 * d_1 - U_hi V_lo := -N*P_2 - U_hi ...FMAs
508 // ...compute the contribution from N_0*d_1 and -N*P_3
511 // t := U_lo + V_lo + w ...any order
513 // ...at this point, the mathematical value
514 // ...s + U_hi + V_hi + t approximates the true reduced argument
515 // ...accurately. Just need to compute this accurately.
517 // ...Calculate U_hi + V_hi accurately:
519 // if |U_hi| >= |V_hi| then
520 // a := (U_hi - A) + V_hi
522 // a := (V_hi - A) + U_hi
524 // ...order in computing "a" must be observed. This branch is
525 // ...best implemented by predicates.
526 // ...A + a is U_hi + V_hi accurately. Moreover, "a" is
527 // ...much smaller than A: |a| <= (1/2)ulp(A).
529 // ...Just need to calculate s + A + a + t
530 // C_hi := s + A t := t + a
531 // C_lo := (s - C_hi) + A
534 // ...Final steps for reduction
536 // c := (C_hi - r) + C_lo
538 // ...At this point, we have r and c
539 // ...And all we need is a couple of terms of the corresponding
543 // poly := c + r*FR_rsq*(S_1 + FR_rsq*S_2)
546 // poly := FR_rsq*(C_1 + FR_rsq*C_2)
550 // If i_0 = 1, result := -result
552 // Last operation. Perform in user-set rounding mode
554 // result := (i_0 == 0? result + poly :
558 // Large Arguments: For arguments above 2**63, a Payne-Hanek
559 // style argument reduction is used and pi_by_2 reduce is called.
566 LOCAL_OBJECT_START(FSINCOSL_CONSTANTS)
569 data8 0xA2F9836E4E44152A, 0x00003FFE // Inv_pi_by_2
570 data8 0xC84D32B0CE81B9F1, 0x00004016 // P_0
571 data8 0xC90FDAA22168C235, 0x00003FFF // P_1
572 data8 0xECE675D1FC8F8CBB, 0x0000BFBD // P_2
573 data8 0xB7ED8FBBACC19C60, 0x0000BF7C // P_3
574 data8 0x8D848E89DBD171A1, 0x0000BFBF // d_1
575 data8 0xD5394C3618A66F8E, 0x0000BF7C // d_2
576 LOCAL_OBJECT_END(FSINCOSL_CONSTANTS)
578 LOCAL_OBJECT_START(sincosl_table_d)
579 data8 0xC90FDAA22168C234, 0x00003FFE // pi_by_4
580 data8 0xA397E5046EC6B45A, 0x00003FE7 // Inv_P_0
581 data4 0x3E000000, 0xBE000000 // 2^-3 and -2^-3
582 data4 0x2F000000, 0xAF000000 // 2^-33 and -2^-33
583 data4 0x9E000000, 0x00000000 // -2^-67
584 data4 0x00000000, 0x00000000 // pad
585 LOCAL_OBJECT_END(sincosl_table_d)
587 LOCAL_OBJECT_START(sincosl_table_pp)
588 data8 0xCC8ABEBCA21C0BC9, 0x00003FCE // PP_8
589 data8 0xD7468A05720221DA, 0x0000BFD6 // PP_7
590 data8 0xB092382F640AD517, 0x00003FDE // PP_6
591 data8 0xD7322B47D1EB75A4, 0x0000BFE5 // PP_5
592 data8 0xFFFFFFFFFFFFFFFE, 0x0000BFFD // C_1
593 data8 0xAAAA000000000000, 0x0000BFFC // PP_1_hi
594 data8 0xB8EF1D2ABAF69EEA, 0x00003FEC // PP_4
595 data8 0xD00D00D00D03BB69, 0x0000BFF2 // PP_3
596 data8 0x8888888888888962, 0x00003FF8 // PP_2
597 data8 0xAAAAAAAAAAAB0000, 0x0000BFEC // PP_1_lo
598 LOCAL_OBJECT_END(sincosl_table_pp)
600 LOCAL_OBJECT_START(sincosl_table_qq)
601 data8 0xD56232EFC2B0FE52, 0x00003FD2 // QQ_8
602 data8 0xC9C99ABA2B48DCA6, 0x0000BFDA // QQ_7
603 data8 0x8F76C6509C716658, 0x00003FE2 // QQ_6
604 data8 0x93F27DBAFDA8D0FC, 0x0000BFE9 // QQ_5
605 data8 0xAAAAAAAAAAAAAAAA, 0x0000BFFC // S_1
606 data8 0x8000000000000000, 0x0000BFFE // QQ_1
607 data8 0xD00D00D00C6E5041, 0x00003FEF // QQ_4
608 data8 0xB60B60B60B607F60, 0x0000BFF5 // QQ_3
609 data8 0xAAAAAAAAAAAAAA9B, 0x00003FFA // QQ_2
610 LOCAL_OBJECT_END(sincosl_table_qq)
612 LOCAL_OBJECT_START(sincosl_table_c)
613 data8 0xFFFFFFFFFFFFFFFE, 0x0000BFFD // C_1
614 data8 0xAAAAAAAAAAAA719F, 0x00003FFA // C_2
615 data8 0xB60B60B60356F994, 0x0000BFF5 // C_3
616 data8 0xD00CFFD5B2385EA9, 0x00003FEF // C_4
617 data8 0x93E4BD18292A14CD, 0x0000BFE9 // C_5
618 LOCAL_OBJECT_END(sincosl_table_c)
620 LOCAL_OBJECT_START(sincosl_table_s)
621 data8 0xAAAAAAAAAAAAAAAA, 0x0000BFFC // S_1
622 data8 0x88888888888868DB, 0x00003FF8 // S_2
623 data8 0xD00D00D0055EFD4B, 0x0000BFF2 // S_3
624 data8 0xB8EF1C5D839730B9, 0x00003FEC // S_4
625 data8 0xD71EA3A4E5B3F492, 0x0000BFE5 // S_5
626 data4 0x38800000, 0xB8800000 // two**-14 and -two**-14
627 LOCAL_OBJECT_END(sincosl_table_s)
636 FR_inv_pi_2to63 = f10
640 FR_N_float_signif = f14
644 FR_Neg_Two_to_M14 = f36
646 FR_Neg_Two_to_M33 = f38
647 FR_Neg_Two_to_M67 = f39
710 FR_Neg_Two_to_M3 = f93
744 // Added for unwind support
753 GLOBAL_IEEE754_ENTRY(sinl)
755 alloc r32 = ar.pfs,0,27,2,0
756 movl GR_sig_inv_pi = 0xa2f9836e4e44152a // significand of 1/pi
759 mov GR_Sin_or_Cos = 0x0
760 movl GR_rshf_2to64 = 0x47e8000000000000 // 1.1000 2^(63+64)
765 addl GR_ad_p = @ltoff(FSINCOSL_CONSTANTS#), gp
766 fclass.m p6, p0 = FR_Input_X, 0x1E3 // Test x natval, nan, inf
767 mov GR_exp_2_to_m3 = 0xffff - 3 // Exponent of 2^-3
771 fnorm.s1 FR_norm_x = FR_Input_X // Normalize x
772 br.cond.sptk SINCOSL_CONTINUE
776 GLOBAL_IEEE754_END(sinl)
778 GLOBAL_IEEE754_ENTRY(cosl)
780 alloc r32 = ar.pfs,0,27,2,0
781 movl GR_sig_inv_pi = 0xa2f9836e4e44152a // significand of 1/pi
784 mov GR_Sin_or_Cos = 0x1
785 movl GR_rshf_2to64 = 0x47e8000000000000 // 1.1000 2^(63+64)
790 addl GR_ad_p = @ltoff(FSINCOSL_CONSTANTS#), gp
791 fclass.m p6, p0 = FR_Input_X, 0x1E3 // Test x natval, nan, inf
792 mov GR_exp_2_to_m3 = 0xffff - 3 // Exponent of 2^-3
796 fnorm.s1 FR_norm_x = FR_Input_X // Normalize x
803 setf.sig FR_inv_pi_2to63 = GR_sig_inv_pi // Form 1/pi * 2^63
805 mov GR_exp_2tom64 = 0xffff - 64 // Scaling constant to compute N
808 setf.d FR_rshf_2to64 = GR_rshf_2to64 // Form const 1.1000 * 2^(63+64)
809 movl GR_rshf = 0x43e8000000000000 // Form const 1.1000 * 2^63
814 ld8 GR_ad_p = [GR_ad_p] // Point to Inv_pi_by_2
815 fclass.m p7, p0 = FR_Input_X, 0x0b // Test x denormal
821 getf.exp GR_signexp_x = FR_Input_X // Get sign and exponent of x
822 fclass.m p10, p0 = FR_Input_X, 0x007 // Test x zero
826 mov GR_exp_mask = 0x1ffff // Exponent mask
828 (p6) br.cond.spnt SINCOSL_SPECIAL // Branch if x natval, nan, inf
833 setf.exp FR_2tom64 = GR_exp_2tom64 // Form 2^-64 for scaling N_float
835 add GR_ad_d = 0x70, GR_ad_p // Point to constant table d
838 setf.d FR_rshf = GR_rshf // Form right shift const 1.1000 * 2^63
839 mov GR_exp_m2_to_m3 = 0x2fffc // Form -(2^-3)
840 (p7) br.cond.spnt SINCOSL_DENORMAL // Branch if x denormal
846 and GR_exp_x = GR_exp_mask, GR_signexp_x // Get exponent of x
847 fclass.nm p8, p0 = FR_Input_X, 0x1FF // Test x unsupported type
848 mov GR_exp_2_to_63 = 0xffff + 63 // Exponent of 2^63
851 add GR_ad_pp = 0x40, GR_ad_d // Point to constant table pp
852 mov GR_exp_2_to_24 = 0xffff + 24 // Exponent of 2^24
853 (p10) br.cond.spnt SINCOSL_ZERO // Branch if x zero
858 ldfe FR_Inv_pi_by_2 = [GR_ad_p], 16 // Load 2/pi
859 fcmp.eq.s0 p15, p0 = FR_Input_X, f0 // Dummy to set denormal
860 add GR_ad_qq = 0xa0, GR_ad_pp // Point to constant table qq
863 ldfe FR_Pi_by_4 = [GR_ad_d], 16 // Load pi/4 for range test
865 cmp.ge p10,p0 = GR_exp_x, GR_exp_2_to_63 // Is |x| >= 2^63
870 ldfe FR_P_0 = [GR_ad_p], 16 // Load P_0 for pi/4 <= |x| < 2^63
871 fmerge.s FR_abs_x = f1, FR_norm_x // |x|
872 add GR_ad_c = 0x90, GR_ad_qq // Point to constant table c
875 ldfe FR_Inv_P_0 = [GR_ad_d], 16 // Load 1/P_0 for pi/4 <= |x| < 2^63
877 cmp.ge p7,p0 = GR_exp_x, GR_exp_2_to_24 // Is |x| >= 2^24
882 ldfe FR_P_1 = [GR_ad_p], 16 // Load P_1 for pi/4 <= |x| < 2^63
884 add GR_ad_s = 0x50, GR_ad_c // Point to constant table s
887 ldfe FR_PP_8 = [GR_ad_pp], 16 // Load PP_8 for 2^-3 < |r| < pi/4
894 ldfe FR_P_2 = [GR_ad_p], 16 // Load P_2 for pi/4 <= |x| < 2^63
896 add GR_ad_ce = 0x40, GR_ad_c // Point to end of constant table c
899 ldfe FR_QQ_8 = [GR_ad_qq], 16 // Load QQ_8 for 2^-3 < |r| < pi/4
906 ldfe FR_QQ_7 = [GR_ad_qq], 16 // Load QQ_7 for 2^-3 < |r| < pi/4
907 fma.s1 FR_N_float_signif = FR_Input_X, FR_inv_pi_2to63, FR_rshf_2to64
908 add GR_ad_se = 0x40, GR_ad_s // Point to end of constant table s
911 ldfe FR_PP_7 = [GR_ad_pp], 16 // Load PP_7 for 2^-3 < |r| < pi/4
912 mov GR_ad_s1 = GR_ad_s // Save pointer to S_1
913 (p10) br.cond.spnt SINCOSL_ARG_TOO_LARGE // Branch if |x| >= 2^63
914 // Use Payne-Hanek Reduction
919 ldfe FR_P_3 = [GR_ad_p], 16 // Load P_3 for pi/4 <= |x| < 2^63
920 fmerge.se FR_r = FR_norm_x, FR_norm_x // r = x, in case |x| < pi/4
921 add GR_ad_m14 = 0x50, GR_ad_s // Point to constant table m14
924 ldfps FR_Two_to_M3, FR_Neg_Two_to_M3 = [GR_ad_d], 8
925 fma.s1 FR_rsq = FR_norm_x, FR_norm_x, f0 // rsq = x*x, in case |x| < pi/4
926 (p7) br.cond.spnt SINCOSL_LARGER_ARG // Branch if 2^24 <= |x| < 2^63
932 ldfe FR_PP_6 = [GR_ad_pp], 16 // Load PP_6 for normal path
933 ldfe FR_QQ_6 = [GR_ad_qq], 16 // Load QQ_6 for normal path
934 fmerge.se FR_c = f0, f0 // c = 0 in case |x| < pi/4
939 ldfe FR_PP_5 = [GR_ad_pp], 16 // Load PP_5 for normal path
940 ldfe FR_QQ_5 = [GR_ad_qq], 16 // Load QQ_5 for normal path
945 // Here if 0 < |x| < 2^24
947 ldfe FR_S_5 = [GR_ad_se], -16 // Load S_5 if i_1=0
948 fcmp.lt.s1 p6, p7 = FR_abs_x, FR_Pi_by_4 // Test |x| < pi/4
952 ldfe FR_C_5 = [GR_ad_ce], -16 // Load C_5 if i_1=1
953 fms.s1 FR_N_float = FR_N_float_signif, FR_2tom64, FR_rshf
959 ldfe FR_S_4 = [GR_ad_se], -16 // Load S_4 if i_1=0
960 ldfe FR_C_4 = [GR_ad_ce], -16 // Load C_4 if i_1=1
967 // Check if Arg < pi/4
970 // Case 2: Convert integer N_fix back to normalized floating-point value.
971 // Case 1: p8 is only affected when p6 is set
974 // Grab the integer part of N and call it N_fix
977 (p7) ldfps FR_Two_to_M33, FR_Neg_Two_to_M33 = [GR_ad_d], 8
978 (p6) fma.s1 FR_r_cubed = FR_r, FR_rsq, f0 // r^3 if |x| < pi/4
979 (p6) mov GR_N_Inc = GR_Sin_or_Cos // N_Inc if |x| < pi/4
983 // If |x| < pi/4, r = x and c = 0
984 // lf |x| < pi/4, is x < 2**(-3).
988 (p7) getf.sig GR_N_Inc = FR_N_float_signif
989 (p6) cmp.lt.unc p8,p0 = GR_exp_x, GR_exp_2_to_m3 // Is |x| < 2^-3
990 (p6) tbit.z p9,p10 = GR_N_Inc, 0 // p9 if i_1=0, N mod 4 = 0,1
991 // p10 if i_1=1, N mod 4 = 2,3
996 // lf |x| < pi/4, is -2**(-3)< x < 2**(-3) - set p8.
998 // Create the right N for |x| < pi/4 and otherwise
999 // Case 2: Place integer part of N in GP register
1005 (p8) br.cond.spnt SINCOSL_SMALL_R_0 // Branch if 0 < |x| < 2^-3
1006 (p6) br.cond.spnt SINCOSL_NORMAL_R_0 // Branch if 2^-3 <= |x| < pi/4
1010 // Here if pi/4 <= |x| < 2^24
1012 ldfs FR_Neg_Two_to_M67 = [GR_ad_d], 8 // Load -2^-67
1013 fnma.s1 FR_s = FR_N_float, FR_P_1, FR_Input_X // s = -N * P_1 + Arg
1014 add GR_N_Inc = GR_N_Inc, GR_Sin_or_Cos // Adjust N_Inc for sin/cos
1018 fma.s1 FR_w = FR_N_float, FR_P_2, f0 // w = N * P_2
1025 fms.s1 FR_r = FR_s, f1, FR_w // r = s - w, assume |s| >= 2^-33
1026 tbit.z p9,p10 = GR_N_Inc, 0 // p9 if i_1=0, N mod 4 = 0,1
1027 // p10 if i_1=1, N mod 4 = 2,3
1033 fcmp.lt.s1 p7, p6 = FR_s, FR_Two_to_M33
1040 (p7) fcmp.gt.s1 p7, p6 = FR_s, FR_Neg_Two_to_M33 // p6 if |s| >= 2^-33, else p7
1047 fms.s1 FR_c = FR_s, f1, FR_r // c = s - r, for |s| >= 2^-33
1052 fma.s1 FR_rsq = FR_r, FR_r, f0 // rsq = r * r, for |s| >= 2^-33
1059 (p7) fma.s1 FR_w = FR_N_float, FR_P_3, f0
1065 (p9) ldfe FR_C_1 = [GR_ad_pp], 16 // Load C_1 if i_1=0
1066 (p10) ldfe FR_S_1 = [GR_ad_qq], 16 // Load S_1 if i_1=1
1067 frcpa.s1 FR_r_hi, p15 = f1, FR_r // r_hi = frcpa(r)
1073 (p6) fcmp.lt.unc.s1 p8, p13 = FR_r, FR_Two_to_M3 // If big s, test r with 2^-3
1080 (p7) fma.s1 FR_U_1 = FR_N_float, FR_P_2, FR_w
1086 // For big s: r = s - w: No futher reduction is necessary
1087 // For small s: w = N * P_3 (change sign) More reduction
1091 (p8) fcmp.gt.s1 p8, p13 = FR_r, FR_Neg_Two_to_M3 // If big s, p8 if |r| < 2^-3
1097 (p9) fma.s1 FR_poly = FR_rsq, FR_PP_8, FR_PP_7 // poly = rsq*PP_8+PP_7 if i_1=0
1102 (p10) fma.s1 FR_poly = FR_rsq, FR_QQ_8, FR_QQ_7 // poly = rsq*QQ_8+QQ_7 if i_1=1
1109 (p7) fms.s1 FR_r = FR_s, f1, FR_U_1
1116 (p6) fma.s1 FR_r_cubed = FR_r, FR_rsq, f0 // rcubed = r * rsq
1123 // For big s: Is |r| < 2**(-3)?
1124 // For big s: c = S - r
1125 // For small s: U_1 = N * P_2 + w
1127 // If p8 is set, prepare to branch to Small_R.
1128 // If p9 is set, prepare to branch to Normal_R.
1129 // For big s, r is complete here.
1132 // For big s: c = c + w (w has not been negated.)
1133 // For small s: r = S - U_1
1136 (p6) fms.s1 FR_c = FR_c, f1, FR_w
1141 (p8) br.cond.spnt SINCOSL_SMALL_R_1 // Branch if |s|>=2^-33, |r| < 2^-3,
1142 // and pi/4 <= |x| < 2^24
1143 (p13) br.cond.sptk SINCOSL_NORMAL_R_1 // Branch if |s|>=2^-33, |r| >= 2^-3,
1144 // and pi/4 <= |x| < 2^24
1150 // Here if |s| < 2^-33, and pi/4 <= |x| < 2^24
1153 fms.s1 FR_U_2 = FR_N_float, FR_P_2, FR_U_1
1165 // Get [i_0,i_1] - two lsb of N_fix_gr.
1166 // Do dummy fmpy so inexact is always set.
1168 tbit.z p9,p10 = GR_N_Inc, 0 // p9 if i_1=0, N mod 4 = 0,1
1169 // p10 if i_1=1, N mod 4 = 2,3
1174 // For small s: U_2 = N * P_2 - U_1
1175 // S_1 stored constant - grab the one stored with the
1179 ldfe FR_S_1 = [GR_ad_s1], 16
1181 // Check if i_1 and i_0 != 0
1183 (p10) fma.s1 FR_poly = f0, f1, FR_Neg_Two_to_M67
1184 tbit.z p11,p12 = GR_N_Inc, 1 // p11 if i_0=0, N mod 4 = 0,2
1185 // p12 if i_0=1, N mod 4 = 1,3
1191 fms.s1 FR_s = FR_s, f1, FR_r
1201 fma.s1 FR_rsq = FR_r, FR_r, f0
1206 fma.s1 FR_U_2 = FR_U_2, f1, FR_w
1211 fmerge.se FR_tmp_result = FR_r, FR_r
1216 (p10) fma.s1 FR_tmp_result = f0, f1, f1
1223 // Save r as the result.
1225 fms.s1 FR_c = FR_s, f1, FR_U_1
1231 // if ( i_1 ==0) poly = c + S_1*r*r*r
1234 (p12) fnma.s1 FR_tmp_result = FR_tmp_result, f1, f0
1239 fma.s1 FR_r = FR_S_1, FR_r, f0
1244 fma.s0 FR_S_1 = FR_S_1, FR_S_1, f0
1250 // If i_1 != 0, poly = 2**(-67)
1252 fms.s1 FR_c = FR_c, f1, FR_U_2
1260 (p9) fma.s1 FR_poly = FR_r, FR_rsq, FR_c
1266 // i_0 != 0, so Result = -Result
1268 (p11) fma.s0 FR_Result = FR_tmp_result, f1, FR_poly
1273 (p12) fms.s0 FR_Result = FR_tmp_result, f1, FR_poly
1275 // if (i_0 == 0), Result = Result + poly
1276 // else Result = Result - poly
1278 br.ret.sptk b0 // Exit if |s| < 2^-33, and pi/4 <= |x| < 2^24
1284 // Here if 2^24 <= |x| < 2^63
1287 ldfe FR_d_1 = [GR_ad_p], 16 // Load d_1 for |x| >= 2^24 path
1288 fma.s1 FR_N_0 = FR_Input_X, FR_Inv_P_0, f0
1294 // N_0 = Arg * Inv_P_0
1296 // Load values 2**(-14) and -2**(-14)
1298 ldfps FR_Two_to_M14, FR_Neg_Two_to_M14 = [GR_ad_m14]
1302 ldfe FR_d_2 = [GR_ad_p], 16 // Load d_2 for |x| >= 2^24 path
1310 fcvt.fx.s1 FR_N_0_fix = FR_N_0
1316 // N_0_fix = integer part of N_0
1318 fcvt.xf FR_N_0 = FR_N_0_fix
1324 // Make N_0 the integer part
1326 fnma.s1 FR_ArgPrime = FR_N_0, FR_P_0, FR_Input_X
1331 fma.s1 FR_w = FR_N_0, FR_d_1, f0
1337 // Arg' = -N_0 * P_0 + Arg
1340 fma.s1 FR_N_float = FR_ArgPrime, FR_Inv_pi_by_2, f0
1348 fcvt.fx.s1 FR_N_fix = FR_N_float
1354 // N_fix is the integer part
1356 fcvt.xf FR_N_float = FR_N_fix
1360 getf.sig GR_N_Inc = FR_N_fix
1367 add GR_N_Inc = GR_N_Inc, GR_Sin_or_Cos ;;
1372 // N is the integer part of the reduced-reduced argument.
1373 // Put the integer in a GP register
1375 fnma.s1 FR_s = FR_N_float, FR_P_1, FR_ArgPrime
1380 fnma.s1 FR_w = FR_N_float, FR_P_2, FR_w
1386 // s = -N*P_1 + Arg'
1388 // N_fix_gr = N_fix_gr + N_inc
1390 fcmp.lt.unc.s1 p9, p8 = FR_s, FR_Two_to_M14
1395 (p9) fcmp.gt.s1 p9, p8 = FR_s, FR_Neg_Two_to_M14 // p9 if |s| < 2^-14
1402 // For |s| > 2**(-14) r = S + w (r complete)
1403 // Else U_hi = N_0 * d_1
1405 (p9) fma.s1 FR_V_hi = FR_N_float, FR_P_2, f0
1410 (p9) fma.s1 FR_U_hi = FR_N_0, FR_d_1, f0
1416 // Either S <= -2**(-14) or S >= 2**(-14)
1417 // or -2**(-14) < s < 2**(-14)
1419 (p8) fma.s1 FR_r = FR_s, f1, FR_w
1424 (p9) fma.s1 FR_w = FR_N_float, FR_P_3, f0
1430 // We need abs of both U_hi and V_hi - don't
1431 // worry about switched sign of V_hi.
1433 (p9) fms.s1 FR_A = FR_U_hi, f1, FR_V_hi
1439 // Big s: finish up c = (S - r) + w (c complete)
1440 // Case 4: A = U_hi + V_hi
1441 // Note: Worry about switched sign of V_hi, so subtract instead of add.
1443 (p9) fnma.s1 FR_V_lo = FR_N_float, FR_P_2, FR_V_hi
1449 (p9) fms.s1 FR_U_lo = FR_N_0, FR_d_1, FR_U_hi
1453 (p9) fmerge.s FR_V_hiabs = f0, FR_V_hi
1457 //(p9) fmerge.s f8= FR_V_lo,FR_V_lo
1458 //(p9) br.ret.sptk b0
1463 // For big s: c = S - r
1464 // For small s do more work: U_lo = N_0 * d_1 - U_hi
1466 (p9) fmerge.s FR_U_hiabs = f0, FR_U_hi
1472 // For big s: Is |r| < 2**(-3)
1473 // For big s: if p12 set, prepare to branch to Small_R.
1474 // For big s: If p13 set, prepare to branch to Normal_R.
1476 (p8) fms.s1 FR_c = FR_s, f1, FR_r
1482 // For small S: V_hi = N * P_2
1484 // Note the product does not include the (-) as in the writeup
1485 // so (-) missing for V_hi and w.
1487 (p8) fcmp.lt.unc.s1 p12, p13 = FR_r, FR_Two_to_M3
1492 (p12) fcmp.gt.s1 p12, p13 = FR_r, FR_Neg_Two_to_M3
1497 (p8) fma.s1 FR_c = FR_c, f1, FR_w
1502 (p9) fms.s1 FR_w = FR_N_0, FR_d_2, FR_w
1503 (p12) br.cond.spnt SINCOSL_SMALL_R // Branch if |r| < 2^-3
1504 // and 2^24 <= |x| < 2^63
1511 (p13) br.cond.sptk SINCOSL_NORMAL_R // Branch if |r| >= 2^-3
1512 // and 2^24 <= |x| < 2^63
1516 SINCOSL_LARGER_S_TINY:
1518 // Here if |s| < 2^-14, and 2^24 <= |x| < 2^63
1523 // Big s: Vector off when |r| < 2**(-3). Recall that p8 will be true.
1524 // The remaining stuff is for Case 4.
1525 // Small s: V_lo = N * P_2 + U_hi (U_hi is in place of V_hi in writeup)
1526 // Note: the (-) is still missing for V_lo.
1527 // Small s: w = w + N_0 * d_2
1528 // Note: the (-) is now incorporated in w.
1530 fcmp.ge.unc.s1 p7, p8 = FR_U_hiabs, FR_V_hiabs
1537 fma.s1 FR_t = FR_U_lo, f1, FR_V_lo
1547 (p7) fms.s1 FR_a = FR_U_hi, f1, FR_A
1552 (p8) fma.s1 FR_a = FR_V_hi, f1, FR_A
1559 // Is U_hiabs >= V_hiabs?
1562 fma.s1 FR_C_hi = FR_s, f1, FR_A
1566 ldfe FR_C_1 = [GR_ad_c], 16 ;;
1567 ldfe FR_C_2 = [GR_ad_c], 64
1571 // c = c + C_lo finished.
1575 ldfe FR_S_1 = [GR_ad_s], 16
1579 fma.s1 FR_t = FR_t, f1, FR_w
1583 // r and c have been computed.
1584 // Make sure ftz mode is set - should be automatic when using wre
1586 // Get [i_0,i_1] - two lsb of N_fix.
1590 ldfe FR_S_2 = [GR_ad_s], 64
1594 (p7) fms.s1 FR_a = FR_a, f1, FR_V_hi
1595 tbit.z p9,p10 = GR_N_Inc, 0 // p9 if i_1=0, N mod 4 = 0,1
1596 // p10 if i_1=1, N mod 4 = 2,3
1602 // For larger u than v: a = U_hi - A
1603 // Else a = V_hi - A (do an add to account for missing (-) on V_hi
1605 fms.s1 FR_C_lo = FR_s, f1, FR_C_hi
1610 (p8) fms.s1 FR_a = FR_U_hi, f1, FR_a
1611 tbit.z p11,p12 = GR_N_Inc, 1 // p11 if i_0=0, N mod 4 = 0,2
1612 // p12 if i_0=1, N mod 4 = 1,3
1619 // If u > v: a = (U_hi - A) + V_hi
1620 // Else a = (V_hi - A) + U_hi
1621 // In each case account for negative missing from V_hi.
1623 fma.s1 FR_C_lo = FR_C_lo, f1, FR_A
1629 // C_lo = (S - C_hi) + A
1631 fma.s1 FR_t = FR_t, f1, FR_a
1639 fma.s1 FR_C_lo = FR_C_lo, f1, FR_t
1647 fma.s1 FR_r = FR_C_hi, f1, FR_C_lo
1655 fma.s1 FR_rsq = FR_r, FR_r, f0
1663 fms.s1 FR_c = FR_C_hi, f1, FR_r
1669 // if i_1 ==0: poly = S_2 * FR_rsq + S_1
1670 // else poly = C_2 * FR_rsq + C_1
1672 (p9) fma.s1 FR_tmp_result = f0, f1, FR_r
1677 (p10) fma.s1 FR_tmp_result = f0, f1, f1
1683 // Compute r_cube = FR_rsq * r
1685 (p9) fma.s1 FR_poly = FR_rsq, FR_S_2, FR_S_1
1690 (p10) fma.s1 FR_poly = FR_rsq, FR_C_2, FR_C_1
1696 // Compute FR_rsq = r * r
1699 fma.s1 FR_r_cubed = FR_rsq, FR_r, f0
1708 fma.s1 FR_c = FR_c, f1, FR_C_lo
1714 // if i_1 ==0: poly = r_cube * poly + c
1715 // else poly = FR_rsq * poly
1717 (p12) fms.s1 FR_tmp_result = f0, f1, FR_tmp_result
1723 // if i_1 ==0: Result = r
1724 // else Result = 1.0
1726 (p9) fma.s1 FR_poly = FR_r_cubed, FR_poly, FR_c
1731 (p10) fma.s1 FR_poly = FR_rsq, FR_poly, f0
1737 // if i_0 !=0: Result = -Result
1739 (p11) fma.s0 FR_Result = FR_tmp_result, f1, FR_poly
1744 (p12) fms.s0 FR_Result = FR_tmp_result, f1, FR_poly
1746 // if i_0 == 0: Result = Result + poly
1747 // else Result = Result - poly
1749 br.ret.sptk b0 // Exit for |s| < 2^-14, and 2^24 <= |x| < 2^63
1756 // Here if |r| < 2^-3
1758 // Enter with r, c, and N_Inc computed
1760 // Compare both i_1 and i_0 with 0.
1761 // if i_1 == 0, set p9.
1762 // if i_0 == 0, set p11.
1767 fma.s1 FR_rsq = FR_r, FR_r, f0 // rsq = r * r
1768 tbit.z p9,p10 = GR_N_Inc, 0 // p9 if i_1=0, N mod 4 = 0,1
1769 // p10 if i_1=1, N mod 4 = 2,3
1774 (p9) ldfe FR_S_5 = [GR_ad_se], -16 // Load S_5 if i_1=0
1775 (p10) ldfe FR_C_5 = [GR_ad_ce], -16 // Load C_5 if i_1=1
1781 (p9) ldfe FR_S_4 = [GR_ad_se], -16 // Load S_4 if i_1=0
1782 (p10) ldfe FR_C_4 = [GR_ad_ce], -16 // Load C_4 if i_1=1
1788 // Entry point for 2^-3 < |x| < pi/4
1789 .pred.rel "mutex",p9,p10
1791 // Entry point for pi/4 < |x| < 2^24 and |r| < 2^-3
1792 .pred.rel "mutex",p9,p10
1794 (p9) ldfe FR_S_3 = [GR_ad_se], -16 // Load S_3 if i_1=0
1795 fma.s1 FR_Z = FR_rsq, FR_rsq, f0 // Z = rsq * rsq
1799 (p10) ldfe FR_C_3 = [GR_ad_ce], -16 // Load C_3 if i_1=1
1800 (p10) fnma.s1 FR_c = FR_c, FR_r, f0 // c = -c * r if i_1=0
1806 (p9) ldfe FR_S_2 = [GR_ad_se], -16 // Load S_2 if i_1=0
1807 (p10) ldfe FR_C_2 = [GR_ad_ce], -16 // Load C_2 if i_1=1
1808 (p10) fmerge.s FR_r = f1, f1
1813 (p9) ldfe FR_S_1 = [GR_ad_se], -16 // Load S_1 if i_1=0
1814 (p10) ldfe FR_C_1 = [GR_ad_ce], -16 // Load C_1 if i_1=1
1821 (p9) fma.s1 FR_Z = FR_Z, FR_r, f0 // Z = Z * r if i_1=0
1828 (p9) fma.s1 FR_poly_lo = FR_rsq, FR_S_5, FR_S_4 // poly_lo=rsq*S_5+S_4 if i_1=0
1833 (p10) fma.s1 FR_poly_lo = FR_rsq, FR_C_5, FR_C_4 // poly_lo=rsq*C_5+C_4 if i_1=1
1840 (p9) fma.s1 FR_poly_hi = FR_rsq, FR_S_2, FR_S_1 // poly_hi=rsq*S_2+S_1 if i_1=0
1845 (p10) fma.s1 FR_poly_hi = FR_rsq, FR_C_2, FR_C_1 // poly_hi=rsq*C_2+C_1 if i_1=1
1852 fma.s1 FR_Z = FR_Z, FR_rsq, f0 // Z = Z * rsq
1859 (p9) fma.s1 FR_poly_lo = FR_rsq, FR_poly_lo, FR_S_3 // p_lo=p_lo*rsq+S_3, i_1=0
1864 (p10) fma.s1 FR_poly_lo = FR_rsq, FR_poly_lo, FR_C_3 // p_lo=p_lo*rsq+C_3, i_1=1
1871 (p9) fma.s0 FR_inexact = FR_S_4, FR_S_4, f0 // Dummy op to set inexact
1872 tbit.z p11,p12 = GR_N_Inc, 1 // p11 if i_0=0, N mod 4 = 0,2
1873 // p12 if i_0=1, N mod 4 = 1,3
1877 (p10) fma.s0 FR_inexact = FR_C_1, FR_C_1, f0 // Dummy op to set inexact
1884 (p9) fma.s1 FR_poly_hi = FR_poly_hi, FR_rsq, f0 // p_hi=p_hi*rsq if i_1=0
1889 (p10) fma.s1 FR_poly_hi = FR_poly_hi, FR_rsq, f0 // p_hi=p_hi*rsq if i_1=1
1896 fma.s1 FR_poly = FR_Z, FR_poly_lo, FR_c // poly=Z*poly_lo+c
1903 (p9) fma.s1 FR_poly_hi = FR_r, FR_poly_hi, f0 // p_hi=r*p_hi if i_1=0
1910 (p12) fms.s1 FR_r = f0, f1, FR_r // r = -r if i_0=1
1917 fma.s1 FR_poly = FR_poly, f1, FR_poly_hi // poly=poly+poly_hi
1923 // if (i_0 == 0) Result = r + poly
1924 // if (i_0 != 0) Result = r - poly
1928 (p11) fma.s0 FR_Result = FR_r, f1, FR_poly
1933 (p12) fms.s0 FR_Result = FR_r, f1, FR_poly
1934 br.ret.sptk b0 // Exit for |r| < 2^-3
1941 // Here if 2^-3 <= |r| < pi/4
1942 // THIS IS THE MAIN PATH
1944 // Enter with r, c, and N_Inc having been computed
1947 ldfe FR_PP_6 = [GR_ad_pp], 16 // Load PP_6
1948 fma.s1 FR_rsq = FR_r, FR_r, f0 // rsq = r * r
1949 tbit.z p9,p10 = GR_N_Inc, 0 // p9 if i_1=0, N mod 4 = 0,1
1950 // p10 if i_1=1, N mod 4 = 2,3
1953 ldfe FR_QQ_6 = [GR_ad_qq], 16 // Load QQ_6
1960 (p9) ldfe FR_PP_5 = [GR_ad_pp], 16 // Load PP_5 if i_1=0
1961 (p10) ldfe FR_QQ_5 = [GR_ad_qq], 16 // Load QQ_5 if i_1=1
1967 // Entry for 2^-3 < |x| < pi/4
1968 .pred.rel "mutex",p9,p10
1970 (p9) ldfe FR_C_1 = [GR_ad_pp], 16 // Load C_1 if i_1=0
1971 (p10) ldfe FR_S_1 = [GR_ad_qq], 16 // Load S_1 if i_1=1
1972 frcpa.s1 FR_r_hi, p6 = f1, FR_r // r_hi = frcpa(r)
1978 (p9) fma.s1 FR_poly = FR_rsq, FR_PP_8, FR_PP_7 // poly = rsq*PP_8+PP_7 if i_1=0
1983 (p10) fma.s1 FR_poly = FR_rsq, FR_QQ_8, FR_QQ_7 // poly = rsq*QQ_8+QQ_7 if i_1=1
1990 fma.s1 FR_r_cubed = FR_r, FR_rsq, f0 // rcubed = r * rsq
1997 // Entry for pi/4 <= |x| < 2^24
1998 .pred.rel "mutex",p9,p10
2000 (p9) ldfe FR_PP_1 = [GR_ad_pp], 16 // Load PP_1_hi if i_1=0
2001 (p10) ldfe FR_QQ_1 = [GR_ad_qq], 16 // Load QQ_1 if i_1=1
2002 frcpa.s1 FR_r_hi, p6 = f1, FR_r_hi // r_hi = frpca(frcpa(r))
2007 (p9) ldfe FR_PP_4 = [GR_ad_pp], 16 // Load PP_4 if i_1=0
2008 (p9) fma.s1 FR_poly = FR_rsq, FR_poly, FR_PP_6 // poly = rsq*poly+PP_6 if i_1=0
2012 (p10) ldfe FR_QQ_4 = [GR_ad_qq], 16 // Load QQ_4 if i_1=1
2013 (p10) fma.s1 FR_poly = FR_rsq, FR_poly, FR_QQ_6 // poly = rsq*poly+QQ_6 if i_1=1
2020 (p9) fma.s1 FR_corr = FR_C_1, FR_rsq, f0 // corr = C_1 * rsq if i_1=0
2025 (p10) fma.s1 FR_corr = FR_S_1, FR_r_cubed, FR_r // corr = S_1 * r^3 + r if i_1=1
2031 (p9) ldfe FR_PP_3 = [GR_ad_pp], 16 // Load PP_3 if i_1=0
2032 fma.s1 FR_r_hi_sq = FR_r_hi, FR_r_hi, f0 // r_hi_sq = r_hi * r_hi
2036 (p10) ldfe FR_QQ_3 = [GR_ad_qq], 16 // Load QQ_3 if i_1=1
2037 fms.s1 FR_r_lo = FR_r, f1, FR_r_hi // r_lo = r - r_hi
2043 (p9) ldfe FR_PP_2 = [GR_ad_pp], 16 // Load PP_2 if i_1=0
2044 (p9) fma.s1 FR_poly = FR_rsq, FR_poly, FR_PP_5 // poly = rsq*poly+PP_5 if i_1=0
2048 (p10) ldfe FR_QQ_2 = [GR_ad_qq], 16 // Load QQ_2 if i_1=1
2049 (p10) fma.s1 FR_poly = FR_rsq, FR_poly, FR_QQ_5 // poly = rsq*poly+QQ_5 if i_1=1
2055 (p9) ldfe FR_PP_1_lo = [GR_ad_pp], 16 // Load PP_1_lo if i_1=0
2056 (p9) fma.s1 FR_corr = FR_corr, FR_c, FR_c // corr = corr * c + c if i_1=0
2061 (p10) fnma.s1 FR_corr = FR_corr, FR_c, f0 // corr = -corr * c if i_1=1
2068 (p9) fma.s1 FR_U_lo = FR_r, FR_r_hi, FR_r_hi_sq // U_lo = r*r_hi+r_hi_sq, i_1=0
2073 (p10) fma.s1 FR_U_lo = FR_r_hi, f1, FR_r // U_lo = r_hi + r if i_1=1
2080 (p9) fma.s1 FR_U_hi = FR_r_hi, FR_r_hi_sq, f0 // U_hi = r_hi*r_hi_sq if i_1=0
2085 (p10) fma.s1 FR_U_hi = FR_QQ_1, FR_r_hi_sq, f1 // U_hi = QQ_1*r_hi_sq+1, i_1=1
2092 (p9) fma.s1 FR_poly = FR_rsq, FR_poly, FR_PP_4 // poly = poly*rsq+PP_4 if i_1=0
2097 (p10) fma.s1 FR_poly = FR_rsq, FR_poly, FR_QQ_4 // poly = poly*rsq+QQ_4 if i_1=1
2104 (p9) fma.s1 FR_U_lo = FR_r, FR_r, FR_U_lo // U_lo = r * r + U_lo if i_1=0
2109 (p10) fma.s1 FR_U_lo = FR_r_lo, FR_U_lo, f0 // U_lo = r_lo * U_lo if i_1=1
2116 (p9) fma.s1 FR_U_hi = FR_PP_1, FR_U_hi, f0 // U_hi = PP_1 * U_hi if i_1=0
2123 (p9) fma.s1 FR_poly = FR_rsq, FR_poly, FR_PP_3 // poly = poly*rsq+PP_3 if i_1=0
2128 (p10) fma.s1 FR_poly = FR_rsq, FR_poly, FR_QQ_3 // poly = poly*rsq+QQ_3 if i_1=1
2135 (p9) fma.s1 FR_U_lo = FR_r_lo, FR_U_lo, f0 // U_lo = r_lo * U_lo if i_1=0
2140 (p10) fma.s1 FR_U_lo = FR_QQ_1,FR_U_lo, f0 // U_lo = QQ_1 * U_lo if i_1=1
2147 (p9) fma.s1 FR_U_hi = FR_r, f1, FR_U_hi // U_hi = r + U_hi if i_1=0
2154 (p9) fma.s1 FR_poly = FR_rsq, FR_poly, FR_PP_2 // poly = poly*rsq+PP_2 if i_1=0
2159 (p10) fma.s1 FR_poly = FR_rsq, FR_poly, FR_QQ_2 // poly = poly*rsq+QQ_2 if i_1=1
2166 (p9) fma.s1 FR_U_lo = FR_PP_1, FR_U_lo, f0 // U_lo = PP_1 * U_lo if i_1=0
2173 (p9) fma.s1 FR_poly = FR_rsq, FR_poly, FR_PP_1_lo // poly =poly*rsq+PP1lo i_1=0
2178 (p10) fma.s1 FR_poly = FR_rsq, FR_poly, f0 // poly = poly*rsq if i_1=1
2185 fma.s1 FR_V = FR_U_lo, f1, FR_corr // V = U_lo + corr
2186 tbit.z p11,p12 = GR_N_Inc, 1 // p11 if i_0=0, N mod 4 = 0,2
2187 // p12 if i_0=1, N mod 4 = 1,3
2193 (p9) fma.s0 FR_inexact = FR_PP_5, FR_PP_4, f0 // Dummy op to set inexact
2198 (p10) fma.s0 FR_inexact = FR_QQ_5, FR_QQ_5, f0 // Dummy op to set inexact
2205 (p9) fma.s1 FR_poly = FR_r_cubed, FR_poly, f0 // poly = poly*r^3 if i_1=0
2210 (p10) fma.s1 FR_poly = FR_rsq, FR_poly, f0 // poly = poly*rsq if i_1=1
2217 (p11) fma.s1 FR_tmp_result = f0, f1, f1// tmp_result=+1.0 if i_0=0
2222 (p12) fms.s1 FR_tmp_result = f0, f1, f1// tmp_result=-1.0 if i_0=1
2229 fma.s1 FR_V = FR_poly, f1, FR_V // V = poly + V
2234 // If i_0 = 0 Result = U_hi + V
2235 // If i_0 = 1 Result = -U_hi - V
2238 (p11) fma.s0 FR_Result = FR_tmp_result, FR_U_hi, FR_V
2243 (p12) fms.s0 FR_Result = FR_tmp_result, FR_U_hi, FR_V
2244 br.ret.sptk b0 // Exit for 2^-3 <= |r| < pi/4
2251 cmp.eq.unc p6, p7 = 0x1, GR_Sin_or_Cos
2259 (p7) fmerge.s FR_Result = FR_Input_X, FR_Input_X // If sin, result = input
2264 (p6) fma.s0 FR_Result = f1, f1, f0 // If cos, result=1.0
2265 br.ret.sptk b0 // Exit for x=0
2272 getf.exp GR_signexp_x = FR_norm_x // Get sign and exponent of x
2274 br.cond.sptk SINCOSL_COMMON // Return to common code
2282 // Path for Arg = +/- QNaN, SNaN, Inf
2283 // Invalid can be raised. SNaNs
2286 fmpy.s0 FR_Result = FR_Input_X, f0
2290 GLOBAL_IEEE754_END(cosl)
2292 // *******************************************************************
2293 // *******************************************************************
2294 // *******************************************************************
2296 // Special Code to handle very large argument case.
2297 // Call int __libm_pi_by_2_reduce(x,r,c) for |arguments| >= 2**63
2298 // The interface is custom:
2300 // (Arg or x) is in f8
2305 // Be sure to allocate at least 2 GP registers as output registers for
2306 // __libm_pi_by_2_reduce. This routine uses r59-60. These are used as
2307 // scratch registers within the __libm_pi_by_2_reduce routine (for speed).
2309 // We know also that __libm_pi_by_2_reduce preserves f10-15, f71-127. We
2310 // use this to eliminate save/restore of key fp registers in this calling
2313 // *******************************************************************
2314 // *******************************************************************
2315 // *******************************************************************
2317 LOCAL_LIBM_ENTRY(__libm_callout)
2318 SINCOSL_ARG_TOO_LARGE:
2322 .save ar.pfs,GR_SAVE_PFS
2323 mov GR_SAVE_PFS=ar.pfs // Save ar.pfs
2327 setf.exp FR_Two_to_M3 = GR_exp_2_to_m3 // Form 2^-3
2328 mov GR_SAVE_GP=gp // Save gp
2329 .save b0, GR_SAVE_B0
2330 mov GR_SAVE_B0=b0 // Save b0
2335 // Call argument reduction with x in f8
2336 // Returns with N in r8, r in f8, c in f9
2337 // Assumes f71-127 are preserved across the call
2340 setf.exp FR_Neg_Two_to_M3 = GR_exp_m2_to_m3 // Form -(2^-3)
2342 br.call.sptk b0=__libm_pi_by_2_reduce#
2346 add GR_N_Inc = GR_Sin_or_Cos,r8
2347 fcmp.lt.unc.s1 p6, p0 = FR_r, FR_Two_to_M3
2348 mov b0 = GR_SAVE_B0 // Restore return address
2352 mov gp = GR_SAVE_GP // Restore gp
2353 (p6) fcmp.gt.unc.s1 p6, p0 = FR_r, FR_Neg_Two_to_M3
2354 mov ar.pfs = GR_SAVE_PFS // Restore ar.pfs
2359 (p6) br.cond.spnt SINCOSL_SMALL_R // Branch if |r|< 2^-3 for |x| >= 2^63
2360 br.cond.sptk SINCOSL_NORMAL_R // Branch if |r|>=2^-3 for |x| >= 2^63
2363 LOCAL_LIBM_END(__libm_callout)
2364 .type __libm_pi_by_2_reduce#,@function
2365 .global __libm_pi_by_2_reduce#