4 // Copyright (c) 2000 - 2004, Intel Corporation
5 // All rights reserved.
7 // Contributed 2000 by the Intel Numerics Group, Intel Corporation
9 // Redistribution and use in source and binary forms, with or without
10 // modification, are permitted provided that the following conditions are
13 // * Redistributions of source code must retain the above copyright
14 // notice, this list of conditions and the following disclaimer.
16 // * Redistributions in binary form must reproduce the above copyright
17 // notice, this list of conditions and the following disclaimer in the
18 // documentation and/or other materials provided with the distribution.
20 // * The name of Intel Corporation may not be used to endorse or promote
21 // products derived from this software without specific prior written
24 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
25 // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
26 // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
27 // A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL INTEL OR ITS
28 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
29 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
30 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
31 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY
32 // OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY OR TORT (INCLUDING
33 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
34 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 // Intel Corporation is the author of this code, and requests that all
37 // problem reports or change requests be submitted to it directly at
38 // http://www.intel.com/software/products/opensource/libraries/num.htm.
40 //*********************************************************************
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
53 // 10/26/04 Avoided using r14-31 as scratch so not clobbered by dynamic loader
55 //*********************************************************************
57 // Functions: tanl(x) = tangent(x), for double-extended precision x values
58 // cotl(x) = cotangent(x), for double-extended precision x values
60 //*********************************************************************
64 // Floating-Point Registers: f8 (Input and Return Value)
68 // General Purpose Registers:
71 // Predicate Registers: p6-p15
73 //*********************************************************************
75 // IEEE Special Conditions for tanl:
77 // Denormal fault raised on denormal inputs
78 // Overflow exceptions do not occur
79 // Underflow exceptions raised when appropriate for tan
80 // (No specialized error handling for this routine)
81 // Inexact raised when appropriate by algorithm
88 //*********************************************************************
90 // IEEE Special Conditions for cotl:
92 // Denormal fault raised on denormal inputs
93 // Overflow exceptions occur at zero and near zero
94 // Underflow exceptions do not occur
95 // Inexact raised when appropriate by algorithm
100 // cotl(+/-0) = +/-Inf and error handling is called
102 //*********************************************************************
104 // Below are mathematical and algorithmic descriptions for tanl.
105 // For cotl we use next identity cot(x) = -tan(x + Pi/2).
106 // So, to compute cot(x) we just need to increment N (N = N + 1)
107 // and invert sign of the computed result.
109 //*********************************************************************
111 // Mathematical Description
113 // We consider the computation of FPTANL of Arg. Now, given
115 // Arg = N pi/2 + alpha, |alpha| <= pi/4,
117 // basic mathematical relationship shows that
119 // tan( Arg ) = tan( alpha ) if N is even;
120 // = -cot( alpha ) otherwise.
122 // The value of alpha is obtained by argument reduction and
123 // represented by two working precision numbers r and c where
125 // alpha = r + c accurately.
127 // The reduction method is described in a previous write up.
128 // The argument reduction scheme identifies 4 cases. For Cases 2
129 // and 4, because |alpha| is small, tan(r+c) and -cot(r+c) can be
130 // computed very easily by 2 or 3 terms of the Taylor series
131 // expansion as follows:
136 // tan(r + c) = r + c + r^3/3 ...accurately
137 // -cot(r + c) = -1/(r+c) + r/3 ...accurately
142 // tan(r + c) = r + c + r^3/3 + 2r^5/15 ...accurately
143 // -cot(r + c) = -1/(r+c) + r/3 + r^3/45 ...accurately
146 // The only cases left are Cases 1 and 3 of the argument reduction
147 // procedure. These two cases will be merged since after the
148 // argument is reduced in either cases, we have the reduced argument
149 // represented as r + c and that the magnitude |r + c| is not small
150 // enough to allow the usage of a very short approximation.
152 // The greatest challenge of this task is that the second terms of
153 // the Taylor series for tan(r) and -cot(r)
155 // r + r^3/3 + 2 r^5/15 + ...
159 // -1/r + r/3 + r^3/45 + ...
161 // are not very small when |r| is close to pi/4 and the rounding
162 // errors will be a concern if simple polynomial accumulation is
163 // used. When |r| < 2^(-2), however, the second terms will be small
164 // enough (5 bits or so of right shift) that a normal Horner
165 // recurrence suffices. Hence there are two cases that we consider
166 // in the accurate computation of tan(r) and cot(r), |r| <= pi/4.
168 // Case small_r: |r| < 2^(-2)
169 // --------------------------
171 // Since Arg = N pi/4 + r + c accurately, we have
173 // tan(Arg) = tan(r+c) for N even,
174 // = -cot(r+c) otherwise.
176 // Here for this case, both tan(r) and -cot(r) can be approximated
177 // by simple polynomials:
179 // tan(r) = r + P1_1 r^3 + P1_2 r^5 + ... + P1_9 r^19
180 // -cot(r) = -1/r + Q1_1 r + Q1_2 r^3 + ... + Q1_7 r^13
182 // accurately. Since |r| is relatively small, tan(r+c) and
183 // -cot(r+c) can be accurately approximated by replacing r with
184 // r+c only in the first two terms of the corresponding polynomials.
186 // Note that P1_1 (and Q1_1 for that matter) approximates 1/3 to
187 // almost 64 sig. bits, thus
189 // P1_1 (r+c)^3 = P1_1 r^3 + c * r^2 accurately.
193 // tan(r+c) = r + P1_1 r^3 + P1_2 r^5 + ... + P1_9 r^19
196 // -cot(r+c) = -1/(r+c) + Q1_1 r + Q1_2 r^3 + ... + Q1_7 r^13
200 // Case normal_r: 2^(-2) <= |r| <= pi/4
201 // ------------------------------------
203 // This case is more likely than the previous one if one considers
204 // r to be uniformly distributed in [-pi/4 pi/4].
206 // The required calculation is either
208 // tan(r + c) = tan(r) + correction, or
209 // -cot(r + c) = -cot(r) + correction.
213 // tan(r + c) = tan(r) + c tan'(r) + O(c^2)
214 // = tan(r) + c sec^2(r) + O(c^2)
215 // = tan(r) + c SEC_sq ...accurately
216 // as long as SEC_sq approximates sec^2(r)
217 // to, say, 5 bits or so.
221 // -cot(r + c) = -cot(r) - c cot'(r) + O(c^2)
222 // = -cot(r) + c csc^2(r) + O(c^2)
223 // = -cot(r) + c CSC_sq ...accurately
224 // as long as CSC_sq approximates csc^2(r)
225 // to, say, 5 bits or so.
227 // We therefore concentrate on accurately calculating tan(r) and
228 // cot(r) for a working-precision number r, |r| <= pi/4 to within
231 // We will employ a table-driven approach. Let
233 // r = sgn_r * 2^k * 1.b_1 b_2 ... b_5 ... b_63
234 // = sgn_r * ( B + x )
238 // B = 2^k * 1.b_1 b_2 ... b_5 1
243 // tan( B + x ) = ------------------------
247 // | tan(B) + tan(x) |
249 // = tan(B) + | ------------------------ - tan(B) |
250 // | 1 - tan(B)*tan(x) |
254 // = tan(B) + ------------------------
257 // (1/[sin(B)*cos(B)]) * tan(x)
258 // = tan(B) + --------------------------------
262 // Clearly, the values of tan(B), cot(B) and 1/(sin(B)*cos(B)) are
263 // calculated beforehand and stored in a table. Since
265 // |x| <= 2^k * 2^(-6) <= 2^(-7) (because k = -1, -2)
267 // a very short polynomial will be sufficient to approximate tan(x)
268 // accurately. The details involved in computing the last expression
269 // will be given in the next section on algorithm description.
272 // Now, we turn to the case where cot( B + x ) is needed.
276 // cot( B + x ) = ------------------------
280 // | 1 - tan(B)*tan(x) |
282 // = cot(B) + | ----------------------- - cot(B) |
283 // | tan(B) + tan(x) |
286 // [tan(B) + cot(B)] * tan(x)
287 // = cot(B) - ----------------------------
290 // (1/[sin(B)*cos(B)]) * tan(x)
291 // = cot(B) - --------------------------------
295 // Note that the values of tan(B), cot(B) and 1/(sin(B)*cos(B)) that
296 // are needed are the same set of values needed in the previous
299 // Finally, we can put all the ingredients together as follows:
301 // Arg = N * pi/2 + r + c ...accurately
303 // tan(Arg) = tan(r) + correction if N is even;
304 // = -cot(r) + correction otherwise.
306 // For Cases 2 and 4,
309 // tan(Arg) = tan(r + c) = r + c + r^3/3 N even
310 // = -cot(r + c) = -1/(r+c) + r/3 N odd
312 // tan(Arg) = tan(r + c) = r + c + r^3/3 + 2r^5/15 N even
313 // = -cot(r + c) = -1/(r+c) + r/3 + r^3/45 N odd
316 // For Cases 1 and 3,
318 // Case small_r: |r| < 2^(-2)
320 // tan(Arg) = r + P1_1 r^3 + P1_2 r^5 + ... + P1_9 r^19
321 // + c*(1 + r^2) N even
323 // = -1/(r+c) + Q1_1 r + Q1_2 r^3 + ... + Q1_7 r^13
326 // Case normal_r: 2^(-2) <= |r| <= pi/4
328 // tan(Arg) = tan(r) + c * sec^2(r) N even
329 // = -cot(r) + c * csc^2(r) otherwise
333 // tan(Arg) = tan(r) + c*sec^2(r)
334 // = tan( sgn_r * (B+x) ) + c * sec^2(|r|)
335 // = sgn_r * ( tan(B+x) + sgn_r*c*sec^2(|r|) )
336 // = sgn_r * ( tan(B+x) + sgn_r*c*sec^2(B) )
338 // since B approximates |r| to 2^(-6) in relative accuracy.
340 // / (1/[sin(B)*cos(B)]) * tan(x)
341 // tan(Arg) = sgn_r * | tan(B) + --------------------------------
349 // CORR = sgn_r*c*tan(B)*SC_inv(B); SC_inv(B) = 1/(sin(B)*cos(B)).
353 // tan(Arg) = -cot(r) + c*csc^2(r)
354 // = -cot( sgn_r * (B+x) ) + c * csc^2(|r|)
355 // = sgn_r * ( -cot(B+x) + sgn_r*c*csc^2(|r|) )
356 // = sgn_r * ( -cot(B+x) + sgn_r*c*csc^2(B) )
358 // since B approximates |r| to 2^(-6) in relative accuracy.
360 // / (1/[sin(B)*cos(B)]) * tan(x)
361 // tan(Arg) = sgn_r * | -cot(B) + --------------------------------
369 // CORR = sgn_r*c*cot(B)*SC_inv(B); SC_inv(B) = 1/(sin(B)*cos(B)).
372 // The actual algorithm prescribes how all the mathematical formulas
376 // 2. Algorithmic Description
377 // ==========================
379 // 2.1 Computation for Cases 2 and 4.
380 // ----------------------------------
382 // For Case 2, we use two-term polynomials.
387 // Poly := c + r * rsq * P1_1
388 // Result := r + Poly ...in user-defined rounding
391 // S_hi := -frcpa(r) ...8 bits
392 // S_hi := S_hi + S_hi*(1 + S_hi*r) ...16 bits
393 // S_hi := S_hi + S_hi*(1 + S_hi*r) ...32 bits
394 // S_hi := S_hi + S_hi*(1 + S_hi*r) ...64 bits
395 // S_lo := S_hi*( (1 + S_hi*r) + S_hi*c )
396 // ...S_hi + S_lo is -1/(r+c) to extra precision
397 // S_lo := S_lo + Q1_1*r
399 // Result := S_hi + S_lo ...in user-defined rounding
401 // For Case 4, we use three-term polynomials
406 // Poly := c + r * rsq * (P1_1 + rsq * P1_2)
407 // Result := r + Poly ...in user-defined rounding
410 // S_hi := -frcpa(r) ...8 bits
411 // S_hi := S_hi + S_hi*(1 + S_hi*r) ...16 bits
412 // S_hi := S_hi + S_hi*(1 + S_hi*r) ...32 bits
413 // S_hi := S_hi + S_hi*(1 + S_hi*r) ...64 bits
414 // S_lo := S_hi*( (1 + S_hi*r) + S_hi*c )
415 // ...S_hi + S_lo is -1/(r+c) to extra precision
417 // P := Q1_1 + rsq*Q1_2
418 // S_lo := S_lo + r*P
420 // Result := S_hi + S_lo ...in user-defined rounding
423 // Note that the coefficients P1_1, P1_2, Q1_1, and Q1_2 are
424 // the same as those used in the small_r case of Cases 1 and 3
428 // 2.2 Computation for Cases 1 and 3.
429 // ----------------------------------
430 // This is further divided into the case of small_r,
431 // where |r| < 2^(-2), and the case of normal_r, where |r| lies between
434 // Algorithm for the case of small_r
435 // ---------------------------------
439 // Poly1 := rsq*(P1_1 + rsq*(P1_2 + rsq*P1_3))
440 // r_to_the_8 := rsq * rsq
441 // r_to_the_8 := r_to_the_8 * r_to_the_8
442 // Poly2 := P1_4 + rsq*(P1_5 + rsq*(P1_6 + ... rsq*P1_9))
443 // CORR := c * ( 1 + rsq )
444 // Poly := Poly1 + r_to_the_8*Poly2
445 // Poly := r*Poly + CORR
446 // Result := r + Poly ...in user-defined rounding
447 // ...note that Poly1 and r_to_the_8 can be computed in parallel
448 // ...with Poly2 (Poly1 is intentionally set to be much
449 // ...shorter than Poly2 so that r_to_the_8 and CORR can be hidden)
452 // S_hi := -frcpa(r) ...8 bits
453 // S_hi := S_hi + S_hi*(1 + S_hi*r) ...16 bits
454 // S_hi := S_hi + S_hi*(1 + S_hi*r) ...32 bits
455 // S_hi := S_hi + S_hi*(1 + S_hi*r) ...64 bits
456 // S_lo := S_hi*( (1 + S_hi*r) + S_hi*c )
457 // ...S_hi + S_lo is -1/(r+c) to extra precision
458 // S_lo := S_lo + Q1_1*c
460 // ...S_hi and S_lo are computed in parallel with
463 // P := Q1_1 + rsq*(Q1_2 + rsq*(Q1_3 + ... + rsq*Q1_7))
465 // Poly := r*P + S_lo
466 // Result := S_hi + Poly ...in user-defined rounding
469 // Algorithm for the case of normal_r
470 // ----------------------------------
472 // Here, we first consider the computation of tan( r + c ). As
473 // presented in the previous section,
475 // tan( r + c ) = tan(r) + c * sec^2(r)
476 // = sgn_r * [ tan(B+x) + CORR ]
477 // CORR = sgn_r * c * tan(B) * 1/[sin(B)*cos(B)]
479 // because sec^2(r) = sec^(|r|), and B approximate |r| to 6.5 bits.
482 // / (1/[sin(B)*cos(B)]) * tan(x)
483 // sgn_r * | tan(B) + -------------------------------- +
490 // The values of tan(B), cot(B) and 1/(sin(B)*cos(B)) are
491 // calculated beforehand and stored in a table. Specifically,
492 // the table values are
494 // tan(B) as T_hi + T_lo;
495 // cot(B) as C_hi + C_lo;
496 // 1/[sin(B)*cos(B)] as SC_inv
498 // T_hi, C_hi are in double-precision memory format;
499 // T_lo, C_lo are in single-precision memory format;
500 // SC_inv is in extended-precision memory format.
502 // The value of tan(x) will be approximated by a short polynomial of
505 // tan(x) as x + x * P, where
506 // P = x^2 * (P2_1 + x^2 * (P2_2 + x^2 * P2_3))
508 // Because |x| <= 2^(-7), cot(B) - x approximates cot(B) - tan(x)
509 // to a relative accuracy better than 2^(-20). Thus, a good
510 // initial guess of 1/( cot(B) - tan(x) ) to initiate the iterative
513 // 1/(cot(B) - tan(x)) is approximately
515 // tan(B)/(1 - x*tan(B)) is approximately
516 // T_hi / ( 1 - T_hi * x ) is approximately
518 // T_hi * [ 1 + (Thi * x) + (T_hi * x)^2 ]
520 // The calculation of tan(r+c) therefore proceed as follows:
525 // V_hi := T_hi*(1 + Tx*(1 + Tx))
526 // P := xsq * (P1_1 + xsq*(P1_2 + xsq*P1_3))
527 // ...V_hi serves as an initial guess of 1/(cot(B) - tan(x))
528 // ...good to about 20 bits of accuracy
532 // ...D is a double precision denominator: cot(B) - tan(x)
534 // V_hi := V_hi + V_hi*(1 - V_hi*D)
535 // ....V_hi approximates 1/(cot(B)-tan(x)) to 40 bits
537 // V_lo := V_hi * ( [ (1 - V_hi*C_hi) + V_hi*tanx ]
538 // - V_hi*C_lo ) ...observe all order
539 // ...V_hi + V_lo approximates 1/(cot(B) - tan(x))
540 // ...to extra accuracy
542 // ... SC_inv(B) * (x + x*P)
543 // ... tan(B) + ------------------------- + CORR
544 // ... cot(B) - (x + x*P)
546 // ... = tan(B) + SC_inv(B)*(x + x*P)*(V_hi + V_lo) + CORR
550 // CORR := sgn_r * c * SC_inv * T_hi
552 // ...put the ingredients together to compute
553 // ... SC_inv(B) * (x + x*P)
554 // ... tan(B) + ------------------------- + CORR
555 // ... cot(B) - (x + x*P)
557 // ... = tan(B) + SC_inv(B)*(x + x*P)*(V_hi + V_lo) + CORR
559 // ... = T_hi + T_lo + CORR +
560 // ... Sx * V_hi + Sx * V_lo + Sx * P *(V_hi + V_lo)
562 // CORR := CORR + T_lo
563 // tail := V_lo + P*(V_hi + V_lo)
564 // tail := Sx * tail + CORR
565 // tail := Sx * V_hi + tail
566 // T_hi := sgn_r * T_hi
568 // ...T_hi + sgn_r*tail now approximate
569 // ...sgn_r*(tan(B+x) + CORR) accurately
571 // Result := T_hi + sgn_r*tail ...in user-defined
572 // ...rounding control
573 // ...It is crucial that independent paths be fully
574 // ...exploited for performance's sake.
577 // Next, we consider the computation of -cot( r + c ). As
578 // presented in the previous section,
580 // -cot( r + c ) = -cot(r) + c * csc^2(r)
581 // = sgn_r * [ -cot(B+x) + CORR ]
582 // CORR = sgn_r * c * cot(B) * 1/[sin(B)*cos(B)]
584 // because csc^2(r) = csc^(|r|), and B approximate |r| to 6.5 bits.
587 // / (1/[sin(B)*cos(B)]) * tan(x)
588 // sgn_r * | -cot(B) + -------------------------------- +
595 // The values of tan(B), cot(B) and 1/(sin(B)*cos(B)) are
596 // calculated beforehand and stored in a table. Specifically,
597 // the table values are
599 // tan(B) as T_hi + T_lo;
600 // cot(B) as C_hi + C_lo;
601 // 1/[sin(B)*cos(B)] as SC_inv
603 // T_hi, C_hi are in double-precision memory format;
604 // T_lo, C_lo are in single-precision memory format;
605 // SC_inv is in extended-precision memory format.
607 // The value of tan(x) will be approximated by a short polynomial of
610 // tan(x) as x + x * P, where
611 // P = x^2 * (P2_1 + x^2 * (P2_2 + x^2 * P2_3))
613 // Because |x| <= 2^(-7), tan(B) + x approximates tan(B) + tan(x)
614 // to a relative accuracy better than 2^(-18). Thus, a good
615 // initial guess of 1/( tan(B) + tan(x) ) to initiate the iterative
618 // 1/(tan(B) + tan(x)) is approximately
620 // cot(B)/(1 + x*cot(B)) is approximately
621 // C_hi / ( 1 + C_hi * x ) is approximately
623 // C_hi * [ 1 - (C_hi * x) + (C_hi * x)^2 ]
625 // The calculation of -cot(r+c) therefore proceed as follows:
630 // V_hi := C_hi*(1 - Cx*(1 - Cx))
631 // P := xsq * (P1_1 + xsq*(P1_2 + xsq*P1_3))
632 // ...V_hi serves as an initial guess of 1/(tan(B) + tan(x))
633 // ...good to about 18 bits of accuracy
637 // ...D is a double precision denominator: tan(B) + tan(x)
639 // V_hi := V_hi + V_hi*(1 - V_hi*D)
640 // ....V_hi approximates 1/(tan(B)+tan(x)) to 40 bits
642 // V_lo := V_hi * ( [ (1 - V_hi*T_hi) - V_hi*tanx ]
643 // - V_hi*T_lo ) ...observe all order
644 // ...V_hi + V_lo approximates 1/(tan(B) + tan(x))
645 // ...to extra accuracy
647 // ... SC_inv(B) * (x + x*P)
648 // ... -cot(B) + ------------------------- + CORR
649 // ... tan(B) + (x + x*P)
651 // ... =-cot(B) + SC_inv(B)*(x + x*P)*(V_hi + V_lo) + CORR
655 // CORR := sgn_r * c * SC_inv * C_hi
657 // ...put the ingredients together to compute
658 // ... SC_inv(B) * (x + x*P)
659 // ... -cot(B) + ------------------------- + CORR
660 // ... tan(B) + (x + x*P)
662 // ... =-cot(B) + SC_inv(B)*(x + x*P)*(V_hi + V_lo) + CORR
664 // ... =-C_hi - C_lo + CORR +
665 // ... Sx * V_hi + Sx * V_lo + Sx * P *(V_hi + V_lo)
667 // CORR := CORR - C_lo
668 // tail := V_lo + P*(V_hi + V_lo)
669 // tail := Sx * tail + CORR
670 // tail := Sx * V_hi + tail
671 // C_hi := -sgn_r * C_hi
673 // ...C_hi + sgn_r*tail now approximates
674 // ...sgn_r*(-cot(B+x) + CORR) accurately
676 // Result := C_hi + sgn_r*tail in user-defined rounding control
677 // ...It is crucial that independent paths be fully
678 // ...exploited for performance's sake.
680 // 3. Implementation Notes
681 // =======================
683 // Table entries T_hi, T_lo; C_hi, C_lo; SC_inv
685 // Recall that 2^(-2) <= |r| <= pi/4;
687 // r = sgn_r * 2^k * 1.b_1 b_2 ... b_63
691 // B = 2^k * 1.b_1 b_2 b_3 b_4 b_5 1
693 // Thus, for k = -2, possible values of B are
695 // B = 2^(-2) * ( 1 + index/32 + 1/64 ),
696 // index ranges from 0 to 31
698 // For k = -1, however, since |r| <= pi/4 = 0.78...
699 // possible values of B are
701 // B = 2^(-1) * ( 1 + index/32 + 1/64 )
702 // index ranges from 0 to 19.
709 LOCAL_OBJECT_START(TANL_BASE_CONSTANTS)
712 data8 0xA2F9836E4E44152A, 0x00003FFE // two_by_pi
713 data8 0xC84D32B0CE81B9F1, 0x00004016 // P_0
714 data8 0xC90FDAA22168C235, 0x00003FFF // P_1
715 data8 0xECE675D1FC8F8CBB, 0x0000BFBD // P_2
716 data8 0xB7ED8FBBACC19C60, 0x0000BF7C // P_3
717 LOCAL_OBJECT_END(TANL_BASE_CONSTANTS)
719 LOCAL_OBJECT_START(tanl_table_2)
720 data8 0xC90FDAA22168C234, 0x00003FFE // PI_BY_4
721 data8 0xA397E5046EC6B45A, 0x00003FE7 // Inv_P_0
722 data8 0x8D848E89DBD171A1, 0x0000BFBF // d_1
723 data8 0xD5394C3618A66F8E, 0x0000BF7C // d_2
724 data4 0x3E800000 // two**-2
725 data4 0xBE800000 // -two**-2
726 data4 0x00000000 // pad
727 data4 0x00000000 // pad
728 LOCAL_OBJECT_END(tanl_table_2)
730 LOCAL_OBJECT_START(tanl_table_p1)
731 data8 0xAAAAAAAAAAAAAABD, 0x00003FFD // P1_1
732 data8 0x8888888888882E6A, 0x00003FFC // P1_2
733 data8 0xDD0DD0DD0F0177B6, 0x00003FFA // P1_3
734 data8 0xB327A440646B8C6D, 0x00003FF9 // P1_4
735 data8 0x91371B251D5F7D20, 0x00003FF8 // P1_5
736 data8 0xEB69A5F161C67914, 0x00003FF6 // P1_6
737 data8 0xBEDD37BE019318D2, 0x00003FF5 // P1_7
738 data8 0x9979B1463C794015, 0x00003FF4 // P1_8
739 data8 0x8EBD21A38C6EB58A, 0x00003FF3 // P1_9
740 LOCAL_OBJECT_END(tanl_table_p1)
742 LOCAL_OBJECT_START(tanl_table_q1)
743 data8 0xAAAAAAAAAAAAAAB4, 0x00003FFD // Q1_1
744 data8 0xB60B60B60B5FC93E, 0x00003FF9 // Q1_2
745 data8 0x8AB355E00C9BBFBF, 0x00003FF6 // Q1_3
746 data8 0xDDEBBC89CBEE3D4C, 0x00003FF2 // Q1_4
747 data8 0xB3548A685F80BBB6, 0x00003FEF // Q1_5
748 data8 0x913625604CED5BF1, 0x00003FEC // Q1_6
749 data8 0xF189D95A8EE92A83, 0x00003FE8 // Q1_7
750 LOCAL_OBJECT_END(tanl_table_q1)
752 LOCAL_OBJECT_START(tanl_table_p2)
753 data8 0xAAAAAAAAAAAB362F, 0x00003FFD // P2_1
754 data8 0x88888886E97A6097, 0x00003FFC // P2_2
755 data8 0xDD108EE025E716A1, 0x00003FFA // P2_3
756 LOCAL_OBJECT_END(tanl_table_p2)
758 LOCAL_OBJECT_START(tanl_table_tm2)
760 // Entries T_hi double-precision memory format
761 // Index = 0,1,...,31 B = 2^(-2)*(1+Index/32+1/64)
762 // Entries T_lo single-precision memory format
763 // Index = 0,1,...,31 B = 2^(-2)*(1+Index/32+1/64)
765 data8 0x3FD09BC362400794
766 data4 0x23A05C32, 0x00000000
767 data8 0x3FD124A9DFFBC074
768 data4 0x240078B2, 0x00000000
769 data8 0x3FD1AE235BD4920F
770 data4 0x23826B8E, 0x00000000
771 data8 0x3FD2383515E2701D
772 data4 0x22D31154, 0x00000000
773 data8 0x3FD2C2E463739C2D
774 data4 0x2265C9E2, 0x00000000
775 data8 0x3FD34E36AFEEA48B
776 data4 0x245C05EB, 0x00000000
777 data8 0x3FD3DA317DBB35D1
778 data4 0x24749F2D, 0x00000000
779 data8 0x3FD466DA67321619
780 data4 0x2462CECE, 0x00000000
781 data8 0x3FD4F4371F94A4D5
782 data4 0x246D0DF1, 0x00000000
783 data8 0x3FD5824D740C3E6D
784 data4 0x240A85B5, 0x00000000
785 data8 0x3FD611234CB1E73D
786 data4 0x23F96E33, 0x00000000
787 data8 0x3FD6A0BEAD9EA64B
788 data4 0x247C5393, 0x00000000
789 data8 0x3FD73125B804FD01
790 data4 0x241F3B29, 0x00000000
791 data8 0x3FD7C25EAB53EE83
792 data4 0x2479989B, 0x00000000
793 data8 0x3FD8546FE6640EED
794 data4 0x23B343BC, 0x00000000
795 data8 0x3FD8E75FE8AF1892
796 data4 0x241454D1, 0x00000000
797 data8 0x3FD97B3553928BDA
798 data4 0x238613D9, 0x00000000
799 data8 0x3FDA0FF6EB9DE4DE
800 data4 0x22859FA7, 0x00000000
801 data8 0x3FDAA5AB99ECF92D
802 data4 0x237A6D06, 0x00000000
803 data8 0x3FDB3C5A6D8F1796
804 data4 0x23952F6C, 0x00000000
805 data8 0x3FDBD40A9CFB8BE4
806 data4 0x2280FC95, 0x00000000
807 data8 0x3FDC6CC387943100
808 data4 0x245D2EC0, 0x00000000
809 data8 0x3FDD068CB736C500
810 data4 0x23C4AD7D, 0x00000000
811 data8 0x3FDDA16DE1DDBC31
812 data4 0x23D076E6, 0x00000000
813 data8 0x3FDE3D6EEB515A93
814 data4 0x244809A6, 0x00000000
815 data8 0x3FDEDA97E6E9E5F1
816 data4 0x220856C8, 0x00000000
817 data8 0x3FDF78F11963CE69
818 data4 0x244BE993, 0x00000000
819 data8 0x3FE00C417D635BCE
820 data4 0x23D21799, 0x00000000
821 data8 0x3FE05CAB1C302CD3
822 data4 0x248A1B1D, 0x00000000
823 data8 0x3FE0ADB9DB6A1FA0
824 data4 0x23D53E33, 0x00000000
825 data8 0x3FE0FF724A20BA81
826 data4 0x24DB9ED5, 0x00000000
827 data8 0x3FE151D9153FA6F5
828 data4 0x24E9E451, 0x00000000
829 LOCAL_OBJECT_END(tanl_table_tm2)
831 LOCAL_OBJECT_START(tanl_table_tm1)
833 // Entries T_hi double-precision memory format
834 // Index = 0,1,...,19 B = 2^(-1)*(1+Index/32+1/64)
835 // Entries T_lo single-precision memory format
836 // Index = 0,1,...,19 B = 2^(-1)*(1+Index/32+1/64)
838 data8 0x3FE1CEC4BA1BE39E
839 data4 0x24B60F9E, 0x00000000
840 data8 0x3FE277E45ABD9B2D
841 data4 0x248C2474, 0x00000000
842 data8 0x3FE324180272B110
843 data4 0x247B8311, 0x00000000
844 data8 0x3FE3D38B890E2DF0
845 data4 0x24C55751, 0x00000000
846 data8 0x3FE4866D46236871
847 data4 0x24E5BC34, 0x00000000
848 data8 0x3FE53CEE45E044B0
849 data4 0x24001BA4, 0x00000000
850 data8 0x3FE5F74282EC06E4
851 data4 0x24B973DC, 0x00000000
852 data8 0x3FE6B5A125DF43F9
853 data4 0x24895440, 0x00000000
854 data8 0x3FE77844CAFD348C
855 data4 0x240021CA, 0x00000000
856 data8 0x3FE83F6BCEED6B92
857 data4 0x24C45372, 0x00000000
858 data8 0x3FE90B58A34F3665
859 data4 0x240DAD33, 0x00000000
860 data8 0x3FE9DC522C1E56B4
861 data4 0x24F846CE, 0x00000000
862 data8 0x3FEAB2A427041578
863 data4 0x2323FB6E, 0x00000000
864 data8 0x3FEB8E9F9DD8C373
865 data4 0x24B3090B, 0x00000000
866 data8 0x3FEC709B65C9AA7B
867 data4 0x2449F611, 0x00000000
868 data8 0x3FED58F4ACCF8435
869 data4 0x23616A7E, 0x00000000
870 data8 0x3FEE480F97635082
871 data4 0x24C2FEAE, 0x00000000
872 data8 0x3FEF3E57F0ACC544
873 data4 0x242CE964, 0x00000000
874 data8 0x3FF01E20F7E06E4B
875 data4 0x2480D3EE, 0x00000000
876 data8 0x3FF0A1258A798A69
877 data4 0x24DB8967, 0x00000000
878 LOCAL_OBJECT_END(tanl_table_tm1)
880 LOCAL_OBJECT_START(tanl_table_cm2)
882 // Entries C_hi double-precision memory format
883 // Index = 0,1,...,31 B = 2^(-2)*(1+Index/32+1/64)
884 // Entries C_lo single-precision memory format
885 // Index = 0,1,...,31 B = 2^(-2)*(1+Index/32+1/64)
887 data8 0x400ED3E2E63EFBD0
888 data4 0x259D94D4, 0x00000000
889 data8 0x400DDDB4C515DAB5
890 data4 0x245F0537, 0x00000000
891 data8 0x400CF57ABE19A79F
892 data4 0x25D4EA9F, 0x00000000
893 data8 0x400C1A06D15298ED
894 data4 0x24AE40A0, 0x00000000
895 data8 0x400B4A4C164B2708
896 data4 0x25A5AAB6, 0x00000000
897 data8 0x400A855A5285B068
898 data4 0x25524F18, 0x00000000
899 data8 0x4009CA5A3FFA549F
900 data4 0x24C999C0, 0x00000000
901 data8 0x4009188A646AF623
902 data4 0x254FD801, 0x00000000
903 data8 0x40086F3C6084D0E7
904 data4 0x2560F5FD, 0x00000000
905 data8 0x4007CDD2A29A76EE
906 data4 0x255B9D19, 0x00000000
907 data8 0x400733BE6C8ECA95
908 data4 0x25CB021B, 0x00000000
909 data8 0x4006A07E1F8DDC52
910 data4 0x24AB4722, 0x00000000
911 data8 0x4006139BC298AD58
912 data4 0x252764E2, 0x00000000
913 data8 0x40058CABBAD7164B
914 data4 0x24DAF5DB, 0x00000000
915 data8 0x40050B4BAE31A5D3
916 data4 0x25EA20F4, 0x00000000
917 data8 0x40048F2189F85A8A
918 data4 0x2583A3E8, 0x00000000
919 data8 0x400417DAA862380D
920 data4 0x25DCC4CC, 0x00000000
921 data8 0x4003A52B1088FCFE
922 data4 0x2430A492, 0x00000000
923 data8 0x400336CCCD3527D5
924 data4 0x255F77CF, 0x00000000
925 data8 0x4002CC7F5760766D
926 data4 0x25DA0BDA, 0x00000000
927 data8 0x4002660711CE02E3
928 data4 0x256FF4A2, 0x00000000
929 data8 0x4002032CD37BBE04
930 data4 0x25208AED, 0x00000000
931 data8 0x4001A3BD7F050775
932 data4 0x24B72DD6, 0x00000000
933 data8 0x40014789A554848A
934 data4 0x24AB4DAA, 0x00000000
935 data8 0x4000EE65323E81B7
936 data4 0x2584C440, 0x00000000
937 data8 0x4000982721CF1293
938 data4 0x25C9428D, 0x00000000
939 data8 0x400044A93D415EEB
940 data4 0x25DC8482, 0x00000000
941 data8 0x3FFFE78FBD72C577
942 data4 0x257F5070, 0x00000000
943 data8 0x3FFF4AC375EFD28E
944 data4 0x23EBBF7A, 0x00000000
945 data8 0x3FFEB2AF60B52DDE
946 data4 0x22EECA07, 0x00000000
947 data8 0x3FFE1F1935204180
948 data4 0x24191079, 0x00000000
949 data8 0x3FFD8FCA54F7E60A
950 data4 0x248D3058, 0x00000000
951 LOCAL_OBJECT_END(tanl_table_cm2)
953 LOCAL_OBJECT_START(tanl_table_cm1)
955 // Entries C_hi double-precision memory format
956 // Index = 0,1,...,19 B = 2^(-1)*(1+Index/32+1/64)
957 // Entries C_lo single-precision memory format
958 // Index = 0,1,...,19 B = 2^(-1)*(1+Index/32+1/64)
960 data8 0x3FFCC06A79F6FADE
961 data4 0x239C7886, 0x00000000
962 data8 0x3FFBB91F891662A6
963 data4 0x250BD191, 0x00000000
964 data8 0x3FFABFB6529F155D
965 data4 0x256CC3E6, 0x00000000
966 data8 0x3FF9D3002E964AE9
967 data4 0x250843E3, 0x00000000
968 data8 0x3FF8F1EF89DCB383
969 data4 0x2277C87E, 0x00000000
970 data8 0x3FF81B937C87DBD6
971 data4 0x256DA6CF, 0x00000000
972 data8 0x3FF74F141042EDE4
973 data4 0x2573D28A, 0x00000000
974 data8 0x3FF68BAF1784B360
975 data4 0x242E489A, 0x00000000
976 data8 0x3FF5D0B57C923C4C
977 data4 0x2532D940, 0x00000000
978 data8 0x3FF51D88F418EF20
979 data4 0x253C7DD6, 0x00000000
980 data8 0x3FF4719A02F88DAE
981 data4 0x23DB59BF, 0x00000000
982 data8 0x3FF3CC6649DA0788
983 data4 0x252B4756, 0x00000000
984 data8 0x3FF32D770B980DB8
985 data4 0x23FE585F, 0x00000000
986 data8 0x3FF2945FE56C987A
987 data4 0x25378A63, 0x00000000
988 data8 0x3FF200BDB16523F6
989 data4 0x247BB2E0, 0x00000000
990 data8 0x3FF172358CE27778
991 data4 0x24446538, 0x00000000
992 data8 0x3FF0E873FDEFE692
993 data4 0x2514638F, 0x00000000
994 data8 0x3FF0632C33154062
995 data4 0x24A7FC27, 0x00000000
996 data8 0x3FEFC42EB3EF115F
997 data4 0x248FD0FE, 0x00000000
998 data8 0x3FEEC9E8135D26F6
999 data4 0x2385C719, 0x00000000
1000 LOCAL_OBJECT_END(tanl_table_cm1)
1002 LOCAL_OBJECT_START(tanl_table_scim2)
1004 // Entries SC_inv in Swapped IEEE format (extended)
1005 // Index = 0,1,...,31 B = 2^(-2)*(1+Index/32+1/64)
1007 data8 0x839D6D4A1BF30C9E, 0x00004001
1008 data8 0x80092804554B0EB0, 0x00004001
1009 data8 0xF959F94CA1CF0DE9, 0x00004000
1010 data8 0xF3086BA077378677, 0x00004000
1011 data8 0xED154515CCD4723C, 0x00004000
1012 data8 0xE77909441C27CF25, 0x00004000
1013 data8 0xE22D037D8DDACB88, 0x00004000
1014 data8 0xDD2B2D8A89C73522, 0x00004000
1015 data8 0xD86E1A23BB2C1171, 0x00004000
1016 data8 0xD3F0E288DFF5E0F9, 0x00004000
1017 data8 0xCFAF16B1283BEBD5, 0x00004000
1018 data8 0xCBA4AFAA0D88DD53, 0x00004000
1019 data8 0xC7CE03CCCA67C43D, 0x00004000
1020 data8 0xC427BC820CA0DDB0, 0x00004000
1021 data8 0xC0AECD57F13D8CAB, 0x00004000
1022 data8 0xBD606C3871ECE6B1, 0x00004000
1023 data8 0xBA3A0A96A44C4929, 0x00004000
1024 data8 0xB7394F6FE5CCCEC1, 0x00004000
1025 data8 0xB45C12039637D8BC, 0x00004000
1026 data8 0xB1A0552892CB051B, 0x00004000
1027 data8 0xAF04432B6BA2FFD0, 0x00004000
1028 data8 0xAC862A237221235F, 0x00004000
1029 data8 0xAA2478AF5F00A9D1, 0x00004000
1030 data8 0xA7DDBB0C81E082BF, 0x00004000
1031 data8 0xA5B0987D45684FEE, 0x00004000
1032 data8 0xA39BD0F5627A8F53, 0x00004000
1033 data8 0xA19E3B036EC5C8B0, 0x00004000
1034 data8 0x9FB6C1F091CD7C66, 0x00004000
1035 data8 0x9DE464101FA3DF8A, 0x00004000
1036 data8 0x9C263139A8F6B888, 0x00004000
1037 data8 0x9A7B4968C27B0450, 0x00004000
1038 data8 0x98E2DB7E5EE614EE, 0x00004000
1039 LOCAL_OBJECT_END(tanl_table_scim2)
1041 LOCAL_OBJECT_START(tanl_table_scim1)
1043 // Entries SC_inv in Swapped IEEE format (extended)
1044 // Index = 0,1,...,19 B = 2^(-1)*(1+Index/32+1/64)
1046 data8 0x969F335C13B2B5BA, 0x00004000
1047 data8 0x93D446D9D4C0F548, 0x00004000
1048 data8 0x9147094F61B798AF, 0x00004000
1049 data8 0x8EF317CC758787AC, 0x00004000
1050 data8 0x8CD498B3B99EEFDB, 0x00004000
1051 data8 0x8AE82A7DDFF8BC37, 0x00004000
1052 data8 0x892AD546E3C55D42, 0x00004000
1053 data8 0x8799FEA9D15573C1, 0x00004000
1054 data8 0x86335F88435A4B4C, 0x00004000
1055 data8 0x84F4FB6E3E93A87B, 0x00004000
1056 data8 0x83DD195280A382FB, 0x00004000
1057 data8 0x82EA3D7FA4CB8C9E, 0x00004000
1058 data8 0x821B247C6861D0A8, 0x00004000
1059 data8 0x816EBED163E8D244, 0x00004000
1060 data8 0x80E42D9127E4CFC6, 0x00004000
1061 data8 0x807ABF8D28E64AFD, 0x00004000
1062 data8 0x8031EF26863B4FD8, 0x00004000
1063 data8 0x800960ADAE8C11FD, 0x00004000
1064 data8 0x8000E1475FDBEC21, 0x00004000
1065 data8 0x80186650A07791FA, 0x00004000
1066 LOCAL_OBJECT_END(tanl_table_scim1)
1069 Save_Norm_Arg = f8 // For input to reduction routine
1071 r = f8 // For output from reduction routine
1072 c = f9 // For output from reduction routine
1127 NEGTWO_TO_NEG14 = f76
1128 NEGTWO_TO_NEG33 = f77
1144 NEGTWO_TO_NEG2 = f93
1165 FR_inv_pi_2to63 = f113
1166 FR_rshf_2to64 = f114
1198 GR_exp_2_to_63 = r55
1199 GR_exp_2_to_24 = r56
1204 GR_exp_m2tom14 = r61
1206 GR_exp_m2tom33 = r63
1212 GR_Parameter_X = r67
1213 GR_Parameter_Y = r68
1214 GR_Parameter_RESULT = r69
1215 GR_Parameter_Tag = r70
1219 .global __libm_tanl#
1220 .global __libm_cotl#
1225 LOCAL_LIBM_ENTRY(cotl)
1228 alloc r32 = ar.pfs, 0,35,4,0
1229 movl GR_sig_inv_pi = 0xa2f9836e4e44152a // significand of 1/pi
1232 mov GR_exp_mask = 0x1ffff // Exponent mask
1233 movl GR_rshf_2to64 = 0x47e8000000000000 // 1.1000 2^(63+64)
1237 // Check for NatVals, Infs , NaNs, and Zeros
1239 getf.exp GR_signexp_x = Arg // Get sign and exponent of x
1240 fclass.m p6,p0 = Arg, 0x1E7 // Test for natval, nan, inf, zero
1244 addl table_base = @ltoff(TANL_BASE_CONSTANTS), gp // Pointer to table ptr
1245 fnorm.s1 Norm_Arg = Arg // Normalize x
1246 br.cond.sptk COMMON_PATH
1249 LOCAL_LIBM_END(cotl)
1255 GLOBAL_IEEE754_ENTRY(tanl)
1258 alloc r32 = ar.pfs, 0,35,4,0
1259 movl GR_sig_inv_pi = 0xa2f9836e4e44152a // significand of 1/pi
1262 mov GR_exp_mask = 0x1ffff // Exponent mask
1263 movl GR_rshf_2to64 = 0x47e8000000000000 // 1.1000 2^(63+64)
1267 // Check for NatVals, Infs , NaNs, and Zeros
1269 getf.exp GR_signexp_x = Arg // Get sign and exponent of x
1270 fclass.m p6,p0 = Arg, 0x1E7 // Test for natval, nan, inf, zero
1274 addl table_base = @ltoff(TANL_BASE_CONSTANTS), gp // Pointer to table ptr
1275 fnorm.s1 Norm_Arg = Arg // Normalize x
1279 // Common path for both tanl and cotl
1282 setf.sig FR_inv_pi_2to63 = GR_sig_inv_pi // Form 1/pi * 2^63
1283 fclass.m p9, p0 = Arg, 0x0b // Test x denormal
1284 mov GR_exp_2tom64 = 0xffff - 64 // Scaling constant to compute N
1287 setf.d FR_rshf_2to64 = GR_rshf_2to64 // Form const 1.1000 * 2^(63+64)
1288 movl GR_rshf = 0x43e8000000000000 // Form const 1.1000 * 2^63
1292 // Check for everything - if false, then must be pseudo-zero or pseudo-nan.
1293 // Branch out to deal with special values.
1296 fclass.nm p7,p0 = Arg, 0x1FF // Test x unsupported
1297 mov GR_exp_2_to_63 = 0xffff + 63 // Exponent of 2^63
1300 ld8 table_base = [table_base] // Get pointer to constant table
1301 fms.s1 mOne = f0, f0, f1
1302 (p6) br.cond.spnt TANL_SPECIAL // Branch if x natval, nan, inf, zero
1307 setf.sig fp_tmp = gr_tmp // Make a constant so fmpy produces inexact
1308 mov GR_exp_2_to_24 = 0xffff + 24 // Exponent of 2^24
1309 (p9) br.cond.spnt TANL_DENORMAL // Branch if x denormal
1314 // Return to here if x denormal
1316 // Do fcmp to generate Denormal exception
1317 // - can't do FNORM (will generate Underflow when U is unmasked!)
1318 // Branch out to deal with unsupporteds values.
1320 setf.exp FR_2tom64 = GR_exp_2tom64 // Form 2^-64 for scaling N_float
1321 fcmp.eq.s0 p0, p6 = Arg, f1 // Dummy to flag denormals
1322 add table_ptr1 = 0, table_base // Point to tanl_table_1
1325 setf.d FR_rshf = GR_rshf // Form right shift const 1.1000 * 2^63
1326 add table_ptr2 = 80, table_base // Point to tanl_table_2
1327 (p7) br.cond.spnt TANL_UNSUPPORTED // Branch if x unsupported type
1332 and GR_exp_x = GR_exp_mask, GR_signexp_x // Get exponent of x
1333 fmpy.s1 Save_Norm_Arg = Norm_Arg, f1 // Save x if large arg reduction
1334 dep.z bmask1 = 0x7c, 56, 8 // Form mask to get 5 msb of r
1335 // bmask1 = 0x7c00000000000000
1340 // Decide about the paths to take:
1341 // Set PR_6 if |Arg| >= 2**63
1342 // Set PR_9 if |Arg| < 2**24 - CASE 1 OR 2
1343 // OTHERWISE Set PR_8 - CASE 3 OR 4
1345 // Branch out if the magnitude of the input argument is >= 2^63
1346 // - do this branch before the next.
1348 ldfe two_by_PI = [table_ptr1],16 // Load 2/pi
1350 dep.z bmask2 = 0x41, 57, 7 // Form mask to OR to produce B
1351 // bmask2 = 0x8200000000000000
1354 ldfe PI_BY_4 = [table_ptr2],16 // Load pi/4
1355 cmp.ge p6,p0 = GR_exp_x, GR_exp_2_to_63 // Is |x| >= 2^63
1356 (p6) br.cond.spnt TANL_ARG_TOO_LARGE // Branch if |x| >= 2^63
1361 ldfe P_0 = [table_ptr1],16 // Load P_0
1362 ldfe Inv_P_0 = [table_ptr2],16 // Load Inv_P_0
1368 ldfe P_1 = [table_ptr1],16 // Load P_1
1369 fmerge.s Abs_Arg = f0, Norm_Arg // Get |x|
1370 mov GR_exp_m2tom33 = 0x2ffff - 33 // Form signexp of -2^-33
1373 ldfe d_1 = [table_ptr2],16 // Load d_1 for 2^24 <= |x| < 2^63
1375 mov GR_exp_2tom33 = 0xffff - 33 // Form signexp of 2^-33
1380 ldfe P_2 = [table_ptr1],16 // Load P_2
1381 ldfe d_2 = [table_ptr2],16 // Load d_2 for 2^24 <= |x| < 2^63
1382 cmp.ge p8,p0 = GR_exp_x, GR_exp_2_to_24 // Is |x| >= 2^24
1386 // Use special scaling to right shift so N=Arg * 2/pi is in rightmost bits
1387 // Branch to Cases 3 or 4 if Arg <= -2**24 or Arg >= 2**24
1389 ldfe P_3 = [table_ptr1],16 // Load P_3
1390 fma.s1 N_fix = Norm_Arg, FR_inv_pi_2to63, FR_rshf_2to64
1391 (p8) br.cond.spnt TANL_LARGER_ARG // Branch if 2^24 <= |x| < 2^63
1395 // Here if 0 < |x| < 2^24
1396 // ARGUMENT REDUCTION CODE - CASE 1 and 2
1399 setf.exp TWO_TO_NEG33 = GR_exp_2tom33 // Form 2^-33
1400 setf.exp NEGTWO_TO_NEG33 = GR_exp_m2tom33 // Form -2^-33
1401 fmerge.s r = Norm_Arg,Norm_Arg // Assume r=x, ok if |x| < pi/4
1406 // If |Arg| < pi/4, set PR_8, else pi/4 <=|Arg| < 2^24 - set PR_9.
1408 // Case 2: Convert integer N_fix back to normalized floating-point value.
1410 getf.sig sig_r = Norm_Arg // Get sig_r if 1/4 <= |x| < pi/4
1411 fcmp.lt.s1 p8,p9= Abs_Arg,PI_BY_4 // Test |x| < pi/4
1412 mov GR_exp_2tom2 = 0xffff - 2 // Form signexp of 2^-2
1415 ldfps TWO_TO_NEG2, NEGTWO_TO_NEG2 = [table_ptr2] // Load 2^-2, -2^-2
1416 fms.s1 N = N_fix, FR_2tom64, FR_rshf // Use scaling to get N floated
1417 mov N_fix_gr = r0 // Assume N=0, ok if |x| < pi/4
1422 // Case 1: Is |r| < 2**(-2).
1423 // Arg is the same as r in this case.
1427 // Case 2: Place integer part of N in GP register.
1429 (p9) getf.sig N_fix_gr = N_fix
1430 fmerge.s c = f0, f0 // Assume c=0, ok if |x| < pi/4
1431 cmp.lt p10, p0 = GR_exp_x, GR_exp_2tom2 // Test if |x| < 1/4
1436 setf.sig B_mask1 = bmask1 // Form mask to get 5 msb of r
1438 mov exp_r = GR_exp_x // Get exp_r if 1/4 <= |x| < pi/4
1441 setf.sig B_mask2 = bmask2 // Form mask to form B from r
1442 (p10) br.cond.spnt TANL_SMALL_R // Branch if 0 < |x| < 1/4
1443 (p8) br.cond.spnt TANL_NORMAL_R // Branch if 1/4 <= |x| < pi/4
1447 // Here if pi/4 <= |x| < 2^24
1449 // Case 1: PR_3 is only affected when PR_1 is set.
1452 // Case 2: w = N * P_2
1453 // Case 2: s_val = -N * P_1 + Arg
1458 fnma.s1 s_val = N, P_1, Norm_Arg
1463 fmpy.s1 w = N, P_2 // w = N * P_2 for |s| >= 2^-33
1468 // Case 2_reduce: w = N * P_3 (change sign)
1471 fmpy.s1 w2 = N, P_3 // w = N * P_3 for |s| < 2^-33
1476 // Case 1_reduce: r = s + w (change sign)
1479 fsub.s1 r = s_val, w // r = s_val - w for |s| >= 2^-33
1484 // Case 2_reduce: U_1 = N * P_2 + w
1487 fma.s1 U_1 = N, P_2, w2 // U_1 = N * P_2 + w for |s| < 2^-33
1493 // Decide between case_1 and case_2 reduce:
1494 // Case 1_reduce: |s| >= 2**(-33)
1495 // Case 2_reduce: |s| < 2**(-33)
1499 fcmp.lt.s1 p9, p8 = s_val, TWO_TO_NEG33
1506 (p9) fcmp.gt.s1 p9, p8 = s_val, NEGTWO_TO_NEG33
1511 // Case 1_reduce: c = s - r
1514 fsub.s1 c = s_val, r // c = s_val - r for |s| >= 2^-33
1519 // Case 2_reduce: r is complete here - continue to calculate c .
1523 (p9) fsub.s1 r = s_val, U_1
1528 (p9) fms.s1 U_2 = N, P_2, U_1
1534 // Case 1_reduce: Is |r| < 2**(-2), if so set PR_10
1545 (p8) fcmp.lt.unc.s1 p10, p13 = r, TWO_TO_NEG2
1551 (p8) getf.sig sig_r = r // Get signif of r if |s| >= 2^-33
1558 (p8) getf.exp exp_r = r // Extract signexp of r if |s| >= 2^-33
1559 (p10) fcmp.gt.s1 p10, p13 = r, NEGTWO_TO_NEG2
1564 // Case 1_reduce: c is complete here.
1565 // Case 1: Branch to SMALL_R or NORMAL_R.
1566 // c = c + w (w has not been negated.)
1569 (p8) fsub.s1 c = c, w // c = c - w for |s| >= 2^-33
1574 (p10) br.cond.spnt TANL_SMALL_R // Branch if pi/4 < |x| < 2^24 and |r|<1/4
1575 (p13) br.cond.sptk TANL_NORMAL_R_A // Branch if pi/4 < |x| < 2^24 and |r|>=1/4
1580 // Here if pi/4 < |x| < 2^24 and |s| < 2^-33
1582 // Is i_1 = lsb of N_fix_gr even or odd?
1583 // if i_1 == 0, set p11, else set p12.
1587 fsub.s1 s_val = s_val, r
1588 add N_fix_gr = N_fix_gr, cot_flag // N = N + 1 (for cotl)
1594 // U_2 = N * P_2 - U_1
1595 // Not needed until later.
1597 fadd.s1 U_2 = U_2, w2
1610 // c is complete here
1611 // Argument reduction ends here.
1616 tbit.z p11, p12 = N_fix_gr, 0 ;; // Set p11 if N even, p12 if odd
1621 (p12) frcpa.s1 S_hi,p0 = f1, r
1626 fsub.s1 c = s_val, U_1
1632 add table_ptr1 = 160, table_base ;; // Point to tanl_table_p1
1633 ldfe P1_1 = [table_ptr1],144
1637 // Load P1_1 and point to Q1_1 .
1640 ldfe Q1_1 = [table_ptr1]
1642 // N even: rsq = r * Z
1643 // N odd: S_hi = frcpa(r)
1645 (p12) fmerge.ns S_hi = S_hi, S_hi
1654 (p9) fsub.s1 c = c, U_2
1659 (p12) fma.s1 poly1 = S_hi, r, f1
1665 // N odd: Change sign of S_hi
1667 (p11) fmpy.s1 rsq = rsq, P1_1
1672 (p12) fma.s1 S_hi = S_hi, poly1, S_hi
1678 // N even: rsq = rsq * P1_1
1679 // N odd: poly1 = 1.0 + S_hi * r 16 bits partial account for necessary
1681 (p11) fma.s1 Poly = r, rsq, c
1687 // N even: Poly = c + r * rsq
1688 // N odd: S_hi = S_hi + S_hi*poly1 16 bits account for necessary
1690 (p12) fma.s1 poly1 = S_hi, r, f1
1691 (p11) tbit.z.unc p14, p15 = cot_flag, 0 ;; // p14=1 for tanl; p15=1 for cotl
1696 // N even: Result = Poly + r
1697 // N odd: poly1 = 1.0 + S_hi * r 32 bits partial
1699 (p14) fadd.s0 Result = r, Poly // for tanl
1704 (p15) fms.s0 Result = r, mOne, Poly // for cotl
1711 (p12) fma.s1 S_hi = S_hi, poly1, S_hi
1717 // N even: Result1 = Result + r
1718 // N odd: S_hi = S_hi * poly1 + S_hi 32 bits
1720 (p12) fma.s1 poly1 = S_hi, r, f1
1726 // N odd: poly1 = S_hi * r + 1.0 64 bits partial
1728 (p12) fma.s1 S_hi = S_hi, poly1, S_hi
1734 // N odd: poly1 = S_hi * poly + 1.0 64 bits
1736 (p12) fma.s1 poly1 = S_hi, r, f1
1742 // N odd: poly1 = S_hi * r + 1.0
1744 (p12) fma.s1 poly1 = S_hi, c, poly1
1750 // N odd: poly1 = S_hi * c + poly1
1752 (p12) fmpy.s1 S_lo = S_hi, poly1
1758 // N odd: S_lo = S_hi * poly1
1760 (p12) fma.s1 S_lo = Q1_1, r, S_lo
1761 (p12) tbit.z.unc p14, p15 = cot_flag, 0 // p14=1 for tanl; p15=1 for cotl
1766 // N odd: Result = S_hi + S_lo
1768 fmpy.s0 fp_tmp = fp_tmp, fp_tmp // Dummy mult to set inexact
1774 // N odd: S_lo = S_lo + Q1_1 * r
1776 (p14) fadd.s0 Result = S_hi, S_lo // for tanl
1781 (p15) fms.s0 Result = S_hi, mOne, S_lo // for cotl
1782 br.ret.sptk b0 ;; // Exit for pi/4 <= |x| < 2^24 and |s| < 2^-33
1787 // Here if 2^24 <= |x| < 2^63
1789 // ARGUMENT REDUCTION CODE - CASE 3 and 4
1793 mov GR_exp_2tom14 = 0xffff - 14 // Form signexp of 2^-14
1794 mov GR_exp_m2tom14 = 0x2ffff - 14 // Form signexp of -2^-14
1795 fmpy.s1 N_0 = Norm_Arg, Inv_P_0
1800 setf.exp TWO_TO_NEG14 = GR_exp_2tom14 // Form 2^-14
1801 setf.exp NEGTWO_TO_NEG14 = GR_exp_m2tom14// Form -2^-14
1808 // Adjust table_ptr1 to beginning of table.
1809 // N_0 = Arg * Inv_P_0
1812 add table_ptr2 = 144, table_base ;; // Point to 2^-2
1813 ldfps TWO_TO_NEG2, NEGTWO_TO_NEG2 = [table_ptr2]
1819 // N_0_fix = integer part of N_0 .
1822 // Make N_0 the integer part.
1826 fcvt.fx.s1 N_0_fix = N_0
1830 setf.sig B_mask1 = bmask1 // Form mask to get 5 msb of r
1831 fcvt.xf N_0 = N_0_fix
1835 setf.sig B_mask2 = bmask2 // Form mask to form B from r
1836 fnma.s1 ArgPrime = N_0, P_0, Norm_Arg
1841 fmpy.s1 w = N_0, d_1
1845 // ArgPrime = -N_0 * P_0 + Arg
1849 // N = ArgPrime * 2/pi
1851 // fcvt.fx.s1 N_fix = N
1852 // Use special scaling to right shift so N=Arg * 2/pi is in rightmost bits
1853 // Branch to Cases 3 or 4 if Arg <= -2**24 or Arg >= 2**24
1856 fma.s1 N_fix = ArgPrime, FR_inv_pi_2to63, FR_rshf_2to64
1860 // Convert integer N_fix back to normalized floating-point value.
1863 fms.s1 N = N_fix, FR_2tom64, FR_rshf // Use scaling to get N floated
1869 // N is the integer part of the reduced-reduced argument.
1870 // Put the integer in a GP register.
1873 getf.sig N_fix_gr = N_fix
1880 // s_val = -N*P_1 + ArgPrime
1885 fnma.s1 s_val = N, P_1, ArgPrime
1890 fnma.s1 w = N, P_2, w
1895 // Case 4: V_hi = N * P_2
1896 // Case 4: U_hi = N_0 * d_1
1899 fmpy.s1 V_hi = N, P_2 // V_hi = N * P_2 for |s| < 2^-14
1904 fmpy.s1 U_hi = N_0, d_1 // U_hi = N_0 * d_1 for |s| < 2^-14
1909 // Case 3: r = s_val + w (Z complete)
1910 // Case 4: w = N * P_3
1913 fadd.s1 r = s_val, w // r = s_val + w for |s| >= 2^-14
1918 fmpy.s1 w2 = N, P_3 // w = N * P_3 for |s| < 2^-14
1923 // Case 4: A = U_hi + V_hi
1924 // Note: Worry about switched sign of V_hi, so subtract instead of add.
1925 // Case 4: V_lo = -N * P_2 - V_hi (U_hi is in place of V_hi in writeup)
1926 // Note: the (-) is still missing for V_hi.
1929 fsub.s1 A = U_hi, V_hi // A = U_hi - V_hi for |s| < 2^-14
1934 fnma.s1 V_lo = N, P_2, V_hi // V_lo = V_hi - N * P_2 for |s| < 2^-14
1939 // Decide between case 3 and 4:
1940 // Case 3: |s| >= 2**(-14) Set p10
1941 // Case 4: |s| < 2**(-14) Set p11
1943 // Case 4: U_lo = N_0 * d_1 - U_hi
1946 fms.s1 U_lo = N_0, d_1, U_hi // U_lo = N_0*d_1 - U_hi for |s| < 2^-14
1951 fcmp.lt.s1 p11, p10 = s_val, TWO_TO_NEG14
1956 // Case 4: We need abs of both U_hi and V_hi - dont
1957 // worry about switched sign of V_hi.
1960 fabs V_hiabs = V_hi // |V_hi| for |s| < 2^-14
1965 (p11) fcmp.gt.s1 p11, p10 = s_val, NEGTWO_TO_NEG14
1970 // Case 3: c = s_val - r
1973 fabs U_hiabs = U_hi // |U_hi| for |s| < 2^-14
1978 fsub.s1 c = s_val, r // c = s_val - r for |s| >= 2^-14
1983 // For Case 3, |s| >= 2^-14, determine if |r| < 1/4
1985 // Case 4: C_hi = s_val + A
1989 (p11) fadd.s1 C_hi = s_val, A // C_hi = s_val + A for |s| < 2^-14
1994 (p10) fcmp.lt.unc.s1 p14, p15 = r, TWO_TO_NEG2
2000 getf.sig sig_r = r // Get signif of r if |s| >= 2^-33
2006 // Case 4: t = U_lo + V_lo
2008 getf.exp exp_r = r // Extract signexp of r if |s| >= 2^-33
2009 (p11) fadd.s1 t = U_lo, V_lo // t = U_lo + V_lo for |s| < 2^-14
2014 (p14) fcmp.gt.s1 p14, p15 = r, NEGTWO_TO_NEG2
2019 // Case 3: c = (s - r) + w (c complete)
2022 (p10) fadd.s1 c = c, w // c = c + w for |s| >= 2^-14
2027 (p14) br.cond.spnt TANL_SMALL_R // Branch if 2^24 <= |x| < 2^63 and |r|< 1/4
2028 (p15) br.cond.sptk TANL_NORMAL_R_A // Branch if 2^24 <= |x| < 2^63 and |r|>=1/4
2033 // Here if 2^24 <= |x| < 2^63 and |s| < 2^-14 >>>>>>> Case 4.
2035 // Case 4: Set P_12 if U_hiabs >= V_hiabs
2036 // Case 4: w = w + N_0 * d_2
2037 // Note: the (-) is now incorporated in w .
2039 add table_ptr1 = 160, table_base // Point to tanl_table_p1
2040 fcmp.ge.unc.s1 p12, p13 = U_hiabs, V_hiabs
2045 fms.s1 w2 = N_0, d_2, w2
2050 // Case 4: C_lo = s_val - C_hi
2052 ldfe P1_1 = [table_ptr1], 16 // Load P1_1
2053 fsub.s1 C_lo = s_val, C_hi
2059 // Case 4: a = U_hi - A
2060 // a = V_hi - A (do an add to account for missing (-) on V_hi
2063 ldfe P1_2 = [table_ptr1], 128 // Load P1_2
2064 (p12) fsub.s1 a = U_hi, A
2069 (p13) fadd.s1 a = V_hi, A
2074 // Case 4: t = U_lo + V_lo + w
2076 ldfe Q1_1 = [table_ptr1], 16 // Load Q1_1
2082 // Case 4: a = (U_hi - A) + V_hi
2083 // a = (V_hi - A) + U_hi
2084 // In each case account for negative missing form V_hi .
2087 ldfe Q1_2 = [table_ptr1], 16 // Load Q1_2
2088 (p12) fsub.s1 a = a, V_hi
2093 (p13) fsub.s1 a = U_hi, a
2099 // Case 4: C_lo = (s_val - C_hi) + A
2103 fadd.s1 C_lo = C_lo, A
2107 // Case 4: t = t + a
2116 // Case 4: C_lo = C_lo + t
2117 // Case 4: r = C_hi + C_lo
2120 fadd.s1 C_lo = C_lo, t
2127 fadd.s1 r = C_hi, C_lo
2133 // Case 4: c = C_hi - r
2143 add N_fix_gr = N_fix_gr, cot_flag // N = N + 1 (for cotl)
2147 // Case 4: c = c + C_lo finished.
2149 // Is i_1 = lsb of N_fix_gr even or odd?
2150 // if i_1 == 0, set PR_11, else set PR_12.
2154 fadd.s1 c = c , C_lo
2155 tbit.z p11, p12 = N_fix_gr, 0
2159 // r and c have been computed.
2162 (p12) frcpa.s1 S_hi, p0 = f1, r
2168 // N odd: Change sign of S_hi
2170 (p11) fma.s1 Poly = rsq, P1_2, P1_1
2175 (p12) fma.s1 P = rsq, Q1_2, Q1_1
2181 // N odd: Result = S_hi + S_lo (User supplied rounding mode for C1)
2183 fmpy.s0 fp_tmp = fp_tmp, fp_tmp // Dummy mult to set inexact
2189 // N even: rsq = r * r
2190 // N odd: S_hi = frcpa(r)
2192 (p12) fmerge.ns S_hi = S_hi, S_hi
2198 // N even: rsq = rsq * P1_2 + P1_1
2199 // N odd: poly1 = 1.0 + S_hi * r 16 bits partial account for necessary
2201 (p11) fmpy.s1 Poly = rsq, Poly
2206 (p12) fma.s1 poly1 = S_hi, r,f1
2207 (p11) tbit.z.unc p14, p15 = cot_flag, 0 // p14=1 for tanl; p15=1 for cotl
2212 // N even: Poly = Poly * rsq
2213 // N odd: S_hi = S_hi + S_hi*poly1 16 bits account for necessary
2215 (p11) fma.s1 Poly = r, Poly, c
2220 (p12) fma.s1 S_hi = S_hi, poly1, S_hi
2226 // N odd: S_hi = S_hi * poly1 + S_hi 32 bits
2228 (p14) fadd.s0 Result = r, Poly // for tanl
2232 .pred.rel "mutex",p15,p12
2235 (p15) fms.s0 Result = r, mOne, Poly // for cotl
2240 (p12) fma.s1 poly1 = S_hi, r, f1
2246 // N even: Poly = Poly * r + c
2247 // N odd: poly1 = 1.0 + S_hi * r 32 bits partial
2249 (p12) fma.s1 S_hi = S_hi, poly1, S_hi
2254 (p12) fma.s1 poly1 = S_hi, r, f1
2260 // N even: Result = Poly + r (Rounding mode S0)
2261 // N odd: poly1 = S_hi * r + 1.0 64 bits partial
2263 (p12) fma.s1 S_hi = S_hi, poly1, S_hi
2269 // N odd: poly1 = S_hi * poly + S_hi 64 bits
2271 (p12) fma.s1 poly1 = S_hi, r, f1
2277 // N odd: poly1 = S_hi * r + 1.0
2279 (p12) fma.s1 poly1 = S_hi, c, poly1
2285 // N odd: poly1 = S_hi * c + poly1
2287 (p12) fmpy.s1 S_lo = S_hi, poly1
2293 // N odd: S_lo = S_hi * poly1
2295 (p12) fma.s1 S_lo = P, r, S_lo
2296 (p12) tbit.z.unc p14, p15 = cot_flag, 0 ;; // p14=1 for tanl; p15=1 for cotl
2301 (p14) fadd.s0 Result = S_hi, S_lo // for tanl
2307 // N odd: S_lo = S_lo + r * P
2309 (p15) fms.s0 Result = S_hi, mOne, S_lo // for cotl
2310 br.ret.sptk b0 ;; // Exit for 2^24 <= |x| < 2^63 and |s| < 2^-14
2315 // Here if |r| < 1/4
2316 // r and c have been computed.
2317 // *****************************************************************
2318 // *****************************************************************
2319 // *****************************************************************
2320 // N odd: S_hi = frcpa(r)
2321 // Get [i_1] - lsb of N_fix_gr. Set p11 if N even, p12 if N odd.
2322 // N even: rsq = r * r
2324 add table_ptr1 = 160, table_base // Point to tanl_table_p1
2325 frcpa.s1 S_hi, p0 = f1, r // S_hi for N odd
2326 add N_fix_gr = N_fix_gr, cot_flag // N = N + 1 (for cotl)
2329 add table_ptr2 = 400, table_base // Point to Q1_7
2336 ldfe P1_1 = [table_ptr1], 16
2338 ldfe P1_2 = [table_ptr1], 16
2339 tbit.z p11, p12 = N_fix_gr, 0
2345 ldfe P1_3 = [table_ptr1], 96
2352 (p11) ldfe P1_9 = [table_ptr1], -16
2353 (p12) fmerge.ns S_hi = S_hi, S_hi
2358 (p11) fmpy.s1 r_to_the_8 = rsq, rsq
2364 // N even: Poly2 = P1_7 + Poly2 * rsq
2365 // N odd: poly2 = Q1_5 + poly2 * rsq
2368 (p11) ldfe P1_8 = [table_ptr1], -16
2369 (p11) fadd.s1 CORR = rsq, f1
2375 // N even: Poly1 = P1_2 + P1_3 * rsq
2376 // N odd: poly1 = 1.0 + S_hi * r
2377 // 16 bits partial account for necessary (-1)
2380 (p11) ldfe P1_7 = [table_ptr1], -16
2382 (p11) ldfe P1_6 = [table_ptr1], -16
2388 // N even: Poly1 = P1_1 + Poly1 * rsq
2389 // N odd: S_hi = S_hi + S_hi * poly1) 16 bits account for necessary
2392 // N even: Poly2 = P1_5 + Poly2 * rsq
2393 // N odd: poly2 = Q1_3 + poly2 * rsq
2396 (p11) ldfe P1_5 = [table_ptr1], -16
2397 (p11) fmpy.s1 r_to_the_8 = r_to_the_8, r_to_the_8
2402 (p12) fma.s1 poly1 = S_hi, r, f1
2408 // N even: Poly1 = Poly1 * rsq
2409 // N odd: poly1 = 1.0 + S_hi * r 32 bits partial
2413 // N even: CORR = CORR * c
2414 // N odd: S_hi = S_hi * poly1 + S_hi 32 bits
2418 // N even: Poly2 = P1_6 + Poly2 * rsq
2419 // N odd: poly2 = Q1_4 + poly2 * rsq
2423 (p11) ldfe P1_4 = [table_ptr1], -16
2425 (p11) fmpy.s1 CORR = CORR, c
2431 (p11) fma.s1 Poly1 = P1_3, rsq, P1_2
2435 (p12) ldfe Q1_7 = [table_ptr2], -16
2436 (p12) fma.s1 S_hi = S_hi, poly1, S_hi
2440 (p12) ldfe Q1_6 = [table_ptr2], -16
2441 (p11) fma.s1 Poly2 = P1_9, rsq, P1_8
2445 (p12) ldfe Q1_5 = [table_ptr2], -16 ;;
2446 (p12) ldfe Q1_4 = [table_ptr2], -16
2450 (p12) ldfe Q1_3 = [table_ptr2], -16
2452 // N even: Poly2 = P1_8 + P1_9 * rsq
2453 // N odd: poly2 = Q1_6 + Q1_7 * rsq
2455 (p11) fma.s1 Poly1 = Poly1, rsq, P1_1
2459 (p12) ldfe Q1_2 = [table_ptr2], -16
2460 (p12) fma.s1 poly1 = S_hi, r, f1
2464 (p12) ldfe Q1_1 = [table_ptr2], -16
2465 (p11) fma.s1 Poly2 = Poly2, rsq, P1_7
2471 // N even: CORR = rsq + 1
2472 // N even: r_to_the_8 = rsq * rsq
2474 (p11) fmpy.s1 Poly1 = Poly1, rsq
2479 (p12) fma.s1 S_hi = S_hi, poly1, S_hi
2484 (p12) fma.s1 poly2 = Q1_7, rsq, Q1_6
2489 (p11) fma.s1 Poly2 = Poly2, rsq, P1_6
2494 (p12) fma.s1 poly1 = S_hi, r, f1
2499 (p12) fma.s1 poly2 = poly2, rsq, Q1_5
2504 (p11) fma.s1 Poly2= Poly2, rsq, P1_5
2509 (p12) fma.s1 S_hi = S_hi, poly1, S_hi
2514 (p12) fma.s1 poly2 = poly2, rsq, Q1_4
2520 // N even: r_to_the_8 = r_to_the_8 * r_to_the_8
2521 // N odd: poly1 = S_hi * r + 1.0 64 bits partial
2523 (p11) fma.s1 Poly2 = Poly2, rsq, P1_4
2529 // N even: Poly = CORR + Poly * r
2530 // N odd: P = Q1_1 + poly2 * rsq
2532 (p12) fma.s1 poly1 = S_hi, r, f1
2537 (p12) fma.s1 poly2 = poly2, rsq, Q1_3
2543 // N even: Poly2 = P1_4 + Poly2 * rsq
2544 // N odd: poly2 = Q1_2 + poly2 * rsq
2546 (p11) fma.s1 Poly = Poly2, r_to_the_8, Poly1
2551 (p12) fma.s1 poly1 = S_hi, c, poly1
2556 (p12) fma.s1 poly2 = poly2, rsq, Q1_2
2563 // N even: Poly = Poly1 + Poly2 * r_to_the_8
2564 // N odd: S_hi = S_hi * poly1 + S_hi 64 bits
2566 (p11) fma.s1 Poly = Poly, r, CORR
2572 // N even: Result = r + Poly (User supplied rounding mode)
2573 // N odd: poly1 = S_hi * c + poly1
2575 (p12) fmpy.s1 S_lo = S_hi, poly1
2576 (p11) tbit.z.unc p14, p15 = cot_flag, 0 // p14=1 for tanl; p15=1 for cotl
2580 (p12) fma.s1 P = poly2, rsq, Q1_1
2586 // N odd: poly1 = S_hi * r + 1.0
2589 // N odd: S_lo = S_hi * poly1
2591 (p14) fadd.s0 Result = Poly, r // for tanl
2596 (p15) fms.s0 Result = Poly, mOne, r // for cotl
2603 // N odd: S_lo = Q1_1 * c + S_lo
2605 (p12) fma.s1 S_lo = Q1_1, c, S_lo
2610 fmpy.s0 fp_tmp = fp_tmp, fp_tmp // Dummy mult to set inexact
2616 // N odd: Result = S_lo + r * P
2618 (p12) fma.s1 Result = P, r, S_lo
2619 (p12) tbit.z.unc p14, p15 = cot_flag, 0 ;; // p14=1 for tanl; p15=1 for cotl
2623 // N odd: Result = Result + S_hi (user supplied rounding mode)
2627 (p14) fadd.s0 Result = Result, S_hi // for tanl
2632 (p15) fms.s0 Result = Result, mOne, S_hi // for cotl
2633 br.ret.sptk b0 ;; // Exit |r| < 1/4 path
2638 // Here if 1/4 <= |x| < pi/4 or if |x| >= 2^63 and |r| >= 1/4
2639 // *******************************************************************
2640 // *******************************************************************
2641 // *******************************************************************
2643 // r and c have been computed.
2653 // Enter here if pi/4 <= |x| < 2^63 and |r| >= 1/4
2654 // Get the 5 bits or r for the lookup. 1.xxxxx ....
2656 add table_ptr1 = 416, table_base // Point to tanl_table_p2
2657 mov GR_exp_2tom65 = 0xffff - 65 // Scaling constant for B
2658 extr.u lookup = sig_r, 58, 5
2663 ldfe P2_1 = [table_ptr1], 16
2664 setf.exp TWO_TO_NEG65 = GR_exp_2tom65 // 2^-65 for scaling B if exp_r=-2
2665 add N_fix_gr = N_fix_gr, cot_flag // N = N + 1 (for cotl)
2669 .pred.rel "mutex",p11,p12
2670 // B = 2^63 * 1.xxxxx 100...0
2672 ldfe P2_2 = [table_ptr1], 16
2674 mov table_offset = 512 // Assume table offset is 512
2679 ldfe P2_3 = [table_ptr1], 16
2680 fmerge.s Pos_r = f1, r
2681 tbit.nz p8,p9 = exp_r, 0
2685 // Is B = 2** -2 or B= 2** -1? If 2**-1, then
2686 // we want an offset of 512 for table addressing.
2688 add table_ptr2 = 1296, table_base // Point to tanl_table_cm2
2689 (p9) shladd table_offset = lookup, 4, table_offset
2690 (p8) shladd table_offset = lookup, 4, r0
2695 add table_ptr1 = table_ptr1, table_offset // Point to T_hi
2696 add table_ptr2 = table_ptr2, table_offset // Point to C_hi
2697 add table_ptr3 = 2128, table_base // Point to tanl_table_scim2
2702 ldfd T_hi = [table_ptr1], 8 // Load T_hi
2704 ldfd C_hi = [table_ptr2], 8 // Load C_hi
2705 add table_ptr3 = table_ptr3, table_offset // Point to SC_inv
2712 // Convert B so it has the same exponent as Pos_r before subtracting
2714 ldfs T_lo = [table_ptr1] // Load T_lo
2715 (p9) fnma.s1 x = B, FR_2tom64, Pos_r
2720 (p8) fnma.s1 x = B, TWO_TO_NEG65, Pos_r
2726 ldfs C_lo = [table_ptr2] // Load C_lo
2733 ldfe SC_inv = [table_ptr3] // Load SC_inv
2734 fmerge.s sgn_r = r, f1
2735 tbit.z p11, p12 = N_fix_gr, 0 // p11 if N even, p12 if odd
2742 // N even: Tx = T_hi * x
2744 // N even: Tx1 = Tx + 1
2745 // N odd: Cx1 = 1 - Cx
2755 (p11) fmpy.s1 Tx = T_hi, x
2761 // N odd: Cx = C_hi * x
2765 (p12) fmpy.s1 Cx = C_hi, x
2770 // N even and odd: P = P2_3 + P2_2 * xsq
2774 fma.s1 P = P2_3, xsq, P2_2
2779 (p11) fadd.s1 Tx1 = Tx, f1
2785 // N even: D = C_hi - tanx
2786 // N odd: D = T_hi + tanx
2788 (p11) fmpy.s1 CORR = SC_inv, T_hi
2793 fmpy.s1 Sx = SC_inv, x
2798 (p12) fmpy.s1 CORR = SC_inv, C_hi
2803 (p12) fsub.s1 V_hi = f1, Cx
2808 fma.s1 P = P, xsq, P2_1
2814 // N even and odd: P = P2_1 + P * xsq
2816 (p11) fma.s1 V_hi = Tx, Tx1, f1
2822 // N even: Result = sgn_r * tail + T_hi (user rounding mode for C1)
2823 // N odd: Result = sgn_r * tail + C_hi (user rounding mode for C1)
2825 fmpy.s0 fp_tmp = fp_tmp, fp_tmp // Dummy mult to set inexact
2830 fmpy.s1 CORR = CORR, c
2835 (p12) fnma.s1 V_hi = Cx,V_hi,f1
2841 // N even: V_hi = Tx * Tx1 + 1
2842 // N odd: Cx1 = 1 - Cx * Cx1
2850 // N even and odd: P = P * xsq
2852 (p11) fmpy.s1 V_hi = V_hi, T_hi
2858 // N even and odd: tail = P * tail + V_lo
2860 (p11) fmpy.s1 T_hi = sgn_r, T_hi
2865 fmpy.s1 CORR = CORR, sgn_r
2870 (p12) fmpy.s1 V_hi = V_hi,C_hi
2876 // N even: V_hi = T_hi * V_hi
2877 // N odd: V_hi = C_hi * V_hi
2879 fma.s1 tanx = P, x, x
2884 (p12) fnmpy.s1 C_hi = sgn_r, C_hi
2890 // N even: V_lo = 1 - V_hi + C_hi
2891 // N odd: V_lo = 1 - V_hi + T_hi
2893 (p11) fadd.s1 CORR = CORR, T_lo
2898 (p12) fsub.s1 CORR = CORR, C_lo
2904 // N even and odd: tanx = x + x * P
2905 // N even and odd: Sx = SC_inv * x
2907 (p11) fsub.s1 D = C_hi, tanx
2912 (p12) fadd.s1 D = T_hi, tanx
2918 // N odd: CORR = SC_inv * C_hi
2919 // N even: CORR = SC_inv * T_hi
2921 fnma.s1 D = V_hi, D, f1
2927 // N even and odd: D = 1 - V_hi * D
2928 // N even and odd: CORR = CORR * c
2930 fma.s1 V_hi = V_hi, D, V_hi
2936 // N even and odd: V_hi = V_hi + V_hi * D
2937 // N even and odd: CORR = sgn_r * CORR
2939 (p11) fnma.s1 V_lo = V_hi, C_hi, f1
2944 (p12) fnma.s1 V_lo = V_hi, T_hi, f1
2950 // N even: CORR = COOR + T_lo
2951 // N odd: CORR = CORR - C_lo
2953 (p11) fma.s1 V_lo = tanx, V_hi, V_lo
2954 tbit.nz p15, p0 = cot_flag, 0 // p15=1 if we compute cotl
2958 (p12) fnma.s1 V_lo = tanx, V_hi, V_lo
2964 (p15) fms.s1 T_hi = f0, f0, T_hi // to correct result's sign for cotl
2969 (p15) fms.s1 C_hi = f0, f0, C_hi // to correct result's sign for cotl
2975 (p15) fms.s1 sgn_r = f0, f0, sgn_r // to correct result's sign for cotl
2982 // N even: V_lo = V_lo + V_hi * tanx
2983 // N odd: V_lo = V_lo - V_hi * tanx
2985 (p11) fnma.s1 V_lo = C_lo, V_hi, V_lo
2990 (p12) fnma.s1 V_lo = T_lo, V_hi, V_lo
2996 // N even: V_lo = V_lo - V_hi * C_lo
2997 // N odd: V_lo = V_lo - V_hi * T_lo
2999 fmpy.s1 V_lo = V_hi, V_lo
3005 // N even and odd: V_lo = V_lo * V_hi
3007 fadd.s1 tail = V_hi, V_lo
3013 // N even and odd: tail = V_hi + V_lo
3015 fma.s1 tail = tail, P, V_lo
3021 // N even: T_hi = sgn_r * T_hi
3022 // N odd : C_hi = -sgn_r * C_hi
3024 fma.s1 tail = tail, Sx, CORR
3030 // N even and odd: tail = Sx * tail + CORR
3032 fma.s1 tail = V_hi, Sx, tail
3038 // N even an odd: tail = Sx * V_hi + tail
3040 (p11) fma.s0 Result = sgn_r, tail, T_hi
3045 (p12) fma.s0 Result = sgn_r, tail, C_hi
3046 br.ret.sptk b0 ;; // Exit for 1/4 <= |r| < pi/4
3050 // Here if x denormal
3052 getf.exp GR_signexp_x = Norm_Arg // Get sign and exponent of x
3054 br.cond.sptk TANL_COMMON // Return to common code
3062 // Code for NaNs, Unsupporteds, Infs, or +/- zero ?
3063 // Invalid raised for Infs and SNaNs.
3068 fmerge.s f10 = f8, f8 // Save input for error call
3069 tbit.nz p6, p7 = cot_flag, 0 // p6=1 if we compute cotl
3075 (p6) fclass.m p6, p7 = f8, 0x7 // Test for zero (cotl only)
3080 .pred.rel "mutex", p6, p7
3082 (p6) mov GR_Parameter_Tag = 225 // (cotl)
3083 (p6) frcpa.s0 f8, p0 = f1, f8 // cotl(+-0) = +-Inf
3088 (p7) fmpy.s0 f8 = f8, f0
3093 GLOBAL_IEEE754_END(tanl)
3096 LOCAL_LIBM_ENTRY(__libm_error_region)
3101 add GR_Parameter_Y=-32,sp // Parameter 2 value
3103 .save ar.pfs,GR_SAVE_PFS
3104 mov GR_SAVE_PFS=ar.pfs // Save ar.pfs
3108 add sp=-64,sp // Create new stack
3110 mov GR_SAVE_GP=gp // Save gp
3115 stfe [GR_Parameter_Y] = f1,16 // STORE Parameter 2 on stack
3116 add GR_Parameter_X = 16,sp // Parameter 1 address
3117 .save b0, GR_SAVE_B0
3118 mov GR_SAVE_B0=b0 // Save b0
3124 stfe [GR_Parameter_X] = f10 // STORE Parameter 1 on stack
3125 add GR_Parameter_RESULT = 0,GR_Parameter_Y // Parameter 3 address
3129 stfe [GR_Parameter_Y] = f8 // STORE Parameter 3 on stack
3130 add GR_Parameter_Y = -16,GR_Parameter_Y
3131 br.call.sptk b0=__libm_error_support# // Call error handling function
3136 add GR_Parameter_RESULT = 48,sp
3141 ldfe f8 = [GR_Parameter_RESULT] // Get return result off stack
3143 add sp = 64,sp // Restore stack pointer
3144 mov b0 = GR_SAVE_B0 // Restore return address
3147 mov gp = GR_SAVE_GP // Restore gp
3148 mov ar.pfs = GR_SAVE_PFS // Restore ar.pfs
3149 br.ret.sptk b0 // Return
3152 LOCAL_LIBM_END(__libm_error_region)
3154 .type __libm_error_support#,@function
3155 .global __libm_error_support#
3158 // *******************************************************************
3159 // *******************************************************************
3160 // *******************************************************************
3162 // Special Code to handle very large argument case.
3163 // Call int __libm_pi_by_2_reduce(x,r,c) for |arguments| >= 2**63
3164 // The interface is custom:
3166 // (Arg or x) is in f8
3171 // We know also that __libm_pi_by_2_reduce preserves f10-15, f71-127. We
3172 // use this to eliminate save/restore of key fp registers in this calling
3175 // *******************************************************************
3176 // *******************************************************************
3177 // *******************************************************************
3179 LOCAL_LIBM_ENTRY(__libm_callout)
3183 add table_ptr2 = 144, table_base // Point to 2^-2
3185 .save ar.pfs,GR_SAVE_PFS
3186 mov GR_SAVE_PFS=ar.pfs // Save ar.pfs
3192 ldfps TWO_TO_NEG2, NEGTWO_TO_NEG2 = [table_ptr2]
3193 setf.sig B_mask1 = bmask1 // Form mask to get 5 msb of r
3194 .save b0, GR_SAVE_B0
3195 mov GR_SAVE_B0=b0 // Save b0
3200 // Call argument reduction with x in f8
3201 // Returns with N in r8, r in f8, c in f9
3202 // Assumes f71-127 are preserved across the call
3205 setf.sig B_mask2 = bmask2 // Form mask to form B from r
3206 mov GR_SAVE_GP=gp // Save gp
3207 br.call.sptk b0=__libm_pi_by_2_reduce#
3215 getf.sig sig_r = r // Extract significand of r
3216 fcmp.lt.s1 p6, p0 = r, TWO_TO_NEG2
3217 mov gp = GR_SAVE_GP // Restore gp
3222 getf.exp exp_r = r // Extract signexp of r
3224 mov b0 = GR_SAVE_B0 // Restore return address
3233 (p6) fcmp.gt.unc.s1 p6, p0 = r, NEGTWO_TO_NEG2
3234 mov ar.pfs = GR_SAVE_PFS // Restore pfs
3240 (p6) br.cond.spnt TANL_SMALL_R // Branch if |r| < 1/4
3241 br.cond.sptk TANL_NORMAL_R // Branch if 1/4 <= |r| < pi/4
3245 LOCAL_LIBM_END(__libm_callout)
3247 .type __libm_pi_by_2_reduce#,@function
3248 .global __libm_pi_by_2_reduce#