3 // Copyright (C) 2000, 2001, Intel Corporation
4 // All rights reserved.
7 // Redistribution and use in source and binary forms, with or without
8 // modification, are permitted provided that the following conditions are
11 // * Redistributions of source code must retain the above copyright
12 // notice, this list of conditions and the following disclaimer.
14 // * Redistributions in binary form must reproduce the above copyright
15 // notice, this list of conditions and the following disclaimer in the
16 // documentation and/or other materials provided with the distribution.
18 // * The name of Intel Corporation may not be used to endorse or promote
19 // products derived from this software without specific prior written
22 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
23 // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
24 // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
25 // A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL INTEL OR ITS
26 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
27 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
28 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
29 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY
30 // OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY OR TORT (INCLUDING
31 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
32 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
34 // Intel Corporation is the author of this code, and requests that all
35 // problem reports or change requests be submitted to it directly at
36 // http://developer.intel.com/opensource.
38 // *********************************************************************
41 // 02/02/00 Initial Version
42 // 4/04/00 Unwind support added
43 // 12/28/00 Fixed false invalid flags
45 // *********************************************************************
47 // Function: tan(x) = tangent(x), for double precision x values
49 // *********************************************************************
51 // Accuracy: Very accurate for double-precision values
53 // *********************************************************************
57 // Floating-Point Registers: f8 (Input and Return Value)
61 // General Purpose Registers:
63 // r49-r50 (Used to pass arguments to pi_by_2 reduce routine)
65 // Predicate Registers: p6-p15
67 // *********************************************************************
69 // IEEE Special Conditions:
71 // Denormal fault raised on denormal inputs
72 // Overflow exceptions do not occur
73 // Underflow exceptions raised when appropriate for tan
74 // (No specialized error handling for this routine)
75 // Inexact raised when appropriate by algorithm
82 // *********************************************************************
84 // Mathematical Description
86 // We consider the computation of FPTAN of Arg. Now, given
88 // Arg = N pi/2 + alpha, |alpha| <= pi/4,
90 // basic mathematical relationship shows that
92 // tan( Arg ) = tan( alpha ) if N is even;
93 // = -cot( alpha ) otherwise.
95 // The value of alpha is obtained by argument reduction and
96 // represented by two working precision numbers r and c where
98 // alpha = r + c accurately.
100 // The reduction method is described in a previous write up.
101 // The argument reduction scheme identifies 4 cases. For Cases 2
102 // and 4, because |alpha| is small, tan(r+c) and -cot(r+c) can be
103 // computed very easily by 2 or 3 terms of the Taylor series
104 // expansion as follows:
109 // tan(r + c) = r + c + r^3/3 ...accurately
110 // -cot(r + c) = -1/(r+c) + r/3 ...accurately
115 // tan(r + c) = r + c + r^3/3 + 2r^5/15 ...accurately
116 // -cot(r + c) = -1/(r+c) + r/3 + r^3/45 ...accurately
119 // The only cases left are Cases 1 and 3 of the argument reduction
120 // procedure. These two cases will be merged since after the
121 // argument is reduced in either cases, we have the reduced argument
122 // represented as r + c and that the magnitude |r + c| is not small
123 // enough to allow the usage of a very short approximation.
125 // The greatest challenge of this task is that the second terms of
126 // the Taylor series for tan(r) and -cot(r)
128 // r + r^3/3 + 2 r^5/15 + ...
132 // -1/r + r/3 + r^3/45 + ...
134 // are not very small when |r| is close to pi/4 and the rounding
135 // errors will be a concern if simple polynomial accumulation is
136 // used. When |r| < 2^(-2), however, the second terms will be small
137 // enough (5 bits or so of right shift) that a normal Horner
138 // recurrence suffices. Hence there are two cases that we consider
139 // in the accurate computation of tan(r) and cot(r), |r| <= pi/4.
141 // Case small_r: |r| < 2^(-2)
142 // --------------------------
144 // Since Arg = N pi/4 + r + c accurately, we have
146 // tan(Arg) = tan(r+c) for N even,
147 // = -cot(r+c) otherwise.
149 // Here for this case, both tan(r) and -cot(r) can be approximated
150 // by simple polynomials:
152 // tan(r) = r + P1_1 r^3 + P1_2 r^5 + ... + P1_9 r^19
153 // -cot(r) = -1/r + Q1_1 r + Q1_2 r^3 + ... + Q1_7 r^13
155 // accurately. Since |r| is relatively small, tan(r+c) and
156 // -cot(r+c) can be accurately approximated by replacing r with
157 // r+c only in the first two terms of the corresponding polynomials.
159 // Note that P1_1 (and Q1_1 for that matter) approximates 1/3 to
160 // almost 64 sig. bits, thus
162 // P1_1 (r+c)^3 = P1_1 r^3 + c * r^2 accurately.
166 // tan(r+c) = r + P1_1 r^3 + P1_2 r^5 + ... + P1_9 r^19
169 // -cot(r+c) = -1/(r+c) + Q1_1 r + Q1_2 r^3 + ... + Q1_7 r^13
173 // Case normal_r: 2^(-2) <= |r| <= pi/4
174 // ------------------------------------
176 // This case is more likely than the previous one if one considers
177 // r to be uniformly distributed in [-pi/4 pi/4].
179 // The required calculation is either
181 // tan(r + c) = tan(r) + correction, or
182 // -cot(r + c) = -cot(r) + correction.
186 // tan(r + c) = tan(r) + c tan'(r) + O(c^2)
187 // = tan(r) + c sec^2(r) + O(c^2)
188 // = tan(r) + c SEC_sq ...accurately
189 // as long as SEC_sq approximates sec^2(r)
190 // to, say, 5 bits or so.
194 // -cot(r + c) = -cot(r) - c cot'(r) + O(c^2)
195 // = -cot(r) + c csc^2(r) + O(c^2)
196 // = -cot(r) + c CSC_sq ...accurately
197 // as long as CSC_sq approximates csc^2(r)
198 // to, say, 5 bits or so.
200 // We therefore concentrate on accurately calculating tan(r) and
201 // cot(r) for a working-precision number r, |r| <= pi/4 to within
204 // We will employ a table-driven approach. Let
206 // r = sgn_r * 2^k * 1.b_1 b_2 ... b_5 ... b_63
207 // = sgn_r * ( B + x )
211 // B = 2^k * 1.b_1 b_2 ... b_5 1
216 // tan( B + x ) = ------------------------
220 // | tan(B) + tan(x) |
222 // = tan(B) + | ------------------------ - tan(B) |
223 // | 1 - tan(B)*tan(x) |
227 // = tan(B) + ------------------------
230 // (1/[sin(B)*cos(B)]) * tan(x)
231 // = tan(B) + --------------------------------
235 // Clearly, the values of tan(B), cot(B) and 1/(sin(B)*cos(B)) are
236 // calculated beforehand and stored in a table. Since
238 // |x| <= 2^k * 2^(-6) <= 2^(-7) (because k = -1, -2)
240 // a very short polynomial will be sufficient to approximate tan(x)
241 // accurately. The details involved in computing the last expression
242 // will be given in the next section on algorithm description.
245 // Now, we turn to the case where cot( B + x ) is needed.
249 // cot( B + x ) = ------------------------
253 // | 1 - tan(B)*tan(x) |
255 // = cot(B) + | ----------------------- - cot(B) |
256 // | tan(B) + tan(x) |
259 // [tan(B) + cot(B)] * tan(x)
260 // = cot(B) - ----------------------------
263 // (1/[sin(B)*cos(B)]) * tan(x)
264 // = cot(B) - --------------------------------
268 // Note that the values of tan(B), cot(B) and 1/(sin(B)*cos(B)) that
269 // are needed are the same set of values needed in the previous
272 // Finally, we can put all the ingredients together as follows:
274 // Arg = N * pi/2 + r + c ...accurately
276 // tan(Arg) = tan(r) + correction if N is even;
277 // = -cot(r) + correction otherwise.
279 // For Cases 2 and 4,
282 // tan(Arg) = tan(r + c) = r + c + r^3/3 N even
283 // = -cot(r + c) = -1/(r+c) + r/3 N odd
285 // tan(Arg) = tan(r + c) = r + c + r^3/3 + 2r^5/15 N even
286 // = -cot(r + c) = -1/(r+c) + r/3 + r^3/45 N odd
289 // For Cases 1 and 3,
291 // Case small_r: |r| < 2^(-2)
293 // tan(Arg) = r + P1_1 r^3 + P1_2 r^5 + ... + P1_9 r^19
294 // + c*(1 + r^2) N even
296 // = -1/(r+c) + Q1_1 r + Q1_2 r^3 + ... + Q1_7 r^13
299 // Case normal_r: 2^(-2) <= |r| <= pi/4
301 // tan(Arg) = tan(r) + c * sec^2(r) N even
302 // = -cot(r) + c * csc^2(r) otherwise
306 // tan(Arg) = tan(r) + c*sec^2(r)
307 // = tan( sgn_r * (B+x) ) + c * sec^2(|r|)
308 // = sgn_r * ( tan(B+x) + sgn_r*c*sec^2(|r|) )
309 // = sgn_r * ( tan(B+x) + sgn_r*c*sec^2(B) )
311 // since B approximates |r| to 2^(-6) in relative accuracy.
313 // / (1/[sin(B)*cos(B)]) * tan(x)
314 // tan(Arg) = sgn_r * | tan(B) + --------------------------------
322 // CORR = sgn_r*c*tan(B)*SC_inv(B); SC_inv(B) = 1/(sin(B)*cos(B)).
326 // tan(Arg) = -cot(r) + c*csc^2(r)
327 // = -cot( sgn_r * (B+x) ) + c * csc^2(|r|)
328 // = sgn_r * ( -cot(B+x) + sgn_r*c*csc^2(|r|) )
329 // = sgn_r * ( -cot(B+x) + sgn_r*c*csc^2(B) )
331 // since B approximates |r| to 2^(-6) in relative accuracy.
333 // / (1/[sin(B)*cos(B)]) * tan(x)
334 // tan(Arg) = sgn_r * | -cot(B) + --------------------------------
342 // CORR = sgn_r*c*cot(B)*SC_inv(B); SC_inv(B) = 1/(sin(B)*cos(B)).
345 // The actual algorithm prescribes how all the mathematical formulas
349 // 2. Algorithmic Description
350 // ==========================
352 // 2.1 Computation for Cases 2 and 4.
353 // ----------------------------------
355 // For Case 2, we use two-term polynomials.
360 // Result := c + r * rsq * P1_1
361 // Result := r + Result ...in user-defined rounding
364 // S_hi := -frcpa(r) ...8 bits
365 // S_hi := S_hi + S_hi*(1 + S_hi*r) ...16 bits
366 // S_hi := S_hi + S_hi*(1 + S_hi*r) ...32 bits
367 // S_hi := S_hi + S_hi*(1 + S_hi*r) ...64 bits
368 // S_lo := S_hi*( (1 + S_hi*r) + S_hi*c )
369 // ...S_hi + S_lo is -1/(r+c) to extra precision
370 // S_lo := S_lo + Q1_1*r
372 // Result := S_hi + S_lo ...in user-defined rounding
374 // For Case 4, we use three-term polynomials
379 // Result := c + r * rsq * (P1_1 + rsq * P1_2)
380 // Result := r + Result ...in user-defined rounding
383 // S_hi := -frcpa(r) ...8 bits
384 // S_hi := S_hi + S_hi*(1 + S_hi*r) ...16 bits
385 // S_hi := S_hi + S_hi*(1 + S_hi*r) ...32 bits
386 // S_hi := S_hi + S_hi*(1 + S_hi*r) ...64 bits
387 // S_lo := S_hi*( (1 + S_hi*r) + S_hi*c )
388 // ...S_hi + S_lo is -1/(r+c) to extra precision
390 // P := Q1_1 + rsq*Q1_2
391 // S_lo := S_lo + r*P
393 // Result := S_hi + S_lo ...in user-defined rounding
396 // Note that the coefficients P1_1, P1_2, Q1_1, and Q1_2 are
397 // the same as those used in the small_r case of Cases 1 and 3
401 // 2.2 Computation for Cases 1 and 3.
402 // ----------------------------------
403 // This is further divided into the case of small_r,
404 // where |r| < 2^(-2), and the case of normal_r, where |r| lies between
407 // Algorithm for the case of small_r
408 // ---------------------------------
412 // Poly1 := rsq*(P1_1 + rsq*(P1_2 + rsq*P1_3))
413 // r_to_the_8 := rsq * rsq
414 // r_to_the_8 := r_to_the_8 * r_to_the_8
415 // Poly2 := P1_4 + rsq*(P1_5 + rsq*(P1_6 + ... rsq*P1_9))
416 // CORR := c * ( 1 + rsq )
417 // Poly := Poly1 + r_to_the_8*Poly2
418 // Result := r*Poly + CORR
419 // Result := r + Result ...in user-defined rounding
420 // ...note that Poly1 and r_to_the_8 can be computed in parallel
421 // ...with Poly2 (Poly1 is intentionally set to be much
422 // ...shorter than Poly2 so that r_to_the_8 and CORR can be hidden)
425 // S_hi := -frcpa(r) ...8 bits
426 // S_hi := S_hi + S_hi*(1 + S_hi*r) ...16 bits
427 // S_hi := S_hi + S_hi*(1 + S_hi*r) ...32 bits
428 // S_hi := S_hi + S_hi*(1 + S_hi*r) ...64 bits
429 // S_lo := S_hi*( (1 + S_hi*r) + S_hi*c )
430 // ...S_hi + S_lo is -1/(r+c) to extra precision
431 // S_lo := S_lo + Q1_1*c
433 // ...S_hi and S_lo are computed in parallel with
436 // P := Q1_1 + rsq*(Q1_2 + rsq*(Q1_3 + ... + rsq*Q1_7))
438 // Result := r*P + S_lo
439 // Result := S_hi + Result ...in user-defined rounding
442 // Algorithm for the case of normal_r
443 // ----------------------------------
445 // Here, we first consider the computation of tan( r + c ). As
446 // presented in the previous section,
448 // tan( r + c ) = tan(r) + c * sec^2(r)
449 // = sgn_r * [ tan(B+x) + CORR ]
450 // CORR = sgn_r * c * tan(B) * 1/[sin(B)*cos(B)]
452 // because sec^2(r) = sec^(|r|), and B approximate |r| to 6.5 bits.
455 // / (1/[sin(B)*cos(B)]) * tan(x)
456 // sgn_r * | tan(B) + -------------------------------- +
463 // The values of tan(B), cot(B) and 1/(sin(B)*cos(B)) are
464 // calculated beforehand and stored in a table. Specifically,
465 // the table values are
467 // tan(B) as T_hi + T_lo;
468 // cot(B) as C_hi + C_lo;
469 // 1/[sin(B)*cos(B)] as SC_inv
471 // T_hi, C_hi are in double-precision memory format;
472 // T_lo, C_lo are in single-precision memory format;
473 // SC_inv is in extended-precision memory format.
475 // The value of tan(x) will be approximated by a short polynomial of
478 // tan(x) as x + x * P, where
479 // P = x^2 * (P2_1 + x^2 * (P2_2 + x^2 * P2_3))
481 // Because |x| <= 2^(-7), cot(B) - x approximates cot(B) - tan(x)
482 // to a relative accuracy better than 2^(-20). Thus, a good
483 // initial guess of 1/( cot(B) - tan(x) ) to initiate the iterative
486 // 1/(cot(B) - tan(x)) is approximately
488 // tan(B)/(1 - x*tan(B)) is approximately
489 // T_hi / ( 1 - T_hi * x ) is approximately
491 // T_hi * [ 1 + (Thi * x) + (T_hi * x)^2 ]
493 // The calculation of tan(r+c) therefore proceed as follows:
498 // V_hi := T_hi*(1 + Tx*(1 + Tx))
499 // P := xsq * (P1_1 + xsq*(P1_2 + xsq*P1_3))
500 // ...V_hi serves as an initial guess of 1/(cot(B) - tan(x))
501 // ...good to about 20 bits of accuracy
505 // ...D is a double precision denominator: cot(B) - tan(x)
507 // V_hi := V_hi + V_hi*(1 - V_hi*D)
508 // ....V_hi approximates 1/(cot(B)-tan(x)) to 40 bits
510 // V_lo := V_hi * ( [ (1 - V_hi*C_hi) + V_hi*tanx ]
511 // - V_hi*C_lo ) ...observe all order
512 // ...V_hi + V_lo approximates 1/(cot(B) - tan(x))
513 // ...to extra accuracy
515 // ... SC_inv(B) * (x + x*P)
516 // ... tan(B) + ------------------------- + CORR
517 // ... cot(B) - (x + x*P)
519 // ... = tan(B) + SC_inv(B)*(x + x*P)*(V_hi + V_lo) + CORR
523 // CORR := sgn_r * c * SC_inv * T_hi
525 // ...put the ingredients together to compute
526 // ... SC_inv(B) * (x + x*P)
527 // ... tan(B) + ------------------------- + CORR
528 // ... cot(B) - (x + x*P)
530 // ... = tan(B) + SC_inv(B)*(x + x*P)*(V_hi + V_lo) + CORR
532 // ... = T_hi + T_lo + CORR +
533 // ... Sx * V_hi + Sx * V_lo + Sx * P *(V_hi + V_lo)
535 // CORR := CORR + T_lo
536 // tail := V_lo + P*(V_hi + V_lo)
537 // tail := Sx * tail + CORR
538 // tail := Sx * V_hi + tail
539 // T_hi := sgn_r * T_hi
541 // ...T_hi + sgn_r*tail now approximate
542 // ...sgn_r*(tan(B+x) + CORR) accurately
544 // Result := T_hi + sgn_r*tail ...in user-defined
545 // ...rounding control
546 // ...It is crucial that independent paths be fully
547 // ...exploited for performance's sake.
550 // Next, we consider the computation of -cot( r + c ). As
551 // presented in the previous section,
553 // -cot( r + c ) = -cot(r) + c * csc^2(r)
554 // = sgn_r * [ -cot(B+x) + CORR ]
555 // CORR = sgn_r * c * cot(B) * 1/[sin(B)*cos(B)]
557 // because csc^2(r) = csc^(|r|), and B approximate |r| to 6.5 bits.
560 // / (1/[sin(B)*cos(B)]) * tan(x)
561 // sgn_r * | -cot(B) + -------------------------------- +
568 // The values of tan(B), cot(B) and 1/(sin(B)*cos(B)) are
569 // calculated beforehand and stored in a table. Specifically,
570 // the table values are
572 // tan(B) as T_hi + T_lo;
573 // cot(B) as C_hi + C_lo;
574 // 1/[sin(B)*cos(B)] as SC_inv
576 // T_hi, C_hi are in double-precision memory format;
577 // T_lo, C_lo are in single-precision memory format;
578 // SC_inv is in extended-precision memory format.
580 // The value of tan(x) will be approximated by a short polynomial of
583 // tan(x) as x + x * P, where
584 // P = x^2 * (P2_1 + x^2 * (P2_2 + x^2 * P2_3))
586 // Because |x| <= 2^(-7), tan(B) + x approximates tan(B) + tan(x)
587 // to a relative accuracy better than 2^(-18). Thus, a good
588 // initial guess of 1/( tan(B) + tan(x) ) to initiate the iterative
591 // 1/(tan(B) + tan(x)) is approximately
593 // cot(B)/(1 + x*cot(B)) is approximately
594 // C_hi / ( 1 + C_hi * x ) is approximately
596 // C_hi * [ 1 - (C_hi * x) + (C_hi * x)^2 ]
598 // The calculation of -cot(r+c) therefore proceed as follows:
603 // V_hi := C_hi*(1 - Cx*(1 - Cx))
604 // P := xsq * (P1_1 + xsq*(P1_2 + xsq*P1_3))
605 // ...V_hi serves as an initial guess of 1/(tan(B) + tan(x))
606 // ...good to about 18 bits of accuracy
610 // ...D is a double precision denominator: tan(B) + tan(x)
612 // V_hi := V_hi + V_hi*(1 - V_hi*D)
613 // ....V_hi approximates 1/(tan(B)+tan(x)) to 40 bits
615 // V_lo := V_hi * ( [ (1 - V_hi*T_hi) - V_hi*tanx ]
616 // - V_hi*T_lo ) ...observe all order
617 // ...V_hi + V_lo approximates 1/(tan(B) + tan(x))
618 // ...to extra accuracy
620 // ... SC_inv(B) * (x + x*P)
621 // ... -cot(B) + ------------------------- + CORR
622 // ... tan(B) + (x + x*P)
624 // ... =-cot(B) + SC_inv(B)*(x + x*P)*(V_hi + V_lo) + CORR
628 // CORR := sgn_r * c * SC_inv * C_hi
630 // ...put the ingredients together to compute
631 // ... SC_inv(B) * (x + x*P)
632 // ... -cot(B) + ------------------------- + CORR
633 // ... tan(B) + (x + x*P)
635 // ... =-cot(B) + SC_inv(B)*(x + x*P)*(V_hi + V_lo) + CORR
637 // ... =-C_hi - C_lo + CORR +
638 // ... Sx * V_hi + Sx * V_lo + Sx * P *(V_hi + V_lo)
640 // CORR := CORR - C_lo
641 // tail := V_lo + P*(V_hi + V_lo)
642 // tail := Sx * tail + CORR
643 // tail := Sx * V_hi + tail
644 // C_hi := -sgn_r * C_hi
646 // ...C_hi + sgn_r*tail now approximates
647 // ...sgn_r*(-cot(B+x) + CORR) accurately
649 // Result := C_hi + sgn_r*tail in user-defined rounding control
650 // ...It is crucial that independent paths be fully
651 // ...exploited for performance's sake.
653 // 3. Implementation Notes
654 // =======================
656 // Table entries T_hi, T_lo; C_hi, C_lo; SC_inv
658 // Recall that 2^(-2) <= |r| <= pi/4;
660 // r = sgn_r * 2^k * 1.b_1 b_2 ... b_63
664 // B = 2^k * 1.b_1 b_2 b_3 b_4 b_5 1
666 // Thus, for k = -2, possible values of B are
668 // B = 2^(-2) * ( 1 + index/32 + 1/64 ),
669 // index ranges from 0 to 31
671 // For k = -1, however, since |r| <= pi/4 = 0.78...
672 // possible values of B are
674 // B = 2^(-1) * ( 1 + index/32 + 1/64 )
675 // index ranges from 0 to 19.
679 #include "libm_support.h"
690 .type TAN_BASE_CONSTANTS, @object
691 data4 0x4B800000, 0xCB800000, 0x38800000, 0xB8800000 // two**24, -two**24
692 // two**-14, -two**-14
693 data4 0x4E44152A, 0xA2F9836E, 0x00003FFE, 0x00000000 // two_by_pi
694 data4 0xCE81B9F1, 0xC84D32B0, 0x00004016, 0x00000000 // P_0
695 data4 0x2168C235, 0xC90FDAA2, 0x00003FFF, 0x00000000 // P_1
696 data4 0xFC8F8CBB, 0xECE675D1, 0x0000BFBD, 0x00000000 // P_2
697 data4 0xACC19C60, 0xB7ED8FBB, 0x0000BF7C, 0x00000000 // P_3
698 data4 0x5F000000, 0xDF000000, 0x00000000, 0x00000000 // two_to_63, -two_to_63
699 data4 0x6EC6B45A, 0xA397E504, 0x00003FE7, 0x00000000 // Inv_P_0
700 data4 0xDBD171A1, 0x8D848E89, 0x0000BFBF, 0x00000000 // d_1
701 data4 0x18A66F8E, 0xD5394C36, 0x0000BF7C, 0x00000000 // d_2
702 data4 0x2168C234, 0xC90FDAA2, 0x00003FFE, 0x00000000 // PI_BY_4
703 data4 0x2168C234, 0xC90FDAA2, 0x0000BFFE, 0x00000000 // MPI_BY_4
704 data4 0x3E800000, 0xBE800000, 0x00000000, 0x00000000 // two**-2, -two**-2
705 data4 0x2F000000, 0xAF000000, 0x00000000, 0x00000000 // two**-33, -two**-33
706 data4 0xAAAAAABD, 0xAAAAAAAA, 0x00003FFD, 0x00000000 // P1_1
707 data4 0x88882E6A, 0x88888888, 0x00003FFC, 0x00000000 // P1_2
708 data4 0x0F0177B6, 0xDD0DD0DD, 0x00003FFA, 0x00000000 // P1_3
709 data4 0x646B8C6D, 0xB327A440, 0x00003FF9, 0x00000000 // P1_4
710 data4 0x1D5F7D20, 0x91371B25, 0x00003FF8, 0x00000000 // P1_5
711 data4 0x61C67914, 0xEB69A5F1, 0x00003FF6, 0x00000000 // P1_6
712 data4 0x019318D2, 0xBEDD37BE, 0x00003FF5, 0x00000000 // P1_7
713 data4 0x3C794015, 0x9979B146, 0x00003FF4, 0x00000000 // P1_8
714 data4 0x8C6EB58A, 0x8EBD21A3, 0x00003FF3, 0x00000000 // P1_9
715 data4 0xAAAAAAB4, 0xAAAAAAAA, 0x00003FFD, 0x00000000 // Q1_1
716 data4 0x0B5FC93E, 0xB60B60B6, 0x00003FF9, 0x00000000 // Q1_2
717 data4 0x0C9BBFBF, 0x8AB355E0, 0x00003FF6, 0x00000000 // Q1_3
718 data4 0xCBEE3D4C, 0xDDEBBC89, 0x00003FF2, 0x00000000 // Q1_4
719 data4 0x5F80BBB6, 0xB3548A68, 0x00003FEF, 0x00000000 // Q1_5
720 data4 0x4CED5BF1, 0x91362560, 0x00003FEC, 0x00000000 // Q1_6
721 data4 0x8EE92A83, 0xF189D95A, 0x00003FE8, 0x00000000 // Q1_7
722 data4 0xAAAB362F, 0xAAAAAAAA, 0x00003FFD, 0x00000000 // P2_1
723 data4 0xE97A6097, 0x88888886, 0x00003FFC, 0x00000000 // P2_2
724 data4 0x25E716A1, 0xDD108EE0, 0x00003FFA, 0x00000000 // P2_3
726 // Entries T_hi double-precision memory format
727 // Index = 0,1,...,31 B = 2^(-2)*(1+Index/32+1/64)
728 // Entries T_lo single-precision memory format
729 // Index = 0,1,...,31 B = 2^(-2)*(1+Index/32+1/64)
731 data4 0x62400794, 0x3FD09BC3, 0x23A05C32, 0x00000000
732 data4 0xDFFBC074, 0x3FD124A9, 0x240078B2, 0x00000000
733 data4 0x5BD4920F, 0x3FD1AE23, 0x23826B8E, 0x00000000
734 data4 0x15E2701D, 0x3FD23835, 0x22D31154, 0x00000000
735 data4 0x63739C2D, 0x3FD2C2E4, 0x2265C9E2, 0x00000000
736 data4 0xAFEEA48B, 0x3FD34E36, 0x245C05EB, 0x00000000
737 data4 0x7DBB35D1, 0x3FD3DA31, 0x24749F2D, 0x00000000
738 data4 0x67321619, 0x3FD466DA, 0x2462CECE, 0x00000000
739 data4 0x1F94A4D5, 0x3FD4F437, 0x246D0DF1, 0x00000000
740 data4 0x740C3E6D, 0x3FD5824D, 0x240A85B5, 0x00000000
741 data4 0x4CB1E73D, 0x3FD61123, 0x23F96E33, 0x00000000
742 data4 0xAD9EA64B, 0x3FD6A0BE, 0x247C5393, 0x00000000
743 data4 0xB804FD01, 0x3FD73125, 0x241F3B29, 0x00000000
744 data4 0xAB53EE83, 0x3FD7C25E, 0x2479989B, 0x00000000
745 data4 0xE6640EED, 0x3FD8546F, 0x23B343BC, 0x00000000
746 data4 0xE8AF1892, 0x3FD8E75F, 0x241454D1, 0x00000000
747 data4 0x53928BDA, 0x3FD97B35, 0x238613D9, 0x00000000
748 data4 0xEB9DE4DE, 0x3FDA0FF6, 0x22859FA7, 0x00000000
749 data4 0x99ECF92D, 0x3FDAA5AB, 0x237A6D06, 0x00000000
750 data4 0x6D8F1796, 0x3FDB3C5A, 0x23952F6C, 0x00000000
751 data4 0x9CFB8BE4, 0x3FDBD40A, 0x2280FC95, 0x00000000
752 data4 0x87943100, 0x3FDC6CC3, 0x245D2EC0, 0x00000000
753 data4 0xB736C500, 0x3FDD068C, 0x23C4AD7D, 0x00000000
754 data4 0xE1DDBC31, 0x3FDDA16D, 0x23D076E6, 0x00000000
755 data4 0xEB515A93, 0x3FDE3D6E, 0x244809A6, 0x00000000
756 data4 0xE6E9E5F1, 0x3FDEDA97, 0x220856C8, 0x00000000
757 data4 0x1963CE69, 0x3FDF78F1, 0x244BE993, 0x00000000
758 data4 0x7D635BCE, 0x3FE00C41, 0x23D21799, 0x00000000
759 data4 0x1C302CD3, 0x3FE05CAB, 0x248A1B1D, 0x00000000
760 data4 0xDB6A1FA0, 0x3FE0ADB9, 0x23D53E33, 0x00000000
761 data4 0x4A20BA81, 0x3FE0FF72, 0x24DB9ED5, 0x00000000
762 data4 0x153FA6F5, 0x3FE151D9, 0x24E9E451, 0x00000000
764 // Entries T_hi double-precision memory format
765 // Index = 0,1,...,19 B = 2^(-1)*(1+Index/32+1/64)
766 // Entries T_lo single-precision memory format
767 // Index = 0,1,...,19 B = 2^(-1)*(1+Index/32+1/64)
769 data4 0xBA1BE39E, 0x3FE1CEC4, 0x24B60F9E, 0x00000000
770 data4 0x5ABD9B2D, 0x3FE277E4, 0x248C2474, 0x00000000
771 data4 0x0272B110, 0x3FE32418, 0x247B8311, 0x00000000
772 data4 0x890E2DF0, 0x3FE3D38B, 0x24C55751, 0x00000000
773 data4 0x46236871, 0x3FE4866D, 0x24E5BC34, 0x00000000
774 data4 0x45E044B0, 0x3FE53CEE, 0x24001BA4, 0x00000000
775 data4 0x82EC06E4, 0x3FE5F742, 0x24B973DC, 0x00000000
776 data4 0x25DF43F9, 0x3FE6B5A1, 0x24895440, 0x00000000
777 data4 0xCAFD348C, 0x3FE77844, 0x240021CA, 0x00000000
778 data4 0xCEED6B92, 0x3FE83F6B, 0x24C45372, 0x00000000
779 data4 0xA34F3665, 0x3FE90B58, 0x240DAD33, 0x00000000
780 data4 0x2C1E56B4, 0x3FE9DC52, 0x24F846CE, 0x00000000
781 data4 0x27041578, 0x3FEAB2A4, 0x2323FB6E, 0x00000000
782 data4 0x9DD8C373, 0x3FEB8E9F, 0x24B3090B, 0x00000000
783 data4 0x65C9AA7B, 0x3FEC709B, 0x2449F611, 0x00000000
784 data4 0xACCF8435, 0x3FED58F4, 0x23616A7E, 0x00000000
785 data4 0x97635082, 0x3FEE480F, 0x24C2FEAE, 0x00000000
786 data4 0xF0ACC544, 0x3FEF3E57, 0x242CE964, 0x00000000
787 data4 0xF7E06E4B, 0x3FF01E20, 0x2480D3EE, 0x00000000
788 data4 0x8A798A69, 0x3FF0A125, 0x24DB8967, 0x00000000
790 // Entries C_hi double-precision memory format
791 // Index = 0,1,...,31 B = 2^(-2)*(1+Index/32+1/64)
792 // Entries C_lo single-precision memory format
793 // Index = 0,1,...,31 B = 2^(-2)*(1+Index/32+1/64)
795 data4 0xE63EFBD0, 0x400ED3E2, 0x259D94D4, 0x00000000
796 data4 0xC515DAB5, 0x400DDDB4, 0x245F0537, 0x00000000
797 data4 0xBE19A79F, 0x400CF57A, 0x25D4EA9F, 0x00000000
798 data4 0xD15298ED, 0x400C1A06, 0x24AE40A0, 0x00000000
799 data4 0x164B2708, 0x400B4A4C, 0x25A5AAB6, 0x00000000
800 data4 0x5285B068, 0x400A855A, 0x25524F18, 0x00000000
801 data4 0x3FFA549F, 0x4009CA5A, 0x24C999C0, 0x00000000
802 data4 0x646AF623, 0x4009188A, 0x254FD801, 0x00000000
803 data4 0x6084D0E7, 0x40086F3C, 0x2560F5FD, 0x00000000
804 data4 0xA29A76EE, 0x4007CDD2, 0x255B9D19, 0x00000000
805 data4 0x6C8ECA95, 0x400733BE, 0x25CB021B, 0x00000000
806 data4 0x1F8DDC52, 0x4006A07E, 0x24AB4722, 0x00000000
807 data4 0xC298AD58, 0x4006139B, 0x252764E2, 0x00000000
808 data4 0xBAD7164B, 0x40058CAB, 0x24DAF5DB, 0x00000000
809 data4 0xAE31A5D3, 0x40050B4B, 0x25EA20F4, 0x00000000
810 data4 0x89F85A8A, 0x40048F21, 0x2583A3E8, 0x00000000
811 data4 0xA862380D, 0x400417DA, 0x25DCC4CC, 0x00000000
812 data4 0x1088FCFE, 0x4003A52B, 0x2430A492, 0x00000000
813 data4 0xCD3527D5, 0x400336CC, 0x255F77CF, 0x00000000
814 data4 0x5760766D, 0x4002CC7F, 0x25DA0BDA, 0x00000000
815 data4 0x11CE02E3, 0x40026607, 0x256FF4A2, 0x00000000
816 data4 0xD37BBE04, 0x4002032C, 0x25208AED, 0x00000000
817 data4 0x7F050775, 0x4001A3BD, 0x24B72DD6, 0x00000000
818 data4 0xA554848A, 0x40014789, 0x24AB4DAA, 0x00000000
819 data4 0x323E81B7, 0x4000EE65, 0x2584C440, 0x00000000
820 data4 0x21CF1293, 0x40009827, 0x25C9428D, 0x00000000
821 data4 0x3D415EEB, 0x400044A9, 0x25DC8482, 0x00000000
822 data4 0xBD72C577, 0x3FFFE78F, 0x257F5070, 0x00000000
823 data4 0x75EFD28E, 0x3FFF4AC3, 0x23EBBF7A, 0x00000000
824 data4 0x60B52DDE, 0x3FFEB2AF, 0x22EECA07, 0x00000000
825 data4 0x35204180, 0x3FFE1F19, 0x24191079, 0x00000000
826 data4 0x54F7E60A, 0x3FFD8FCA, 0x248D3058, 0x00000000
828 // Entries C_hi double-precision memory format
829 // Index = 0,1,...,19 B = 2^(-1)*(1+Index/32+1/64)
830 // Entries C_lo single-precision memory format
831 // Index = 0,1,...,19 B = 2^(-1)*(1+Index/32+1/64)
833 data4 0x79F6FADE, 0x3FFCC06A, 0x239C7886, 0x00000000
834 data4 0x891662A6, 0x3FFBB91F, 0x250BD191, 0x00000000
835 data4 0x529F155D, 0x3FFABFB6, 0x256CC3E6, 0x00000000
836 data4 0x2E964AE9, 0x3FF9D300, 0x250843E3, 0x00000000
837 data4 0x89DCB383, 0x3FF8F1EF, 0x2277C87E, 0x00000000
838 data4 0x7C87DBD6, 0x3FF81B93, 0x256DA6CF, 0x00000000
839 data4 0x1042EDE4, 0x3FF74F14, 0x2573D28A, 0x00000000
840 data4 0x1784B360, 0x3FF68BAF, 0x242E489A, 0x00000000
841 data4 0x7C923C4C, 0x3FF5D0B5, 0x2532D940, 0x00000000
842 data4 0xF418EF20, 0x3FF51D88, 0x253C7DD6, 0x00000000
843 data4 0x02F88DAE, 0x3FF4719A, 0x23DB59BF, 0x00000000
844 data4 0x49DA0788, 0x3FF3CC66, 0x252B4756, 0x00000000
845 data4 0x0B980DB8, 0x3FF32D77, 0x23FE585F, 0x00000000
846 data4 0xE56C987A, 0x3FF2945F, 0x25378A63, 0x00000000
847 data4 0xB16523F6, 0x3FF200BD, 0x247BB2E0, 0x00000000
848 data4 0x8CE27778, 0x3FF17235, 0x24446538, 0x00000000
849 data4 0xFDEFE692, 0x3FF0E873, 0x2514638F, 0x00000000
850 data4 0x33154062, 0x3FF0632C, 0x24A7FC27, 0x00000000
851 data4 0xB3EF115F, 0x3FEFC42E, 0x248FD0FE, 0x00000000
852 data4 0x135D26F6, 0x3FEEC9E8, 0x2385C719, 0x00000000
854 // Entries SC_inv in Swapped IEEE format (extended)
855 // Index = 0,1,...,31 B = 2^(-2)*(1+Index/32+1/64)
857 data4 0x1BF30C9E, 0x839D6D4A, 0x00004001, 0x00000000
858 data4 0x554B0EB0, 0x80092804, 0x00004001, 0x00000000
859 data4 0xA1CF0DE9, 0xF959F94C, 0x00004000, 0x00000000
860 data4 0x77378677, 0xF3086BA0, 0x00004000, 0x00000000
861 data4 0xCCD4723C, 0xED154515, 0x00004000, 0x00000000
862 data4 0x1C27CF25, 0xE7790944, 0x00004000, 0x00000000
863 data4 0x8DDACB88, 0xE22D037D, 0x00004000, 0x00000000
864 data4 0x89C73522, 0xDD2B2D8A, 0x00004000, 0x00000000
865 data4 0xBB2C1171, 0xD86E1A23, 0x00004000, 0x00000000
866 data4 0xDFF5E0F9, 0xD3F0E288, 0x00004000, 0x00000000
867 data4 0x283BEBD5, 0xCFAF16B1, 0x00004000, 0x00000000
868 data4 0x0D88DD53, 0xCBA4AFAA, 0x00004000, 0x00000000
869 data4 0xCA67C43D, 0xC7CE03CC, 0x00004000, 0x00000000
870 data4 0x0CA0DDB0, 0xC427BC82, 0x00004000, 0x00000000
871 data4 0xF13D8CAB, 0xC0AECD57, 0x00004000, 0x00000000
872 data4 0x71ECE6B1, 0xBD606C38, 0x00004000, 0x00000000
873 data4 0xA44C4929, 0xBA3A0A96, 0x00004000, 0x00000000
874 data4 0xE5CCCEC1, 0xB7394F6F, 0x00004000, 0x00000000
875 data4 0x9637D8BC, 0xB45C1203, 0x00004000, 0x00000000
876 data4 0x92CB051B, 0xB1A05528, 0x00004000, 0x00000000
877 data4 0x6BA2FFD0, 0xAF04432B, 0x00004000, 0x00000000
878 data4 0x7221235F, 0xAC862A23, 0x00004000, 0x00000000
879 data4 0x5F00A9D1, 0xAA2478AF, 0x00004000, 0x00000000
880 data4 0x81E082BF, 0xA7DDBB0C, 0x00004000, 0x00000000
881 data4 0x45684FEE, 0xA5B0987D, 0x00004000, 0x00000000
882 data4 0x627A8F53, 0xA39BD0F5, 0x00004000, 0x00000000
883 data4 0x6EC5C8B0, 0xA19E3B03, 0x00004000, 0x00000000
884 data4 0x91CD7C66, 0x9FB6C1F0, 0x00004000, 0x00000000
885 data4 0x1FA3DF8A, 0x9DE46410, 0x00004000, 0x00000000
886 data4 0xA8F6B888, 0x9C263139, 0x00004000, 0x00000000
887 data4 0xC27B0450, 0x9A7B4968, 0x00004000, 0x00000000
888 data4 0x5EE614EE, 0x98E2DB7E, 0x00004000, 0x00000000
890 // Entries SC_inv in Swapped IEEE format (extended)
891 // Index = 0,1,...,19 B = 2^(-1)*(1+Index/32+1/64)
893 data4 0x13B2B5BA, 0x969F335C, 0x00004000, 0x00000000
894 data4 0xD4C0F548, 0x93D446D9, 0x00004000, 0x00000000
895 data4 0x61B798AF, 0x9147094F, 0x00004000, 0x00000000
896 data4 0x758787AC, 0x8EF317CC, 0x00004000, 0x00000000
897 data4 0xB99EEFDB, 0x8CD498B3, 0x00004000, 0x00000000
898 data4 0xDFF8BC37, 0x8AE82A7D, 0x00004000, 0x00000000
899 data4 0xE3C55D42, 0x892AD546, 0x00004000, 0x00000000
900 data4 0xD15573C1, 0x8799FEA9, 0x00004000, 0x00000000
901 data4 0x435A4B4C, 0x86335F88, 0x00004000, 0x00000000
902 data4 0x3E93A87B, 0x84F4FB6E, 0x00004000, 0x00000000
903 data4 0x80A382FB, 0x83DD1952, 0x00004000, 0x00000000
904 data4 0xA4CB8C9E, 0x82EA3D7F, 0x00004000, 0x00000000
905 data4 0x6861D0A8, 0x821B247C, 0x00004000, 0x00000000
906 data4 0x63E8D244, 0x816EBED1, 0x00004000, 0x00000000
907 data4 0x27E4CFC6, 0x80E42D91, 0x00004000, 0x00000000
908 data4 0x28E64AFD, 0x807ABF8D, 0x00004000, 0x00000000
909 data4 0x863B4FD8, 0x8031EF26, 0x00004000, 0x00000000
910 data4 0xAE8C11FD, 0x800960AD, 0x00004000, 0x00000000
911 data4 0x5FDBEC21, 0x8000E147, 0x00004000, 0x00000000
912 data4 0xA07791FA, 0x80186650, 0x00004000, 0x00000000
976 NEGTWO_TO_NEG14 = f76
977 NEGTWO_TO_NEG33 = f77
1033 GR_Parameter_X = r49
1034 GR_Parameter_r = r50
1046 alloc r32 = ar.pfs, 0,17,2,0
1047 (p0) fclass.m.unc p6,p0 = Arg, 0x1E7
1054 (p0) fclass.nm.unc p7,p0 = Arg, 0x1FF
1060 (p0) addl table_ptr1 = @ltoff(TAN_BASE_CONSTANTS), gp
1067 ld8 table_ptr1 = [table_ptr1]
1068 setf.sig fp_tmp = gr_tmp // Make a constant so fmpy produces inexact
1074 // Check for NatVals, Infs , NaNs, and Zeros
1075 // Check for everything - if false, then must be pseudo-zero
1077 // Local table pointer
1081 (p0) add table_ptr2 = 96, table_ptr1
1082 (p6) br.cond.spnt __libm_TAN_SPECIAL
1083 (p7) br.cond.spnt __libm_TAN_SPECIAL ;;
1087 // Branch out to deal with unsupporteds and special values.
1091 (p0) ldfs TWO_TO_24 = [table_ptr1],4
1092 (p0) ldfs TWO_TO_63 = [table_ptr2],4
1094 // Load -2**24, load -2**63.
1096 (p0) fcmp.eq.s0 p0, p6 = Arg, f1 ;;
1100 (p0) ldfs NEGTWO_TO_63 = [table_ptr2],12
1101 (p0) fnorm.s1 Arg = Arg
1105 // Load 2**24, Load 2**63.
1109 (p0) ldfs NEGTWO_TO_24 = [table_ptr1],12 ;;
1111 // Do fcmp to generate Denormal exception
1112 // - can't do FNORM (will generate Underflow when U is unmasked!)
1113 // Normalize input argument.
1115 (p0) ldfe two_by_PI = [table_ptr1],16
1120 (p0) ldfe Inv_P_0 = [table_ptr2],16 ;;
1121 (p0) ldfe d_1 = [table_ptr2],16
1125 // Decide about the paths to take:
1126 // PR_1 and PR_3 set if -2**24 < Arg < 2**24 - CASE 1 OR 2
1127 // OTHERWISE - CASE 3 OR 4
1128 // Load inverse of P_0 .
1129 // Set PR_6 if Arg <= -2**63
1130 // Are there any Infs, NaNs, or zeros?
1134 (p0) ldfe P_0 = [table_ptr1],16 ;;
1135 (p0) ldfe d_2 = [table_ptr2],16
1139 // Set PR_8 if Arg <= -2**24
1140 // Set PR_6 if Arg >= 2**63
1144 (p0) ldfe P_1 = [table_ptr1],16 ;;
1145 (p0) ldfe PI_BY_4 = [table_ptr2],16
1149 // Set PR_8 if Arg >= 2**24
1153 (p0) ldfe P_2 = [table_ptr1],16 ;;
1154 (p0) ldfe MPI_BY_4 = [table_ptr2],16
1158 // Load P_2 and PI_BY_4
1162 (p0) ldfe P_3 = [table_ptr1],16
1169 (p0) fcmp.le.unc.s1 p6,p7 = Arg,NEGTWO_TO_63
1175 (p0) fcmp.le.unc.s1 p8,p9 = Arg,NEGTWO_TO_24
1181 (p7) fcmp.ge.s1 p6,p0 = Arg,TWO_TO_63
1187 (p9) fcmp.ge.s1 p8,p0 = Arg,TWO_TO_24
1195 // Load P_3 and -PI_BY_4
1197 (p6) br.cond.spnt TAN_ARG_TOO_LARGE ;;
1206 // Branch out if we have a special argument.
1207 // Branch out if the magnitude of the input argument is too large
1208 // - do this branch before the next.
1210 (p8) br.cond.spnt TAN_LARGER_ARG ;;
1213 // Branch to Cases 3 or 4 if Arg <= -2**24 or Arg >= 2**24
1217 (p0) ldfs TWO_TO_NEG2 = [table_ptr2],4
1218 // ARGUMENT REDUCTION CODE - CASE 1 and 2
1221 (p0) fmpy.s1 N = Arg,two_by_PI
1226 (p0) ldfs NEGTWO_TO_NEG2 = [table_ptr2],12
1230 (p0) fcmp.lt.unc.s1 p8,p9= Arg,PI_BY_4
1237 // if Arg < pi/4, set PR_8.
1239 (p8) fcmp.gt.s1 p8,p9= Arg,MPI_BY_4
1243 // Case 1: Is |r| < 2**(-2).
1244 // Arg is the same as r in this case.
1250 (p8) mov N_fix_gr = r0
1252 // if Arg > -pi/4, reset PR_8.
1253 // Select the case when |Arg| < pi/4 - set PR[8] = true.
1254 // Else Select the case when |Arg| >= pi/4 - set PR[9] = true.
1256 (p0) fcvt.fx.s1 N_fix = N
1263 // Grab the integer part of N .
1277 (p8) fcmp.lt.unc.s1 p10, p11 = Arg, TWO_TO_NEG2
1283 (p10) fcmp.gt.s1 p10,p0 = Arg, NEGTWO_TO_NEG2
1290 // Case 2: Place integer part of N in GP register.
1292 (p9) fcvt.xf N = N_fix
1297 (p9) getf.sig N_fix_gr = N_fix
1300 // Case 2: Convert integer N_fix back to normalized floating-point value.
1302 (p10) br.cond.spnt TAN_SMALL_R ;;
1308 (p8) br.cond.sptk TAN_NORMAL_R ;;
1311 // Case 1: PR_3 is only affected when PR_1 is set.
1315 (p9) ldfs TWO_TO_NEG33 = [table_ptr2], 4 ;;
1317 // Case 2: Load 2**(-33).
1319 (p9) ldfs NEGTWO_TO_NEG33 = [table_ptr2], 4
1326 // Case 2: Load -2**(-33).
1328 (p9) fnma.s1 s_val = N, P_1, Arg
1334 (p9) fmpy.s1 w = N, P_2
1341 // Case 2: w = N * P_2
1342 // Case 2: s_val = -N * P_1 + Arg
1344 (p0) fcmp.lt.unc.s1 p9,p8 = s_val, TWO_TO_NEG33
1351 // Decide between case_1 and case_2 reduce:
1353 (p9) fcmp.gt.s1 p9, p8 = s_val, NEGTWO_TO_NEG33
1360 // Case 1_reduce: s <= -2**(-33) or s >= 2**(-33)
1361 // Case 2_reduce: -2**(-33) < s < 2**(-33)
1363 (p8) fsub.s1 r = s_val, w
1369 (p9) fmpy.s1 w = N, P_3
1375 (p9) fma.s1 U_1 = N, P_2, w
1382 // Case 1_reduce: Is |r| < 2**(-2), if so set PR_10
1385 (p8) fsub.s1 c = s_val, r
1392 // Case 1_reduce: r = s + w (change sign)
1393 // Case 2_reduce: w = N * P_3 (change sign)
1395 (p8) fcmp.lt.unc.s1 p10, p11 = r, TWO_TO_NEG2
1401 (p10) fcmp.gt.s1 p10, p11 = r, NEGTWO_TO_NEG2
1407 (p9) fsub.s1 r = s_val, U_1
1414 // Case 1_reduce: c is complete here.
1415 // c = c + w (w has not been negated.)
1416 // Case 2_reduce: r is complete here - continue to calculate c .
1419 (p9) fms.s1 U_2 = N, P_2, U_1
1426 // Case 1_reduce: c = s - r
1427 // Case 2_reduce: U_1 = N * P_2 + w
1429 (p8) fsub.s1 c = c, w
1435 (p9) fsub.s1 s_val = s_val, r
1443 // U_2 = N * P_2 - U_1
1444 // Not needed until later.
1446 (p9) fadd.s1 U_2 = U_2, w
1452 (p10) br.cond.spnt TAN_SMALL_R ;;
1458 (p11) br.cond.sptk TAN_NORMAL_R ;;
1466 // c is complete here
1467 // Argument reduction ends here.
1469 (p9) extr.u i_1 = N_fix_gr, 0, 1 ;;
1470 (p9) cmp.eq.unc p11, p12 = 0x0000,i_1 ;;
1476 // Is i_1 even or odd?
1477 // if i_1 == 0, set p11, else set p12.
1479 (p11) fmpy.s1 rsq = r, Z
1485 (p12) frcpa.s1 S_hi,p0 = f1, r
1490 // Case 1: Branch to SMALL_R or NORMAL_R.
1491 // Case 1 is done now.
1495 (p9) addl table_ptr1 = @ltoff(TAN_BASE_CONSTANTS), gp
1496 (p9) fsub.s1 c = s_val, U_1
1502 (p9) ld8 table_ptr1 = [table_ptr1]
1509 (p9) add table_ptr1 = 224, table_ptr1 ;;
1510 (p9) ldfe P1_1 = [table_ptr1],144
1514 // Get [i_1] - lsb of N_fix_gr .
1515 // Load P1_1 and point to Q1_1 .
1519 (p9) ldfe Q1_1 = [table_ptr1] , 0
1521 // N even: rsq = r * Z
1522 // N odd: S_hi = frcpa(r)
1524 (p12) fmerge.ns S_hi = S_hi, S_hi
1534 (p9) fsub.s1 c = c, U_2
1540 (p12) fma.s1 poly1 = S_hi, r, f1
1547 // N odd: Change sign of S_hi
1549 (p11) fmpy.s1 rsq = rsq, P1_1
1555 (p12) fma.s1 S_hi = S_hi, poly1, S_hi
1562 // N even: rsq = rsq * P1_1
1563 // N odd: poly1 = 1.0 + S_hi * r 16 bits partial account for necessary
1565 (p11) fma.s1 Result = r, rsq, c
1572 // N even: Result = c + r * rsq
1573 // N odd: S_hi = S_hi + S_hi*poly1 16 bits account for necessary
1575 (p12) fma.s1 poly1 = S_hi, r, f1
1582 // N even: Result = Result + r
1583 // N odd: poly1 = 1.0 + S_hi * r 32 bits partial
1585 (p11) fadd.s0 Result = r, Result
1591 (p12) fma.s1 S_hi = S_hi, poly1, S_hi
1598 // N even: Result1 = Result + r
1599 // N odd: S_hi = S_hi * poly1 + S_hi 32 bits
1601 (p12) fma.s1 poly1 = S_hi, r, f1
1608 // N odd: poly1 = S_hi * r + 1.0 64 bits partial
1610 (p12) fma.s1 S_hi = S_hi, poly1, S_hi
1617 // N odd: poly1 = S_hi * poly + 1.0 64 bits
1619 (p12) fma.s1 poly1 = S_hi, r, f1
1626 // N odd: poly1 = S_hi * r + 1.0
1628 (p12) fma.s1 poly1 = S_hi, c, poly1
1635 // N odd: poly1 = S_hi * c + poly1
1637 (p12) fmpy.s1 S_lo = S_hi, poly1
1644 // N odd: S_lo = S_hi * poly1
1646 (p12) fma.s1 S_lo = Q1_1, r, S_lo
1653 // N odd: Result = S_hi + S_lo
1655 (p0) fmpy.s0 fp_tmp = fp_tmp, fp_tmp // Dummy mult to set inexact
1662 // N odd: S_lo = S_lo + Q1_1 * r
1664 (p12) fadd.s0 Result = S_hi, S_lo
1665 (p0) br.ret.sptk b0 ;;
1672 (p0) addl table_ptr1 = @ltoff(TAN_BASE_CONSTANTS), gp
1674 (p0) fmpy.s1 N_0 = Arg, Inv_P_0
1679 // ARGUMENT REDUCTION CODE - CASE 3 and 4
1682 // Adjust table_ptr1 to beginning of table.
1683 // N_0 = Arg * Inv_P_0
1688 (p0) ld8 table_ptr1 = [table_ptr1]
1696 (p0) add table_ptr1 = 8, table_ptr1 ;;
1700 (p0) ldfs TWO_TO_NEG14 = [table_ptr1], 4
1708 (p0) ldfs NEGTWO_TO_NEG14 = [table_ptr1], 180 ;;
1710 // N_0_fix = integer part of N_0 .
1711 // Adjust table_ptr1 to beginning of table.
1713 (p0) ldfs TWO_TO_NEG2 = [table_ptr1], 4
1717 // Make N_0 the integer part.
1721 (p0) ldfs NEGTWO_TO_NEG2 = [table_ptr1]
1725 (p0) fcvt.fx.s1 N_0_fix = N_0
1731 (p0) fcvt.xf N_0 = N_0_fix
1737 (p0) fnma.s1 ArgPrime = N_0, P_0, Arg
1743 (p0) fmpy.s1 w = N_0, d_1
1750 // ArgPrime = -N_0 * P_0 + Arg
1753 (p0) fmpy.s1 N = ArgPrime, two_by_PI
1760 // N = ArgPrime * 2/pi
1762 (p0) fcvt.fx.s1 N_fix = N
1769 // N_fix is the integer part.
1771 (p0) fcvt.xf N = N_fix
1776 (p0) getf.sig N_fix_gr = N_fix
1784 // N is the integer part of the reduced-reduced argument.
1785 // Put the integer in a GP register.
1787 (p0) fnma.s1 s_val = N, P_1, ArgPrime
1793 (p0) fnma.s1 w = N, P_2, w
1800 // s_val = -N*P_1 + ArgPrime
1803 (p0) fcmp.lt.unc.s1 p11, p10 = s_val, TWO_TO_NEG14
1809 (p11) fcmp.gt.s1 p11, p10 = s_val, NEGTWO_TO_NEG14
1816 // Case 3: r = s_val + w (Z complete)
1817 // Case 4: U_hi = N_0 * d_1
1819 (p10) fmpy.s1 V_hi = N, P_2
1825 (p11) fmpy.s1 U_hi = N_0, d_1
1832 // Case 3: r = s_val + w (Z complete)
1833 // Case 4: U_hi = N_0 * d_1
1835 (p11) fmpy.s1 V_hi = N, P_2
1841 (p11) fmpy.s1 U_hi = N_0, d_1
1848 // Decide between case 3 and 4:
1849 // Case 3: s <= -2**(-14) or s >= 2**(-14)
1850 // Case 4: -2**(-14) < s < 2**(-14)
1852 (p10) fadd.s1 r = s_val, w
1858 (p11) fmpy.s1 w = N, P_3
1865 // Case 4: We need abs of both U_hi and V_hi - dont
1866 // worry about switched sign of V_hi .
1868 (p11) fsub.s1 A = U_hi, V_hi
1875 // Case 4: A = U_hi + V_hi
1876 // Note: Worry about switched sign of V_hi, so subtract instead of add.
1878 (p11) fnma.s1 V_lo = N, P_2, V_hi
1884 (p11) fms.s1 U_lo = N_0, d_1, U_hi
1890 (p11) fabs V_hiabs = V_hi
1897 // Case 4: V_hi = N * P_2
1899 // Note the product does not include the (-) as in the writeup
1900 // so (-) missing for V_hi and w .
1901 (p10) fadd.s1 r = s_val, w
1908 // Case 3: c = s_val - r
1909 // Case 4: U_lo = N_0 * d_1 - U_hi
1911 (p11) fabs U_hiabs = U_hi
1917 (p11) fmpy.s1 w = N, P_3
1924 // Case 4: Set P_12 if U_hiabs >= V_hiabs
1926 (p11) fadd.s1 C_hi = s_val, A
1933 // Case 4: C_hi = s_val + A
1935 (p11) fadd.s1 t = U_lo, V_lo
1942 // Case 3: Is |r| < 2**(-2), if so set PR_7
1944 // Case 3: If PR_7 is set, prepare to branch to Small_R.
1945 // Case 3: If PR_8 is set, prepare to branch to Normal_R.
1947 (p10) fsub.s1 c = s_val, r
1954 // Case 3: c = (s - r) + w (c complete)
1956 (p11) fcmp.ge.unc.s1 p12, p13 = U_hiabs, V_hiabs
1962 (p11) fms.s1 w = N_0, d_2, w
1969 // Case 4: V_hi = N * P_2
1971 // Note the product does not include the (-) as in the writeup
1972 // so (-) missing for V_hi and w .
1974 (p10) fcmp.lt.unc.s1 p14, p15 = r, TWO_TO_NEG2
1980 (p14) fcmp.gt.s1 p14, p15 = r, NEGTWO_TO_NEG2
1987 // Case 4: V_lo = -N * P_2 - V_hi (U_hi is in place of V_hi in writeup)
1988 // Note: the (-) is still missing for V_hi .
1989 // Case 4: w = w + N_0 * d_2
1990 // Note: the (-) is now incorporated in w .
1992 (p10) fadd.s1 c = c, w
1994 // Case 4: t = U_lo + V_lo
1995 // Note: remember V_lo should be (-), subtract instead of add. NO
1997 (p14) br.cond.spnt TAN_SMALL_R ;;
2003 (p15) br.cond.spnt TAN_NORMAL_R ;;
2009 // Case 3: Vector off when |r| < 2**(-2). Recall that PR_3 will be true.
2010 // The remaining stuff is for Case 4.
2012 (p12) fsub.s1 a = U_hi, A
2013 (p11) extr.u i_1 = N_fix_gr, 0, 1 ;;
2019 // Case 4: C_lo = s_val - C_hi
2021 (p11) fadd.s1 t = t, w
2027 (p13) fadd.s1 a = V_hi, A
2032 // Case 4: a = U_hi - A
2033 // a = V_hi - A (do an add to account for missing (-) on V_hi
2037 (p11) addl table_ptr1 = @ltoff(TAN_BASE_CONSTANTS), gp
2038 (p11) fsub.s1 C_lo = s_val, C_hi
2044 (p11) ld8 table_ptr1 = [table_ptr1]
2051 // Case 4: a = (U_hi - A) + V_hi
2052 // a = (V_hi - A) + U_hi
2053 // In each case account for negative missing form V_hi .
2056 // Case 4: C_lo = (s_val - C_hi) + A
2060 (p11) add table_ptr1 = 224, table_ptr1 ;;
2061 (p11) ldfe P1_1 = [table_ptr1], 16
2066 (p11) ldfe P1_2 = [table_ptr1], 128
2068 // Case 4: w = U_lo + V_lo + w
2070 (p12) fsub.s1 a = a, V_hi
2074 // Case 4: r = C_hi + C_lo
2078 (p11) ldfe Q1_1 = [table_ptr1], 16
2079 (p11) fadd.s1 C_lo = C_lo, A
2083 // Case 4: c = C_hi - r
2084 // Get [i_1] - lsb of N_fix_gr.
2088 (p11) ldfe Q1_2 = [table_ptr1], 16
2095 (p13) fsub.s1 a = U_hi, a
2101 (p11) fadd.s1 t = t, a
2108 // Case 4: t = t + a
2110 (p11) fadd.s1 C_lo = C_lo, t
2117 // Case 4: C_lo = C_lo + t
2119 (p11) fadd.s1 r = C_hi, C_lo
2125 (p11) fsub.s1 c = C_hi, r
2132 // Case 4: c = c + C_lo finished.
2133 // Is i_1 even or odd?
2134 // if i_1 == 0, set PR_4, else set PR_5.
2136 // r and c have been computed.
2137 // We known whether this is the sine or cosine routine.
2138 // Make sure ftz mode is set - should be automatic when using wre
2139 (p0) fmpy.s1 rsq = r, r
2145 (p11) fadd.s1 c = c , C_lo
2146 (p11) cmp.eq.unc p11, p12 = 0x0000, i_1 ;;
2151 (p12) frcpa.s1 S_hi, p0 = f1, r
2158 // N odd: Change sign of S_hi
2160 (p11) fma.s1 Result = rsq, P1_2, P1_1
2166 (p12) fma.s1 P = rsq, Q1_2, Q1_1
2173 // N odd: Result = S_hi + S_lo (User supplied rounding mode for C1)
2175 (p0) fmpy.s0 fp_tmp = fp_tmp, fp_tmp // Dummy mult to set inexact
2182 // N even: rsq = r * r
2183 // N odd: S_hi = frcpa(r)
2185 (p12) fmerge.ns S_hi = S_hi, S_hi
2192 // N even: rsq = rsq * P1_2 + P1_1
2193 // N odd: poly1 = 1.0 + S_hi * r 16 bits partial account for necessary
2195 (p11) fmpy.s1 Result = rsq, Result
2201 (p12) fma.s1 poly1 = S_hi, r,f1
2208 // N even: Result = Result * rsq
2209 // N odd: S_hi = S_hi + S_hi*poly1 16 bits account for necessary
2211 (p11) fma.s1 Result = r, Result, c
2217 (p12) fma.s1 S_hi = S_hi, poly1, S_hi
2224 // N odd: S_hi = S_hi * poly1 + S_hi 32 bits
2226 (p11) fadd.s0 Result= r, Result
2232 (p12) fma.s1 poly1 = S_hi, r, f1
2239 // N even: Result = Result * r + c
2240 // N odd: poly1 = 1.0 + S_hi * r 32 bits partial
2242 (p12) fma.s1 S_hi = S_hi, poly1, S_hi
2248 (p12) fma.s1 poly1 = S_hi, r, f1
2255 // N even: Result1 = Result + r (Rounding mode S0)
2256 // N odd: poly1 = S_hi * r + 1.0 64 bits partial
2258 (p12) fma.s1 S_hi = S_hi, poly1, S_hi
2265 // N odd: poly1 = S_hi * poly + S_hi 64 bits
2267 (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
2283 // N odd: poly1 = S_hi * c + poly1
2285 (p12) fmpy.s1 S_lo = S_hi, poly1
2292 // N odd: S_lo = S_hi * poly1
2294 (p12) fma.s1 S_lo = P, r, S_lo
2301 // N odd: S_lo = S_lo + r * P
2303 (p12) fadd.s0 Result = S_hi, S_lo
2304 (p0) br.ret.sptk b0 ;;
2312 (p0) extr.u i_1 = N_fix_gr, 0, 1 ;;
2313 (p0) cmp.eq.unc p11, p12 = 0x0000, i_1
2318 (p0) fmpy.s1 rsq = r, r
2324 (p12) frcpa.s1 S_hi, p0 = f1, r
2329 (p0) addl table_ptr1 = @ltoff(TAN_BASE_CONSTANTS), gp
2336 (p0) ld8 table_ptr1 = [table_ptr1]
2342 // *****************************************************************
2343 // *****************************************************************
2344 // *****************************************************************
2347 (p0) add table_ptr1 = 224, table_ptr1 ;;
2348 (p0) ldfe P1_1 = [table_ptr1], 16
2351 // r and c have been computed.
2352 // We known whether this is the sine or cosine routine.
2353 // Make sure ftz mode is set - should be automatic when using wre
2357 (p0) ldfe P1_2 = [table_ptr1], 16
2358 (p11) fmpy.s1 r_to_the_8 = rsq, rsq
2362 // Set table_ptr1 to beginning of constant table.
2363 // Get [i_1] - lsb of N_fix_gr.
2367 (p0) ldfe P1_3 = [table_ptr1], 96
2369 // N even: rsq = r * r
2370 // N odd: S_hi = frcpa(r)
2372 (p12) fmerge.ns S_hi = S_hi, S_hi
2376 // Is i_1 even or odd?
2377 // if i_1 == 0, set PR_11.
2378 // if i_1 != 0, set PR_12.
2382 (p11) ldfe P1_9 = [table_ptr1], -16
2384 // N even: Poly2 = P1_7 + Poly2 * rsq
2385 // N odd: poly2 = Q1_5 + poly2 * rsq
2387 (p11) fadd.s1 CORR = rsq, f1
2392 (p11) ldfe P1_8 = [table_ptr1], -16 ;;
2394 // N even: Poly1 = P1_2 + P1_3 * rsq
2395 // N odd: poly1 = 1.0 + S_hi * r
2396 // 16 bits partial account for necessary (-1)
2398 (p11) ldfe P1_7 = [table_ptr1], -16
2402 // N even: Poly1 = P1_1 + Poly1 * rsq
2403 // N odd: S_hi = S_hi + S_hi * poly1) 16 bits account for necessary
2407 (p11) ldfe P1_6 = [table_ptr1], -16
2409 // N even: Poly2 = P1_5 + Poly2 * rsq
2410 // N odd: poly2 = Q1_3 + poly2 * rsq
2412 (p11) fmpy.s1 r_to_the_8 = r_to_the_8, r_to_the_8
2416 // N even: Poly1 = Poly1 * rsq
2417 // N odd: poly1 = 1.0 + S_hi * r 32 bits partial
2421 (p11) ldfe P1_5 = [table_ptr1], -16
2422 (p12) fma.s1 poly1 = S_hi, r, f1
2426 // N even: CORR = CORR * c
2427 // N odd: S_hi = S_hi * poly1 + S_hi 32 bits
2431 // N even: Poly2 = P1_6 + Poly2 * rsq
2432 // N odd: poly2 = Q1_4 + poly2 * rsq
2435 (p0) addl table_ptr2 = @ltoff(TAN_BASE_CONSTANTS), gp
2436 (p11) ldfe P1_4 = [table_ptr1], -16
2437 (p11) fmpy.s1 CORR = CORR, c
2443 (p0) ld8 table_ptr2 = [table_ptr2]
2451 (p0) add table_ptr2 = 464, table_ptr2
2458 (p11) fma.s1 Poly1 = P1_3, rsq, P1_2
2463 (p0) ldfe Q1_7 = [table_ptr2], -16
2464 (p12) fma.s1 S_hi = S_hi, poly1, S_hi
2469 (p0) ldfe Q1_6 = [table_ptr2], -16
2470 (p11) fma.s1 Poly2 = P1_9, rsq, P1_8
2475 (p0) ldfe Q1_5 = [table_ptr2], -16 ;;
2476 (p12) ldfe Q1_4 = [table_ptr2], -16
2481 (p12) ldfe Q1_3 = [table_ptr2], -16
2483 // N even: Poly2 = P1_8 + P1_9 * rsq
2484 // N odd: poly2 = Q1_6 + Q1_7 * rsq
2486 (p11) fma.s1 Poly1 = Poly1, rsq, P1_1
2491 (p12) ldfe Q1_2 = [table_ptr2], -16
2492 (p12) fma.s1 poly1 = S_hi, r, f1
2497 (p12) ldfe Q1_1 = [table_ptr2], -16
2498 (p11) fma.s1 Poly2 = Poly2, rsq, P1_7
2505 // N even: CORR = rsq + 1
2506 // N even: r_to_the_8 = rsq * rsq
2508 (p11) fmpy.s1 Poly1 = Poly1, rsq
2514 (p12) fma.s1 S_hi = S_hi, poly1, S_hi
2520 (p12) fma.s1 poly2 = Q1_7, rsq, Q1_6
2526 (p11) fma.s1 Poly2 = Poly2, rsq, P1_6
2532 (p12) fma.s1 poly1 = S_hi, r, f1
2538 (p12) fma.s1 poly2 = poly2, rsq, Q1_5
2544 (p11) fma.s1 Poly2= Poly2, rsq, P1_5
2550 (p12) fma.s1 S_hi = S_hi, poly1, S_hi
2556 (p12) fma.s1 poly2 = poly2, rsq, Q1_4
2563 // N even: r_to_the_8 = r_to_the_8 * r_to_the_8
2564 // N odd: poly1 = S_hi * r + 1.0 64 bits partial
2566 (p11) fma.s1 Poly2 = Poly2, rsq, P1_4
2573 // N even: Result = CORR + Poly * r
2574 // N odd: P = Q1_1 + poly2 * rsq
2576 (p12) fma.s1 poly1 = S_hi, r, f1
2582 (p12) fma.s1 poly2 = poly2, rsq, Q1_3
2589 // N even: Poly2 = P1_4 + Poly2 * rsq
2590 // N odd: poly2 = Q1_2 + poly2 * rsq
2592 (p11) fma.s1 Poly = Poly2, r_to_the_8, Poly1
2598 (p12) fma.s1 poly1 = S_hi, c, poly1
2604 (p12) fma.s1 poly2 = poly2, rsq, Q1_2
2611 // N even: Poly = Poly1 + Poly2 * r_to_the_8
2612 // N odd: S_hi = S_hi * poly1 + S_hi 64 bits
2614 (p11) fma.s1 Result = Poly, r, CORR
2621 // N even: Result = r + Result (User supplied rounding mode)
2622 // N odd: poly1 = S_hi * c + poly1
2624 (p12) fmpy.s1 S_lo = S_hi, poly1
2630 (p12) fma.s1 P = poly2, rsq, Q1_1
2637 // N odd: poly1 = S_hi * r + 1.0
2639 (p11) fadd.s0 Result = Result, r
2646 // N odd: S_lo = S_hi * poly1
2648 (p12) fma.s1 S_lo = Q1_1, c, S_lo
2655 // N odd: Result = Result + S_hi (user supplied rounding mode)
2657 (p0) fmpy.s0 fp_tmp = fp_tmp, fp_tmp // Dummy mult to set inexact
2664 // N odd: S_lo = Q1_1 * c + S_lo
2666 (p12) fma.s1 Result = P, r, S_lo
2673 // N odd: Result = S_lo + r * P
2675 (p12) fadd.s0 Result = Result, S_hi
2676 (p0) br.ret.sptk b0 ;;
2683 (p0) getf.sig sig_r = r
2684 // *******************************************************************
2685 // *******************************************************************
2686 // *******************************************************************
2688 // r and c have been computed.
2689 // Make sure ftz mode is set - should be automatic when using wre
2692 // Get [i_1] - lsb of N_fix_gr alone.
2694 (p0) fmerge.s Pos_r = f1, r
2695 (p0) extr.u i_1 = N_fix_gr, 0, 1 ;;
2700 (p0) fmerge.s sgn_r = r, f1
2701 (p0) cmp.eq.unc p11, p12 = 0x0000, i_1 ;;
2707 (p0) extr.u lookup = sig_r, 58, 5
2712 (p0) movl Create_B = 0x8200000000000000 ;;
2716 (p0) addl table_ptr1 = @ltoff(TAN_BASE_CONSTANTS), gp
2718 (p0) dep Create_B = lookup, Create_B, 58, 5
2723 // Get [i_1] - lsb of N_fix_gr alone.
2729 ld8 table_ptr1 = [table_ptr1]
2738 (p0) setf.sig B = Create_B
2740 // Set table_ptr1 and table_ptr2 to base address of
2743 (p0) add table_ptr1 = 480, table_ptr1 ;;
2749 // Is i_1 or i_0 == 0 ?
2750 // Create the constant 1 00000 1000000000000000000000...
2752 (p0) ldfe P2_1 = [table_ptr1], 16
2758 (p0) getf.exp exp_r = Pos_r
2763 // Get r's significand
2767 (p0) ldfe P2_2 = [table_ptr1], 16 ;;
2769 // Get the 5 bits or r for the lookup. 1.xxxxx ....
2771 // Grab lsb of exp of B
2773 (p0) ldfe P2_3 = [table_ptr1], 16
2779 (p0) andcm table_offset = 0x0001, exp_r ;;
2780 (p0) shl table_offset = table_offset, 9 ;;
2786 // Deposit 0 00000 1000000000000000000000... on
2787 // 1 xxxxx yyyyyyyyyyyyyyyyyyyyyy...,
2788 // getting rid of the ys.
2789 // Is B = 2** -2 or B= 2** -1? If 2**-1, then
2790 // we want an offset of 512 for table addressing.
2792 (p0) shladd table_offset = lookup, 4, table_offset ;;
2794 // B = ........ 1xxxxx 1000000000000000000...
2796 (p0) add table_ptr1 = table_ptr1, table_offset ;;
2802 // B = ........ 1xxxxx 1000000000000000000...
2803 // Convert B so it has the same exponent as Pos_r
2805 (p0) ldfd T_hi = [table_ptr1], 8
2816 (p0) addl table_ptr2 = @ltoff(TAN_BASE_CONSTANTS), gp
2817 (p0) ldfs T_lo = [table_ptr1]
2818 (p0) fmerge.se B = Pos_r, B
2823 ld8 table_ptr2 = [table_ptr2]
2830 (p0) add table_ptr2 = 1360, table_ptr2
2832 (p0) add table_ptr2 = table_ptr2, table_offset ;;
2836 (p0) ldfd C_hi = [table_ptr2], 8
2837 (p0) fsub.s1 x = Pos_r, B
2842 (p0) ldfs C_lo = [table_ptr2],255
2846 // N even: Tx = T_hi * x
2848 // Load C_lo - increment pointer to get SC_inv
2849 // - cant get all the way, do an add later.
2851 (p0) add table_ptr2 = 569, table_ptr2 ;;
2854 // N even: Tx1 = Tx + 1
2855 // N odd: Cx1 = 1 - Cx
2859 (p0) ldfe SC_inv = [table_ptr2], 0
2866 (p0) fmpy.s1 xsq = x, x
2872 (p11) fmpy.s1 Tx = T_hi, x
2878 (p12) fmpy.s1 Cx = C_hi, x
2885 // N odd: Cx = C_hi * x
2887 (p0) fma.s1 P = P2_3, xsq, P2_2
2894 // N even and odd: P = P2_3 + P2_2 * xsq
2896 (p11) fadd.s1 Tx1 = Tx, f1
2903 // N even: D = C_hi - tanx
2904 // N odd: D = T_hi + tanx
2906 (p11) fmpy.s1 CORR = SC_inv, T_hi
2912 (p0) fmpy.s1 Sx = SC_inv, x
2918 (p12) fmpy.s1 CORR = SC_inv, C_hi
2924 (p12) fsub.s1 V_hi = f1, Cx
2930 (p0) fma.s1 P = P, xsq, P2_1
2937 // N even and odd: P = P2_1 + P * xsq
2939 (p11) fma.s1 V_hi = Tx, Tx1, f1
2946 // N even: Result = sgn_r * tail + T_hi (user rounding mode for C1)
2947 // N odd: Result = sgn_r * tail + C_hi (user rounding mode for C1)
2949 (p0) fmpy.s0 fp_tmp = fp_tmp, fp_tmp // Dummy mult to set inexact
2955 (p0) fmpy.s1 CORR = CORR, c
2961 (p12) fnma.s1 V_hi = Cx,V_hi,f1
2968 // N even: V_hi = Tx * Tx1 + 1
2969 // N odd: Cx1 = 1 - Cx * Cx1
2971 (p0) fmpy.s1 P = P, xsq
2978 // N even and odd: P = P * xsq
2980 (p11) fmpy.s1 V_hi = V_hi, T_hi
2987 // N even and odd: tail = P * tail + V_lo
2989 (p11) fmpy.s1 T_hi = sgn_r, T_hi
2995 (p0) fmpy.s1 CORR = CORR, sgn_r
3001 (p12) fmpy.s1 V_hi = V_hi,C_hi
3008 // N even: V_hi = T_hi * V_hi
3009 // N odd: V_hi = C_hi * V_hi
3011 (p0) fma.s1 tanx = P, x, x
3017 (p12) fnmpy.s1 C_hi = sgn_r, C_hi
3024 // N even: V_lo = 1 - V_hi + C_hi
3025 // N odd: V_lo = 1 - V_hi + T_hi
3027 (p11) fadd.s1 CORR = CORR, T_lo
3033 (p12) fsub.s1 CORR = CORR, C_lo
3040 // N even and odd: tanx = x + x * P
3041 // N even and odd: Sx = SC_inv * x
3043 (p11) fsub.s1 D = C_hi, tanx
3049 (p12) fadd.s1 D = T_hi, tanx
3056 // N odd: CORR = SC_inv * C_hi
3057 // N even: CORR = SC_inv * T_hi
3059 (p0) fnma.s1 D = V_hi, D, f1
3066 // N even and odd: D = 1 - V_hi * D
3067 // N even and odd: CORR = CORR * c
3069 (p0) fma.s1 V_hi = V_hi, D, V_hi
3076 // N even and odd: V_hi = V_hi + V_hi * D
3077 // N even and odd: CORR = sgn_r * CORR
3079 (p11) fnma.s1 V_lo = V_hi, C_hi, f1
3085 (p12) fnma.s1 V_lo = V_hi, T_hi, f1
3092 // N even: CORR = COOR + T_lo
3093 // N odd: CORR = CORR - C_lo
3095 (p11) fma.s1 V_lo = tanx, V_hi, V_lo
3101 (p12) fnma.s1 V_lo = tanx, V_hi, V_lo
3108 // N even: V_lo = V_lo + V_hi * tanx
3109 // N odd: V_lo = V_lo - V_hi * tanx
3111 (p11) fnma.s1 V_lo = C_lo, V_hi, V_lo
3117 (p12) fnma.s1 V_lo = T_lo, V_hi, V_lo
3124 // N even: V_lo = V_lo - V_hi * C_lo
3125 // N odd: V_lo = V_lo - V_hi * T_lo
3127 (p0) fmpy.s1 V_lo = V_hi, V_lo
3134 // N even and odd: V_lo = V_lo * V_hi
3136 (p0) fadd.s1 tail = V_hi, V_lo
3143 // N even and odd: tail = V_hi + V_lo
3145 (p0) fma.s1 tail = tail, P, V_lo
3152 // N even: T_hi = sgn_r * T_hi
3153 // N odd : C_hi = -sgn_r * C_hi
3155 (p0) fma.s1 tail = tail, Sx, CORR
3162 // N even and odd: tail = Sx * tail + CORR
3164 (p0) fma.s1 tail = V_hi, Sx, tail
3171 // N even an odd: tail = Sx * V_hi + tail
3173 (p11) fma.s0 Result = sgn_r, tail, T_hi
3179 (p12) fma.s0 Result = sgn_r, tail, C_hi
3180 (p0) br.ret.sptk b0 ;;
3184 ASM_SIZE_DIRECTIVE(__libm_tan)
3188 // *******************************************************************
3189 // *******************************************************************
3190 // *******************************************************************
3192 // Special Code to handle very large argument case.
3193 // Call int pi_by_2_reduce(&x,&r)
3194 // for |arguments| >= 2**63
3195 // (Arg or x) is in f8
3196 // Address to save r and c as double
3198 // (1) (2) (3) (call) (4)
3199 // sp -> + psp -> + psp -> + sp -> +
3201 // | r50 ->| <- r50 f0 ->| r50 -> | -> c
3203 // sp-32 -> | <- r50 f0 ->| f0 ->| <- r50 r49 -> | -> r
3205 // | r49 ->| <- r49 Arg ->| <- r49 | -> x
3207 // sp -64 ->| sp -64 ->| sp -64 ->| |
3209 // save pfs save b0 restore gp
3210 // save gp restore b0
3215 .proc __libm_callout
3221 add GR_Parameter_r =-32,sp // Parameter: r address
3223 .save ar.pfs,GR_SAVE_PFS
3224 mov GR_SAVE_PFS=ar.pfs // Save ar.pfs
3228 add sp=-64,sp // Create new stack
3230 mov GR_SAVE_GP=gp // Save gp
3235 stfe [GR_Parameter_r ] = f0,16 // Clear Parameter r on stack
3236 add GR_Parameter_X = 16,sp // Parameter x address
3237 .save b0, GR_SAVE_B0
3238 mov GR_SAVE_B0=b0 // Save b0
3244 stfe [GR_Parameter_r ] = f0,-16 // Clear Parameter c on stack
3249 stfe [GR_Parameter_X] = Arg // Store Parameter x on stack
3251 (p0) br.call.sptk b0=__libm_pi_by_2_reduce#
3258 mov gp = GR_SAVE_GP // Restore gp
3259 (p0) mov N_fix_gr = r8
3265 (p0) ldfe Arg =[GR_Parameter_X],16
3266 (p0) ldfs TWO_TO_NEG2 = [table_ptr2],4
3273 (p0) ldfe r =[GR_Parameter_r ],16
3274 (p0) ldfs NEGTWO_TO_NEG2 = [table_ptr2],4
3279 (p0) ldfe c =[GR_Parameter_r ]
3289 (p0) fcmp.lt.unc.s1 p6, p0 = r, TWO_TO_NEG2
3290 mov b0 = GR_SAVE_B0 // Restore return address
3296 (p6) fcmp.gt.unc.s1 p6, p0 = r, NEGTWO_TO_NEG2
3297 mov ar.pfs = GR_SAVE_PFS // Restore ar.pfs
3303 add sp = 64,sp // Restore stack pointer
3304 (p6) br.cond.spnt TAN_SMALL_R
3305 (p0) br.cond.sptk TAN_NORMAL_R
3308 .endp __libm_callout
3309 ASM_SIZE_DIRECTIVE(__libm_callout)
3312 .proc __libm_TAN_SPECIAL
3316 // Code for NaNs, Unsupporteds, Infs, or +/- zero ?
3317 // Invalid raised for Infs and SNaNs.
3322 (p0) fmpy.s0 Arg = Arg, f0
3325 .endp __libm_TAN_SPECIAL
3326 ASM_SIZE_DIRECTIVE(__libm_TAN_SPECIAL)
3329 .type __libm_pi_by_2_reduce#,@function
3330 .global __libm_pi_by_2_reduce#