1 .file "libm_sincos_large.s"
4 // Copyright (c) 2002 - 2003, Intel Corporation
5 // All rights reserved.
7 // Contributed 2002 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.
41 //==============================================================
42 // 02/15/02 Initial version
43 // 05/13/02 Changed interface to __libm_pi_by_2_reduce
44 // 02/10/03 Reordered header: .section, .global, .proc, .align;
45 // used data8 for long double table values
46 // 05/15/03 Reformatted data tables
49 // Overview of operation
50 //==============================================================
52 // These functions calculate the sin and cos for inputs
57 // They accept argument in f8
58 // and return result in f8 without final rounding
60 // __libm_sincos_large#
61 // It accepts argument in f8
62 // and returns cos in f8 and sin in f9 without final rounding
65 //*********************************************************************
67 // Accuracy: Within .7 ulps for 80-bit floating point values
68 // Very accurate for double precision values
70 //*********************************************************************
74 // Floating-Point Registers: f8 as Input Value, f8 and f9 as Return Values
77 // General Purpose Registers:
79 // r44-r45 (Used to pass arguments to pi_by_2 reduce routine)
81 // Predicate Registers: p6-p13
83 //*********************************************************************
85 // IEEE Special Conditions:
87 // Denormal fault raised on denormal inputs
88 // Overflow exceptions do not occur
89 // Underflow exceptions raised when appropriate for sin
90 // (No specialized error handling for this routine)
91 // Inexact raised when appropriate by algorithm
102 //*********************************************************************
104 // Mathematical Description
105 // ========================
107 // The computation of FSIN and FCOS is best handled in one piece of
108 // code. The main reason is that given any argument Arg, computation
109 // of trigonometric functions first calculate N and an approximation
112 // Arg = N pi/2 + alpha, |alpha| <= pi/4.
116 // cos( Arg ) = sin( (N+1) pi/2 + alpha ),
118 // therefore, the code for computing sine will produce cosine as long
119 // as 1 is added to N immediately after the argument reduction
127 // Arg = M pi/2 + alpha, |alpha| <= pi/4,
129 // let I = M mod 4, or I be the two lsb of M when M is represented
130 // as 2's complement. I = [i_0 i_1]. Then
132 // sin( Arg ) = (-1)^i_0 sin( alpha ) if i_1 = 0,
133 // = (-1)^i_0 cos( alpha ) if i_1 = 1.
137 // sin ((-pi/2 + alpha) = (-1) cos (alpha)
139 // sin (alpha) = sin (alpha)
141 // sin (pi/2 + alpha) = cos (alpha)
143 // sin (pi + alpha) = (-1) sin (alpha)
145 // sin ((3/2)pi + alpha) = (-1) cos (alpha)
147 // The value of alpha is obtained by argument reduction and
148 // represented by two working precision numbers r and c where
150 // alpha = r + c accurately.
152 // The reduction method is described in a previous write up.
153 // The argument reduction scheme identifies 4 cases. For Cases 2
154 // and 4, because |alpha| is small, sin(r+c) and cos(r+c) can be
155 // computed very easily by 2 or 3 terms of the Taylor series
156 // expansion as follows:
161 // sin(r + c) = r + c - r^3/6 accurately
162 // cos(r + c) = 1 - 2^(-67) accurately
167 // sin(r + c) = r + c - r^3/6 + r^5/120 accurately
168 // cos(r + c) = 1 - r^2/2 + r^4/24 accurately
170 // The only cases left are Cases 1 and 3 of the argument reduction
171 // procedure. These two cases will be merged since after the
172 // argument is reduced in either cases, we have the reduced argument
173 // represented as r + c and that the magnitude |r + c| is not small
174 // enough to allow the usage of a very short approximation.
176 // The required calculation is either
178 // sin(r + c) = sin(r) + correction, or
179 // cos(r + c) = cos(r) + correction.
183 // sin(r + c) = sin(r) + c sin'(r) + O(c^2)
184 // = sin(r) + c cos (r) + O(c^2)
185 // = sin(r) + c(1 - r^2/2) accurately.
188 // cos(r + c) = cos(r) - c sin(r) + O(c^2)
189 // = cos(r) - c(r - r^3/6) accurately.
191 // We therefore concentrate on accurately calculating sin(r) and
192 // cos(r) for a working-precision number r, |r| <= pi/4 to within
195 // The greatest challenge of this task is that the second terms of
198 // r - r^3/3! + r^r/5! - ...
202 // 1 - r^2/2! + r^4/4! - ...
204 // are not very small when |r| is close to pi/4 and the rounding
205 // errors will be a concern if simple polynomial accumulation is
206 // used. When |r| < 2^-3, however, the second terms will be small
207 // enough (6 bits or so of right shift) that a normal Horner
208 // recurrence suffices. Hence there are two cases that we consider
209 // in the accurate computation of sin(r) and cos(r), |r| <= pi/4.
211 // Case small_r: |r| < 2^(-3)
212 // --------------------------
214 // Since Arg = M pi/4 + r + c accurately, and M mod 4 is [i_0 i_1],
217 // sin(Arg) = (-1)^i_0 * sin(r + c) if i_1 = 0
218 // = (-1)^i_0 * cos(r + c) if i_1 = 1
220 // can be accurately approximated by
222 // sin(Arg) = (-1)^i_0 * [sin(r) + c] if i_1 = 0
223 // = (-1)^i_0 * [cos(r) - c*r] if i_1 = 1
225 // because |r| is small and thus the second terms in the correction
228 // Finally, sin(r) and cos(r) are approximated by polynomials of
231 // sin(r) = r + S_1 r^3 + S_2 r^5 + ... + S_5 r^11
232 // cos(r) = 1 + C_1 r^2 + C_2 r^4 + ... + C_5 r^10
234 // We can make use of predicates to selectively calculate
235 // sin(r) or cos(r) based on i_1.
237 // Case normal_r: 2^(-3) <= |r| <= pi/4
238 // ------------------------------------
240 // This case is more likely than the previous one if one considers
241 // r to be uniformly distributed in [-pi/4 pi/4]. Again,
243 // sin(Arg) = (-1)^i_0 * sin(r + c) if i_1 = 0
244 // = (-1)^i_0 * cos(r + c) if i_1 = 1.
246 // Because |r| is now larger, we need one extra term in the
247 // correction. sin(Arg) can be accurately approximated by
249 // sin(Arg) = (-1)^i_0 * [sin(r) + c(1-r^2/2)] if i_1 = 0
250 // = (-1)^i_0 * [cos(r) - c*r*(1 - r^2/6)] i_1 = 1.
252 // Finally, sin(r) and cos(r) are approximated by polynomials of
255 // sin(r) = r + PP_1_hi r^3 + PP_1_lo r^3 +
256 // PP_2 r^5 + ... + PP_8 r^17
258 // cos(r) = 1 + QQ_1 r^2 + QQ_2 r^4 + ... + QQ_8 r^16
260 // where PP_1_hi is only about 16 bits long and QQ_1 is -1/2.
261 // The crux in accurate computation is to calculate
263 // r + PP_1_hi r^3 or 1 + QQ_1 r^2
265 // accurately as two pieces: U_hi and U_lo. The way to achieve this
266 // is to obtain r_hi as a 10 sig. bit number that approximates r to
267 // roughly 8 bits or so of accuracy. (One convenient way is
269 // r_hi := frcpa( frcpa( r ) ).)
273 // r + PP_1_hi r^3 = r + PP_1_hi r_hi^3 +
274 // PP_1_hi (r^3 - r_hi^3)
275 // = [r + PP_1_hi r_hi^3] +
276 // [PP_1_hi (r - r_hi)
277 // (r^2 + r_hi r + r_hi^2) ]
280 // Since r_hi is only 10 bit long and PP_1_hi is only 16 bit long,
281 // PP_1_hi * r_hi^3 is only at most 46 bit long and thus computed
282 // exactly. Furthermore, r and PP_1_hi r_hi^3 are of opposite sign
283 // and that there is no more than 8 bit shift off between r and
284 // PP_1_hi * r_hi^3. Hence the sum, U_hi, is representable and thus
285 // calculated without any error. Finally, the fact that
287 // |U_lo| <= 2^(-8) |U_hi|
289 // says that U_hi + U_lo is approximating r + PP_1_hi r^3 to roughly
290 // 8 extra bits of accuracy.
294 // 1 + QQ_1 r^2 = [1 + QQ_1 r_hi^2] +
295 // [QQ_1 (r - r_hi)(r + r_hi)]
298 // Summarizing, we calculate r_hi = frcpa( frcpa( r ) ).
302 // U_hi := r + PP_1_hi * r_hi^3
303 // U_lo := PP_1_hi * (r - r_hi) * (r^2 + r*r_hi + r_hi^2)
304 // poly := PP_1_lo r^3 + PP_2 r^5 + ... + PP_8 r^17
305 // correction := c * ( 1 + C_1 r^2 )
309 // U_hi := 1 + QQ_1 * r_hi * r_hi
310 // U_lo := QQ_1 * (r - r_hi) * (r + r_hi)
311 // poly := QQ_2 * r^4 + QQ_3 * r^6 + ... + QQ_8 r^16
312 // correction := -c * r * (1 + S_1 * r^2)
318 // V := poly + ( U_lo + correction )
320 // / U_hi + V if i_0 = 0
322 // \ (-U_hi) - V if i_0 = 1
324 // It is important that in the last step, negation of U_hi is
325 // performed prior to the subtraction which is to be performed in
326 // the user-set rounding mode.
329 // Algorithmic Description
330 // =======================
332 // The argument reduction algorithm is tightly integrated into FSIN
333 // and FCOS which share the same code. The following is complete and
334 // self-contained. The argument reduction description given
335 // previously is repeated below.
338 // Step 0. Initialization.
340 // If FSIN is invoked, set N_inc := 0; else if FCOS is invoked,
343 // Step 1. Check for exceptional and special cases.
345 // * If Arg is +-0, +-inf, NaN, NaT, go to Step 10 for special
347 // * If |Arg| < 2^24, go to Step 2 for reduction of moderate
348 // arguments. This is the most likely case.
349 // * If |Arg| < 2^63, go to Step 8 for pre-reduction of large
351 // * If |Arg| >= 2^63, go to Step 10 for special handling.
353 // Step 2. Reduction of moderate arguments.
355 // If |Arg| < pi/4 ...quick branch
356 // N_fix := N_inc (integer)
359 // Branch to Step 4, Case_1_complete
360 // Else ...cf. argument reduction
361 // N := Arg * two_by_PI (fp)
362 // N_fix := fcvt.fx( N ) (int)
363 // N := fcvt.xf( N_fix )
364 // N_fix := N_fix + N_inc
365 // s := Arg - N * P_1 (first piece of pi/2)
366 // w := -N * P_2 (second piece of pi/2)
369 // go to Step 3, Case_1_reduce
371 // go to Step 7, Case_2_reduce
375 // Step 3. Case_1_reduce.
378 // c := (s - r) + w ...observe order
380 // Step 4. Case_1_complete
382 // ...At this point, the reduced argument alpha is
383 // ...accurately represented as r + c.
384 // If |r| < 2^(-3), go to Step 6, small_r.
388 // Let [i_0 i_1] by the 2 lsb of N_fix.
390 // r_hi := frcpa( frcpa( r ) )
394 // poly := r*FR_rsq*(PP_1_lo + FR_rsq*(PP_2 + ... FR_rsq*PP_8))
395 // U_hi := r + PP_1_hi*r_hi*r_hi*r_hi ...any order
396 // U_lo := PP_1_hi*r_lo*(r*r + r*r_hi + r_hi*r_hi)
397 // correction := c + c*C_1*FR_rsq ...any order
399 // poly := FR_rsq*FR_rsq*(QQ_2 + FR_rsq*(QQ_3 + ... + FR_rsq*QQ_8))
400 // U_hi := 1 + QQ_1 * r_hi * r_hi ...any order
401 // U_lo := QQ_1 * r_lo * (r + r_hi)
402 // correction := -c*(r + S_1*FR_rsq*r) ...any order
405 // V := poly + (U_lo + correction) ...observe order
407 // result := (i_0 == 0? 1.0 : -1.0)
409 // Last instruction in user-set rounding mode
411 // result := (i_0 == 0? result*U_hi + V :
418 // ...Use flush to zero mode without causing exception
419 // Let [i_0 i_1] be the two lsb of N_fix.
424 // z := FR_rsq*FR_rsq; z := FR_rsq*z *r
425 // poly_lo := S_3 + FR_rsq*(S_4 + FR_rsq*S_5)
426 // poly_hi := r*FR_rsq*(S_1 + FR_rsq*S_2)
430 // z := FR_rsq*FR_rsq; z := FR_rsq*z
431 // poly_lo := C_3 + FR_rsq*(C_4 + FR_rsq*C_5)
432 // poly_hi := FR_rsq*(C_1 + FR_rsq*C_2)
433 // correction := -c*r
437 // poly := poly_hi + (z * poly_lo + correction)
439 // If i_0 = 1, result := -result
441 // Last operation. Perform in user-set rounding mode
443 // result := (i_0 == 0? result + poly :
447 // Step 7. Case_2_reduce.
449 // ...Refer to the write up for argument reduction for
450 // ...rationale. The reduction algorithm below is taken from
451 // ...argument reduction description and integrated this.
454 // U_1 := N*P_2 + w ...FMA
455 // U_2 := (N*P_2 - U_1) + w ...2 FMA
456 // ...U_1 + U_2 is N*(P_2+P_3) accurately
459 // c := ( (s - r) - U_1 ) - U_2
461 // ...The mathematical sum r + c approximates the reduced
462 // ...argument accurately. Note that although compared to
463 // ...Case 1, this case requires much more work to reduce
464 // ...the argument, the subsequent calculation needed for
465 // ...any of the trigonometric function is very little because
466 // ...|alpha| < 1.01*2^(-33) and thus two terms of the
467 // ...Taylor series expansion suffices.
470 // poly := c + S_1 * r * r * r ...any order
477 // If i_0 = 1, result := -result
479 // Last operation. Perform in user-set rounding mode
481 // result := (i_0 == 0? result + poly :
487 // Step 8. Pre-reduction of large arguments.
489 // ...Again, the following reduction procedure was described
490 // ...in the separate write up for argument reduction, which
491 // ...is tightly integrated here.
493 // N_0 := Arg * Inv_P_0
494 // N_0_fix := fcvt.fx( N_0 )
495 // N_0 := fcvt.xf( N_0_fix)
497 // Arg' := Arg - N_0 * P_0
499 // N := Arg' * two_by_PI
500 // N_fix := fcvt.fx( N )
501 // N := fcvt.xf( N_fix )
502 // N_fix := N_fix + N_inc
504 // s := Arg' - N * P_1
513 // Step 9. Case_4_reduce.
515 // ...first obtain N_0*d_1 and -N*P_2 accurately
516 // U_hi := N_0 * d_1 V_hi := -N*P_2
517 // U_lo := N_0 * d_1 - U_hi V_lo := -N*P_2 - U_hi ...FMAs
519 // ...compute the contribution from N_0*d_1 and -N*P_3
522 // t := U_lo + V_lo + w ...any order
524 // ...at this point, the mathematical value
525 // ...s + U_hi + V_hi + t approximates the true reduced argument
526 // ...accurately. Just need to compute this accurately.
528 // ...Calculate U_hi + V_hi accurately:
530 // if |U_hi| >= |V_hi| then
531 // a := (U_hi - A) + V_hi
533 // a := (V_hi - A) + U_hi
535 // ...order in computing "a" must be observed. This branch is
536 // ...best implemented by predicates.
537 // ...A + a is U_hi + V_hi accurately. Moreover, "a" is
538 // ...much smaller than A: |a| <= (1/2)ulp(A).
540 // ...Just need to calculate s + A + a + t
541 // C_hi := s + A t := t + a
542 // C_lo := (s - C_hi) + A
545 // ...Final steps for reduction
547 // c := (C_hi - r) + C_lo
549 // ...At this point, we have r and c
550 // ...And all we need is a couple of terms of the corresponding
554 // poly := c + r*FR_rsq*(S_1 + FR_rsq*S_2)
557 // poly := FR_rsq*(C_1 + FR_rsq*C_2)
561 // If i_0 = 1, result := -result
563 // Last operation. Perform in user-set rounding mode
565 // result := (i_0 == 0? result + poly :
569 // Large Arguments: For arguments above 2**63, a Payne-Hanek
570 // style argument reduction is used and pi_by_2 reduce is called.
577 LOCAL_OBJECT_START(FSINCOS_CONSTANTS)
579 data4 0x4B800000 // two**24
580 data4 0xCB800000 // -two**24
581 data4 0x00000000 // pad
582 data4 0x00000000 // pad
583 data8 0xA2F9836E4E44152A, 0x00003FFE // Inv_pi_by_2
584 data8 0xC84D32B0CE81B9F1, 0x00004016 // P_0
585 data8 0xC90FDAA22168C235, 0x00003FFF // P_1
586 data8 0xECE675D1FC8F8CBB, 0x0000BFBD // P_2
587 data8 0xB7ED8FBBACC19C60, 0x0000BF7C // P_3
588 data4 0x5F000000 // two**63
589 data4 0xDF000000 // -two**63
590 data4 0x00000000 // pad
591 data4 0x00000000 // pad
592 data8 0xA397E5046EC6B45A, 0x00003FE7 // Inv_P_0
593 data8 0x8D848E89DBD171A1, 0x0000BFBF // d_1
594 data8 0xD5394C3618A66F8E, 0x0000BF7C // d_2
595 data8 0xC90FDAA22168C234, 0x00003FFE // pi_by_4
596 data8 0xC90FDAA22168C234, 0x0000BFFE // neg_pi_by_4
597 data4 0x3E000000 // two**-3
598 data4 0xBE000000 // -two**-3
599 data4 0x00000000 // pad
600 data4 0x00000000 // pad
601 data4 0x2F000000 // two**-33
602 data4 0xAF000000 // -two**-33
603 data4 0x9E000000 // -two**-67
604 data4 0x00000000 // pad
605 data8 0xCC8ABEBCA21C0BC9, 0x00003FCE // PP_8
606 data8 0xD7468A05720221DA, 0x0000BFD6 // PP_7
607 data8 0xB092382F640AD517, 0x00003FDE // PP_6
608 data8 0xD7322B47D1EB75A4, 0x0000BFE5 // PP_5
609 data8 0xFFFFFFFFFFFFFFFE, 0x0000BFFD // C_1
610 data8 0xAAAA000000000000, 0x0000BFFC // PP_1_hi
611 data8 0xB8EF1D2ABAF69EEA, 0x00003FEC // PP_4
612 data8 0xD00D00D00D03BB69, 0x0000BFF2 // PP_3
613 data8 0x8888888888888962, 0x00003FF8 // PP_2
614 data8 0xAAAAAAAAAAAB0000, 0x0000BFEC // PP_1_lo
615 data8 0xD56232EFC2B0FE52, 0x00003FD2 // QQ_8
616 data8 0xC9C99ABA2B48DCA6, 0x0000BFDA // QQ_7
617 data8 0x8F76C6509C716658, 0x00003FE2 // QQ_6
618 data8 0x93F27DBAFDA8D0FC, 0x0000BFE9 // QQ_5
619 data8 0xAAAAAAAAAAAAAAAA, 0x0000BFFC // S_1
620 data8 0x8000000000000000, 0x0000BFFE // QQ_1
621 data8 0xD00D00D00C6E5041, 0x00003FEF // QQ_4
622 data8 0xB60B60B60B607F60, 0x0000BFF5 // QQ_3
623 data8 0xAAAAAAAAAAAAAA9B, 0x00003FFA // QQ_2
624 data8 0xFFFFFFFFFFFFFFFE, 0x0000BFFD // C_1
625 data8 0xAAAAAAAAAAAA719F, 0x00003FFA // C_2
626 data8 0xB60B60B60356F994, 0x0000BFF5 // C_3
627 data8 0xD00CFFD5B2385EA9, 0x00003FEF // C_4
628 data8 0x93E4BD18292A14CD, 0x0000BFE9 // C_5
629 data8 0xAAAAAAAAAAAAAAAA, 0x0000BFFC // S_1
630 data8 0x88888888888868DB, 0x00003FF8 // S_2
631 data8 0xD00D00D0055EFD4B, 0x0000BFF2 // S_3
632 data8 0xB8EF1C5D839730B9, 0x00003FEC // S_4
633 data8 0xD71EA3A4E5B3F492, 0x0000BFE5 // S_5
634 data4 0x38800000 // two**-14
635 data4 0xB8800000 // -two**-14
636 LOCAL_OBJECT_END(FSINCOS_CONSTANTS)
638 // sin and cos registers
651 FR_Neg_Two_to_24 = f36
653 FR_Neg_Two_to_M14 = f37
654 FR_Neg_Two_to_M33 = f38
655 FR_Neg_Two_to_M67 = f39
717 FR_Neg_Two_to_M3 = f93
719 FR_Neg_Two_to_63 = f94
738 // sincos combined routine registers
741 GR_SINCOS_SAVE_PFS = r32
742 GR_SINCOS_SAVE_B0 = r33
743 GR_SINCOS_SAVE_GP = r34
747 FR_SINCOS_RES_SIN = f101
753 GLOBAL_LIBM_ENTRY(__libm_sincos_large)
756 alloc GR_SINCOS_SAVE_PFS = ar.pfs,0,3,0,0
757 fma.s1 FR_SINCOS_ARG = f8, f1, f0 // Save argument for sin and cos
758 mov GR_SINCOS_SAVE_B0 = b0
762 mov GR_SINCOS_SAVE_GP = gp
764 br.call.sptk b0 = __libm_sin_large // Call sin
769 fma.s1 FR_SINCOS_RES_SIN = f8, f1, f0 // Save sin result
775 fma.s1 f8 = FR_SINCOS_ARG, f1, f0 // Arg for cos
776 br.call.sptk b0 = __libm_cos_large // Call cos
780 mov gp = GR_SINCOS_SAVE_GP
781 fma.s1 f9 = FR_SINCOS_RES_SIN, f1, f0 // Out sin result
782 mov b0 = GR_SINCOS_SAVE_B0
787 mov ar.pfs = GR_SINCOS_SAVE_PFS
788 br.ret.sptk b0 // sincos_large exit
791 GLOBAL_LIBM_END(__libm_sincos_large)
796 GLOBAL_LIBM_ENTRY(__libm_sin_large)
799 alloc GR_Table_Base = ar.pfs,0,12,2,0
800 movl GR_Sin_or_Cos = 0x0 ;;
805 addl GR_Table_Base = @ltoff(FSINCOS_CONSTANTS#), gp
811 ld8 GR_Table_Base = [GR_Table_Base]
821 br.cond.sptk SINCOS_CONTINUE ;;
824 GLOBAL_LIBM_END(__libm_sin_large)
826 GLOBAL_LIBM_ENTRY(__libm_cos_large)
829 alloc GR_Table_Base= ar.pfs,0,12,2,0
830 movl GR_Sin_or_Cos = 0x1 ;;
835 addl GR_Table_Base = @ltoff(FSINCOS_CONSTANTS#), gp
841 ld8 GR_Table_Base = [GR_Table_Base]
848 // Load Table Address
853 add GR_Table_Base1 = 96, GR_Table_Base
854 ldfs FR_Two_to_24 = [GR_Table_Base], 4
862 // Load 2**24, load 2**63.
864 ldfs FR_Neg_Two_to_24 = [GR_Table_Base], 12
869 ldfs FR_Two_to_63 = [GR_Table_Base1], 4
871 // Check for unnormals - unsupported operands. We do not want
872 // to generate denormal exception
873 // Check for NatVals, QNaNs, SNaNs, +/-Infs
874 // Check for EM unsupporteds
877 fclass.m.unc p6, p8 = FR_Input_X, 0x1E3
883 fclass.nm.unc p8, p0 = FR_Input_X, 0x1FF
884 // GR_Sin_or_Cos denotes
889 ldfs FR_Neg_Two_to_63 = [GR_Table_Base1], 12
890 fclass.m.unc p10, p0 = FR_Input_X, 0x007
891 (p6) br.cond.spnt SINCOS_SPECIAL ;;
897 (p8) br.cond.spnt SINCOS_SPECIAL ;;
904 // Branch if +/- NaN, Inf.
905 // Load -2**24, load -2**63.
907 (p10) br.cond.spnt SINCOS_ZERO ;;
911 ldfe FR_Inv_pi_by_2 = [GR_Table_Base], 16
912 ldfe FR_Inv_P_0 = [GR_Table_Base1], 16
918 ldfe FR_d_1 = [GR_Table_Base1], 16
922 // Raise possible denormal operand flag with useful fcmp
924 // Load Inv_P_0 for pre-reduction
929 ldfe FR_P_0 = [GR_Table_Base], 16
930 ldfe FR_d_2 = [GR_Table_Base1], 16
941 ldfe FR_P_1 = [GR_Table_Base], 16 ;;
947 ldfe FR_P_2 = [GR_Table_Base], 16
953 ldfe FR_P_3 = [GR_Table_Base], 16
954 fcmp.le.unc.s1 p7, p8 = FR_Input_X, FR_Neg_Two_to_24
960 // Branch if +/- zero.
961 // Decide about the paths to take:
962 // If -2**24 < FR_Input_X < 2**24 - CASE 1 OR 2
963 // OTHERWISE - CASE 3 OR 4
965 fcmp.le.unc.s1 p10, p11 = FR_Input_X, FR_Neg_Two_to_63
971 (p8) fcmp.ge.s1 p7, p0 = FR_Input_X, FR_Two_to_24
976 ldfe FR_Pi_by_4 = [GR_Table_Base1], 16
977 (p11) fcmp.ge.s1 p10, p0 = FR_Input_X, FR_Two_to_63
982 ldfe FR_Neg_Pi_by_4 = [GR_Table_Base1], 16 ;;
983 ldfs FR_Two_to_M3 = [GR_Table_Base1], 4
988 ldfs FR_Neg_Two_to_M3 = [GR_Table_Base1], 12
998 (p10) br.cond.spnt SINCOS_ARG_TOO_LARGE ;;
1005 // Branch out if x >= 2**63. Use Payne-Hanek Reduction
1007 (p7) br.cond.spnt SINCOS_LARGER_ARG ;;
1013 // Branch if Arg <= -2**24 or Arg >= 2**24 and use pre-reduction.
1015 fma.s1 FR_N_float = FR_Input_X, FR_Inv_pi_by_2, f0
1021 fcmp.lt.unc.s1 p6, p7 = FR_Input_X, FR_Pi_by_4
1028 // Select the case when |Arg| < pi/4
1029 // Else Select the case when |Arg| >= pi/4
1031 fcvt.fx.s1 FR_N_fix = FR_N_float
1039 // Check if Arg < pi/4
1041 (p6) fcmp.gt.s1 p6, p7 = FR_Input_X, FR_Neg_Pi_by_4
1045 // Case 2: Convert integer N_fix back to normalized floating-point value.
1046 // Case 1: p8 is only affected when p6 is set
1050 (p7) ldfs FR_Two_to_M33 = [GR_Table_Base1], 4
1052 // Grab the integer part of N and call it N_fix
1054 (p6) fmerge.se FR_r = FR_Input_X, FR_Input_X
1055 // If |x| < pi/4, r = x and c = 0
1056 // lf |x| < pi/4, is x < 2**(-3).
1059 (p6) mov GR_N_Inc = GR_Sin_or_Cos ;;
1064 (p7) ldfs FR_Neg_Two_to_M33 = [GR_Table_Base1], 4
1065 (p6) fmerge.se FR_c = f0, f0
1070 (p6) fcmp.lt.unc.s1 p8, p9 = FR_Input_X, FR_Two_to_M3
1077 // lf |x| < pi/4, is -2**(-3)< x < 2**(-3) - set p8.
1079 // Create the right N for |x| < pi/4 and otherwise
1080 // Case 2: Place integer part of N in GP register
1082 (p7) fcvt.xf FR_N_float = FR_N_fix
1088 (p7) getf.sig GR_N_Inc = FR_N_fix
1089 (p8) fcmp.gt.s1 p8, p0 = FR_Input_X, FR_Neg_Two_to_M3 ;;
1096 // Load 2**(-33), -2**(-33)
1098 (p8) br.cond.spnt SINCOS_SMALL_R ;;
1104 (p6) br.cond.sptk SINCOS_NORMAL_R ;;
1107 // if |x| < pi/4, branch based on |x| < 2**(-3) or otherwise.
1110 // In this branch, |x| >= pi/4.
1114 ldfs FR_Neg_Two_to_M67 = [GR_Table_Base1], 8
1118 fnma.s1 FR_s = FR_N_float, FR_P_1, FR_Input_X
1121 // s = -N * P_1 + Arg
1123 add GR_N_Inc = GR_N_Inc, GR_Sin_or_Cos
1128 fma.s1 FR_w = FR_N_float, FR_P_2, f0
1135 // Adjust N_fix by N_inc to determine whether sine or
1136 // cosine is being calculated
1138 fcmp.lt.unc.s1 p7, p6 = FR_s, FR_Two_to_M33
1144 (p7) fcmp.gt.s1 p7, p6 = FR_s, FR_Neg_Two_to_M33
1150 // Remember x >= pi/4.
1151 // Is s <= -2**(-33) or s >= 2**(-33) (p6)
1152 // or -2**(-33) < s < 2**(-33) (p7)
1153 (p6) fms.s1 FR_r = FR_s, f1, FR_w
1159 (p7) fma.s1 FR_w = FR_N_float, FR_P_3, f0
1165 (p7) fma.s1 FR_U_1 = FR_N_float, FR_P_2, FR_w
1171 (p6) fms.s1 FR_c = FR_s, f1, FR_r
1178 // For big s: r = s - w: No futher reduction is necessary
1179 // For small s: w = N * P_3 (change sign) More reduction
1181 (p6) fcmp.lt.unc.s1 p8, p9 = FR_r, FR_Two_to_M3
1187 (p8) fcmp.gt.s1 p8, p9 = FR_r, FR_Neg_Two_to_M3
1193 (p7) fms.s1 FR_r = FR_s, f1, FR_U_1
1200 // For big s: Is |r| < 2**(-3)?
1201 // For big s: c = S - r
1202 // For small s: U_1 = N * P_2 + w
1204 // If p8 is set, prepare to branch to Small_R.
1205 // If p9 is set, prepare to branch to Normal_R.
1206 // For big s, r is complete here.
1208 (p6) fms.s1 FR_c = FR_c, f1, FR_w
1210 // For big s: c = c + w (w has not been negated.)
1211 // For small s: r = S - U_1
1213 (p8) br.cond.spnt SINCOS_SMALL_R ;;
1219 (p9) br.cond.sptk SINCOS_NORMAL_R ;;
1223 (p7) add GR_Table_Base1 = 224, GR_Table_Base1
1225 // Branch to SINCOS_SMALL_R or SINCOS_NORMAL_R
1227 (p7) fms.s1 FR_U_2 = FR_N_float, FR_P_2, FR_U_1
1233 (p7) extr.u GR_i_1 = GR_N_Inc, 0, 1
1239 // Get [i_0,i_1] - two lsb of N_fix_gr.
1240 // Do dummy fmpy so inexact is always set.
1242 (p7) cmp.eq.unc p9, p10 = 0x0, GR_i_1
1243 (p7) extr.u GR_i_0 = GR_N_Inc, 1, 1 ;;
1246 // For small s: U_2 = N * P_2 - U_1
1247 // S_1 stored constant - grab the one stored with the
1252 (p7) ldfe FR_S_1 = [GR_Table_Base1], 16
1254 // Check if i_1 and i_0 != 0
1256 (p10) fma.s1 FR_poly = f0, f1, FR_Neg_Two_to_M67
1257 (p7) cmp.eq.unc p11, p12 = 0x0, GR_i_0 ;;
1262 (p7) fms.s1 FR_s = FR_s, f1, FR_r
1273 (p7) fma.s1 FR_rsq = FR_r, FR_r, f0
1279 (p7) fma.s1 FR_U_2 = FR_U_2, f1, FR_w
1285 //(p7) fmerge.se FR_Input_X = FR_r, FR_r
1286 (p7) fmerge.se FR_prelim = FR_r, FR_r
1292 //(p10) fma.s1 FR_Input_X = f0, f1, f1
1293 (p10) fma.s1 FR_prelim = f0, f1, f1
1301 // Save r as the result.
1303 (p7) fms.s1 FR_c = FR_s, f1, FR_U_1
1310 // if ( i_1 ==0) poly = c + S_1*r*r*r
1313 //(p12) fnma.s1 FR_Input_X = FR_Input_X, f1, f0
1314 (p12) fnma.s1 FR_prelim = FR_prelim, f1, f0
1320 (p7) fma.s1 FR_r = FR_S_1, FR_r, f0
1326 (p7) fma.d.s1 FR_S_1 = FR_S_1, FR_S_1, f0
1333 // If i_1 != 0, poly = 2**(-67)
1335 (p7) fms.s1 FR_c = FR_c, f1, FR_U_2
1344 (p9) fma.s1 FR_poly = FR_r, FR_rsq, FR_c
1351 // i_0 != 0, so Result = -Result
1353 (p11) fma.s1 FR_Input_X = FR_prelim, f1, FR_poly
1359 (p12) fms.s1 FR_Input_X = FR_prelim, f1, FR_poly
1361 // if (i_0 == 0), Result = Result + poly
1362 // else Result = Result - poly
1370 fma.s1 FR_N_0 = FR_Input_X, FR_Inv_P_0, f0
1375 // This path for argument > 2*24
1376 // Adjust table_ptr1 to beginning of table.
1381 addl GR_Table_Base = @ltoff(FSINCOS_CONSTANTS#), gp
1387 ld8 GR_Table_Base = [GR_Table_Base]
1396 // N_0 = Arg * Inv_P_0
1400 add GR_Table_Base = 688, GR_Table_Base ;;
1401 ldfs FR_Two_to_M14 = [GR_Table_Base], 4
1406 ldfs FR_Neg_Two_to_M14 = [GR_Table_Base], 0
1414 // Load values 2**(-14) and -2**(-14)
1416 fcvt.fx.s1 FR_N_0_fix = FR_N_0
1423 // N_0_fix = integer part of N_0
1425 fcvt.xf FR_N_0 = FR_N_0_fix
1432 // Make N_0 the integer part
1434 fnma.s1 FR_ArgPrime = FR_N_0, FR_P_0, FR_Input_X
1440 fma.s1 FR_w = FR_N_0, FR_d_1, f0
1447 // Arg' = -N_0 * P_0 + Arg
1450 fma.s1 FR_N_float = FR_ArgPrime, FR_Inv_pi_by_2, f0
1459 fcvt.fx.s1 FR_N_fix = FR_N_float
1466 // N_fix is the integer part
1468 fcvt.xf FR_N_float = FR_N_fix
1473 getf.sig GR_N_Inc = FR_N_fix
1481 add GR_N_Inc = GR_N_Inc, GR_Sin_or_Cos ;;
1487 // N is the integer part of the reduced-reduced argument.
1488 // Put the integer in a GP register
1490 fnma.s1 FR_s = FR_N_float, FR_P_1, FR_ArgPrime
1496 fnma.s1 FR_w = FR_N_float, FR_P_2, FR_w
1503 // s = -N*P_1 + Arg'
1505 // N_fix_gr = N_fix_gr + N_inc
1507 fcmp.lt.unc.s1 p9, p8 = FR_s, FR_Two_to_M14
1513 (p9) fcmp.gt.s1 p9, p8 = FR_s, FR_Neg_Two_to_M14
1520 // For |s| > 2**(-14) r = S + w (r complete)
1521 // Else U_hi = N_0 * d_1
1523 (p9) fma.s1 FR_V_hi = FR_N_float, FR_P_2, f0
1529 (p9) fma.s1 FR_U_hi = FR_N_0, FR_d_1, f0
1536 // Either S <= -2**(-14) or S >= 2**(-14)
1537 // or -2**(-14) < s < 2**(-14)
1539 (p8) fma.s1 FR_r = FR_s, f1, FR_w
1545 (p9) fma.s1 FR_w = FR_N_float, FR_P_3, f0
1552 // We need abs of both U_hi and V_hi - don't
1553 // worry about switched sign of V_hi.
1555 (p9) fms.s1 FR_A = FR_U_hi, f1, FR_V_hi
1562 // Big s: finish up c = (S - r) + w (c complete)
1563 // Case 4: A = U_hi + V_hi
1564 // Note: Worry about switched sign of V_hi, so subtract instead of add.
1566 (p9) fnma.s1 FR_V_lo = FR_N_float, FR_P_2, FR_V_hi
1572 (p9) fms.s1 FR_U_lo = FR_N_0, FR_d_1, FR_U_hi
1578 (p9) fmerge.s FR_V_hiabs = f0, FR_V_hi
1584 // For big s: c = S - r
1585 // For small s do more work: U_lo = N_0 * d_1 - U_hi
1587 (p9) fmerge.s FR_U_hiabs = f0, FR_U_hi
1594 // For big s: Is |r| < 2**(-3)
1595 // For big s: if p12 set, prepare to branch to Small_R.
1596 // For big s: If p13 set, prepare to branch to Normal_R.
1598 (p8) fms.s1 FR_c = FR_s, f1, FR_r
1605 // For small S: V_hi = N * P_2
1607 // Note the product does not include the (-) as in the writeup
1608 // so (-) missing for V_hi and w.
1610 (p8) fcmp.lt.unc.s1 p12, p13 = FR_r, FR_Two_to_M3
1616 (p12) fcmp.gt.s1 p12, p13 = FR_r, FR_Neg_Two_to_M3
1622 (p8) fma.s1 FR_c = FR_c, f1, FR_w
1628 (p9) fms.s1 FR_w = FR_N_0, FR_d_2, FR_w
1629 (p12) br.cond.spnt SINCOS_SMALL_R ;;
1635 (p13) br.cond.sptk SINCOS_NORMAL_R ;;
1641 // Big s: Vector off when |r| < 2**(-3). Recall that p8 will be true.
1642 // The remaining stuff is for Case 4.
1643 // Small s: V_lo = N * P_2 + U_hi (U_hi is in place of V_hi in writeup)
1644 // Note: the (-) is still missing for V_lo.
1645 // Small s: w = w + N_0 * d_2
1646 // Note: the (-) is now incorporated in w.
1648 (p9) fcmp.ge.unc.s1 p10, p11 = FR_U_hiabs, FR_V_hiabs
1649 extr.u GR_i_1 = GR_N_Inc, 0, 1 ;;
1657 (p9) fma.s1 FR_t = FR_U_lo, f1, FR_V_lo
1658 extr.u GR_i_0 = GR_N_Inc, 1, 1 ;;
1667 (p10) fms.s1 FR_a = FR_U_hi, f1, FR_A
1673 (p11) fma.s1 FR_a = FR_V_hi, f1, FR_A
1680 addl GR_Table_Base = @ltoff(FSINCOS_CONSTANTS#), gp
1686 ld8 GR_Table_Base = [GR_Table_Base]
1694 add GR_Table_Base = 528, GR_Table_Base
1696 // Is U_hiabs >= V_hiabs?
1698 (p9) fma.s1 FR_C_hi = FR_s, f1, FR_A
1703 ldfe FR_C_1 = [GR_Table_Base], 16 ;;
1704 ldfe FR_C_2 = [GR_Table_Base], 64
1711 // c = c + C_lo finished.
1714 ldfe FR_S_1 = [GR_Table_Base], 16
1718 fma.s1 FR_t = FR_t, f1, FR_w ;;
1721 // r and c have been computed.
1722 // Make sure ftz mode is set - should be automatic when using wre
1724 // Get [i_0,i_1] - two lsb of N_fix.
1729 ldfe FR_S_2 = [GR_Table_Base], 64
1733 (p10) fms.s1 FR_a = FR_a, f1, FR_V_hi
1734 cmp.eq.unc p9, p10 = 0x0, GR_i_0
1740 // For larger u than v: a = U_hi - A
1741 // Else a = V_hi - A (do an add to account for missing (-) on V_hi
1743 fms.s1 FR_C_lo = FR_s, f1, FR_C_hi
1749 (p11) fms.s1 FR_a = FR_U_hi, f1, FR_a
1750 cmp.eq.unc p11, p12 = 0x0, GR_i_1
1756 // If u > v: a = (U_hi - A) + V_hi
1757 // Else a = (V_hi - A) + U_hi
1758 // In each case account for negative missing from V_hi.
1760 fma.s1 FR_C_lo = FR_C_lo, f1, FR_A
1767 // C_lo = (S - C_hi) + A
1769 fma.s1 FR_t = FR_t, f1, FR_a
1778 fma.s1 FR_C_lo = FR_C_lo, f1, FR_t
1786 // Adjust Table_Base to beginning of table
1788 fma.s1 FR_r = FR_C_hi, f1, FR_C_lo
1797 fma.s1 FR_rsq = FR_r, FR_r, f0
1804 // Table_Base points to C_1
1807 fms.s1 FR_c = FR_C_hi, f1, FR_r
1814 // if i_1 ==0: poly = S_2 * FR_rsq + S_1
1815 // else poly = C_2 * FR_rsq + C_1
1817 //(p11) fma.s1 FR_Input_X = f0, f1, FR_r
1818 (p11) fma.s1 FR_prelim = f0, f1, FR_r
1824 //(p12) fma.s1 FR_Input_X = f0, f1, f1
1825 (p12) fma.s1 FR_prelim = f0, f1, f1
1832 // Compute r_cube = FR_rsq * r
1834 (p11) fma.s1 FR_poly = FR_rsq, FR_S_2, FR_S_1
1840 (p12) fma.s1 FR_poly = FR_rsq, FR_C_2, FR_C_1
1847 // Compute FR_rsq = r * r
1850 fma.s1 FR_r_cubed = FR_rsq, FR_r, f0
1860 fma.s1 FR_c = FR_c, f1, FR_C_lo
1867 // if i_1 ==0: poly = r_cube * poly + c
1868 // else poly = FR_rsq * poly
1870 //(p10) fms.s1 FR_Input_X = f0, f1, FR_Input_X
1871 (p10) fms.s1 FR_prelim = f0, f1, FR_prelim
1878 // if i_1 ==0: Result = r
1879 // else Result = 1.0
1881 (p11) fma.s1 FR_poly = FR_r_cubed, FR_poly, FR_c
1887 (p12) fma.s1 FR_poly = FR_rsq, FR_poly, f0
1894 // if i_0 !=0: Result = -Result
1896 (p9) fma.s1 FR_Input_X = FR_prelim, f1, FR_poly
1902 (p10) fms.s1 FR_Input_X = FR_prelim, f1, FR_poly
1904 // if i_0 == 0: Result = Result + poly
1905 // else Result = Result - poly
1913 extr.u GR_i_1 = GR_N_Inc, 0, 1 ;;
1916 // Compare both i_1 and i_0 with 0.
1917 // if i_1 == 0, set p9.
1918 // if i_0 == 0, set p11.
1920 cmp.eq.unc p9, p10 = 0x0, GR_i_1 ;;
1925 fma.s1 FR_rsq = FR_r, FR_r, f0
1926 extr.u GR_i_0 = GR_N_Inc, 1, 1 ;;
1934 (p10) fnma.s1 FR_c = FR_c, FR_r, f0
1935 cmp.eq.unc p11, p12 = 0x0, GR_i_0
1939 // ******************************************************************
1940 // ******************************************************************
1941 // ******************************************************************
1942 // r and c have been computed.
1943 // We know whether this is the sine or cosine routine.
1944 // Make sure ftz mode is set - should be automatic when using wre
1947 // Set table_ptr1 to beginning of constant table.
1948 // Get [i_0,i_1] - two lsb of N_fix_gr.
1953 addl GR_Table_Base = @ltoff(FSINCOS_CONSTANTS#), gp
1959 ld8 GR_Table_Base = [GR_Table_Base]
1967 // Set table_ptr1 to point to S_5.
1968 // Set table_ptr1 to point to C_5.
1969 // Compute FR_rsq = r * r
1973 (p9) add GR_Table_Base = 672, GR_Table_Base
1974 (p10) fmerge.s FR_r = f1, f1
1975 (p10) add GR_Table_Base = 592, GR_Table_Base ;;
1978 // Set table_ptr1 to point to S_5.
1979 // Set table_ptr1 to point to C_5.
1983 (p9) ldfe FR_S_5 = [GR_Table_Base], -16 ;;
1985 // if (i_1 == 0) load S_5
1986 // if (i_1 != 0) load C_5
1988 (p9) ldfe FR_S_4 = [GR_Table_Base], -16
1993 (p10) ldfe FR_C_5 = [GR_Table_Base], -16
1995 // Z = FR_rsq * FR_rsq
1997 (p9) ldfe FR_S_3 = [GR_Table_Base], -16
1999 // Compute FR_rsq = r * r
2000 // if (i_1 == 0) load S_4
2001 // if (i_1 != 0) load C_4
2003 fma.s1 FR_Z = FR_rsq, FR_rsq, f0 ;;
2006 // if (i_1 == 0) load S_3
2007 // if (i_1 != 0) load C_3
2011 (p9) ldfe FR_S_2 = [GR_Table_Base], -16 ;;
2013 // if (i_1 == 0) load S_2
2014 // if (i_1 != 0) load C_2
2016 (p9) ldfe FR_S_1 = [GR_Table_Base], -16
2021 (p10) ldfe FR_C_4 = [GR_Table_Base], -16 ;;
2022 (p10) ldfe FR_C_3 = [GR_Table_Base], -16
2027 (p10) ldfe FR_C_2 = [GR_Table_Base], -16 ;;
2028 (p10) ldfe FR_C_1 = [GR_Table_Base], -16
2036 // poly_lo = FR_rsq * C_5 + C_4
2037 // poly_hi = FR_rsq * C_2 + C_1
2039 (p9) fma.s1 FR_Z = FR_Z, FR_r, f0
2046 // if (i_1 == 0) load S_1
2047 // if (i_1 != 0) load C_1
2049 (p9) fma.s1 FR_poly_lo = FR_rsq, FR_S_5, FR_S_4
2057 // dummy fmpy's to flag inexact.
2059 (p9) fma.d.s1 FR_S_4 = FR_S_4, FR_S_4, f0
2066 // poly_lo = FR_rsq * poly_lo + C_3
2067 // poly_hi = FR_rsq * poly_hi
2069 fma.s1 FR_Z = FR_Z, FR_rsq, f0
2075 (p9) fma.s1 FR_poly_hi = FR_rsq, FR_S_2, FR_S_1
2083 // poly_lo = FR_rsq * S_5 + S_4
2084 // poly_hi = FR_rsq * S_2 + S_1
2086 (p10) fma.s1 FR_poly_lo = FR_rsq, FR_C_5, FR_C_4
2094 // Z = Z * r for only one of the small r cases - not there
2095 // in original implementation notes.
2097 (p9) fma.s1 FR_poly_lo = FR_rsq, FR_poly_lo, FR_S_3
2103 (p10) fma.s1 FR_poly_hi = FR_rsq, FR_C_2, FR_C_1
2109 (p10) fma.d.s1 FR_C_1 = FR_C_1, FR_C_1, f0
2115 (p9) fma.s1 FR_poly_hi = FR_poly_hi, FR_rsq, f0
2122 // poly_lo = FR_rsq * poly_lo + S_3
2123 // poly_hi = FR_rsq * poly_hi
2125 (p10) fma.s1 FR_poly_lo = FR_rsq, FR_poly_lo, FR_C_3
2131 (p10) fma.s1 FR_poly_hi = FR_poly_hi, FR_rsq, f0
2138 // if (i_1 == 0): dummy fmpy's to flag inexact
2141 (p9) fma.s1 FR_poly_hi = FR_r, FR_poly_hi, f0
2148 // poly_hi = r * poly_hi
2150 fma.s1 FR_poly = FR_Z, FR_poly_lo, FR_c
2156 (p12) fms.s1 FR_r = f0, f1, FR_r
2163 // poly_hi = Z * poly_lo + c
2164 // if i_0 == 1: r = -r
2166 fma.s1 FR_poly = FR_poly, f1, FR_poly_hi
2172 (p12) fms.s1 FR_Input_X = FR_r, f1, FR_poly
2179 // poly = poly + poly_hi
2181 (p11) fma.s1 FR_Input_X = FR_r, f1, FR_poly
2183 // if (i_0 == 0) Result = r + poly
2184 // if (i_0 != 0) Result = r - poly
2192 extr.u GR_i_1 = GR_N_Inc, 0, 1 ;;
2194 // Set table_ptr1 and table_ptr2 to base address of
2196 cmp.eq.unc p9, p10 = 0x0, GR_i_1 ;;
2201 fma.s1 FR_rsq = FR_r, FR_r, f0
2202 extr.u GR_i_0 = GR_N_Inc, 1, 1 ;;
2207 frcpa.s1 FR_r_hi, p6 = f1, FR_r
2208 cmp.eq.unc p11, p12 = 0x0, GR_i_0
2212 // ******************************************************************
2213 // ******************************************************************
2214 // ******************************************************************
2216 // r and c have been computed.
2217 // We known whether this is the sine or cosine routine.
2218 // Make sure ftz mode is set - should be automatic when using wre
2219 // Get [i_0,i_1] - two lsb of N_fix_gr alone.
2224 addl GR_Table_Base = @ltoff(FSINCOS_CONSTANTS#), gp
2230 ld8 GR_Table_Base = [GR_Table_Base]
2238 (p10) add GR_Table_Base = 384, GR_Table_Base
2239 //(p12) fms.s1 FR_Input_X = f0, f1, f1
2240 (p12) fms.s1 FR_prelim = f0, f1, f1
2241 (p9) add GR_Table_Base = 224, GR_Table_Base ;;
2246 (p10) ldfe FR_QQ_8 = [GR_Table_Base], 16
2248 // if (i_1==0) poly = poly * FR_rsq + PP_1_lo
2249 // else poly = FR_rsq * poly
2251 //(p11) fma.s1 FR_Input_X = f0, f1, f1 ;;
2252 (p11) fma.s1 FR_prelim = f0, f1, f1 ;;
2256 (p10) ldfe FR_QQ_7 = [GR_Table_Base], 16
2258 // Adjust table pointers based on i_0
2259 // Compute rsq = r * r
2261 (p9) ldfe FR_PP_8 = [GR_Table_Base], 16
2262 fma.s1 FR_r_cubed = FR_r, FR_rsq, f0 ;;
2266 (p9) ldfe FR_PP_7 = [GR_Table_Base], 16
2267 (p10) ldfe FR_QQ_6 = [GR_Table_Base], 16
2269 // Load PP_8 and QQ_8; PP_7 and QQ_7
2271 frcpa.s1 FR_r_hi, p6 = f1, FR_r_hi ;;
2274 // if (i_1==0) poly = PP_7 + FR_rsq * PP_8.
2275 // else poly = QQ_7 + FR_rsq * QQ_8.
2279 (p9) ldfe FR_PP_6 = [GR_Table_Base], 16
2280 (p10) ldfe FR_QQ_5 = [GR_Table_Base], 16
2285 (p9) ldfe FR_PP_5 = [GR_Table_Base], 16
2286 (p10) ldfe FR_S_1 = [GR_Table_Base], 16
2291 (p10) ldfe FR_QQ_1 = [GR_Table_Base], 16
2292 (p9) ldfe FR_C_1 = [GR_Table_Base], 16
2297 (p10) ldfe FR_QQ_4 = [GR_Table_Base], 16 ;;
2298 (p9) ldfe FR_PP_1 = [GR_Table_Base], 16
2303 (p10) ldfe FR_QQ_3 = [GR_Table_Base], 16
2305 // if (i_1=0) corr = corr + c*c
2306 // else corr = corr * c
2308 (p9) ldfe FR_PP_4 = [GR_Table_Base], 16
2309 (p10) fma.s1 FR_poly = FR_rsq, FR_QQ_8, FR_QQ_7 ;;
2312 // if (i_1=0) poly = rsq * poly + PP_5
2313 // else poly = rsq * poly + QQ_5
2314 // Load PP_4 or QQ_4
2318 (p9) ldfe FR_PP_3 = [GR_Table_Base], 16
2319 (p10) ldfe FR_QQ_2 = [GR_Table_Base], 16
2321 // r_hi = frcpa(frcpa(r)).
2322 // r_cube = r * FR_rsq.
2324 (p9) fma.s1 FR_poly = FR_rsq, FR_PP_8, FR_PP_7 ;;
2327 // Do dummy multiplies so inexact is always set.
2331 (p9) ldfe FR_PP_2 = [GR_Table_Base], 16
2335 (p9) fma.s1 FR_U_lo = FR_r_hi, FR_r_hi, f0
2341 (p9) ldfe FR_PP_1_lo = [GR_Table_Base], 16
2342 (p10) fma.s1 FR_corr = FR_S_1, FR_r_cubed, FR_r
2347 (p10) fma.s1 FR_poly = FR_rsq, FR_poly, FR_QQ_6
2354 // if (i_1=0) U_lo = r_hi * r_hi
2355 // else U_lo = r_hi + r
2357 (p9) fma.s1 FR_corr = FR_C_1, FR_rsq, f0
2364 // if (i_1=0) corr = C_1 * rsq
2365 // else corr = S_1 * r_cubed + r
2367 (p9) fma.s1 FR_poly = FR_rsq, FR_poly, FR_PP_6
2373 (p10) fma.s1 FR_U_lo = FR_r_hi, f1, FR_r
2380 // if (i_1=0) U_hi = r_hi + U_hi
2381 // else U_hi = QQ_1 * U_hi + 1
2383 (p9) fma.s1 FR_U_lo = FR_r, FR_r_hi, FR_U_lo
2390 // U_hi = r_hi * r_hi
2392 fms.s1 FR_r_lo = FR_r, f1, FR_r_hi
2399 // Load PP_1, PP_6, PP_5, and C_1
2400 // Load QQ_1, QQ_6, QQ_5, and S_1
2402 fma.s1 FR_U_hi = FR_r_hi, FR_r_hi, f0
2408 (p10) fma.s1 FR_poly = FR_rsq, FR_poly, FR_QQ_5
2414 (p10) fnma.s1 FR_corr = FR_corr, FR_c, f0
2421 // if (i_1=0) U_lo = r * r_hi + U_lo
2422 // else U_lo = r_lo * U_lo
2424 (p9) fma.s1 FR_corr = FR_corr, FR_c, FR_c
2430 (p9) fma.s1 FR_poly = FR_rsq, FR_poly, FR_PP_5
2437 // if (i_1 =0) U_hi = r + U_hi
2438 // if (i_1 =0) U_lo = r_lo * U_lo
2441 (p9) fma.d.s1 FR_PP_5 = FR_PP_5, FR_PP_4, f0
2447 (p9) fma.s1 FR_U_lo = FR_r, FR_r, FR_U_lo
2453 (p10) fma.s1 FR_U_lo = FR_r_lo, FR_U_lo, f0
2460 // if (i_1=0) poly = poly * rsq + PP_6
2461 // else poly = poly * rsq + QQ_6
2463 (p9) fma.s1 FR_U_hi = FR_r_hi, FR_U_hi, f0
2469 (p10) fma.s1 FR_poly = FR_rsq, FR_poly, FR_QQ_4
2475 (p10) fma.s1 FR_U_hi = FR_QQ_1, FR_U_hi, f1
2481 (p10) fma.d.s1 FR_QQ_5 = FR_QQ_5, FR_QQ_5, f0
2488 // if (i_1!=0) U_hi = PP_1 * U_hi
2489 // if (i_1!=0) U_lo = r * r + U_lo
2490 // Load PP_3 or QQ_3
2492 (p9) fma.s1 FR_poly = FR_rsq, FR_poly, FR_PP_4
2498 (p9) fma.s1 FR_U_lo = FR_r_lo, FR_U_lo, f0
2504 (p10) fma.s1 FR_U_lo = FR_QQ_1,FR_U_lo, f0
2510 (p9) fma.s1 FR_U_hi = FR_PP_1, FR_U_hi, f0
2516 (p10) fma.s1 FR_poly = FR_rsq, FR_poly, FR_QQ_3
2525 (p9) fma.s1 FR_poly = FR_rsq, FR_poly, FR_PP_3
2532 // if (i_1==0) poly = FR_rsq * poly + PP_3
2533 // else poly = FR_rsq * poly + QQ_3
2536 (p9) fma.s1 FR_U_lo = FR_PP_1, FR_U_lo, f0
2543 // if (i_1 =0) poly = poly * rsq + pp_r4
2544 // else poly = poly * rsq + qq_r4
2546 (p9) fma.s1 FR_U_hi = FR_r, f1, FR_U_hi
2552 (p10) fma.s1 FR_poly = FR_rsq, FR_poly, FR_QQ_2
2559 // if (i_1==0) U_lo = PP_1_hi * U_lo
2560 // else U_lo = QQ_1 * U_lo
2562 (p9) fma.s1 FR_poly = FR_rsq, FR_poly, FR_PP_2
2569 // if (i_0==0) Result = 1
2572 fma.s1 FR_V = FR_U_lo, f1, FR_corr
2578 (p10) fma.s1 FR_poly = FR_rsq, FR_poly, f0
2585 // if (i_1==0) poly = FR_rsq * poly + PP_2
2586 // else poly = FR_rsq * poly + QQ_2
2588 (p9) fma.s1 FR_poly = FR_rsq, FR_poly, FR_PP_1_lo
2594 (p10) fma.s1 FR_poly = FR_rsq, FR_poly, f0
2603 (p9) fma.s1 FR_poly = FR_r_cubed, FR_poly, f0
2610 // if (i_1==0) poly = r_cube * poly
2611 // else poly = FR_rsq * poly
2613 fma.s1 FR_V = FR_poly, f1, FR_V
2619 //(p12) fms.s1 FR_Input_X = FR_Input_X, FR_U_hi, FR_V
2620 (p12) fms.s1 FR_Input_X = FR_prelim, FR_U_hi, FR_V
2629 //(p11) fma.s1 FR_Input_X = FR_Input_X, FR_U_hi, FR_V
2630 (p11) fma.s1 FR_Input_X = FR_prelim, FR_U_hi, FR_V
2632 // if (i_0==0) Result = Result * U_hi + V
2633 // else Result = Result * U_hi - V
2639 // If cosine, FR_Input_X = 1
2640 // If sine, FR_Input_X = +/-Zero (Input FR_Input_X)
2641 // Results are exact, no exceptions
2646 cmp.eq.unc p6, p7 = 0x1, GR_Sin_or_Cos
2653 (p7) fmerge.s FR_Input_X = FR_Input_X, FR_Input_X
2659 (p6) fmerge.s FR_Input_X = f1, f1
2666 // Path for Arg = +/- QNaN, SNaN, Inf
2667 // Invalid can be raised. SNaNs
2673 fmpy.s1 FR_Input_X = FR_Input_X, f0
2676 GLOBAL_LIBM_END(__libm_cos_large)
2679 // *******************************************************************
2680 // *******************************************************************
2681 // *******************************************************************
2683 // Special Code to handle very large argument case.
2684 // Call int __libm_pi_by_2_reduce(x,r,c) for |arguments| >= 2**63
2685 // The interface is custom:
2687 // (Arg or x) is in f8
2692 // Be sure to allocate at least 2 GP registers as output registers for
2693 // __libm_pi_by_2_reduce. This routine uses r49-50. These are used as
2694 // scratch registers within the __libm_pi_by_2_reduce routine (for speed).
2696 // We know also that __libm_pi_by_2_reduce preserves f10-15, f71-127. We
2697 // use this to eliminate save/restore of key fp registers in this calling
2700 // *******************************************************************
2701 // *******************************************************************
2702 // *******************************************************************
2704 LOCAL_LIBM_ENTRY(__libm_callout_2)
2705 SINCOS_ARG_TOO_LARGE:
2708 // Readjust Table ptr
2710 adds GR_Table_Base1 = -16, GR_Table_Base1
2712 .save ar.pfs,GR_SAVE_PFS
2713 mov GR_SAVE_PFS=ar.pfs // Save ar.pfs
2717 ldfs FR_Two_to_M3 = [GR_Table_Base1],4
2718 mov GR_SAVE_GP=gp // Save gp
2719 .save b0, GR_SAVE_B0
2720 mov GR_SAVE_B0=b0 // Save b0
2725 // Call argument reduction with x in f8
2726 // Returns with N in r8, r in f8, c in f9
2727 // Assumes f71-127 are preserved across the call
2730 ldfs FR_Neg_Two_to_M3 = [GR_Table_Base1],0
2732 br.call.sptk b0=__libm_pi_by_2_reduce#
2736 add GR_N_Inc = GR_Sin_or_Cos,r8
2737 fcmp.lt.unc.s1 p6, p0 = FR_r, FR_Two_to_M3
2738 mov b0 = GR_SAVE_B0 // Restore return address
2742 mov gp = GR_SAVE_GP // Restore gp
2743 (p6) fcmp.gt.unc.s1 p6, p0 = FR_r, FR_Neg_Two_to_M3
2744 mov ar.pfs = GR_SAVE_PFS // Restore ar.pfs
2749 (p6) br.cond.spnt SINCOS_SMALL_R // Branch if |r| < 1/4
2750 br.cond.sptk SINCOS_NORMAL_R ;; // Branch if 1/4 <= |r| < pi/4
2753 LOCAL_LIBM_END(__libm_callout_2)
2755 .type __libm_pi_by_2_reduce#,@function
2756 .global __libm_pi_by_2_reduce#