4 // Copyright (c) 2000 - 2003, 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 //*********************************************************************
44 // 02/02/00 (hand-optimized)
45 // 04/04/00 Unwind support added
46 // 12/28/00 Fixed false invalid flags
47 // 02/06/02 Improved speed
48 // 05/07/02 Changed interface to __libm_pi_by_2_reduce
49 // 05/30/02 Added cotl
50 // 02/10/03 Reordered header: .section, .global, .proc, .align;
51 // used data8 for long double table values
52 // 05/15/03 Reformatted data tables
54 //*********************************************************************
56 // Functions: tanl(x) = tangent(x), for double-extended precision x values
57 // cotl(x) = cotangent(x), for double-extended precision x values
59 //*********************************************************************
63 // Floating-Point Registers: f8 (Input and Return Value)
67 // General Purpose Registers:
70 // Predicate Registers: p6-p15
72 //*********************************************************************
74 // IEEE Special Conditions for tanl:
76 // Denormal fault raised on denormal inputs
77 // Overflow exceptions do not occur
78 // Underflow exceptions raised when appropriate for tan
79 // (No specialized error handling for this routine)
80 // Inexact raised when appropriate by algorithm
87 //*********************************************************************
89 // IEEE Special Conditions for cotl:
91 // Denormal fault raised on denormal inputs
92 // Overflow exceptions occur at zero and near zero
93 // Underflow exceptions do not occur
94 // Inexact raised when appropriate by algorithm
99 // cotl(+/-0) = +/-Inf and error handling is called
101 //*********************************************************************
103 // Below are mathematical and algorithmic descriptions for tanl.
104 // For cotl we use next identity cot(x) = -tan(x + Pi/2).
105 // So, to compute cot(x) we just need to increment N (N = N + 1)
106 // and invert sign of the computed result.
108 //*********************************************************************
110 // Mathematical Description
112 // We consider the computation of FPTANL of Arg. Now, given
114 // Arg = N pi/2 + alpha, |alpha| <= pi/4,
116 // basic mathematical relationship shows that
118 // tan( Arg ) = tan( alpha ) if N is even;
119 // = -cot( alpha ) otherwise.
121 // The value of alpha is obtained by argument reduction and
122 // represented by two working precision numbers r and c where
124 // alpha = r + c accurately.
126 // The reduction method is described in a previous write up.
127 // The argument reduction scheme identifies 4 cases. For Cases 2
128 // and 4, because |alpha| is small, tan(r+c) and -cot(r+c) can be
129 // computed very easily by 2 or 3 terms of the Taylor series
130 // expansion as follows:
135 // tan(r + c) = r + c + r^3/3 ...accurately
136 // -cot(r + c) = -1/(r+c) + r/3 ...accurately
141 // tan(r + c) = r + c + r^3/3 + 2r^5/15 ...accurately
142 // -cot(r + c) = -1/(r+c) + r/3 + r^3/45 ...accurately
145 // The only cases left are Cases 1 and 3 of the argument reduction
146 // procedure. These two cases will be merged since after the
147 // argument is reduced in either cases, we have the reduced argument
148 // represented as r + c and that the magnitude |r + c| is not small
149 // enough to allow the usage of a very short approximation.
151 // The greatest challenge of this task is that the second terms of
152 // the Taylor series for tan(r) and -cot(r)
154 // r + r^3/3 + 2 r^5/15 + ...
158 // -1/r + r/3 + r^3/45 + ...
160 // are not very small when |r| is close to pi/4 and the rounding
161 // errors will be a concern if simple polynomial accumulation is
162 // used. When |r| < 2^(-2), however, the second terms will be small
163 // enough (5 bits or so of right shift) that a normal Horner
164 // recurrence suffices. Hence there are two cases that we consider
165 // in the accurate computation of tan(r) and cot(r), |r| <= pi/4.
167 // Case small_r: |r| < 2^(-2)
168 // --------------------------
170 // Since Arg = N pi/4 + r + c accurately, we have
172 // tan(Arg) = tan(r+c) for N even,
173 // = -cot(r+c) otherwise.
175 // Here for this case, both tan(r) and -cot(r) can be approximated
176 // by simple polynomials:
178 // tan(r) = r + P1_1 r^3 + P1_2 r^5 + ... + P1_9 r^19
179 // -cot(r) = -1/r + Q1_1 r + Q1_2 r^3 + ... + Q1_7 r^13
181 // accurately. Since |r| is relatively small, tan(r+c) and
182 // -cot(r+c) can be accurately approximated by replacing r with
183 // r+c only in the first two terms of the corresponding polynomials.
185 // Note that P1_1 (and Q1_1 for that matter) approximates 1/3 to
186 // almost 64 sig. bits, thus
188 // P1_1 (r+c)^3 = P1_1 r^3 + c * r^2 accurately.
192 // tan(r+c) = r + P1_1 r^3 + P1_2 r^5 + ... + P1_9 r^19
195 // -cot(r+c) = -1/(r+c) + Q1_1 r + Q1_2 r^3 + ... + Q1_7 r^13
199 // Case normal_r: 2^(-2) <= |r| <= pi/4
200 // ------------------------------------
202 // This case is more likely than the previous one if one considers
203 // r to be uniformly distributed in [-pi/4 pi/4].
205 // The required calculation is either
207 // tan(r + c) = tan(r) + correction, or
208 // -cot(r + c) = -cot(r) + correction.
212 // tan(r + c) = tan(r) + c tan'(r) + O(c^2)
213 // = tan(r) + c sec^2(r) + O(c^2)
214 // = tan(r) + c SEC_sq ...accurately
215 // as long as SEC_sq approximates sec^2(r)
216 // to, say, 5 bits or so.
220 // -cot(r + c) = -cot(r) - c cot'(r) + O(c^2)
221 // = -cot(r) + c csc^2(r) + O(c^2)
222 // = -cot(r) + c CSC_sq ...accurately
223 // as long as CSC_sq approximates csc^2(r)
224 // to, say, 5 bits or so.
226 // We therefore concentrate on accurately calculating tan(r) and
227 // cot(r) for a working-precision number r, |r| <= pi/4 to within
230 // We will employ a table-driven approach. Let
232 // r = sgn_r * 2^k * 1.b_1 b_2 ... b_5 ... b_63
233 // = sgn_r * ( B + x )
237 // B = 2^k * 1.b_1 b_2 ... b_5 1
242 // tan( B + x ) = ------------------------
246 // | tan(B) + tan(x) |
248 // = tan(B) + | ------------------------ - tan(B) |
249 // | 1 - tan(B)*tan(x) |
253 // = tan(B) + ------------------------
256 // (1/[sin(B)*cos(B)]) * tan(x)
257 // = tan(B) + --------------------------------
261 // Clearly, the values of tan(B), cot(B) and 1/(sin(B)*cos(B)) are
262 // calculated beforehand and stored in a table. Since
264 // |x| <= 2^k * 2^(-6) <= 2^(-7) (because k = -1, -2)
266 // a very short polynomial will be sufficient to approximate tan(x)
267 // accurately. The details involved in computing the last expression
268 // will be given in the next section on algorithm description.
271 // Now, we turn to the case where cot( B + x ) is needed.
275 // cot( B + x ) = ------------------------
279 // | 1 - tan(B)*tan(x) |
281 // = cot(B) + | ----------------------- - cot(B) |
282 // | tan(B) + tan(x) |
285 // [tan(B) + cot(B)] * tan(x)
286 // = cot(B) - ----------------------------
289 // (1/[sin(B)*cos(B)]) * tan(x)
290 // = cot(B) - --------------------------------
294 // Note that the values of tan(B), cot(B) and 1/(sin(B)*cos(B)) that
295 // are needed are the same set of values needed in the previous
298 // Finally, we can put all the ingredients together as follows:
300 // Arg = N * pi/2 + r + c ...accurately
302 // tan(Arg) = tan(r) + correction if N is even;
303 // = -cot(r) + correction otherwise.
305 // For Cases 2 and 4,
308 // tan(Arg) = tan(r + c) = r + c + r^3/3 N even
309 // = -cot(r + c) = -1/(r+c) + r/3 N odd
311 // tan(Arg) = tan(r + c) = r + c + r^3/3 + 2r^5/15 N even
312 // = -cot(r + c) = -1/(r+c) + r/3 + r^3/45 N odd
315 // For Cases 1 and 3,
317 // Case small_r: |r| < 2^(-2)
319 // tan(Arg) = r + P1_1 r^3 + P1_2 r^5 + ... + P1_9 r^19
320 // + c*(1 + r^2) N even
322 // = -1/(r+c) + Q1_1 r + Q1_2 r^3 + ... + Q1_7 r^13
325 // Case normal_r: 2^(-2) <= |r| <= pi/4
327 // tan(Arg) = tan(r) + c * sec^2(r) N even
328 // = -cot(r) + c * csc^2(r) otherwise
332 // tan(Arg) = tan(r) + c*sec^2(r)
333 // = tan( sgn_r * (B+x) ) + c * sec^2(|r|)
334 // = sgn_r * ( tan(B+x) + sgn_r*c*sec^2(|r|) )
335 // = sgn_r * ( tan(B+x) + sgn_r*c*sec^2(B) )
337 // since B approximates |r| to 2^(-6) in relative accuracy.
339 // / (1/[sin(B)*cos(B)]) * tan(x)
340 // tan(Arg) = sgn_r * | tan(B) + --------------------------------
348 // CORR = sgn_r*c*tan(B)*SC_inv(B); SC_inv(B) = 1/(sin(B)*cos(B)).
352 // tan(Arg) = -cot(r) + c*csc^2(r)
353 // = -cot( sgn_r * (B+x) ) + c * csc^2(|r|)
354 // = sgn_r * ( -cot(B+x) + sgn_r*c*csc^2(|r|) )
355 // = sgn_r * ( -cot(B+x) + sgn_r*c*csc^2(B) )
357 // since B approximates |r| to 2^(-6) in relative accuracy.
359 // / (1/[sin(B)*cos(B)]) * tan(x)
360 // tan(Arg) = sgn_r * | -cot(B) + --------------------------------
368 // CORR = sgn_r*c*cot(B)*SC_inv(B); SC_inv(B) = 1/(sin(B)*cos(B)).
371 // The actual algorithm prescribes how all the mathematical formulas
375 // 2. Algorithmic Description
376 // ==========================
378 // 2.1 Computation for Cases 2 and 4.
379 // ----------------------------------
381 // For Case 2, we use two-term polynomials.
386 // Poly := c + r * rsq * P1_1
387 // Result := r + Poly ...in user-defined rounding
390 // S_hi := -frcpa(r) ...8 bits
391 // S_hi := S_hi + S_hi*(1 + S_hi*r) ...16 bits
392 // S_hi := S_hi + S_hi*(1 + S_hi*r) ...32 bits
393 // S_hi := S_hi + S_hi*(1 + S_hi*r) ...64 bits
394 // S_lo := S_hi*( (1 + S_hi*r) + S_hi*c )
395 // ...S_hi + S_lo is -1/(r+c) to extra precision
396 // S_lo := S_lo + Q1_1*r
398 // Result := S_hi + S_lo ...in user-defined rounding
400 // For Case 4, we use three-term polynomials
405 // Poly := c + r * rsq * (P1_1 + rsq * P1_2)
406 // Result := r + Poly ...in user-defined rounding
409 // S_hi := -frcpa(r) ...8 bits
410 // S_hi := S_hi + S_hi*(1 + S_hi*r) ...16 bits
411 // S_hi := S_hi + S_hi*(1 + S_hi*r) ...32 bits
412 // S_hi := S_hi + S_hi*(1 + S_hi*r) ...64 bits
413 // S_lo := S_hi*( (1 + S_hi*r) + S_hi*c )
414 // ...S_hi + S_lo is -1/(r+c) to extra precision
416 // P := Q1_1 + rsq*Q1_2
417 // S_lo := S_lo + r*P
419 // Result := S_hi + S_lo ...in user-defined rounding
422 // Note that the coefficients P1_1, P1_2, Q1_1, and Q1_2 are
423 // the same as those used in the small_r case of Cases 1 and 3
427 // 2.2 Computation for Cases 1 and 3.
428 // ----------------------------------
429 // This is further divided into the case of small_r,
430 // where |r| < 2^(-2), and the case of normal_r, where |r| lies between
433 // Algorithm for the case of small_r
434 // ---------------------------------
438 // Poly1 := rsq*(P1_1 + rsq*(P1_2 + rsq*P1_3))
439 // r_to_the_8 := rsq * rsq
440 // r_to_the_8 := r_to_the_8 * r_to_the_8
441 // Poly2 := P1_4 + rsq*(P1_5 + rsq*(P1_6 + ... rsq*P1_9))
442 // CORR := c * ( 1 + rsq )
443 // Poly := Poly1 + r_to_the_8*Poly2
444 // Poly := r*Poly + CORR
445 // Result := r + Poly ...in user-defined rounding
446 // ...note that Poly1 and r_to_the_8 can be computed in parallel
447 // ...with Poly2 (Poly1 is intentionally set to be much
448 // ...shorter than Poly2 so that r_to_the_8 and CORR can be hidden)
451 // S_hi := -frcpa(r) ...8 bits
452 // S_hi := S_hi + S_hi*(1 + S_hi*r) ...16 bits
453 // S_hi := S_hi + S_hi*(1 + S_hi*r) ...32 bits
454 // S_hi := S_hi + S_hi*(1 + S_hi*r) ...64 bits
455 // S_lo := S_hi*( (1 + S_hi*r) + S_hi*c )
456 // ...S_hi + S_lo is -1/(r+c) to extra precision
457 // S_lo := S_lo + Q1_1*c
459 // ...S_hi and S_lo are computed in parallel with
462 // P := Q1_1 + rsq*(Q1_2 + rsq*(Q1_3 + ... + rsq*Q1_7))
464 // Poly := r*P + S_lo
465 // Result := S_hi + Poly ...in user-defined rounding
468 // Algorithm for the case of normal_r
469 // ----------------------------------
471 // Here, we first consider the computation of tan( r + c ). As
472 // presented in the previous section,
474 // tan( r + c ) = tan(r) + c * sec^2(r)
475 // = sgn_r * [ tan(B+x) + CORR ]
476 // CORR = sgn_r * c * tan(B) * 1/[sin(B)*cos(B)]
478 // because sec^2(r) = sec^(|r|), and B approximate |r| to 6.5 bits.
481 // / (1/[sin(B)*cos(B)]) * tan(x)
482 // sgn_r * | tan(B) + -------------------------------- +
489 // The values of tan(B), cot(B) and 1/(sin(B)*cos(B)) are
490 // calculated beforehand and stored in a table. Specifically,
491 // the table values are
493 // tan(B) as T_hi + T_lo;
494 // cot(B) as C_hi + C_lo;
495 // 1/[sin(B)*cos(B)] as SC_inv
497 // T_hi, C_hi are in double-precision memory format;
498 // T_lo, C_lo are in single-precision memory format;
499 // SC_inv is in extended-precision memory format.
501 // The value of tan(x) will be approximated by a short polynomial of
504 // tan(x) as x + x * P, where
505 // P = x^2 * (P2_1 + x^2 * (P2_2 + x^2 * P2_3))
507 // Because |x| <= 2^(-7), cot(B) - x approximates cot(B) - tan(x)
508 // to a relative accuracy better than 2^(-20). Thus, a good
509 // initial guess of 1/( cot(B) - tan(x) ) to initiate the iterative
512 // 1/(cot(B) - tan(x)) is approximately
514 // tan(B)/(1 - x*tan(B)) is approximately
515 // T_hi / ( 1 - T_hi * x ) is approximately
517 // T_hi * [ 1 + (Thi * x) + (T_hi * x)^2 ]
519 // The calculation of tan(r+c) therefore proceed as follows:
524 // V_hi := T_hi*(1 + Tx*(1 + Tx))
525 // P := xsq * (P1_1 + xsq*(P1_2 + xsq*P1_3))
526 // ...V_hi serves as an initial guess of 1/(cot(B) - tan(x))
527 // ...good to about 20 bits of accuracy
531 // ...D is a double precision denominator: cot(B) - tan(x)
533 // V_hi := V_hi + V_hi*(1 - V_hi*D)
534 // ....V_hi approximates 1/(cot(B)-tan(x)) to 40 bits
536 // V_lo := V_hi * ( [ (1 - V_hi*C_hi) + V_hi*tanx ]
537 // - V_hi*C_lo ) ...observe all order
538 // ...V_hi + V_lo approximates 1/(cot(B) - tan(x))
539 // ...to extra accuracy
541 // ... SC_inv(B) * (x + x*P)
542 // ... tan(B) + ------------------------- + CORR
543 // ... cot(B) - (x + x*P)
545 // ... = tan(B) + SC_inv(B)*(x + x*P)*(V_hi + V_lo) + CORR
549 // CORR := sgn_r * c * SC_inv * T_hi
551 // ...put the ingredients together to compute
552 // ... SC_inv(B) * (x + x*P)
553 // ... tan(B) + ------------------------- + CORR
554 // ... cot(B) - (x + x*P)
556 // ... = tan(B) + SC_inv(B)*(x + x*P)*(V_hi + V_lo) + CORR
558 // ... = T_hi + T_lo + CORR +
559 // ... Sx * V_hi + Sx * V_lo + Sx * P *(V_hi + V_lo)
561 // CORR := CORR + T_lo
562 // tail := V_lo + P*(V_hi + V_lo)
563 // tail := Sx * tail + CORR
564 // tail := Sx * V_hi + tail
565 // T_hi := sgn_r * T_hi
567 // ...T_hi + sgn_r*tail now approximate
568 // ...sgn_r*(tan(B+x) + CORR) accurately
570 // Result := T_hi + sgn_r*tail ...in user-defined
571 // ...rounding control
572 // ...It is crucial that independent paths be fully
573 // ...exploited for performance's sake.
576 // Next, we consider the computation of -cot( r + c ). As
577 // presented in the previous section,
579 // -cot( r + c ) = -cot(r) + c * csc^2(r)
580 // = sgn_r * [ -cot(B+x) + CORR ]
581 // CORR = sgn_r * c * cot(B) * 1/[sin(B)*cos(B)]
583 // because csc^2(r) = csc^(|r|), and B approximate |r| to 6.5 bits.
586 // / (1/[sin(B)*cos(B)]) * tan(x)
587 // sgn_r * | -cot(B) + -------------------------------- +
594 // The values of tan(B), cot(B) and 1/(sin(B)*cos(B)) are
595 // calculated beforehand and stored in a table. Specifically,
596 // the table values are
598 // tan(B) as T_hi + T_lo;
599 // cot(B) as C_hi + C_lo;
600 // 1/[sin(B)*cos(B)] as SC_inv
602 // T_hi, C_hi are in double-precision memory format;
603 // T_lo, C_lo are in single-precision memory format;
604 // SC_inv is in extended-precision memory format.
606 // The value of tan(x) will be approximated by a short polynomial of
609 // tan(x) as x + x * P, where
610 // P = x^2 * (P2_1 + x^2 * (P2_2 + x^2 * P2_3))
612 // Because |x| <= 2^(-7), tan(B) + x approximates tan(B) + tan(x)
613 // to a relative accuracy better than 2^(-18). Thus, a good
614 // initial guess of 1/( tan(B) + tan(x) ) to initiate the iterative
617 // 1/(tan(B) + tan(x)) is approximately
619 // cot(B)/(1 + x*cot(B)) is approximately
620 // C_hi / ( 1 + C_hi * x ) is approximately
622 // C_hi * [ 1 - (C_hi * x) + (C_hi * x)^2 ]
624 // The calculation of -cot(r+c) therefore proceed as follows:
629 // V_hi := C_hi*(1 - Cx*(1 - Cx))
630 // P := xsq * (P1_1 + xsq*(P1_2 + xsq*P1_3))
631 // ...V_hi serves as an initial guess of 1/(tan(B) + tan(x))
632 // ...good to about 18 bits of accuracy
636 // ...D is a double precision denominator: tan(B) + tan(x)
638 // V_hi := V_hi + V_hi*(1 - V_hi*D)
639 // ....V_hi approximates 1/(tan(B)+tan(x)) to 40 bits
641 // V_lo := V_hi * ( [ (1 - V_hi*T_hi) - V_hi*tanx ]
642 // - V_hi*T_lo ) ...observe all order
643 // ...V_hi + V_lo approximates 1/(tan(B) + tan(x))
644 // ...to extra accuracy
646 // ... SC_inv(B) * (x + x*P)
647 // ... -cot(B) + ------------------------- + CORR
648 // ... tan(B) + (x + x*P)
650 // ... =-cot(B) + SC_inv(B)*(x + x*P)*(V_hi + V_lo) + CORR
654 // CORR := sgn_r * c * SC_inv * C_hi
656 // ...put the ingredients together to compute
657 // ... SC_inv(B) * (x + x*P)
658 // ... -cot(B) + ------------------------- + CORR
659 // ... tan(B) + (x + x*P)
661 // ... =-cot(B) + SC_inv(B)*(x + x*P)*(V_hi + V_lo) + CORR
663 // ... =-C_hi - C_lo + CORR +
664 // ... Sx * V_hi + Sx * V_lo + Sx * P *(V_hi + V_lo)
666 // CORR := CORR - C_lo
667 // tail := V_lo + P*(V_hi + V_lo)
668 // tail := Sx * tail + CORR
669 // tail := Sx * V_hi + tail
670 // C_hi := -sgn_r * C_hi
672 // ...C_hi + sgn_r*tail now approximates
673 // ...sgn_r*(-cot(B+x) + CORR) accurately
675 // Result := C_hi + sgn_r*tail in user-defined rounding control
676 // ...It is crucial that independent paths be fully
677 // ...exploited for performance's sake.
679 // 3. Implementation Notes
680 // =======================
682 // Table entries T_hi, T_lo; C_hi, C_lo; SC_inv
684 // Recall that 2^(-2) <= |r| <= pi/4;
686 // r = sgn_r * 2^k * 1.b_1 b_2 ... b_63
690 // B = 2^k * 1.b_1 b_2 b_3 b_4 b_5 1
692 // Thus, for k = -2, possible values of B are
694 // B = 2^(-2) * ( 1 + index/32 + 1/64 ),
695 // index ranges from 0 to 31
697 // For k = -1, however, since |r| <= pi/4 = 0.78...
698 // possible values of B are
700 // B = 2^(-1) * ( 1 + index/32 + 1/64 )
701 // index ranges from 0 to 19.
708 LOCAL_OBJECT_START(TANL_BASE_CONSTANTS)
711 data8 0xA2F9836E4E44152A, 0x00003FFE // two_by_pi
712 data8 0xC84D32B0CE81B9F1, 0x00004016 // P_0
713 data8 0xC90FDAA22168C235, 0x00003FFF // P_1
714 data8 0xECE675D1FC8F8CBB, 0x0000BFBD // P_2
715 data8 0xB7ED8FBBACC19C60, 0x0000BF7C // P_3
716 LOCAL_OBJECT_END(TANL_BASE_CONSTANTS)
718 LOCAL_OBJECT_START(tanl_table_2)
719 data8 0xC90FDAA22168C234, 0x00003FFE // PI_BY_4
720 data8 0xA397E5046EC6B45A, 0x00003FE7 // Inv_P_0
721 data8 0x8D848E89DBD171A1, 0x0000BFBF // d_1
722 data8 0xD5394C3618A66F8E, 0x0000BF7C // d_2
723 data4 0x3E800000 // two**-2
724 data4 0xBE800000 // -two**-2
725 data4 0x00000000 // pad
726 data4 0x00000000 // pad
727 LOCAL_OBJECT_END(tanl_table_2)
729 LOCAL_OBJECT_START(tanl_table_p1)
730 data8 0xAAAAAAAAAAAAAABD, 0x00003FFD // P1_1
731 data8 0x8888888888882E6A, 0x00003FFC // P1_2
732 data8 0xDD0DD0DD0F0177B6, 0x00003FFA // P1_3
733 data8 0xB327A440646B8C6D, 0x00003FF9 // P1_4
734 data8 0x91371B251D5F7D20, 0x00003FF8 // P1_5
735 data8 0xEB69A5F161C67914, 0x00003FF6 // P1_6
736 data8 0xBEDD37BE019318D2, 0x00003FF5 // P1_7
737 data8 0x9979B1463C794015, 0x00003FF4 // P1_8
738 data8 0x8EBD21A38C6EB58A, 0x00003FF3 // P1_9
739 LOCAL_OBJECT_END(tanl_table_p1)
741 LOCAL_OBJECT_START(tanl_table_q1)
742 data8 0xAAAAAAAAAAAAAAB4, 0x00003FFD // Q1_1
743 data8 0xB60B60B60B5FC93E, 0x00003FF9 // Q1_2
744 data8 0x8AB355E00C9BBFBF, 0x00003FF6 // Q1_3
745 data8 0xDDEBBC89CBEE3D4C, 0x00003FF2 // Q1_4
746 data8 0xB3548A685F80BBB6, 0x00003FEF // Q1_5
747 data8 0x913625604CED5BF1, 0x00003FEC // Q1_6
748 data8 0xF189D95A8EE92A83, 0x00003FE8 // Q1_7
749 LOCAL_OBJECT_END(tanl_table_q1)
751 LOCAL_OBJECT_START(tanl_table_p2)
752 data8 0xAAAAAAAAAAAB362F, 0x00003FFD // P2_1
753 data8 0x88888886E97A6097, 0x00003FFC // P2_2
754 data8 0xDD108EE025E716A1, 0x00003FFA // P2_3
755 LOCAL_OBJECT_END(tanl_table_p2)
757 LOCAL_OBJECT_START(tanl_table_tm2)
759 // Entries T_hi double-precision memory format
760 // Index = 0,1,...,31 B = 2^(-2)*(1+Index/32+1/64)
761 // Entries T_lo single-precision memory format
762 // Index = 0,1,...,31 B = 2^(-2)*(1+Index/32+1/64)
764 data8 0x3FD09BC362400794
765 data4 0x23A05C32, 0x00000000
766 data8 0x3FD124A9DFFBC074
767 data4 0x240078B2, 0x00000000
768 data8 0x3FD1AE235BD4920F
769 data4 0x23826B8E, 0x00000000
770 data8 0x3FD2383515E2701D
771 data4 0x22D31154, 0x00000000
772 data8 0x3FD2C2E463739C2D
773 data4 0x2265C9E2, 0x00000000
774 data8 0x3FD34E36AFEEA48B
775 data4 0x245C05EB, 0x00000000
776 data8 0x3FD3DA317DBB35D1
777 data4 0x24749F2D, 0x00000000
778 data8 0x3FD466DA67321619
779 data4 0x2462CECE, 0x00000000
780 data8 0x3FD4F4371F94A4D5
781 data4 0x246D0DF1, 0x00000000
782 data8 0x3FD5824D740C3E6D
783 data4 0x240A85B5, 0x00000000
784 data8 0x3FD611234CB1E73D
785 data4 0x23F96E33, 0x00000000
786 data8 0x3FD6A0BEAD9EA64B
787 data4 0x247C5393, 0x00000000
788 data8 0x3FD73125B804FD01
789 data4 0x241F3B29, 0x00000000
790 data8 0x3FD7C25EAB53EE83
791 data4 0x2479989B, 0x00000000
792 data8 0x3FD8546FE6640EED
793 data4 0x23B343BC, 0x00000000
794 data8 0x3FD8E75FE8AF1892
795 data4 0x241454D1, 0x00000000
796 data8 0x3FD97B3553928BDA
797 data4 0x238613D9, 0x00000000
798 data8 0x3FDA0FF6EB9DE4DE
799 data4 0x22859FA7, 0x00000000
800 data8 0x3FDAA5AB99ECF92D
801 data4 0x237A6D06, 0x00000000
802 data8 0x3FDB3C5A6D8F1796
803 data4 0x23952F6C, 0x00000000
804 data8 0x3FDBD40A9CFB8BE4
805 data4 0x2280FC95, 0x00000000
806 data8 0x3FDC6CC387943100
807 data4 0x245D2EC0, 0x00000000
808 data8 0x3FDD068CB736C500
809 data4 0x23C4AD7D, 0x00000000
810 data8 0x3FDDA16DE1DDBC31
811 data4 0x23D076E6, 0x00000000
812 data8 0x3FDE3D6EEB515A93
813 data4 0x244809A6, 0x00000000
814 data8 0x3FDEDA97E6E9E5F1
815 data4 0x220856C8, 0x00000000
816 data8 0x3FDF78F11963CE69
817 data4 0x244BE993, 0x00000000
818 data8 0x3FE00C417D635BCE
819 data4 0x23D21799, 0x00000000
820 data8 0x3FE05CAB1C302CD3
821 data4 0x248A1B1D, 0x00000000
822 data8 0x3FE0ADB9DB6A1FA0
823 data4 0x23D53E33, 0x00000000
824 data8 0x3FE0FF724A20BA81
825 data4 0x24DB9ED5, 0x00000000
826 data8 0x3FE151D9153FA6F5
827 data4 0x24E9E451, 0x00000000
828 LOCAL_OBJECT_END(tanl_table_tm2)
830 LOCAL_OBJECT_START(tanl_table_tm1)
832 // Entries T_hi double-precision memory format
833 // Index = 0,1,...,19 B = 2^(-1)*(1+Index/32+1/64)
834 // Entries T_lo single-precision memory format
835 // Index = 0,1,...,19 B = 2^(-1)*(1+Index/32+1/64)
837 data8 0x3FE1CEC4BA1BE39E
838 data4 0x24B60F9E, 0x00000000
839 data8 0x3FE277E45ABD9B2D
840 data4 0x248C2474, 0x00000000
841 data8 0x3FE324180272B110
842 data4 0x247B8311, 0x00000000
843 data8 0x3FE3D38B890E2DF0
844 data4 0x24C55751, 0x00000000
845 data8 0x3FE4866D46236871
846 data4 0x24E5BC34, 0x00000000
847 data8 0x3FE53CEE45E044B0
848 data4 0x24001BA4, 0x00000000
849 data8 0x3FE5F74282EC06E4
850 data4 0x24B973DC, 0x00000000
851 data8 0x3FE6B5A125DF43F9
852 data4 0x24895440, 0x00000000
853 data8 0x3FE77844CAFD348C
854 data4 0x240021CA, 0x00000000
855 data8 0x3FE83F6BCEED6B92
856 data4 0x24C45372, 0x00000000
857 data8 0x3FE90B58A34F3665
858 data4 0x240DAD33, 0x00000000
859 data8 0x3FE9DC522C1E56B4
860 data4 0x24F846CE, 0x00000000
861 data8 0x3FEAB2A427041578
862 data4 0x2323FB6E, 0x00000000
863 data8 0x3FEB8E9F9DD8C373
864 data4 0x24B3090B, 0x00000000
865 data8 0x3FEC709B65C9AA7B
866 data4 0x2449F611, 0x00000000
867 data8 0x3FED58F4ACCF8435
868 data4 0x23616A7E, 0x00000000
869 data8 0x3FEE480F97635082
870 data4 0x24C2FEAE, 0x00000000
871 data8 0x3FEF3E57F0ACC544
872 data4 0x242CE964, 0x00000000
873 data8 0x3FF01E20F7E06E4B
874 data4 0x2480D3EE, 0x00000000
875 data8 0x3FF0A1258A798A69
876 data4 0x24DB8967, 0x00000000
877 LOCAL_OBJECT_END(tanl_table_tm1)
879 LOCAL_OBJECT_START(tanl_table_cm2)
881 // Entries C_hi double-precision memory format
882 // Index = 0,1,...,31 B = 2^(-2)*(1+Index/32+1/64)
883 // Entries C_lo single-precision memory format
884 // Index = 0,1,...,31 B = 2^(-2)*(1+Index/32+1/64)
886 data8 0x400ED3E2E63EFBD0
887 data4 0x259D94D4, 0x00000000
888 data8 0x400DDDB4C515DAB5
889 data4 0x245F0537, 0x00000000
890 data8 0x400CF57ABE19A79F
891 data4 0x25D4EA9F, 0x00000000
892 data8 0x400C1A06D15298ED
893 data4 0x24AE40A0, 0x00000000
894 data8 0x400B4A4C164B2708
895 data4 0x25A5AAB6, 0x00000000
896 data8 0x400A855A5285B068
897 data4 0x25524F18, 0x00000000
898 data8 0x4009CA5A3FFA549F
899 data4 0x24C999C0, 0x00000000
900 data8 0x4009188A646AF623
901 data4 0x254FD801, 0x00000000
902 data8 0x40086F3C6084D0E7
903 data4 0x2560F5FD, 0x00000000
904 data8 0x4007CDD2A29A76EE
905 data4 0x255B9D19, 0x00000000
906 data8 0x400733BE6C8ECA95
907 data4 0x25CB021B, 0x00000000
908 data8 0x4006A07E1F8DDC52
909 data4 0x24AB4722, 0x00000000
910 data8 0x4006139BC298AD58
911 data4 0x252764E2, 0x00000000
912 data8 0x40058CABBAD7164B
913 data4 0x24DAF5DB, 0x00000000
914 data8 0x40050B4BAE31A5D3
915 data4 0x25EA20F4, 0x00000000
916 data8 0x40048F2189F85A8A
917 data4 0x2583A3E8, 0x00000000
918 data8 0x400417DAA862380D
919 data4 0x25DCC4CC, 0x00000000
920 data8 0x4003A52B1088FCFE
921 data4 0x2430A492, 0x00000000
922 data8 0x400336CCCD3527D5
923 data4 0x255F77CF, 0x00000000
924 data8 0x4002CC7F5760766D
925 data4 0x25DA0BDA, 0x00000000
926 data8 0x4002660711CE02E3
927 data4 0x256FF4A2, 0x00000000
928 data8 0x4002032CD37BBE04
929 data4 0x25208AED, 0x00000000
930 data8 0x4001A3BD7F050775
931 data4 0x24B72DD6, 0x00000000
932 data8 0x40014789A554848A
933 data4 0x24AB4DAA, 0x00000000
934 data8 0x4000EE65323E81B7
935 data4 0x2584C440, 0x00000000
936 data8 0x4000982721CF1293
937 data4 0x25C9428D, 0x00000000
938 data8 0x400044A93D415EEB
939 data4 0x25DC8482, 0x00000000
940 data8 0x3FFFE78FBD72C577
941 data4 0x257F5070, 0x00000000
942 data8 0x3FFF4AC375EFD28E
943 data4 0x23EBBF7A, 0x00000000
944 data8 0x3FFEB2AF60B52DDE
945 data4 0x22EECA07, 0x00000000
946 data8 0x3FFE1F1935204180
947 data4 0x24191079, 0x00000000
948 data8 0x3FFD8FCA54F7E60A
949 data4 0x248D3058, 0x00000000
950 LOCAL_OBJECT_END(tanl_table_cm2)
952 LOCAL_OBJECT_START(tanl_table_cm1)
954 // Entries C_hi double-precision memory format
955 // Index = 0,1,...,19 B = 2^(-1)*(1+Index/32+1/64)
956 // Entries C_lo single-precision memory format
957 // Index = 0,1,...,19 B = 2^(-1)*(1+Index/32+1/64)
959 data8 0x3FFCC06A79F6FADE
960 data4 0x239C7886, 0x00000000
961 data8 0x3FFBB91F891662A6
962 data4 0x250BD191, 0x00000000
963 data8 0x3FFABFB6529F155D
964 data4 0x256CC3E6, 0x00000000
965 data8 0x3FF9D3002E964AE9
966 data4 0x250843E3, 0x00000000
967 data8 0x3FF8F1EF89DCB383
968 data4 0x2277C87E, 0x00000000
969 data8 0x3FF81B937C87DBD6
970 data4 0x256DA6CF, 0x00000000
971 data8 0x3FF74F141042EDE4
972 data4 0x2573D28A, 0x00000000
973 data8 0x3FF68BAF1784B360
974 data4 0x242E489A, 0x00000000
975 data8 0x3FF5D0B57C923C4C
976 data4 0x2532D940, 0x00000000
977 data8 0x3FF51D88F418EF20
978 data4 0x253C7DD6, 0x00000000
979 data8 0x3FF4719A02F88DAE
980 data4 0x23DB59BF, 0x00000000
981 data8 0x3FF3CC6649DA0788
982 data4 0x252B4756, 0x00000000
983 data8 0x3FF32D770B980DB8
984 data4 0x23FE585F, 0x00000000
985 data8 0x3FF2945FE56C987A
986 data4 0x25378A63, 0x00000000
987 data8 0x3FF200BDB16523F6
988 data4 0x247BB2E0, 0x00000000
989 data8 0x3FF172358CE27778
990 data4 0x24446538, 0x00000000
991 data8 0x3FF0E873FDEFE692
992 data4 0x2514638F, 0x00000000
993 data8 0x3FF0632C33154062
994 data4 0x24A7FC27, 0x00000000
995 data8 0x3FEFC42EB3EF115F
996 data4 0x248FD0FE, 0x00000000
997 data8 0x3FEEC9E8135D26F6
998 data4 0x2385C719, 0x00000000
999 LOCAL_OBJECT_END(tanl_table_cm1)
1001 LOCAL_OBJECT_START(tanl_table_scim2)
1003 // Entries SC_inv in Swapped IEEE format (extended)
1004 // Index = 0,1,...,31 B = 2^(-2)*(1+Index/32+1/64)
1006 data8 0x839D6D4A1BF30C9E, 0x00004001
1007 data8 0x80092804554B0EB0, 0x00004001
1008 data8 0xF959F94CA1CF0DE9, 0x00004000
1009 data8 0xF3086BA077378677, 0x00004000
1010 data8 0xED154515CCD4723C, 0x00004000
1011 data8 0xE77909441C27CF25, 0x00004000
1012 data8 0xE22D037D8DDACB88, 0x00004000
1013 data8 0xDD2B2D8A89C73522, 0x00004000
1014 data8 0xD86E1A23BB2C1171, 0x00004000
1015 data8 0xD3F0E288DFF5E0F9, 0x00004000
1016 data8 0xCFAF16B1283BEBD5, 0x00004000
1017 data8 0xCBA4AFAA0D88DD53, 0x00004000
1018 data8 0xC7CE03CCCA67C43D, 0x00004000
1019 data8 0xC427BC820CA0DDB0, 0x00004000
1020 data8 0xC0AECD57F13D8CAB, 0x00004000
1021 data8 0xBD606C3871ECE6B1, 0x00004000
1022 data8 0xBA3A0A96A44C4929, 0x00004000
1023 data8 0xB7394F6FE5CCCEC1, 0x00004000
1024 data8 0xB45C12039637D8BC, 0x00004000
1025 data8 0xB1A0552892CB051B, 0x00004000
1026 data8 0xAF04432B6BA2FFD0, 0x00004000
1027 data8 0xAC862A237221235F, 0x00004000
1028 data8 0xAA2478AF5F00A9D1, 0x00004000
1029 data8 0xA7DDBB0C81E082BF, 0x00004000
1030 data8 0xA5B0987D45684FEE, 0x00004000
1031 data8 0xA39BD0F5627A8F53, 0x00004000
1032 data8 0xA19E3B036EC5C8B0, 0x00004000
1033 data8 0x9FB6C1F091CD7C66, 0x00004000
1034 data8 0x9DE464101FA3DF8A, 0x00004000
1035 data8 0x9C263139A8F6B888, 0x00004000
1036 data8 0x9A7B4968C27B0450, 0x00004000
1037 data8 0x98E2DB7E5EE614EE, 0x00004000
1038 LOCAL_OBJECT_END(tanl_table_scim2)
1040 LOCAL_OBJECT_START(tanl_table_scim1)
1042 // Entries SC_inv in Swapped IEEE format (extended)
1043 // Index = 0,1,...,19 B = 2^(-1)*(1+Index/32+1/64)
1045 data8 0x969F335C13B2B5BA, 0x00004000
1046 data8 0x93D446D9D4C0F548, 0x00004000
1047 data8 0x9147094F61B798AF, 0x00004000
1048 data8 0x8EF317CC758787AC, 0x00004000
1049 data8 0x8CD498B3B99EEFDB, 0x00004000
1050 data8 0x8AE82A7DDFF8BC37, 0x00004000
1051 data8 0x892AD546E3C55D42, 0x00004000
1052 data8 0x8799FEA9D15573C1, 0x00004000
1053 data8 0x86335F88435A4B4C, 0x00004000
1054 data8 0x84F4FB6E3E93A87B, 0x00004000
1055 data8 0x83DD195280A382FB, 0x00004000
1056 data8 0x82EA3D7FA4CB8C9E, 0x00004000
1057 data8 0x821B247C6861D0A8, 0x00004000
1058 data8 0x816EBED163E8D244, 0x00004000
1059 data8 0x80E42D9127E4CFC6, 0x00004000
1060 data8 0x807ABF8D28E64AFD, 0x00004000
1061 data8 0x8031EF26863B4FD8, 0x00004000
1062 data8 0x800960ADAE8C11FD, 0x00004000
1063 data8 0x8000E1475FDBEC21, 0x00004000
1064 data8 0x80186650A07791FA, 0x00004000
1065 LOCAL_OBJECT_END(tanl_table_scim1)
1068 Save_Norm_Arg = f8 // For input to reduction routine
1070 r = f8 // For output from reduction routine
1071 c = f9 // For output from reduction routine
1126 NEGTWO_TO_NEG14 = f76
1127 NEGTWO_TO_NEG33 = f77
1143 NEGTWO_TO_NEG2 = f93
1164 FR_inv_pi_2to63 = f113
1165 FR_rshf_2to64 = f114
1178 GR_exp_2_to_63 = r18
1179 GR_exp_2_to_24 = r19
1184 GR_exp_m2tom14 = r24
1186 GR_exp_m2tom33 = r26
1210 GR_Parameter_X = r54
1211 GR_Parameter_Y = r55
1212 GR_Parameter_RESULT = r56
1213 GR_Parameter_Tag = r57
1217 .global __libm_tanl#
1218 .global __libm_cotl#
1223 LOCAL_LIBM_ENTRY(cotl)
1226 alloc r32 = ar.pfs, 0,22,4,0
1227 movl GR_sig_inv_pi = 0xa2f9836e4e44152a // significand of 1/pi
1230 mov GR_exp_mask = 0x1ffff // Exponent mask
1231 movl GR_rshf_2to64 = 0x47e8000000000000 // 1.1000 2^(63+64)
1235 // Check for NatVals, Infs , NaNs, and Zeros
1237 getf.exp GR_signexp_x = Arg // Get sign and exponent of x
1238 fclass.m p6,p0 = Arg, 0x1E7 // Test for natval, nan, inf, zero
1242 addl table_base = @ltoff(TANL_BASE_CONSTANTS), gp // Pointer to table ptr
1243 fnorm.s1 Norm_Arg = Arg // Normalize x
1244 br.cond.sptk COMMON_PATH
1247 LOCAL_LIBM_END(cotl)
1252 GLOBAL_IEEE754_ENTRY(tanl)
1255 alloc r32 = ar.pfs, 0,22,4,0
1256 movl GR_sig_inv_pi = 0xa2f9836e4e44152a // significand of 1/pi
1259 mov GR_exp_mask = 0x1ffff // Exponent mask
1260 movl GR_rshf_2to64 = 0x47e8000000000000 // 1.1000 2^(63+64)
1264 // Check for NatVals, Infs , NaNs, and Zeros
1266 getf.exp GR_signexp_x = Arg // Get sign and exponent of x
1267 fclass.m p6,p0 = Arg, 0x1E7 // Test for natval, nan, inf, zero
1271 addl table_base = @ltoff(TANL_BASE_CONSTANTS), gp // Pointer to table ptr
1272 fnorm.s1 Norm_Arg = Arg // Normalize x
1276 // Common path for both tanl and cotl
1279 setf.sig FR_inv_pi_2to63 = GR_sig_inv_pi // Form 1/pi * 2^63
1280 fclass.m p9, p0 = Arg, 0x0b // Test x denormal
1281 mov GR_exp_2tom64 = 0xffff - 64 // Scaling constant to compute N
1284 setf.d FR_rshf_2to64 = GR_rshf_2to64 // Form const 1.1000 * 2^(63+64)
1285 movl GR_rshf = 0x43e8000000000000 // Form const 1.1000 * 2^63
1289 // Check for everything - if false, then must be pseudo-zero or pseudo-nan.
1290 // Branch out to deal with special values.
1293 fclass.nm p7,p0 = Arg, 0x1FF // Test x unsupported
1294 mov GR_exp_2_to_63 = 0xffff + 63 // Exponent of 2^63
1297 ld8 table_base = [table_base] // Get pointer to constant table
1298 fms.s1 mOne = f0, f0, f1
1299 (p6) br.cond.spnt TANL_SPECIAL // Branch if x natval, nan, inf, zero
1304 setf.sig fp_tmp = gr_tmp // Make a constant so fmpy produces inexact
1305 mov GR_exp_2_to_24 = 0xffff + 24 // Exponent of 2^24
1306 (p9) br.cond.spnt TANL_DENORMAL // Branch if x denormal
1311 // Return to here if x denormal
1313 // Do fcmp to generate Denormal exception
1314 // - can't do FNORM (will generate Underflow when U is unmasked!)
1315 // Branch out to deal with unsupporteds values.
1317 setf.exp FR_2tom64 = GR_exp_2tom64 // Form 2^-64 for scaling N_float
1318 fcmp.eq.s0 p0, p6 = Arg, f1 // Dummy to flag denormals
1319 add table_ptr1 = 0, table_base // Point to tanl_table_1
1322 setf.d FR_rshf = GR_rshf // Form right shift const 1.1000 * 2^63
1323 add table_ptr2 = 80, table_base // Point to tanl_table_2
1324 (p7) br.cond.spnt TANL_UNSUPPORTED // Branch if x unsupported type
1329 and GR_exp_x = GR_exp_mask, GR_signexp_x // Get exponent of x
1330 fmpy.s1 Save_Norm_Arg = Norm_Arg, f1 // Save x if large arg reduction
1331 dep.z bmask1 = 0x7c, 56, 8 // Form mask to get 5 msb of r
1332 // bmask1 = 0x7c00000000000000
1337 // Decide about the paths to take:
1338 // Set PR_6 if |Arg| >= 2**63
1339 // Set PR_9 if |Arg| < 2**24 - CASE 1 OR 2
1340 // OTHERWISE Set PR_8 - CASE 3 OR 4
1342 // Branch out if the magnitude of the input argument is >= 2^63
1343 // - do this branch before the next.
1345 ldfe two_by_PI = [table_ptr1],16 // Load 2/pi
1347 dep.z bmask2 = 0x41, 57, 7 // Form mask to OR to produce B
1348 // bmask2 = 0x8200000000000000
1351 ldfe PI_BY_4 = [table_ptr2],16 // Load pi/4
1352 cmp.ge p6,p0 = GR_exp_x, GR_exp_2_to_63 // Is |x| >= 2^63
1353 (p6) br.cond.spnt TANL_ARG_TOO_LARGE // Branch if |x| >= 2^63
1358 ldfe P_0 = [table_ptr1],16 // Load P_0
1359 ldfe Inv_P_0 = [table_ptr2],16 // Load Inv_P_0
1365 ldfe P_1 = [table_ptr1],16 // Load P_1
1366 fmerge.s Abs_Arg = f0, Norm_Arg // Get |x|
1367 mov GR_exp_m2tom33 = 0x2ffff - 33 // Form signexp of -2^-33
1370 ldfe d_1 = [table_ptr2],16 // Load d_1 for 2^24 <= |x| < 2^63
1372 mov GR_exp_2tom33 = 0xffff - 33 // Form signexp of 2^-33
1377 ldfe P_2 = [table_ptr1],16 // Load P_2
1378 ldfe d_2 = [table_ptr2],16 // Load d_2 for 2^24 <= |x| < 2^63
1379 cmp.ge p8,p0 = GR_exp_x, GR_exp_2_to_24 // Is |x| >= 2^24
1383 // Use special scaling to right shift so N=Arg * 2/pi is in rightmost bits
1384 // Branch to Cases 3 or 4 if Arg <= -2**24 or Arg >= 2**24
1386 ldfe P_3 = [table_ptr1],16 // Load P_3
1387 fma.s1 N_fix = Norm_Arg, FR_inv_pi_2to63, FR_rshf_2to64
1388 (p8) br.cond.spnt TANL_LARGER_ARG // Branch if 2^24 <= |x| < 2^63
1392 // Here if 0 < |x| < 2^24
1393 // ARGUMENT REDUCTION CODE - CASE 1 and 2
1396 setf.exp TWO_TO_NEG33 = GR_exp_2tom33 // Form 2^-33
1397 setf.exp NEGTWO_TO_NEG33 = GR_exp_m2tom33 // Form -2^-33
1398 fmerge.s r = Norm_Arg,Norm_Arg // Assume r=x, ok if |x| < pi/4
1403 // If |Arg| < pi/4, set PR_8, else pi/4 <=|Arg| < 2^24 - set PR_9.
1405 // Case 2: Convert integer N_fix back to normalized floating-point value.
1407 getf.sig sig_r = Norm_Arg // Get sig_r if 1/4 <= |x| < pi/4
1408 fcmp.lt.s1 p8,p9= Abs_Arg,PI_BY_4 // Test |x| < pi/4
1409 mov GR_exp_2tom2 = 0xffff - 2 // Form signexp of 2^-2
1412 ldfps TWO_TO_NEG2, NEGTWO_TO_NEG2 = [table_ptr2] // Load 2^-2, -2^-2
1413 fms.s1 N = N_fix, FR_2tom64, FR_rshf // Use scaling to get N floated
1414 mov N_fix_gr = r0 // Assume N=0, ok if |x| < pi/4
1419 // Case 1: Is |r| < 2**(-2).
1420 // Arg is the same as r in this case.
1424 // Case 2: Place integer part of N in GP register.
1426 (p9) getf.sig N_fix_gr = N_fix
1427 fmerge.s c = f0, f0 // Assume c=0, ok if |x| < pi/4
1428 cmp.lt p10, p0 = GR_exp_x, GR_exp_2tom2 // Test if |x| < 1/4
1433 setf.sig B_mask1 = bmask1 // Form mask to get 5 msb of r
1435 mov exp_r = GR_exp_x // Get exp_r if 1/4 <= |x| < pi/4
1438 setf.sig B_mask2 = bmask2 // Form mask to form B from r
1439 (p10) br.cond.spnt TANL_SMALL_R // Branch if 0 < |x| < 1/4
1440 (p8) br.cond.spnt TANL_NORMAL_R // Branch if 1/4 <= |x| < pi/4
1444 // Here if pi/4 <= |x| < 2^24
1446 // Case 1: PR_3 is only affected when PR_1 is set.
1449 // Case 2: w = N * P_2
1450 // Case 2: s_val = -N * P_1 + Arg
1455 fnma.s1 s_val = N, P_1, Norm_Arg
1460 fmpy.s1 w = N, P_2 // w = N * P_2 for |s| >= 2^-33
1465 // Case 2_reduce: w = N * P_3 (change sign)
1468 fmpy.s1 w2 = N, P_3 // w = N * P_3 for |s| < 2^-33
1473 // Case 1_reduce: r = s + w (change sign)
1476 fsub.s1 r = s_val, w // r = s_val - w for |s| >= 2^-33
1481 // Case 2_reduce: U_1 = N * P_2 + w
1484 fma.s1 U_1 = N, P_2, w2 // U_1 = N * P_2 + w for |s| < 2^-33
1490 // Decide between case_1 and case_2 reduce:
1491 // Case 1_reduce: |s| >= 2**(-33)
1492 // Case 2_reduce: |s| < 2**(-33)
1496 fcmp.lt.s1 p9, p8 = s_val, TWO_TO_NEG33
1503 (p9) fcmp.gt.s1 p9, p8 = s_val, NEGTWO_TO_NEG33
1508 // Case 1_reduce: c = s - r
1511 fsub.s1 c = s_val, r // c = s_val - r for |s| >= 2^-33
1516 // Case 2_reduce: r is complete here - continue to calculate c .
1520 (p9) fsub.s1 r = s_val, U_1
1525 (p9) fms.s1 U_2 = N, P_2, U_1
1531 // Case 1_reduce: Is |r| < 2**(-2), if so set PR_10
1542 (p8) fcmp.lt.unc.s1 p10, p13 = r, TWO_TO_NEG2
1548 (p8) getf.sig sig_r = r // Get signif of r if |s| >= 2^-33
1555 (p8) getf.exp exp_r = r // Extract signexp of r if |s| >= 2^-33
1556 (p10) fcmp.gt.s1 p10, p13 = r, NEGTWO_TO_NEG2
1561 // Case 1_reduce: c is complete here.
1562 // Case 1: Branch to SMALL_R or NORMAL_R.
1563 // c = c + w (w has not been negated.)
1566 (p8) fsub.s1 c = c, w // c = c - w for |s| >= 2^-33
1571 (p10) br.cond.spnt TANL_SMALL_R // Branch if pi/4 < |x| < 2^24 and |r|<1/4
1572 (p13) br.cond.sptk TANL_NORMAL_R_A // Branch if pi/4 < |x| < 2^24 and |r|>=1/4
1577 // Here if pi/4 < |x| < 2^24 and |s| < 2^-33
1579 // Is i_1 = lsb of N_fix_gr even or odd?
1580 // if i_1 == 0, set p11, else set p12.
1584 fsub.s1 s_val = s_val, r
1585 add N_fix_gr = N_fix_gr, cot_flag // N = N + 1 (for cotl)
1591 // U_2 = N * P_2 - U_1
1592 // Not needed until later.
1594 fadd.s1 U_2 = U_2, w2
1607 // c is complete here
1608 // Argument reduction ends here.
1613 tbit.z p11, p12 = N_fix_gr, 0 ;; // Set p11 if N even, p12 if odd
1618 (p12) frcpa.s1 S_hi,p0 = f1, r
1623 fsub.s1 c = s_val, U_1
1629 add table_ptr1 = 160, table_base ;; // Point to tanl_table_p1
1630 ldfe P1_1 = [table_ptr1],144
1634 // Load P1_1 and point to Q1_1 .
1637 ldfe Q1_1 = [table_ptr1]
1639 // N even: rsq = r * Z
1640 // N odd: S_hi = frcpa(r)
1642 (p12) fmerge.ns S_hi = S_hi, S_hi
1651 (p9) fsub.s1 c = c, U_2
1656 (p12) fma.s1 poly1 = S_hi, r, f1
1662 // N odd: Change sign of S_hi
1664 (p11) fmpy.s1 rsq = rsq, P1_1
1669 (p12) fma.s1 S_hi = S_hi, poly1, S_hi
1675 // N even: rsq = rsq * P1_1
1676 // N odd: poly1 = 1.0 + S_hi * r 16 bits partial account for necessary
1678 (p11) fma.s1 Poly = r, rsq, c
1684 // N even: Poly = c + r * rsq
1685 // N odd: S_hi = S_hi + S_hi*poly1 16 bits account for necessary
1687 (p12) fma.s1 poly1 = S_hi, r, f1
1688 (p11) tbit.z.unc p14, p15 = cot_flag, 0 ;; // p14=1 for tanl; p15=1 for cotl
1693 // N even: Result = Poly + r
1694 // N odd: poly1 = 1.0 + S_hi * r 32 bits partial
1696 (p14) fadd.s0 Result = r, Poly // for tanl
1701 (p15) fms.s0 Result = r, mOne, Poly // for cotl
1708 (p12) fma.s1 S_hi = S_hi, poly1, S_hi
1714 // N even: Result1 = Result + r
1715 // N odd: S_hi = S_hi * poly1 + S_hi 32 bits
1717 (p12) fma.s1 poly1 = S_hi, r, f1
1723 // N odd: poly1 = S_hi * r + 1.0 64 bits partial
1725 (p12) fma.s1 S_hi = S_hi, poly1, S_hi
1731 // N odd: poly1 = S_hi * poly + 1.0 64 bits
1733 (p12) fma.s1 poly1 = S_hi, r, f1
1739 // N odd: poly1 = S_hi * r + 1.0
1741 (p12) fma.s1 poly1 = S_hi, c, poly1
1747 // N odd: poly1 = S_hi * c + poly1
1749 (p12) fmpy.s1 S_lo = S_hi, poly1
1755 // N odd: S_lo = S_hi * poly1
1757 (p12) fma.s1 S_lo = Q1_1, r, S_lo
1758 (p12) tbit.z.unc p14, p15 = cot_flag, 0 // p14=1 for tanl; p15=1 for cotl
1763 // N odd: Result = S_hi + S_lo
1765 fmpy.s0 fp_tmp = fp_tmp, fp_tmp // Dummy mult to set inexact
1771 // N odd: S_lo = S_lo + Q1_1 * r
1773 (p14) fadd.s0 Result = S_hi, S_lo // for tanl
1778 (p15) fms.s0 Result = S_hi, mOne, S_lo // for cotl
1779 br.ret.sptk b0 ;; // Exit for pi/4 <= |x| < 2^24 and |s| < 2^-33
1784 // Here if 2^24 <= |x| < 2^63
1786 // ARGUMENT REDUCTION CODE - CASE 3 and 4
1790 mov GR_exp_2tom14 = 0xffff - 14 // Form signexp of 2^-14
1791 mov GR_exp_m2tom14 = 0x2ffff - 14 // Form signexp of -2^-14
1792 fmpy.s1 N_0 = Norm_Arg, Inv_P_0
1797 setf.exp TWO_TO_NEG14 = GR_exp_2tom14 // Form 2^-14
1798 setf.exp NEGTWO_TO_NEG14 = GR_exp_m2tom14// Form -2^-14
1805 // Adjust table_ptr1 to beginning of table.
1806 // N_0 = Arg * Inv_P_0
1809 add table_ptr2 = 144, table_base ;; // Point to 2^-2
1810 ldfps TWO_TO_NEG2, NEGTWO_TO_NEG2 = [table_ptr2]
1816 // N_0_fix = integer part of N_0 .
1819 // Make N_0 the integer part.
1823 fcvt.fx.s1 N_0_fix = N_0
1827 setf.sig B_mask1 = bmask1 // Form mask to get 5 msb of r
1828 fcvt.xf N_0 = N_0_fix
1832 setf.sig B_mask2 = bmask2 // Form mask to form B from r
1833 fnma.s1 ArgPrime = N_0, P_0, Norm_Arg
1838 fmpy.s1 w = N_0, d_1
1842 // ArgPrime = -N_0 * P_0 + Arg
1846 // N = ArgPrime * 2/pi
1848 // fcvt.fx.s1 N_fix = N
1849 // Use special scaling to right shift so N=Arg * 2/pi is in rightmost bits
1850 // Branch to Cases 3 or 4 if Arg <= -2**24 or Arg >= 2**24
1853 fma.s1 N_fix = ArgPrime, FR_inv_pi_2to63, FR_rshf_2to64
1857 // Convert integer N_fix back to normalized floating-point value.
1860 fms.s1 N = N_fix, FR_2tom64, FR_rshf // Use scaling to get N floated
1866 // N is the integer part of the reduced-reduced argument.
1867 // Put the integer in a GP register.
1870 getf.sig N_fix_gr = N_fix
1877 // s_val = -N*P_1 + ArgPrime
1882 fnma.s1 s_val = N, P_1, ArgPrime
1887 fnma.s1 w = N, P_2, w
1892 // Case 4: V_hi = N * P_2
1893 // Case 4: U_hi = N_0 * d_1
1896 fmpy.s1 V_hi = N, P_2 // V_hi = N * P_2 for |s| < 2^-14
1901 fmpy.s1 U_hi = N_0, d_1 // U_hi = N_0 * d_1 for |s| < 2^-14
1906 // Case 3: r = s_val + w (Z complete)
1907 // Case 4: w = N * P_3
1910 fadd.s1 r = s_val, w // r = s_val + w for |s| >= 2^-14
1915 fmpy.s1 w2 = N, P_3 // w = N * P_3 for |s| < 2^-14
1920 // Case 4: A = U_hi + V_hi
1921 // Note: Worry about switched sign of V_hi, so subtract instead of add.
1922 // Case 4: V_lo = -N * P_2 - V_hi (U_hi is in place of V_hi in writeup)
1923 // Note: the (-) is still missing for V_hi.
1926 fsub.s1 A = U_hi, V_hi // A = U_hi - V_hi for |s| < 2^-14
1931 fnma.s1 V_lo = N, P_2, V_hi // V_lo = V_hi - N * P_2 for |s| < 2^-14
1936 // Decide between case 3 and 4:
1937 // Case 3: |s| >= 2**(-14) Set p10
1938 // Case 4: |s| < 2**(-14) Set p11
1940 // Case 4: U_lo = N_0 * d_1 - U_hi
1943 fms.s1 U_lo = N_0, d_1, U_hi // U_lo = N_0*d_1 - U_hi for |s| < 2^-14
1948 fcmp.lt.s1 p11, p10 = s_val, TWO_TO_NEG14
1953 // Case 4: We need abs of both U_hi and V_hi - dont
1954 // worry about switched sign of V_hi.
1957 fabs V_hiabs = V_hi // |V_hi| for |s| < 2^-14
1962 (p11) fcmp.gt.s1 p11, p10 = s_val, NEGTWO_TO_NEG14
1967 // Case 3: c = s_val - r
1970 fabs U_hiabs = U_hi // |U_hi| for |s| < 2^-14
1975 fsub.s1 c = s_val, r // c = s_val - r for |s| >= 2^-14
1980 // For Case 3, |s| >= 2^-14, determine if |r| < 1/4
1982 // Case 4: C_hi = s_val + A
1986 (p11) fadd.s1 C_hi = s_val, A // C_hi = s_val + A for |s| < 2^-14
1991 (p10) fcmp.lt.unc.s1 p14, p15 = r, TWO_TO_NEG2
1997 getf.sig sig_r = r // Get signif of r if |s| >= 2^-33
2003 // Case 4: t = U_lo + V_lo
2005 getf.exp exp_r = r // Extract signexp of r if |s| >= 2^-33
2006 (p11) fadd.s1 t = U_lo, V_lo // t = U_lo + V_lo for |s| < 2^-14
2011 (p14) fcmp.gt.s1 p14, p15 = r, NEGTWO_TO_NEG2
2016 // Case 3: c = (s - r) + w (c complete)
2019 (p10) fadd.s1 c = c, w // c = c + w for |s| >= 2^-14
2024 (p14) br.cond.spnt TANL_SMALL_R // Branch if 2^24 <= |x| < 2^63 and |r|< 1/4
2025 (p15) br.cond.sptk TANL_NORMAL_R_A // Branch if 2^24 <= |x| < 2^63 and |r|>=1/4
2030 // Here if 2^24 <= |x| < 2^63 and |s| < 2^-14 >>>>>>> Case 4.
2032 // Case 4: Set P_12 if U_hiabs >= V_hiabs
2033 // Case 4: w = w + N_0 * d_2
2034 // Note: the (-) is now incorporated in w .
2036 add table_ptr1 = 160, table_base // Point to tanl_table_p1
2037 fcmp.ge.unc.s1 p12, p13 = U_hiabs, V_hiabs
2042 fms.s1 w2 = N_0, d_2, w2
2047 // Case 4: C_lo = s_val - C_hi
2049 ldfe P1_1 = [table_ptr1], 16 // Load P1_1
2050 fsub.s1 C_lo = s_val, C_hi
2056 // Case 4: a = U_hi - A
2057 // a = V_hi - A (do an add to account for missing (-) on V_hi
2060 ldfe P1_2 = [table_ptr1], 128 // Load P1_2
2061 (p12) fsub.s1 a = U_hi, A
2066 (p13) fadd.s1 a = V_hi, A
2071 // Case 4: t = U_lo + V_lo + w
2073 ldfe Q1_1 = [table_ptr1], 16 // Load Q1_1
2079 // Case 4: a = (U_hi - A) + V_hi
2080 // a = (V_hi - A) + U_hi
2081 // In each case account for negative missing form V_hi .
2084 ldfe Q1_2 = [table_ptr1], 16 // Load Q1_2
2085 (p12) fsub.s1 a = a, V_hi
2090 (p13) fsub.s1 a = U_hi, a
2096 // Case 4: C_lo = (s_val - C_hi) + A
2100 fadd.s1 C_lo = C_lo, A
2104 // Case 4: t = t + a
2113 // Case 4: C_lo = C_lo + t
2114 // Case 4: r = C_hi + C_lo
2117 fadd.s1 C_lo = C_lo, t
2124 fadd.s1 r = C_hi, C_lo
2130 // Case 4: c = C_hi - r
2140 add N_fix_gr = N_fix_gr, cot_flag // N = N + 1 (for cotl)
2144 // Case 4: c = c + C_lo finished.
2146 // Is i_1 = lsb of N_fix_gr even or odd?
2147 // if i_1 == 0, set PR_11, else set PR_12.
2151 fadd.s1 c = c , C_lo
2152 tbit.z p11, p12 = N_fix_gr, 0
2156 // r and c have been computed.
2159 (p12) frcpa.s1 S_hi, p0 = f1, r
2165 // N odd: Change sign of S_hi
2167 (p11) fma.s1 Poly = rsq, P1_2, P1_1
2172 (p12) fma.s1 P = rsq, Q1_2, Q1_1
2178 // N odd: Result = S_hi + S_lo (User supplied rounding mode for C1)
2180 fmpy.s0 fp_tmp = fp_tmp, fp_tmp // Dummy mult to set inexact
2186 // N even: rsq = r * r
2187 // N odd: S_hi = frcpa(r)
2189 (p12) fmerge.ns S_hi = S_hi, S_hi
2195 // N even: rsq = rsq * P1_2 + P1_1
2196 // N odd: poly1 = 1.0 + S_hi * r 16 bits partial account for necessary
2198 (p11) fmpy.s1 Poly = rsq, Poly
2203 (p12) fma.s1 poly1 = S_hi, r,f1
2204 (p11) tbit.z.unc p14, p15 = cot_flag, 0 // p14=1 for tanl; p15=1 for cotl
2209 // N even: Poly = Poly * rsq
2210 // N odd: S_hi = S_hi + S_hi*poly1 16 bits account for necessary
2212 (p11) fma.s1 Poly = r, Poly, c
2217 (p12) fma.s1 S_hi = S_hi, poly1, S_hi
2223 // N odd: S_hi = S_hi * poly1 + S_hi 32 bits
2225 (p14) fadd.s0 Result = r, Poly // for tanl
2229 .pred.rel "mutex",p15,p12
2232 (p15) fms.s0 Result = r, mOne, Poly // for cotl
2237 (p12) fma.s1 poly1 = S_hi, r, f1
2243 // N even: Poly = Poly * r + c
2244 // N odd: poly1 = 1.0 + S_hi * r 32 bits partial
2246 (p12) fma.s1 S_hi = S_hi, poly1, S_hi
2251 (p12) fma.s1 poly1 = S_hi, r, f1
2257 // N even: Result = Poly + r (Rounding mode S0)
2258 // N odd: poly1 = S_hi * r + 1.0 64 bits partial
2260 (p12) fma.s1 S_hi = S_hi, poly1, S_hi
2266 // N odd: poly1 = S_hi * poly + S_hi 64 bits
2268 (p12) fma.s1 poly1 = S_hi, r, f1
2274 // N odd: poly1 = S_hi * r + 1.0
2276 (p12) fma.s1 poly1 = S_hi, c, poly1
2282 // N odd: poly1 = S_hi * c + poly1
2284 (p12) fmpy.s1 S_lo = S_hi, poly1
2290 // N odd: S_lo = S_hi * poly1
2292 (p12) fma.s1 S_lo = P, r, S_lo
2293 (p12) tbit.z.unc p14, p15 = cot_flag, 0 ;; // p14=1 for tanl; p15=1 for cotl
2298 (p14) fadd.s0 Result = S_hi, S_lo // for tanl
2304 // N odd: S_lo = S_lo + r * P
2306 (p15) fms.s0 Result = S_hi, mOne, S_lo // for cotl
2307 br.ret.sptk b0 ;; // Exit for 2^24 <= |x| < 2^63 and |s| < 2^-14
2312 // Here if |r| < 1/4
2313 // r and c have been computed.
2314 // *****************************************************************
2315 // *****************************************************************
2316 // *****************************************************************
2317 // N odd: S_hi = frcpa(r)
2318 // Get [i_1] - lsb of N_fix_gr. Set p11 if N even, p12 if N odd.
2319 // N even: rsq = r * r
2321 add table_ptr1 = 160, table_base // Point to tanl_table_p1
2322 frcpa.s1 S_hi, p0 = f1, r // S_hi for N odd
2323 add N_fix_gr = N_fix_gr, cot_flag // N = N + 1 (for cotl)
2326 add table_ptr2 = 400, table_base // Point to Q1_7
2333 ldfe P1_1 = [table_ptr1], 16
2335 ldfe P1_2 = [table_ptr1], 16
2336 tbit.z p11, p12 = N_fix_gr, 0
2342 ldfe P1_3 = [table_ptr1], 96
2349 (p11) ldfe P1_9 = [table_ptr1], -16
2350 (p12) fmerge.ns S_hi = S_hi, S_hi
2355 (p11) fmpy.s1 r_to_the_8 = rsq, rsq
2361 // N even: Poly2 = P1_7 + Poly2 * rsq
2362 // N odd: poly2 = Q1_5 + poly2 * rsq
2365 (p11) ldfe P1_8 = [table_ptr1], -16
2366 (p11) fadd.s1 CORR = rsq, f1
2372 // N even: Poly1 = P1_2 + P1_3 * rsq
2373 // N odd: poly1 = 1.0 + S_hi * r
2374 // 16 bits partial account for necessary (-1)
2377 (p11) ldfe P1_7 = [table_ptr1], -16
2379 (p11) ldfe P1_6 = [table_ptr1], -16
2385 // N even: Poly1 = P1_1 + Poly1 * rsq
2386 // N odd: S_hi = S_hi + S_hi * poly1) 16 bits account for necessary
2389 // N even: Poly2 = P1_5 + Poly2 * rsq
2390 // N odd: poly2 = Q1_3 + poly2 * rsq
2393 (p11) ldfe P1_5 = [table_ptr1], -16
2394 (p11) fmpy.s1 r_to_the_8 = r_to_the_8, r_to_the_8
2399 (p12) fma.s1 poly1 = S_hi, r, f1
2405 // N even: Poly1 = Poly1 * rsq
2406 // N odd: poly1 = 1.0 + S_hi * r 32 bits partial
2410 // N even: CORR = CORR * c
2411 // N odd: S_hi = S_hi * poly1 + S_hi 32 bits
2415 // N even: Poly2 = P1_6 + Poly2 * rsq
2416 // N odd: poly2 = Q1_4 + poly2 * rsq
2420 (p11) ldfe P1_4 = [table_ptr1], -16
2422 (p11) fmpy.s1 CORR = CORR, c
2428 (p11) fma.s1 Poly1 = P1_3, rsq, P1_2
2432 (p12) ldfe Q1_7 = [table_ptr2], -16
2433 (p12) fma.s1 S_hi = S_hi, poly1, S_hi
2437 (p12) ldfe Q1_6 = [table_ptr2], -16
2438 (p11) fma.s1 Poly2 = P1_9, rsq, P1_8
2442 (p12) ldfe Q1_5 = [table_ptr2], -16 ;;
2443 (p12) ldfe Q1_4 = [table_ptr2], -16
2447 (p12) ldfe Q1_3 = [table_ptr2], -16
2449 // N even: Poly2 = P1_8 + P1_9 * rsq
2450 // N odd: poly2 = Q1_6 + Q1_7 * rsq
2452 (p11) fma.s1 Poly1 = Poly1, rsq, P1_1
2456 (p12) ldfe Q1_2 = [table_ptr2], -16
2457 (p12) fma.s1 poly1 = S_hi, r, f1
2461 (p12) ldfe Q1_1 = [table_ptr2], -16
2462 (p11) fma.s1 Poly2 = Poly2, rsq, P1_7
2468 // N even: CORR = rsq + 1
2469 // N even: r_to_the_8 = rsq * rsq
2471 (p11) fmpy.s1 Poly1 = Poly1, rsq
2476 (p12) fma.s1 S_hi = S_hi, poly1, S_hi
2481 (p12) fma.s1 poly2 = Q1_7, rsq, Q1_6
2486 (p11) fma.s1 Poly2 = Poly2, rsq, P1_6
2491 (p12) fma.s1 poly1 = S_hi, r, f1
2496 (p12) fma.s1 poly2 = poly2, rsq, Q1_5
2501 (p11) fma.s1 Poly2= Poly2, rsq, P1_5
2506 (p12) fma.s1 S_hi = S_hi, poly1, S_hi
2511 (p12) fma.s1 poly2 = poly2, rsq, Q1_4
2517 // N even: r_to_the_8 = r_to_the_8 * r_to_the_8
2518 // N odd: poly1 = S_hi * r + 1.0 64 bits partial
2520 (p11) fma.s1 Poly2 = Poly2, rsq, P1_4
2526 // N even: Poly = CORR + Poly * r
2527 // N odd: P = Q1_1 + poly2 * rsq
2529 (p12) fma.s1 poly1 = S_hi, r, f1
2534 (p12) fma.s1 poly2 = poly2, rsq, Q1_3
2540 // N even: Poly2 = P1_4 + Poly2 * rsq
2541 // N odd: poly2 = Q1_2 + poly2 * rsq
2543 (p11) fma.s1 Poly = Poly2, r_to_the_8, Poly1
2548 (p12) fma.s1 poly1 = S_hi, c, poly1
2553 (p12) fma.s1 poly2 = poly2, rsq, Q1_2
2560 // N even: Poly = Poly1 + Poly2 * r_to_the_8
2561 // N odd: S_hi = S_hi * poly1 + S_hi 64 bits
2563 (p11) fma.s1 Poly = Poly, r, CORR
2569 // N even: Result = r + Poly (User supplied rounding mode)
2570 // N odd: poly1 = S_hi * c + poly1
2572 (p12) fmpy.s1 S_lo = S_hi, poly1
2573 (p11) tbit.z.unc p14, p15 = cot_flag, 0 // p14=1 for tanl; p15=1 for cotl
2577 (p12) fma.s1 P = poly2, rsq, Q1_1
2583 // N odd: poly1 = S_hi * r + 1.0
2586 // N odd: S_lo = S_hi * poly1
2588 (p14) fadd.s0 Result = Poly, r // for tanl
2593 (p15) fms.s0 Result = Poly, mOne, r // for cotl
2600 // N odd: S_lo = Q1_1 * c + S_lo
2602 (p12) fma.s1 S_lo = Q1_1, c, S_lo
2607 fmpy.s0 fp_tmp = fp_tmp, fp_tmp // Dummy mult to set inexact
2613 // N odd: Result = S_lo + r * P
2615 (p12) fma.s1 Result = P, r, S_lo
2616 (p12) tbit.z.unc p14, p15 = cot_flag, 0 ;; // p14=1 for tanl; p15=1 for cotl
2620 // N odd: Result = Result + S_hi (user supplied rounding mode)
2624 (p14) fadd.s0 Result = Result, S_hi // for tanl
2629 (p15) fms.s0 Result = Result, mOne, S_hi // for cotl
2630 br.ret.sptk b0 ;; // Exit |r| < 1/4 path
2635 // Here if 1/4 <= |x| < pi/4 or if |x| >= 2^63 and |r| >= 1/4
2636 // *******************************************************************
2637 // *******************************************************************
2638 // *******************************************************************
2640 // r and c have been computed.
2650 // Enter here if pi/4 <= |x| < 2^63 and |r| >= 1/4
2651 // Get the 5 bits or r for the lookup. 1.xxxxx ....
2653 add table_ptr1 = 416, table_base // Point to tanl_table_p2
2654 mov GR_exp_2tom65 = 0xffff - 65 // Scaling constant for B
2655 extr.u lookup = sig_r, 58, 5
2660 ldfe P2_1 = [table_ptr1], 16
2661 setf.exp TWO_TO_NEG65 = GR_exp_2tom65 // 2^-65 for scaling B if exp_r=-2
2662 add N_fix_gr = N_fix_gr, cot_flag // N = N + 1 (for cotl)
2666 .pred.rel "mutex",p11,p12
2667 // B = 2^63 * 1.xxxxx 100...0
2669 ldfe P2_2 = [table_ptr1], 16
2671 mov table_offset = 512 // Assume table offset is 512
2676 ldfe P2_3 = [table_ptr1], 16
2677 fmerge.s Pos_r = f1, r
2678 tbit.nz p8,p9 = exp_r, 0
2682 // Is B = 2** -2 or B= 2** -1? If 2**-1, then
2683 // we want an offset of 512 for table addressing.
2685 add table_ptr2 = 1296, table_base // Point to tanl_table_cm2
2686 (p9) shladd table_offset = lookup, 4, table_offset
2687 (p8) shladd table_offset = lookup, 4, r0
2692 add table_ptr1 = table_ptr1, table_offset // Point to T_hi
2693 add table_ptr2 = table_ptr2, table_offset // Point to C_hi
2694 add table_ptr3 = 2128, table_base // Point to tanl_table_scim2
2699 ldfd T_hi = [table_ptr1], 8 // Load T_hi
2701 ldfd C_hi = [table_ptr2], 8 // Load C_hi
2702 add table_ptr3 = table_ptr3, table_offset // Point to SC_inv
2709 // Convert B so it has the same exponent as Pos_r before subtracting
2711 ldfs T_lo = [table_ptr1] // Load T_lo
2712 (p9) fnma.s1 x = B, FR_2tom64, Pos_r
2717 (p8) fnma.s1 x = B, TWO_TO_NEG65, Pos_r
2723 ldfs C_lo = [table_ptr2] // Load C_lo
2730 ldfe SC_inv = [table_ptr3] // Load SC_inv
2731 fmerge.s sgn_r = r, f1
2732 tbit.z p11, p12 = N_fix_gr, 0 // p11 if N even, p12 if odd
2739 // N even: Tx = T_hi * x
2741 // N even: Tx1 = Tx + 1
2742 // N odd: Cx1 = 1 - Cx
2752 (p11) fmpy.s1 Tx = T_hi, x
2758 // N odd: Cx = C_hi * x
2762 (p12) fmpy.s1 Cx = C_hi, x
2767 // N even and odd: P = P2_3 + P2_2 * xsq
2771 fma.s1 P = P2_3, xsq, P2_2
2776 (p11) fadd.s1 Tx1 = Tx, f1
2782 // N even: D = C_hi - tanx
2783 // N odd: D = T_hi + tanx
2785 (p11) fmpy.s1 CORR = SC_inv, T_hi
2790 fmpy.s1 Sx = SC_inv, x
2795 (p12) fmpy.s1 CORR = SC_inv, C_hi
2800 (p12) fsub.s1 V_hi = f1, Cx
2805 fma.s1 P = P, xsq, P2_1
2811 // N even and odd: P = P2_1 + P * xsq
2813 (p11) fma.s1 V_hi = Tx, Tx1, f1
2819 // N even: Result = sgn_r * tail + T_hi (user rounding mode for C1)
2820 // N odd: Result = sgn_r * tail + C_hi (user rounding mode for C1)
2822 fmpy.s0 fp_tmp = fp_tmp, fp_tmp // Dummy mult to set inexact
2827 fmpy.s1 CORR = CORR, c
2832 (p12) fnma.s1 V_hi = Cx,V_hi,f1
2838 // N even: V_hi = Tx * Tx1 + 1
2839 // N odd: Cx1 = 1 - Cx * Cx1
2847 // N even and odd: P = P * xsq
2849 (p11) fmpy.s1 V_hi = V_hi, T_hi
2855 // N even and odd: tail = P * tail + V_lo
2857 (p11) fmpy.s1 T_hi = sgn_r, T_hi
2862 fmpy.s1 CORR = CORR, sgn_r
2867 (p12) fmpy.s1 V_hi = V_hi,C_hi
2873 // N even: V_hi = T_hi * V_hi
2874 // N odd: V_hi = C_hi * V_hi
2876 fma.s1 tanx = P, x, x
2881 (p12) fnmpy.s1 C_hi = sgn_r, C_hi
2887 // N even: V_lo = 1 - V_hi + C_hi
2888 // N odd: V_lo = 1 - V_hi + T_hi
2890 (p11) fadd.s1 CORR = CORR, T_lo
2895 (p12) fsub.s1 CORR = CORR, C_lo
2901 // N even and odd: tanx = x + x * P
2902 // N even and odd: Sx = SC_inv * x
2904 (p11) fsub.s1 D = C_hi, tanx
2909 (p12) fadd.s1 D = T_hi, tanx
2915 // N odd: CORR = SC_inv * C_hi
2916 // N even: CORR = SC_inv * T_hi
2918 fnma.s1 D = V_hi, D, f1
2924 // N even and odd: D = 1 - V_hi * D
2925 // N even and odd: CORR = CORR * c
2927 fma.s1 V_hi = V_hi, D, V_hi
2933 // N even and odd: V_hi = V_hi + V_hi * D
2934 // N even and odd: CORR = sgn_r * CORR
2936 (p11) fnma.s1 V_lo = V_hi, C_hi, f1
2941 (p12) fnma.s1 V_lo = V_hi, T_hi, f1
2947 // N even: CORR = COOR + T_lo
2948 // N odd: CORR = CORR - C_lo
2950 (p11) fma.s1 V_lo = tanx, V_hi, V_lo
2951 tbit.nz p15, p0 = cot_flag, 0 // p15=1 if we compute cotl
2955 (p12) fnma.s1 V_lo = tanx, V_hi, V_lo
2961 (p15) fms.s1 T_hi = f0, f0, T_hi // to correct result's sign for cotl
2966 (p15) fms.s1 C_hi = f0, f0, C_hi // to correct result's sign for cotl
2972 (p15) fms.s1 sgn_r = f0, f0, sgn_r // to correct result's sign for cotl
2979 // N even: V_lo = V_lo + V_hi * tanx
2980 // N odd: V_lo = V_lo - V_hi * tanx
2982 (p11) fnma.s1 V_lo = C_lo, V_hi, V_lo
2987 (p12) fnma.s1 V_lo = T_lo, V_hi, V_lo
2993 // N even: V_lo = V_lo - V_hi * C_lo
2994 // N odd: V_lo = V_lo - V_hi * T_lo
2996 fmpy.s1 V_lo = V_hi, V_lo
3002 // N even and odd: V_lo = V_lo * V_hi
3004 fadd.s1 tail = V_hi, V_lo
3010 // N even and odd: tail = V_hi + V_lo
3012 fma.s1 tail = tail, P, V_lo
3018 // N even: T_hi = sgn_r * T_hi
3019 // N odd : C_hi = -sgn_r * C_hi
3021 fma.s1 tail = tail, Sx, CORR
3027 // N even and odd: tail = Sx * tail + CORR
3029 fma.s1 tail = V_hi, Sx, tail
3035 // N even an odd: tail = Sx * V_hi + tail
3037 (p11) fma.s0 Result = sgn_r, tail, T_hi
3042 (p12) fma.s0 Result = sgn_r, tail, C_hi
3043 br.ret.sptk b0 ;; // Exit for 1/4 <= |r| < pi/4
3047 // Here if x denormal
3049 getf.exp GR_signexp_x = Norm_Arg // Get sign and exponent of x
3051 br.cond.sptk TANL_COMMON // Return to common code
3059 // Code for NaNs, Unsupporteds, Infs, or +/- zero ?
3060 // Invalid raised for Infs and SNaNs.
3065 fmerge.s f10 = f8, f8 // Save input for error call
3066 tbit.nz p6, p7 = cot_flag, 0 // p6=1 if we compute cotl
3072 (p6) fclass.m p6, p7 = f8, 0x7 // Test for zero (cotl only)
3077 .pred.rel "mutex", p6, p7
3079 (p6) mov GR_Parameter_Tag = 225 // (cotl)
3080 (p6) frcpa.s0 f8, p0 = f1, f8 // cotl(+-0) = +-Inf
3085 (p7) fmpy.s0 f8 = f8, f0
3090 GLOBAL_IEEE754_END(tanl)
3092 LOCAL_LIBM_ENTRY(__libm_error_region)
3097 add GR_Parameter_Y=-32,sp // Parameter 2 value
3099 .save ar.pfs,GR_SAVE_PFS
3100 mov GR_SAVE_PFS=ar.pfs // Save ar.pfs
3104 add sp=-64,sp // Create new stack
3106 mov GR_SAVE_GP=gp // Save gp
3111 stfe [GR_Parameter_Y] = f1,16 // STORE Parameter 2 on stack
3112 add GR_Parameter_X = 16,sp // Parameter 1 address
3113 .save b0, GR_SAVE_B0
3114 mov GR_SAVE_B0=b0 // Save b0
3120 stfe [GR_Parameter_X] = f10 // STORE Parameter 1 on stack
3121 add GR_Parameter_RESULT = 0,GR_Parameter_Y // Parameter 3 address
3125 stfe [GR_Parameter_Y] = f8 // STORE Parameter 3 on stack
3126 add GR_Parameter_Y = -16,GR_Parameter_Y
3127 br.call.sptk b0=__libm_error_support# // Call error handling function
3132 add GR_Parameter_RESULT = 48,sp
3137 ldfe f8 = [GR_Parameter_RESULT] // Get return result off stack
3139 add sp = 64,sp // Restore stack pointer
3140 mov b0 = GR_SAVE_B0 // Restore return address
3143 mov gp = GR_SAVE_GP // Restore gp
3144 mov ar.pfs = GR_SAVE_PFS // Restore ar.pfs
3145 br.ret.sptk b0 // Return
3148 LOCAL_LIBM_END(__libm_error_region)
3150 .type __libm_error_support#,@function
3151 .global __libm_error_support#
3154 // *******************************************************************
3155 // *******************************************************************
3156 // *******************************************************************
3158 // Special Code to handle very large argument case.
3159 // Call int __libm_pi_by_2_reduce(x,r,c) for |arguments| >= 2**63
3160 // The interface is custom:
3162 // (Arg or x) is in f8
3167 // We know also that __libm_pi_by_2_reduce preserves f10-15, f71-127. We
3168 // use this to eliminate save/restore of key fp registers in this calling
3171 // *******************************************************************
3172 // *******************************************************************
3173 // *******************************************************************
3175 LOCAL_LIBM_ENTRY(__libm_callout)
3179 add table_ptr2 = 144, table_base // Point to 2^-2
3181 .save ar.pfs,GR_SAVE_PFS
3182 mov GR_SAVE_PFS=ar.pfs // Save ar.pfs
3188 ldfps TWO_TO_NEG2, NEGTWO_TO_NEG2 = [table_ptr2]
3189 setf.sig B_mask1 = bmask1 // Form mask to get 5 msb of r
3190 .save b0, GR_SAVE_B0
3191 mov GR_SAVE_B0=b0 // Save b0
3196 // Call argument reduction with x in f8
3197 // Returns with N in r8, r in f8, c in f9
3198 // Assumes f71-127 are preserved across the call
3201 setf.sig B_mask2 = bmask2 // Form mask to form B from r
3202 mov GR_SAVE_GP=gp // Save gp
3203 br.call.sptk b0=__libm_pi_by_2_reduce#
3211 getf.sig sig_r = r // Extract significand of r
3212 fcmp.lt.s1 p6, p0 = r, TWO_TO_NEG2
3213 mov gp = GR_SAVE_GP // Restore gp
3218 getf.exp exp_r = r // Extract signexp of r
3220 mov b0 = GR_SAVE_B0 // Restore return address
3229 (p6) fcmp.gt.unc.s1 p6, p0 = r, NEGTWO_TO_NEG2
3230 mov ar.pfs = GR_SAVE_PFS // Restore pfs
3236 (p6) br.cond.spnt TANL_SMALL_R // Branch if |r| < 1/4
3237 br.cond.sptk TANL_NORMAL_R // Branch if 1/4 <= |r| < pi/4
3241 LOCAL_LIBM_END(__libm_callout)
3243 .type __libm_pi_by_2_reduce#,@function
3244 .global __libm_pi_by_2_reduce#