sysdeps/ia64/fpu/s_tanl.S

   1 .file "tancotl.s"
   2
   3
   4 // Copyright (c) 2000 - 2004, Intel Corporation
   5 // All rights reserved.
   6 //
   7 //
   8 // Redistribution and use in source and binary forms, with or without
   9 // modification, are permitted provided that the following conditions are
  10 // met:
  11 //
  12 // * Redistributions of source code must retain the above copyright
  13 // notice, this list of conditions and the following disclaimer.
  14 //
  15 // * Redistributions in binary form must reproduce the above copyright
  16 // notice, this list of conditions and the following disclaimer in the
  17 // documentation and/or other materials provided with the distribution.
  18 //
  19 // * The name of Intel Corporation may not be used to endorse or promote
  20 // products derived from this software without specific prior written
  21 // permission.
  22
  23 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
  24 // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
  25 // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
  26 // A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL INTEL OR ITS
  27 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
  28 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
  29 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
  30 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY
  31 // OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY OR TORT (INCLUDING
  32 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
  33 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  34 //
  35 // Intel Corporation is the author of this code, and requests that all
  36 // problem reports or change requests be submitted to it directly at
  37 // http://www.intel.com/software/products/opensource/libraries/num.htm.
  38 //
  39 //*********************************************************************
  40 //
  41 // History:
  42 //
  43 // 02/02/00 (hand-optimized)
  44 // 04/04/00 Unwind support added
  45 // 12/28/00 Fixed false invalid flags
  46 // 02/06/02 Improved speed
  47 // 05/07/02 Changed interface to __libm_pi_by_2_reduce
  48 // 05/30/02 Added cotl
  49 // 02/10/03 Reordered header: .section, .global, .proc, .align;
  50 //          used data8 for long double table values
  51 // 05/15/03 Reformatted data tables
  52 // 10/26/04 Avoided using r14-31 as scratch so not clobbered by dynamic loader
  53 //
  54 //*********************************************************************
  55 //
  56 // Functions:   tanl(x) = tangent(x), for double-extended precision x values
  57 //              cotl(x) = cotangent(x), for double-extended precision x values
  58 //
  59 //*********************************************************************
  60 //
  61 // Resources Used:
  62 //
  63 //    Floating-Point Registers: f8 (Input and Return Value)
  64 //                              f9-f15
  65 //                              f32-f121
  66 //
  67 //    General Purpose Registers:
  68 //      r32-r70
  69 //
  70 //    Predicate Registers:      p6-p15
  71 //
  72 //*********************************************************************
  73 //
  74 // IEEE Special Conditions for tanl:
  75 //
  76 //    Denormal  fault raised on denormal inputs
  77 //    Overflow exceptions do not occur
  78 //    Underflow exceptions raised when appropriate for tan
  79 //    (No specialized error handling for this routine)
  80 //    Inexact raised when appropriate by algorithm
  81 //
  82 //    tanl(SNaN) = QNaN
  83 //    tanl(QNaN) = QNaN
  84 //    tanl(inf) = QNaN
  85 //    tanl(+/-0) = +/-0
  86 //
  87 //*********************************************************************
  88 //
  89 // IEEE Special Conditions for cotl:
  90 //
  91 //    Denormal  fault raised on denormal inputs
  92 //    Overflow exceptions occur at zero and near zero
  93 //    Underflow exceptions do not occur
  94 //    Inexact raised when appropriate by algorithm
  95 //
  96 //    cotl(SNaN) = QNaN
  97 //    cotl(QNaN) = QNaN
  98 //    cotl(inf) = QNaN
  99 //    cotl(+/-0) = +/-Inf and error handling is called
 100 //
 101 //*********************************************************************
 102 //
 103 //    Below are mathematical and algorithmic descriptions for tanl.
 104 //    For cotl we use next identity cot(x) = -tan(x + Pi/2).
 105 //    So, to compute cot(x) we just need to increment N (N = N + 1)
 106 //    and invert sign of the computed result.
 107 //
 108 //*********************************************************************
 109 //
 110 // Mathematical Description
 111 //
 112 // We consider the computation of FPTANL of Arg. Now, given
 113 //
 114 //      Arg = N pi/2  + alpha,          |alpha| <= pi/4,
 115 //
 116 // basic mathematical relationship shows that
 117 //
 118 //      tan( Arg ) =  tan( alpha )     if N is even;
 119 //                 = -cot( alpha )      otherwise.
 120 //
 121 // The value of alpha is obtained by argument reduction and
 122 // represented by two working precision numbers r and c where
 123 //
 124 //      alpha =  r  +  c     accurately.
 125 //
 126 // The reduction method is described in a previous write up.
 127 // The argument reduction scheme identifies 4 cases. For Cases 2
 128 // and 4, because |alpha| is small, tan(r+c) and -cot(r+c) can be
 129 // computed very easily by 2 or 3 terms of the Taylor series
 130 // expansion as follows:
 131 //
 132 // Case 2:
 133 // -------
 134 //
 135 //      tan(r + c) = r + c + r^3/3          ...accurately
 136 //     -cot(r + c) = -1/(r+c) + r/3          ...accurately
 137 //
 138 // Case 4:
 139 // -------
 140 //
 141 //      tan(r + c) = r + c + r^3/3 + 2r^5/15     ...accurately
 142 //     -cot(r + c) = -1/(r+c) + r/3 + r^3/45     ...accurately
 143 //
 144 //
 145 // The only cases left are Cases 1 and 3 of the argument reduction
 146 // procedure. These two cases will be merged since after the
 147 // argument is reduced in either cases, we have the reduced argument
 148 // represented as r + c and that the magnitude |r + c| is not small
 149 // enough to allow the usage of a very short approximation.
 150 //
 151 // The greatest challenge of this task is that the second terms of
 152 // the Taylor series for tan(r) and -cot(r)
 153 //
 154 //      r + r^3/3 + 2 r^5/15 + ...
 155 //
 156 // and
 157 //
 158 //      -1/r + r/3 + r^3/45 + ...
 159 //
 160 // are not very small when |r| is close to pi/4 and the rounding
 161 // errors will be a concern if simple polynomial accumulation is
 162 // used. When |r| < 2^(-2), however, the second terms will be small
 163 // enough (5 bits or so of right shift) that a normal Horner
 164 // recurrence suffices. Hence there are two cases that we consider
 165 // in the accurate computation of tan(r) and cot(r), |r| <= pi/4.
 166 //
 167 // Case small_r: |r| < 2^(-2)
 168 // --------------------------
 169 //
 170 // Since Arg = N pi/4 + r + c accurately, we have
 171 //
 172 //      tan(Arg) =  tan(r+c)            for N even,
 173 //               = -cot(r+c)            otherwise.
 174 //
 175 // Here for this case, both tan(r) and -cot(r) can be approximated
 176 // by simple polynomials:
 177 //
 178 //      tan(r) =    r + P1_1 r^3 + P1_2 r^5 + ... + P1_9 r^19
 179 //     -cot(r) = -1/r + Q1_1 r   + Q1_2 r^3 + ... + Q1_7 r^13
 180 //
 181 // accurately. Since |r| is relatively small, tan(r+c) and
 182 // -cot(r+c) can be accurately approximated by replacing r with
 183 // r+c only in the first two terms of the corresponding polynomials.
 184 //
 185 // Note that P1_1 (and Q1_1 for that matter) approximates 1/3 to
 186 // almost 64 sig. bits, thus
 187 //
 188 //      P1_1 (r+c)^3 =  P1_1 r^3 + c * r^2     accurately.
 189 //
 190 // Hence,
 191 //
 192 //      tan(r+c) =    r + P1_1 r^3 + P1_2 r^5 + ... + P1_9 r^19
 193 //                     + c*(1 + r^2)
 194 //
 195 //        -cot(r+c) = -1/(r+c) + Q1_1 r   + Q1_2 r^3 + ... + Q1_7 r^13
 196 //               + Q1_1*c
 197 //
 198 //
 199 // Case normal_r: 2^(-2) <= |r| <= pi/4
 200 // ------------------------------------
 201 //
 202 // This case is more likely than the previous one if one considers
 203 // r to be uniformly distributed in [-pi/4 pi/4].
 204 //
 205 // The required calculation is either
 206 //
 207 //      tan(r + c)  =  tan(r)  +  correction,  or
 208 //     -cot(r + c)  = -cot(r)  +  correction.
 209 //
 210 // Specifically,
 211 //
 212 //      tan(r + c) =  tan(r) + c tan'(r)  + O(c^2)
 213 //                 =  tan(r) + c sec^2(r) + O(c^2)
 214 //                 =  tan(r) + c SEC_sq     ...accurately
 215 //                as long as SEC_sq approximates sec^2(r)
 216 //                to, say, 5 bits or so.
 217 //
 218 // Similarly,
 219 //
 220 //     -cot(r + c) = -cot(r) - c cot'(r)  + O(c^2)
 221 //                 = -cot(r) + c csc^2(r) + O(c^2)
 222 //                 = -cot(r) + c CSC_sq     ...accurately
 223 //                as long as CSC_sq approximates csc^2(r)
 224 //                to, say, 5 bits or so.
 225 //
 226 // We therefore concentrate on accurately calculating tan(r) and
 227 // cot(r) for a working-precision number r, |r| <= pi/4 to within
 228 // 0.1% or so.
 229 //
 230 // We will employ a table-driven approach. Let
 231 //
 232 //      r = sgn_r * 2^k * 1.b_1 b_2 ... b_5 ... b_63
 233 //        = sgn_r * ( B + x )
 234 //
 235 // where
 236 //
 237 //      B = 2^k * 1.b_1 b_2 ... b_5 1
 238 //      x = |r| - B
 239 //
 240 // Now,
 241 //                   tan(B)  +   tan(x)
 242 //      tan( B + x ) =  ------------------------
 243 //                   1 -  tan(B)*tan(x)
 244 //
 245 //               /                         \
 246 //               |   tan(B)  +   tan(x)          |
 247
 248 //      = tan(B) +  | ------------------------ - tan(B) |
 249 //               |     1 -  tan(B)*tan(x)          |
 250 //               \                         /
 251 //
 252 //                 sec^2(B) * tan(x)
 253 //      = tan(B) + ------------------------
 254 //                 1 -  tan(B)*tan(x)
 255 //
 256 //                (1/[sin(B)*cos(B)]) * tan(x)
 257 //      = tan(B) + --------------------------------
 258 //                      cot(B)  -  tan(x)
 259 //
 260 //
 261 // Clearly, the values of tan(B), cot(B) and 1/(sin(B)*cos(B)) are
 262 // calculated beforehand and stored in a table. Since
 263 //
 264 //      |x| <= 2^k * 2^(-6)  <= 2^(-7)  (because k = -1, -2)
 265 //
 266 // a very short polynomial will be sufficient to approximate tan(x)
 267 // accurately. The details involved in computing the last expression
 268 // will be given in the next section on algorithm description.
 269 //
 270 //
 271 // Now, we turn to the case where cot( B + x ) is needed.
 272 //
 273 //
 274 //                   1 - tan(B)*tan(x)
 275 //      cot( B + x ) =  ------------------------
 276 //                   tan(B)  +  tan(x)
 277 //
 278 //               /                           \
 279 //               |   1 - tan(B)*tan(x)              |
 280
 281 //      = cot(B) +  | ----------------------- - cot(B) |
 282 //               |     tan(B)  +  tan(x)            |
 283 //               \                           /
 284 //
 285 //               [tan(B) + cot(B)] * tan(x)
 286 //      = cot(B) - ----------------------------
 287 //                   tan(B)  +  tan(x)
 288 //
 289 //                (1/[sin(B)*cos(B)]) * tan(x)
 290 //      = cot(B) - --------------------------------
 291 //                      tan(B)  +  tan(x)
 292 //
 293 //
 294 // Note that the values of tan(B), cot(B) and 1/(sin(B)*cos(B)) that
 295 // are needed are the same set of values needed in the previous
 296 // case.
 297 //
 298 // Finally, we can put all the ingredients together as follows:
 299 //
 300 //      Arg = N * pi/2 +  r + c          ...accurately
 301 //
 302 //      tan(Arg) =  tan(r) + correction    if N is even;
 303 //               = -cot(r) + correction    otherwise.
 304 //
 305 // For Cases 2 and 4,
 306 //
 307 //     Case 2:
 308 //     tan(Arg) =  tan(r + c) = r + c + r^3/3           N even
 309 //              = -cot(r + c) = -1/(r+c) + r/3           N odd
 310 //     Case 4:
 311 //     tan(Arg) =  tan(r + c) = r + c + r^3/3 + 2r^5/15  N even
 312 //              = -cot(r + c) = -1/(r+c) + r/3 + r^3/45  N odd
 313 //
 314 //
 315 // For Cases 1 and 3,
 316 //
 317 //     Case small_r: |r| < 2^(-2)
 318 //
 319 //      tan(Arg) =  r + P1_1 r^3 + P1_2 r^5 + ... + P1_9 r^19
 320 //                     + c*(1 + r^2)               N even
 321 //
 322 //               = -1/(r+c) + Q1_1 r   + Q1_2 r^3 + ... + Q1_7 r^13
 323 //                     + Q1_1*c                    N odd
 324 //
 325 //     Case normal_r: 2^(-2) <= |r| <= pi/4
 326 //
 327 //      tan(Arg) =  tan(r) + c * sec^2(r)     N even
 328 //               = -cot(r) + c * csc^2(r)     otherwise
 329 //
 330 //     For N even,
 331 //
 332 //      tan(Arg) = tan(r) + c*sec^2(r)
 333 //               = tan( sgn_r * (B+x) ) + c * sec^2(|r|)
 334 //               = sgn_r * ( tan(B+x)  + sgn_r*c*sec^2(|r|) )
 335 //               = sgn_r * ( tan(B+x)  + sgn_r*c*sec^2(B) )
 336 //
 337 // since B approximates |r| to 2^(-6) in relative accuracy.
 338 //
 339 //                 /            (1/[sin(B)*cos(B)]) * tan(x)
 340 //    tan(Arg) = sgn_r * | tan(B) + --------------------------------
 341 //                 \                     cot(B)  -  tan(x)
 342 //                                        \
 343 //                       + CORR  |
 344
 345 //                                     /
 346 // where
 347 //
 348 //    CORR = sgn_r*c*tan(B)*SC_inv(B);  SC_inv(B) = 1/(sin(B)*cos(B)).
 349 //
 350 // For N odd,
 351 //
 352 //      tan(Arg) = -cot(r) + c*csc^2(r)
 353 //               = -cot( sgn_r * (B+x) ) + c * csc^2(|r|)
 354 //               = sgn_r * ( -cot(B+x)  + sgn_r*c*csc^2(|r|) )
 355 //               = sgn_r * ( -cot(B+x)  + sgn_r*c*csc^2(B) )
 356 //
 357 // since B approximates |r| to 2^(-6) in relative accuracy.
 358 //
 359 //                 /            (1/[sin(B)*cos(B)]) * tan(x)
 360 //    tan(Arg) = sgn_r * | -cot(B) + --------------------------------
 361 //                 \                     tan(B)  +  tan(x)
 362 //                                        \
 363 //                       + CORR  |
 364
 365 //                                     /
 366 // where
 367 //
 368 //    CORR = sgn_r*c*cot(B)*SC_inv(B);  SC_inv(B) = 1/(sin(B)*cos(B)).
 369 //
 370 //
 371 // The actual algorithm prescribes how all the mathematical formulas
 372 // are calculated.
 373 //
 374 //
 375 // 2. Algorithmic Description
 376 // ==========================
 377 //
 378 // 2.1 Computation for Cases 2 and 4.
 379 // ----------------------------------
 380 //
 381 // For Case 2, we use two-term polynomials.
 382 //
 383 //    For N even,
 384 //
 385 //    rsq := r * r
 386 //    Poly := c + r * rsq * P1_1
 387 //    Result := r + Poly          ...in user-defined rounding
 388 //
 389 //    For N odd,
 390 //    S_hi  := -frcpa(r)               ...8 bits
 391 //    S_hi  := S_hi + S_hi*(1 + S_hi*r)     ...16 bits
 392 //    S_hi  := S_hi + S_hi*(1 + S_hi*r)     ...32 bits
 393 //    S_hi  := S_hi + S_hi*(1 + S_hi*r)     ...64 bits
 394 //    S_lo  := S_hi*( (1 + S_hi*r) + S_hi*c )
 395 //    ...S_hi + S_lo is -1/(r+c) to extra precision
 396 //    S_lo  := S_lo + Q1_1*r
 397 //
 398 //    Result := S_hi + S_lo     ...in user-defined rounding
 399 //
 400 // For Case 4, we use three-term polynomials
 401 //
 402 //    For N even,
 403 //
 404 //    rsq := r * r
 405 //    Poly := c + r * rsq * (P1_1 + rsq * P1_2)
 406 //    Result := r + Poly          ...in user-defined rounding
 407 //
 408 //    For N odd,
 409 //    S_hi  := -frcpa(r)               ...8 bits
 410 //    S_hi  := S_hi + S_hi*(1 + S_hi*r)     ...16 bits
 411 //    S_hi  := S_hi + S_hi*(1 + S_hi*r)     ...32 bits
 412 //    S_hi  := S_hi + S_hi*(1 + S_hi*r)     ...64 bits
 413 //    S_lo  := S_hi*( (1 + S_hi*r) + S_hi*c )
 414 //    ...S_hi + S_lo is -1/(r+c) to extra precision
 415 //    rsq   := r * r
 416 //    P      := Q1_1 + rsq*Q1_2
 417 //    S_lo  := S_lo + r*P
 418 //
 419 //    Result := S_hi + S_lo     ...in user-defined rounding
 420 //
 421 //
 422 // Note that the coefficients P1_1, P1_2, Q1_1, and Q1_2 are
 423 // the same as those used in the small_r case of Cases 1 and 3
 424 // below.
 425 //
 426 //
 427 // 2.2 Computation for Cases 1 and 3.
 428 // ----------------------------------
 429 // This is further divided into the case of small_r,
 430 // where |r| < 2^(-2), and the case of normal_r, where |r| lies between
 431 // 2^(-2) and pi/4.
 432 //
 433 // Algorithm for the case of small_r
 434 // ---------------------------------
 435 //
 436 // For N even,
 437 //      rsq   := r * r
 438 //      Poly1 := rsq*(P1_1 + rsq*(P1_2 + rsq*P1_3))
 439 //      r_to_the_8    := rsq * rsq
 440 //      r_to_the_8    := r_to_the_8 * r_to_the_8
 441 //      Poly2 := P1_4 + rsq*(P1_5 + rsq*(P1_6 + ... rsq*P1_9))
 442 //      CORR  := c * ( 1 + rsq )
 443 //      Poly  := Poly1 + r_to_the_8*Poly2
 444 //      Poly := r*Poly + CORR
 445 //      Result := r + Poly     ...in user-defined rounding
 446 //      ...note that Poly1 and r_to_the_8 can be computed in parallel
 447 //      ...with Poly2 (Poly1 is intentionally set to be much
 448 //      ...shorter than Poly2 so that r_to_the_8 and CORR can be hidden)
 449 //
 450 // For N odd,
 451 //      S_hi  := -frcpa(r)               ...8 bits
 452 //      S_hi  := S_hi + S_hi*(1 + S_hi*r)     ...16 bits
 453 //      S_hi  := S_hi + S_hi*(1 + S_hi*r)     ...32 bits
 454 //      S_hi  := S_hi + S_hi*(1 + S_hi*r)     ...64 bits
 455 //      S_lo  := S_hi*( (1 + S_hi*r) + S_hi*c )
 456 //      ...S_hi + S_lo is -1/(r+c) to extra precision
 457 //      S_lo  := S_lo + Q1_1*c
 458 //
 459 //      ...S_hi and S_lo are computed in parallel with
 460 //      ...the following
 461 //      rsq := r*r
 462 //      P   := Q1_1 + rsq*(Q1_2 + rsq*(Q1_3 + ... + rsq*Q1_7))
 463 //
 464 //      Poly :=  r*P + S_lo
 465 //      Result :=  S_hi  +  Poly      ...in user-defined rounding
 466 //
 467 //
 468 // Algorithm for the case of normal_r
 469 // ----------------------------------
 470 //
 471 // Here, we first consider the computation of tan( r + c ). As
 472 // presented in the previous section,
 473 //
 474 //      tan( r + c )  =  tan(r) + c * sec^2(r)
 475 //                 =  sgn_r * [ tan(B+x) + CORR ]
 476 //      CORR = sgn_r * c * tan(B) * 1/[sin(B)*cos(B)]
 477 //
 478 // because sec^2(r) = sec^(|r|), and B approximate |r| to 6.5 bits.
 479 //
 480 //      tan( r + c ) =
 481 //           /           (1/[sin(B)*cos(B)]) * tan(x)
 482 //      sgn_r * | tan(B) + --------------------------------  +
 483 //           \                     cot(B)  -  tan(x)
 484 //                                \
 485 //                          CORR  |
 486
 487 //                                /
 488 //
 489 // The values of tan(B), cot(B) and 1/(sin(B)*cos(B)) are
 490 // calculated beforehand and stored in a table. Specifically,
 491 // the table values are
 492 //
 493 //      tan(B)             as  T_hi  +  T_lo;
 494 //      cot(B)             as  C_hi  +  C_lo;
 495 //      1/[sin(B)*cos(B)]  as  SC_inv
 496 //
 497 // T_hi, C_hi are in  double-precision  memory format;
 498 // T_lo, C_lo are in  single-precision  memory format;
 499 // SC_inv     is  in extended-precision memory format.
 500 //
 501 // The value of tan(x) will be approximated by a short polynomial of
 502 // the form
 503 //
 504 //      tan(x)  as  x  +  x * P, where
 505 //           P  =   x^2 * (P2_1 + x^2 * (P2_2 + x^2 * P2_3))
 506 //
 507 // Because |x| <= 2^(-7), cot(B) - x approximates cot(B) - tan(x)
 508 // to a relative accuracy better than 2^(-20). Thus, a good
 509 // initial guess of 1/( cot(B) - tan(x) ) to initiate the iterative
 510 // division is:
 511 //
 512 //      1/(cot(B) - tan(x))      is approximately
 513 //      1/(cot(B) -   x)         is
 514 //      tan(B)/(1 - x*tan(B))    is approximately
 515 //      T_hi / ( 1 - T_hi * x )  is approximately
 516 //
 517 //      T_hi * [ 1 + (Thi * x) + (T_hi * x)^2 ]
 518 //
 519 // The calculation of tan(r+c) therefore proceed as follows:
 520 //
 521 //      Tx     := T_hi * x
 522 //      xsq     := x * x
 523 //
 524 //      V_hi     := T_hi*(1 + Tx*(1 + Tx))
 525 //      P     := xsq * (P1_1 + xsq*(P1_2 + xsq*P1_3))
 526 //      ...V_hi serves as an initial guess of 1/(cot(B) - tan(x))
 527 //         ...good to about 20 bits of accuracy
 528 //
 529 //      tanx     := x + x*P
 530 //      D     := C_hi - tanx
 531 //      ...D is a double precision denominator: cot(B) - tan(x)
 532 //
 533 //      V_hi     := V_hi + V_hi*(1 - V_hi*D)
 534 //      ....V_hi approximates 1/(cot(B)-tan(x)) to 40 bits
 535 //
 536 //      V_lo     := V_hi * ( [ (1 - V_hi*C_hi) + V_hi*tanx ]
 537 //                           - V_hi*C_lo )   ...observe all order
 538 //         ...V_hi + V_lo approximates 1/(cot(B) - tan(x))
 539 //      ...to extra accuracy
 540 //
 541 //      ...               SC_inv(B) * (x + x*P)
 542 //      ...   tan(B) +      ------------------------- + CORR
 543 //         ...                cot(B) - (x + x*P)
 544 //      ...
 545 //      ... = tan(B) + SC_inv(B)*(x + x*P)*(V_hi + V_lo) + CORR
 546 //      ...
 547 //
 548 //      Sx     := SC_inv * x
 549 //      CORR     := sgn_r * c * SC_inv * T_hi
 550 //
 551 //      ...put the ingredients together to compute
 552 //      ...               SC_inv(B) * (x + x*P)
 553 //      ...   tan(B) +      ------------------------- + CORR
 554 //         ...                cot(B) - (x + x*P)
 555 //      ...
 556 //      ... = tan(B) + SC_inv(B)*(x + x*P)*(V_hi + V_lo) + CORR
 557 //      ...
 558 //      ... = T_hi + T_lo + CORR +
 559 //      ...    Sx * V_hi + Sx * V_lo + Sx * P *(V_hi + V_lo)
 560 //
 561 //      CORR := CORR + T_lo
 562 //      tail := V_lo + P*(V_hi + V_lo)
 563 //         tail := Sx * tail  +  CORR
 564 //      tail := Sx * V_hi  +  tail
 565 //         T_hi := sgn_r * T_hi
 566 //
 567 //         ...T_hi + sgn_r*tail  now approximate
 568 //      ...sgn_r*(tan(B+x) + CORR) accurately
 569 //
 570 //      Result :=  T_hi + sgn_r*tail  ...in user-defined
 571 //                           ...rounding control
 572 //      ...It is crucial that independent paths be fully
 573 //      ...exploited for performance's sake.
 574 //
 575 //
 576 // Next, we consider the computation of -cot( r + c ). As
 577 // presented in the previous section,
 578 //
 579 //        -cot( r + c )  =  -cot(r) + c * csc^2(r)
 580 //                 =  sgn_r * [ -cot(B+x) + CORR ]
 581 //      CORR = sgn_r * c * cot(B) * 1/[sin(B)*cos(B)]
 582 //
 583 // because csc^2(r) = csc^(|r|), and B approximate |r| to 6.5 bits.
 584 //
 585 //        -cot( r + c ) =
 586 //           /             (1/[sin(B)*cos(B)]) * tan(x)
 587 //      sgn_r * | -cot(B) + --------------------------------  +
 588 //           \                     tan(B)  +  tan(x)
 589 //                                \
 590 //                          CORR  |
 591
 592 //                                /
 593 //
 594 // The values of tan(B), cot(B) and 1/(sin(B)*cos(B)) are
 595 // calculated beforehand and stored in a table. Specifically,
 596 // the table values are
 597 //
 598 //      tan(B)             as  T_hi  +  T_lo;
 599 //      cot(B)             as  C_hi  +  C_lo;
 600 //      1/[sin(B)*cos(B)]  as  SC_inv
 601 //
 602 // T_hi, C_hi are in  double-precision  memory format;
 603 // T_lo, C_lo are in  single-precision  memory format;
 604 // SC_inv     is  in extended-precision memory format.
 605 //
 606 // The value of tan(x) will be approximated by a short polynomial of
 607 // the form
 608 //
 609 //      tan(x)  as  x  +  x * P, where
 610 //           P  =   x^2 * (P2_1 + x^2 * (P2_2 + x^2 * P2_3))
 611 //
 612 // Because |x| <= 2^(-7), tan(B) + x approximates tan(B) + tan(x)
 613 // to a relative accuracy better than 2^(-18). Thus, a good
 614 // initial guess of 1/( tan(B) + tan(x) ) to initiate the iterative
 615 // division is:
 616 //
 617 //      1/(tan(B) + tan(x))      is approximately
 618 //      1/(tan(B) +   x)         is
 619 //      cot(B)/(1 + x*cot(B))    is approximately
 620 //      C_hi / ( 1 + C_hi * x )  is approximately
 621 //
 622 //      C_hi * [ 1 - (C_hi * x) + (C_hi * x)^2 ]
 623 //
 624 // The calculation of -cot(r+c) therefore proceed as follows:
 625 //
 626 //      Cx     := C_hi * x
 627 //      xsq     := x * x
 628 //
 629 //      V_hi     := C_hi*(1 - Cx*(1 - Cx))
 630 //      P     := xsq * (P1_1 + xsq*(P1_2 + xsq*P1_3))
 631 //      ...V_hi serves as an initial guess of 1/(tan(B) + tan(x))
 632 //         ...good to about 18 bits of accuracy
 633 //
 634 //      tanx     := x + x*P
 635 //      D     := T_hi + tanx
 636 //      ...D is a double precision denominator: tan(B) + tan(x)
 637 //
 638 //      V_hi     := V_hi + V_hi*(1 - V_hi*D)
 639 //      ....V_hi approximates 1/(tan(B)+tan(x)) to 40 bits
 640 //
 641 //      V_lo     := V_hi * ( [ (1 - V_hi*T_hi) - V_hi*tanx ]
 642 //                           - V_hi*T_lo )   ...observe all order
 643 //         ...V_hi + V_lo approximates 1/(tan(B) + tan(x))
 644 //      ...to extra accuracy
 645 //
 646 //      ...               SC_inv(B) * (x + x*P)
 647 //      ...  -cot(B) +      ------------------------- + CORR
 648 //         ...                tan(B) + (x + x*P)
 649 //      ...
 650 //      ... =-cot(B) + SC_inv(B)*(x + x*P)*(V_hi + V_lo) + CORR
 651 //      ...
 652 //
 653 //      Sx     := SC_inv * x
 654 //      CORR     := sgn_r * c * SC_inv * C_hi
 655 //
 656 //      ...put the ingredients together to compute
 657 //      ...               SC_inv(B) * (x + x*P)
 658 //      ...  -cot(B) +      ------------------------- + CORR
 659 //         ...                tan(B) + (x + x*P)
 660 //      ...
 661 //      ... =-cot(B) + SC_inv(B)*(x + x*P)*(V_hi + V_lo) + CORR
 662 //      ...
 663 //      ... =-C_hi - C_lo + CORR +
 664 //      ...    Sx * V_hi + Sx * V_lo + Sx * P *(V_hi + V_lo)
 665 //
 666 //      CORR := CORR - C_lo
 667 //      tail := V_lo + P*(V_hi + V_lo)
 668 //         tail := Sx * tail  +  CORR
 669 //      tail := Sx * V_hi  +  tail
 670 //         C_hi := -sgn_r * C_hi
 671 //
 672 //         ...C_hi + sgn_r*tail now approximates
 673 //      ...sgn_r*(-cot(B+x) + CORR) accurately
 674 //
 675 //      Result :=  C_hi + sgn_r*tail   in user-defined rounding control
 676 //      ...It is crucial that independent paths be fully
 677 //      ...exploited for performance's sake.
 678 //
 679 // 3. Implementation Notes
 680 // =======================
 681 //
 682 //   Table entries T_hi, T_lo; C_hi, C_lo; SC_inv
 683 //
 684 //   Recall that 2^(-2) <= |r| <= pi/4;
 685 //
 686 //      r = sgn_r * 2^k * 1.b_1 b_2 ... b_63
 687 //
 688 //   and
 689 //
 690 //        B = 2^k * 1.b_1 b_2 b_3 b_4 b_5 1
 691 //
 692 //   Thus, for k = -2, possible values of B are
 693 //
 694 //          B = 2^(-2) * ( 1 + index/32  +  1/64 ),
 695 //      index ranges from 0 to 31
 696 //
 697 //   For k = -1, however, since |r| <= pi/4 = 0.78...
 698 //   possible values of B are
 699 //
 700 //        B = 2^(-1) * ( 1 + index/32  +  1/64 )
 701 //      index ranges from 0 to 19.
 702 //
 703 //
 704
 705 RODATA
 706 .align 16
 707
 708 LOCAL_OBJECT_START(TANL_BASE_CONSTANTS)
 709
 710 tanl_table_1:
 711 data8    0xA2F9836E4E44152A, 0x00003FFE // two_by_pi
 712 data8    0xC84D32B0CE81B9F1, 0x00004016 // P_0
 713 data8    0xC90FDAA22168C235, 0x00003FFF // P_1
 714 data8    0xECE675D1FC8F8CBB, 0x0000BFBD // P_2
 715 data8    0xB7ED8FBBACC19C60, 0x0000BF7C // P_3
 716 LOCAL_OBJECT_END(TANL_BASE_CONSTANTS)
 717
 718 LOCAL_OBJECT_START(tanl_table_2)
 719 data8    0xC90FDAA22168C234, 0x00003FFE // PI_BY_4
 720 data8    0xA397E5046EC6B45A, 0x00003FE7 // Inv_P_0
 721 data8    0x8D848E89DBD171A1, 0x0000BFBF // d_1
 722 data8    0xD5394C3618A66F8E, 0x0000BF7C // d_2
 723 data4    0x3E800000 // two**-2
 724 data4    0xBE800000 // -two**-2
 725 data4    0x00000000 // pad
 726 data4    0x00000000 // pad
 727 LOCAL_OBJECT_END(tanl_table_2)
 728
 729 LOCAL_OBJECT_START(tanl_table_p1)
 730 data8    0xAAAAAAAAAAAAAABD, 0x00003FFD // P1_1
 731 data8    0x8888888888882E6A, 0x00003FFC // P1_2
 732 data8    0xDD0DD0DD0F0177B6, 0x00003FFA // P1_3
 733 data8    0xB327A440646B8C6D, 0x00003FF9 // P1_4
 734 data8    0x91371B251D5F7D20, 0x00003FF8 // P1_5
 735 data8    0xEB69A5F161C67914, 0x00003FF6 // P1_6
 736 data8    0xBEDD37BE019318D2, 0x00003FF5 // P1_7
 737 data8    0x9979B1463C794015, 0x00003FF4 // P1_8
 738 data8    0x8EBD21A38C6EB58A, 0x00003FF3 // P1_9
 739 LOCAL_OBJECT_END(tanl_table_p1)
 740
 741 LOCAL_OBJECT_START(tanl_table_q1)
 742 data8    0xAAAAAAAAAAAAAAB4, 0x00003FFD // Q1_1
 743 data8    0xB60B60B60B5FC93E, 0x00003FF9 // Q1_2
 744 data8    0x8AB355E00C9BBFBF, 0x00003FF6 // Q1_3
 745 data8    0xDDEBBC89CBEE3D4C, 0x00003FF2 // Q1_4
 746 data8    0xB3548A685F80BBB6, 0x00003FEF // Q1_5
 747 data8    0x913625604CED5BF1, 0x00003FEC // Q1_6
 748 data8    0xF189D95A8EE92A83, 0x00003FE8 // Q1_7
 749 LOCAL_OBJECT_END(tanl_table_q1)
 750
 751 LOCAL_OBJECT_START(tanl_table_p2)
 752 data8    0xAAAAAAAAAAAB362F, 0x00003FFD // P2_1
 753 data8    0x88888886E97A6097, 0x00003FFC // P2_2
 754 data8    0xDD108EE025E716A1, 0x00003FFA // P2_3
 755 LOCAL_OBJECT_END(tanl_table_p2)
 756
 757 LOCAL_OBJECT_START(tanl_table_tm2)
 758 //
 759 //  Entries T_hi   double-precision memory format
 760 //  Index = 0,1,...,31  B = 2^(-2)*(1+Index/32+1/64)
 761 //  Entries T_lo  single-precision memory format
 762 //  Index = 0,1,...,31  B = 2^(-2)*(1+Index/32+1/64)
 763 //
 764 data8 0x3FD09BC362400794
 765 data4 0x23A05C32, 0x00000000
 766 data8 0x3FD124A9DFFBC074
 767 data4 0x240078B2, 0x00000000
 768 data8 0x3FD1AE235BD4920F
 769 data4 0x23826B8E, 0x00000000
 770 data8 0x3FD2383515E2701D
 771 data4 0x22D31154, 0x00000000
 772 data8 0x3FD2C2E463739C2D
 773 data4 0x2265C9E2, 0x00000000
 774 data8 0x3FD34E36AFEEA48B
 775 data4 0x245C05EB, 0x00000000
 776 data8 0x3FD3DA317DBB35D1
 777 data4 0x24749F2D, 0x00000000
 778 data8 0x3FD466DA67321619
 779 data4 0x2462CECE, 0x00000000
 780 data8 0x3FD4F4371F94A4D5
 781 data4 0x246D0DF1, 0x00000000
 782 data8 0x3FD5824D740C3E6D
 783 data4 0x240A85B5, 0x00000000
 784 data8 0x3FD611234CB1E73D
 785 data4 0x23F96E33, 0x00000000
 786 data8 0x3FD6A0BEAD9EA64B
 787 data4 0x247C5393, 0x00000000
 788 data8 0x3FD73125B804FD01
 789 data4 0x241F3B29, 0x00000000
 790 data8 0x3FD7C25EAB53EE83
 791 data4 0x2479989B, 0x00000000
 792 data8 0x3FD8546FE6640EED
 793 data4 0x23B343BC, 0x00000000
 794 data8 0x3FD8E75FE8AF1892
 795 data4 0x241454D1, 0x00000000
 796 data8 0x3FD97B3553928BDA
 797 data4 0x238613D9, 0x00000000
 798 data8 0x3FDA0FF6EB9DE4DE
 799 data4 0x22859FA7, 0x00000000
 800 data8 0x3FDAA5AB99ECF92D
 801 data4 0x237A6D06, 0x00000000
 802 data8 0x3FDB3C5A6D8F1796
 803 data4 0x23952F6C, 0x00000000
 804 data8 0x3FDBD40A9CFB8BE4
 805 data4 0x2280FC95, 0x00000000
 806 data8 0x3FDC6CC387943100
 807 data4 0x245D2EC0, 0x00000000
 808 data8 0x3FDD068CB736C500
 809 data4 0x23C4AD7D, 0x00000000
 810 data8 0x3FDDA16DE1DDBC31
 811 data4 0x23D076E6, 0x00000000
 812 data8 0x3FDE3D6EEB515A93
 813 data4 0x244809A6, 0x00000000
 814 data8 0x3FDEDA97E6E9E5F1
 815 data4 0x220856C8, 0x00000000
 816 data8 0x3FDF78F11963CE69
 817 data4 0x244BE993, 0x00000000
 818 data8 0x3FE00C417D635BCE
 819 data4 0x23D21799, 0x00000000
 820 data8 0x3FE05CAB1C302CD3
 821 data4 0x248A1B1D, 0x00000000
 822 data8 0x3FE0ADB9DB6A1FA0
 823 data4 0x23D53E33, 0x00000000
 824 data8 0x3FE0FF724A20BA81
 825 data4 0x24DB9ED5, 0x00000000
 826 data8 0x3FE151D9153FA6F5
 827 data4 0x24E9E451, 0x00000000
 828 LOCAL_OBJECT_END(tanl_table_tm2)
 829
 830 LOCAL_OBJECT_START(tanl_table_tm1)
 831 //
 832 //  Entries T_hi   double-precision memory format
 833 //  Index = 0,1,...,19  B = 2^(-1)*(1+Index/32+1/64)
 834 //  Entries T_lo  single-precision memory format
 835 //  Index = 0,1,...,19  B = 2^(-1)*(1+Index/32+1/64)
 836 //
 837 data8 0x3FE1CEC4BA1BE39E
 838 data4 0x24B60F9E, 0x00000000
 839 data8 0x3FE277E45ABD9B2D
 840 data4 0x248C2474, 0x00000000
 841 data8 0x3FE324180272B110
 842 data4 0x247B8311, 0x00000000
 843 data8 0x3FE3D38B890E2DF0
 844 data4 0x24C55751, 0x00000000
 845 data8 0x3FE4866D46236871
 846 data4 0x24E5BC34, 0x00000000
 847 data8 0x3FE53CEE45E044B0
 848 data4 0x24001BA4, 0x00000000
 849 data8 0x3FE5F74282EC06E4
 850 data4 0x24B973DC, 0x00000000
 851 data8 0x3FE6B5A125DF43F9
 852 data4 0x24895440, 0x00000000
 853 data8 0x3FE77844CAFD348C
 854 data4 0x240021CA, 0x00000000
 855 data8 0x3FE83F6BCEED6B92
 856 data4 0x24C45372, 0x00000000
 857 data8 0x3FE90B58A34F3665
 858 data4 0x240DAD33, 0x00000000
 859 data8 0x3FE9DC522C1E56B4
 860 data4 0x24F846CE, 0x00000000
 861 data8 0x3FEAB2A427041578
 862 data4 0x2323FB6E, 0x00000000
 863 data8 0x3FEB8E9F9DD8C373
 864 data4 0x24B3090B, 0x00000000
 865 data8 0x3FEC709B65C9AA7B
 866 data4 0x2449F611, 0x00000000
 867 data8 0x3FED58F4ACCF8435
 868 data4 0x23616A7E, 0x00000000
 869 data8 0x3FEE480F97635082
 870 data4 0x24C2FEAE, 0x00000000
 871 data8 0x3FEF3E57F0ACC544
 872 data4 0x242CE964, 0x00000000
 873 data8 0x3FF01E20F7E06E4B
 874 data4 0x2480D3EE, 0x00000000
 875 data8 0x3FF0A1258A798A69
 876 data4 0x24DB8967, 0x00000000
 877 LOCAL_OBJECT_END(tanl_table_tm1)
 878
 879 LOCAL_OBJECT_START(tanl_table_cm2)
 880 //
 881 //  Entries C_hi   double-precision memory format
 882 //  Index = 0,1,...,31  B = 2^(-2)*(1+Index/32+1/64)
 883 //  Entries C_lo  single-precision memory format
 884 //  Index = 0,1,...,31  B = 2^(-2)*(1+Index/32+1/64)
 885 //
 886 data8 0x400ED3E2E63EFBD0
 887 data4 0x259D94D4, 0x00000000
 888 data8 0x400DDDB4C515DAB5
 889 data4 0x245F0537, 0x00000000
 890 data8 0x400CF57ABE19A79F
 891 data4 0x25D4EA9F, 0x00000000
 892 data8 0x400C1A06D15298ED
 893 data4 0x24AE40A0, 0x00000000
 894 data8 0x400B4A4C164B2708
 895 data4 0x25A5AAB6, 0x00000000
 896 data8 0x400A855A5285B068
 897 data4 0x25524F18, 0x00000000
 898 data8 0x4009CA5A3FFA549F
 899 data4 0x24C999C0, 0x00000000
 900 data8 0x4009188A646AF623
 901 data4 0x254FD801, 0x00000000
 902 data8 0x40086F3C6084D0E7
 903 data4 0x2560F5FD, 0x00000000
 904 data8 0x4007CDD2A29A76EE
 905 data4 0x255B9D19, 0x00000000
 906 data8 0x400733BE6C8ECA95
 907 data4 0x25CB021B, 0x00000000
 908 data8 0x4006A07E1F8DDC52
 909 data4 0x24AB4722, 0x00000000
 910 data8 0x4006139BC298AD58
 911 data4 0x252764E2, 0x00000000
 912 data8 0x40058CABBAD7164B
 913 data4 0x24DAF5DB, 0x00000000
 914 data8 0x40050B4BAE31A5D3
 915 data4 0x25EA20F4, 0x00000000
 916 data8 0x40048F2189F85A8A
 917 data4 0x2583A3E8, 0x00000000
 918 data8 0x400417DAA862380D
 919 data4 0x25DCC4CC, 0x00000000
 920 data8 0x4003A52B1088FCFE
 921 data4 0x2430A492, 0x00000000
 922 data8 0x400336CCCD3527D5
 923 data4 0x255F77CF, 0x00000000
 924 data8 0x4002CC7F5760766D
 925 data4 0x25DA0BDA, 0x00000000
 926 data8 0x4002660711CE02E3
 927 data4 0x256FF4A2, 0x00000000
 928 data8 0x4002032CD37BBE04
 929 data4 0x25208AED, 0x00000000
 930 data8 0x4001A3BD7F050775
 931 data4 0x24B72DD6, 0x00000000
 932 data8 0x40014789A554848A
 933 data4 0x24AB4DAA, 0x00000000
 934 data8 0x4000EE65323E81B7
 935 data4 0x2584C440, 0x00000000
 936 data8 0x4000982721CF1293
 937 data4 0x25C9428D, 0x00000000
 938 data8 0x400044A93D415EEB
 939 data4 0x25DC8482, 0x00000000
 940 data8 0x3FFFE78FBD72C577
 941 data4 0x257F5070, 0x00000000
 942 data8 0x3FFF4AC375EFD28E
 943 data4 0x23EBBF7A, 0x00000000
 944 data8 0x3FFEB2AF60B52DDE
 945 data4 0x22EECA07, 0x00000000
 946 data8 0x3FFE1F1935204180
 947 data4 0x24191079, 0x00000000
 948 data8 0x3FFD8FCA54F7E60A
 949 data4 0x248D3058, 0x00000000
 950 LOCAL_OBJECT_END(tanl_table_cm2)
 951
 952 LOCAL_OBJECT_START(tanl_table_cm1)
 953 //
 954 //  Entries C_hi   double-precision memory format
 955 //  Index = 0,1,...,19  B = 2^(-1)*(1+Index/32+1/64)
 956 //  Entries C_lo  single-precision memory format
 957 //  Index = 0,1,...,19  B = 2^(-1)*(1+Index/32+1/64)
 958 //
 959 data8 0x3FFCC06A79F6FADE
 960 data4 0x239C7886, 0x00000000
 961 data8 0x3FFBB91F891662A6
 962 data4 0x250BD191, 0x00000000
 963 data8 0x3FFABFB6529F155D
 964 data4 0x256CC3E6, 0x00000000
 965 data8 0x3FF9D3002E964AE9
 966 data4 0x250843E3, 0x00000000
 967 data8 0x3FF8F1EF89DCB383
 968 data4 0x2277C87E, 0x00000000
 969 data8 0x3FF81B937C87DBD6
 970 data4 0x256DA6CF, 0x00000000
 971 data8 0x3FF74F141042EDE4
 972 data4 0x2573D28A, 0x00000000
 973 data8 0x3FF68BAF1784B360
 974 data4 0x242E489A, 0x00000000
 975 data8 0x3FF5D0B57C923C4C
 976 data4 0x2532D940, 0x00000000
 977 data8 0x3FF51D88F418EF20
 978 data4 0x253C7DD6, 0x00000000
 979 data8 0x3FF4719A02F88DAE
 980 data4 0x23DB59BF, 0x00000000
 981 data8 0x3FF3CC6649DA0788
 982 data4 0x252B4756, 0x00000000
 983 data8 0x3FF32D770B980DB8
 984 data4 0x23FE585F, 0x00000000
 985 data8 0x3FF2945FE56C987A
 986 data4 0x25378A63, 0x00000000
 987 data8 0x3FF200BDB16523F6
 988 data4 0x247BB2E0, 0x00000000
 989 data8 0x3FF172358CE27778
 990 data4 0x24446538, 0x00000000
 991 data8 0x3FF0E873FDEFE692
 992 data4 0x2514638F, 0x00000000
 993 data8 0x3FF0632C33154062
 994 data4 0x24A7FC27, 0x00000000
 995 data8 0x3FEFC42EB3EF115F
 996 data4 0x248FD0FE, 0x00000000
 997 data8 0x3FEEC9E8135D26F6
 998 data4 0x2385C719, 0x00000000
 999 LOCAL_OBJECT_END(tanl_table_cm1)
1000
1001 LOCAL_OBJECT_START(tanl_table_scim2)
1002 //
1003 //  Entries SC_inv in Swapped IEEE format (extended)
1004 //  Index = 0,1,...,31  B = 2^(-2)*(1+Index/32+1/64)
1005 //
1006 data8    0x839D6D4A1BF30C9E, 0x00004001
1007 data8    0x80092804554B0EB0, 0x00004001
1008 data8    0xF959F94CA1CF0DE9, 0x00004000
1009 data8    0xF3086BA077378677, 0x00004000
1010 data8    0xED154515CCD4723C, 0x00004000
1011 data8    0xE77909441C27CF25, 0x00004000
1012 data8    0xE22D037D8DDACB88, 0x00004000
1013 data8    0xDD2B2D8A89C73522, 0x00004000
1014 data8    0xD86E1A23BB2C1171, 0x00004000
1015 data8    0xD3F0E288DFF5E0F9, 0x00004000
1016 data8    0xCFAF16B1283BEBD5, 0x00004000
1017 data8    0xCBA4AFAA0D88DD53, 0x00004000
1018 data8    0xC7CE03CCCA67C43D, 0x00004000
1019 data8    0xC427BC820CA0DDB0, 0x00004000
1020 data8    0xC0AECD57F13D8CAB, 0x00004000
1021 data8    0xBD606C3871ECE6B1, 0x00004000
1022 data8    0xBA3A0A96A44C4929, 0x00004000
1023 data8    0xB7394F6FE5CCCEC1, 0x00004000
1024 data8    0xB45C12039637D8BC, 0x00004000
1025 data8    0xB1A0552892CB051B, 0x00004000
1026 data8    0xAF04432B6BA2FFD0, 0x00004000
1027 data8    0xAC862A237221235F, 0x00004000
1028 data8    0xAA2478AF5F00A9D1, 0x00004000
1029 data8    0xA7DDBB0C81E082BF, 0x00004000
1030 data8    0xA5B0987D45684FEE, 0x00004000
1031 data8    0xA39BD0F5627A8F53, 0x00004000
1032 data8    0xA19E3B036EC5C8B0, 0x00004000
1033 data8    0x9FB6C1F091CD7C66, 0x00004000
1034 data8    0x9DE464101FA3DF8A, 0x00004000
1035 data8    0x9C263139A8F6B888, 0x00004000
1036 data8    0x9A7B4968C27B0450, 0x00004000
1037 data8    0x98E2DB7E5EE614EE, 0x00004000
1038 LOCAL_OBJECT_END(tanl_table_scim2)
1039
1040 LOCAL_OBJECT_START(tanl_table_scim1)
1041 //
1042 //  Entries SC_inv in Swapped IEEE format (extended)
1043 //  Index = 0,1,...,19  B = 2^(-1)*(1+Index/32+1/64)
1044 //
1045 data8    0x969F335C13B2B5BA, 0x00004000
1046 data8    0x93D446D9D4C0F548, 0x00004000
1047 data8    0x9147094F61B798AF, 0x00004000
1048 data8    0x8EF317CC758787AC, 0x00004000
1049 data8    0x8CD498B3B99EEFDB, 0x00004000
1050 data8    0x8AE82A7DDFF8BC37, 0x00004000
1051 data8    0x892AD546E3C55D42, 0x00004000
1052 data8    0x8799FEA9D15573C1, 0x00004000
1053 data8    0x86335F88435A4B4C, 0x00004000
1054 data8    0x84F4FB6E3E93A87B, 0x00004000
1055 data8    0x83DD195280A382FB, 0x00004000
1056 data8    0x82EA3D7FA4CB8C9E, 0x00004000
1057 data8    0x821B247C6861D0A8, 0x00004000
1058 data8    0x816EBED163E8D244, 0x00004000
1059 data8    0x80E42D9127E4CFC6, 0x00004000
1060 data8    0x807ABF8D28E64AFD, 0x00004000
1061 data8    0x8031EF26863B4FD8, 0x00004000
1062 data8    0x800960ADAE8C11FD, 0x00004000
1063 data8    0x8000E1475FDBEC21, 0x00004000
1064 data8    0x80186650A07791FA, 0x00004000
1065 LOCAL_OBJECT_END(tanl_table_scim1)
1066
1067 Arg                 = f8
1068 Save_Norm_Arg       = f8        // For input to reduction routine
1069 Result              = f8
1070 r                   = f8        // For output from reduction routine
1071 c                   = f9        // For output from reduction routine
1072 U_2                 = f10
1073 rsq                 = f11
1074 C_hi                = f12
1075 C_lo                = f13
1076 T_hi                = f14
1077 T_lo                = f15
1078
1079 d_1                 = f33
1080 N_0                 = f34
1081 tail                = f35
1082 tanx                = f36
1083 Cx                  = f37
1084 Sx                  = f38
1085 sgn_r               = f39
1086 CORR                = f40
1087 P                   = f41
1088 D                   = f42
1089 ArgPrime            = f43
1090 P_0                 = f44
1091
1092 P2_1                = f45
1093 P2_2                = f46
1094 P2_3                = f47
1095
1096 P1_1                = f45
1097 P1_2                = f46
1098 P1_3                = f47
1099
1100 P1_4                = f48
1101 P1_5                = f49
1102 P1_6                = f50
1103 P1_7                = f51
1104 P1_8                = f52
1105 P1_9                = f53
1106
1107 x                   = f56
1108 xsq                 = f57
1109 Tx                  = f58
1110 Tx1                 = f59
1111 Set                 = f60
1112 poly1               = f61
1113 poly2               = f62
1114 Poly                = f63
1115 Poly1               = f64
1116 Poly2               = f65
1117 r_to_the_8          = f66
1118 B                   = f67
1119 SC_inv              = f68
1120 Pos_r               = f69
1121 N_0_fix             = f70
1122 d_2                 = f71
1123 PI_BY_4             = f72
1124 TWO_TO_NEG14        = f74
1125 TWO_TO_NEG33        = f75
1126 NEGTWO_TO_NEG14     = f76
1127 NEGTWO_TO_NEG33     = f77
1128 two_by_PI           = f78
1129 N                   = f79
1130 N_fix               = f80
1131 P_1                 = f81
1132 P_2                 = f82
1133 P_3                 = f83
1134 s_val               = f84
1135 w                   = f85
1136 B_mask1             = f86
1137 B_mask2             = f87
1138 w2                  = f88
1139 A                   = f89
1140 a                   = f90
1141 t                   = f91
1142 U_1                 = f92
1143 NEGTWO_TO_NEG2      = f93
1144 TWO_TO_NEG2         = f94
1145 Q1_1                = f95
1146 Q1_2                = f96
1147 Q1_3                = f97
1148 Q1_4                = f98
1149 Q1_5                = f99
1150 Q1_6                = f100
1151 Q1_7                = f101
1152 Q1_8                = f102
1153 S_hi                = f103
1154 S_lo                = f104
1155 V_hi                = f105
1156 V_lo                = f106
1157 U_hi                = f107
1158 U_lo                = f108
1159 U_hiabs             = f109
1160 V_hiabs             = f110
1161 V                   = f111
1162 Inv_P_0             = f112
1163
1164 FR_inv_pi_2to63     = f113
1165 FR_rshf_2to64       = f114
1166 FR_2tom64           = f115
1167 FR_rshf             = f116
1168 Norm_Arg            = f117
1169 Abs_Arg             = f118
1170 TWO_TO_NEG65        = f119
1171 fp_tmp              = f120
1172 mOne                = f121
1173
1174 GR_SAVE_B0     = r33
1175 GR_SAVE_GP     = r34
1176 GR_SAVE_PFS    = r35
1177 table_base     = r36
1178 table_ptr1     = r37
1179 table_ptr2     = r38
1180 table_ptr3     = r39
1181 lookup         = r40
1182 N_fix_gr       = r41
1183 GR_exp_2tom2   = r42
1184 GR_exp_2tom65  = r43
1185 exp_r          = r44
1186 sig_r          = r45
1187 bmask1         = r46
1188 table_offset   = r47
1189 bmask2         = r48
1190 gr_tmp         = r49
1191 cot_flag       = r50
1192
1193 GR_sig_inv_pi  = r51
1194 GR_rshf_2to64  = r52
1195 GR_exp_2tom64  = r53
1196 GR_rshf        = r54
1197 GR_exp_2_to_63 = r55
1198 GR_exp_2_to_24 = r56
1199 GR_signexp_x   = r57
1200 GR_exp_x       = r58
1201 GR_exp_mask    = r59
1202 GR_exp_2tom14  = r60
1203 GR_exp_m2tom14 = r61
1204 GR_exp_2tom33  = r62
1205 GR_exp_m2tom33 = r63
1206
1207 GR_SAVE_B0                  = r64
1208 GR_SAVE_PFS                 = r65
1209 GR_SAVE_GP                  = r66
1210
1211 GR_Parameter_X              = r67
1212 GR_Parameter_Y              = r68
1213 GR_Parameter_RESULT         = r69
1214 GR_Parameter_Tag            = r70
1215
1216
1217 .section .text
1218 .global __libm_tanl#
1219 .global __libm_cotl#
1220
1221 .proc __libm_cotl#
1222 __libm_cotl:
1223 .endp __libm_cotl#
1224 LOCAL_LIBM_ENTRY(cotl)
1225
1226 { .mlx
1227       alloc r32 = ar.pfs, 0,35,4,0
1228       movl GR_sig_inv_pi = 0xa2f9836e4e44152a // significand of 1/pi
1229 }
1230 { .mlx
1231       mov GR_exp_mask = 0x1ffff            // Exponent mask
1232       movl GR_rshf_2to64 = 0x47e8000000000000 // 1.1000 2^(63+64)
1233 }
1234 ;;
1235
1236 //     Check for NatVals, Infs , NaNs, and Zeros
1237 { .mfi
1238       getf.exp GR_signexp_x = Arg          // Get sign and exponent of x
1239       fclass.m  p6,p0 = Arg, 0x1E7         // Test for natval, nan, inf, zero
1240       mov cot_flag = 0x1
1241 }
1242 { .mfb
1243       addl table_base = @ltoff(TANL_BASE_CONSTANTS), gp // Pointer to table ptr
1244       fnorm.s1 Norm_Arg = Arg              // Normalize x
1245       br.cond.sptk COMMON_PATH
1246 };;
1247
1248 LOCAL_LIBM_END(cotl)
1249
1250
1251 .proc __libm_tanl#
1252 __libm_tanl:
1253 .endp __libm_tanl#
1254 GLOBAL_IEEE754_ENTRY(tanl)
1255
1256 { .mlx
1257       alloc r32 = ar.pfs, 0,35,4,0
1258       movl GR_sig_inv_pi = 0xa2f9836e4e44152a // significand of 1/pi
1259 }
1260 { .mlx
1261       mov GR_exp_mask = 0x1ffff            // Exponent mask
1262       movl GR_rshf_2to64 = 0x47e8000000000000 // 1.1000 2^(63+64)
1263 }
1264 ;;
1265
1266 //     Check for NatVals, Infs , NaNs, and Zeros
1267 { .mfi
1268       getf.exp GR_signexp_x = Arg          // Get sign and exponent of x
1269       fclass.m  p6,p0 = Arg, 0x1E7         // Test for natval, nan, inf, zero
1270       mov cot_flag = 0x0
1271 }
1272 { .mfi
1273       addl table_base = @ltoff(TANL_BASE_CONSTANTS), gp // Pointer to table ptr
1274       fnorm.s1 Norm_Arg = Arg              // Normalize x
1275       nop.i 0
1276 };;
1277
1278 // Common path for both tanl and cotl
1279 COMMON_PATH:
1280 { .mfi
1281       setf.sig FR_inv_pi_2to63 = GR_sig_inv_pi // Form 1/pi * 2^63
1282       fclass.m p9, p0 = Arg, 0x0b          // Test x denormal
1283       mov GR_exp_2tom64 = 0xffff - 64      // Scaling constant to compute N
1284 }
1285 { .mlx
1286       setf.d FR_rshf_2to64 = GR_rshf_2to64 // Form const 1.1000 * 2^(63+64)
1287       movl GR_rshf = 0x43e8000000000000    // Form const 1.1000 * 2^63
1288 }
1289 ;;
1290
1291 // Check for everything - if false, then must be pseudo-zero or pseudo-nan.
1292 // Branch out to deal with special values.
1293 { .mfi
1294       addl gr_tmp = -1,r0
1295       fclass.nm  p7,p0 = Arg, 0x1FF        // Test x unsupported
1296       mov GR_exp_2_to_63 = 0xffff + 63     // Exponent of 2^63
1297 }
1298 { .mfb
1299       ld8 table_base = [table_base]        // Get pointer to constant table
1300       fms.s1 mOne = f0, f0, f1
1301 (p6)  br.cond.spnt TANL_SPECIAL            // Branch if x natval, nan, inf, zero
1302 }
1303 ;;
1304
1305 { .mmb
1306       setf.sig fp_tmp = gr_tmp   // Make a constant so fmpy produces inexact
1307       mov GR_exp_2_to_24 = 0xffff + 24     // Exponent of 2^24
1308 (p9)  br.cond.spnt TANL_DENORMAL           // Branch if x denormal
1309 }
1310 ;;
1311
1312 TANL_COMMON:
1313 // Return to here if x denormal
1314 //
1315 // Do fcmp to generate Denormal exception
1316 //  - can't do FNORM (will generate Underflow when U is unmasked!)
1317 // Branch out to deal with unsupporteds values.
1318 { .mfi
1319       setf.exp FR_2tom64 = GR_exp_2tom64 // Form 2^-64 for scaling N_float
1320       fcmp.eq.s0 p0, p6 = Arg, f1        // Dummy to flag denormals
1321       add table_ptr1 = 0, table_base     // Point to tanl_table_1
1322 }
1323 { .mib
1324       setf.d FR_rshf = GR_rshf           // Form right shift const 1.1000 * 2^63
1325       add table_ptr2 = 80, table_base    // Point to tanl_table_2
1326 (p7)  br.cond.spnt TANL_UNSUPPORTED      // Branch if x unsupported type
1327 }
1328 ;;
1329
1330 { .mfi
1331       and GR_exp_x = GR_exp_mask, GR_signexp_x // Get exponent of x
1332       fmpy.s1 Save_Norm_Arg = Norm_Arg, f1     // Save x if large arg reduction
1333       dep.z bmask1 = 0x7c, 56, 8               // Form mask to get 5 msb of r
1334                                                // bmask1 = 0x7c00000000000000
1335 }
1336 ;;
1337
1338 //
1339 //     Decide about the paths to take:
1340 //     Set PR_6 if |Arg| >= 2**63
1341 //     Set PR_9 if |Arg| < 2**24 - CASE 1 OR 2
1342 //     OTHERWISE Set PR_8 - CASE 3 OR 4
1343 //
1344 //     Branch out if the magnitude of the input argument is >= 2^63
1345 //     - do this branch before the next.
1346 { .mfi
1347       ldfe two_by_PI = [table_ptr1],16        // Load 2/pi
1348       nop.f 999
1349       dep.z bmask2 = 0x41, 57, 7              // Form mask to OR to produce B
1350                                               // bmask2 = 0x8200000000000000
1351 }
1352 { .mib
1353       ldfe PI_BY_4 = [table_ptr2],16          // Load pi/4
1354       cmp.ge p6,p0 = GR_exp_x, GR_exp_2_to_63 // Is |x| >= 2^63
1355 (p6)  br.cond.spnt TANL_ARG_TOO_LARGE         // Branch if |x| >= 2^63
1356 }
1357 ;;
1358
1359 { .mmi
1360       ldfe P_0 = [table_ptr1],16              // Load P_0
1361       ldfe Inv_P_0 = [table_ptr2],16          // Load Inv_P_0
1362       nop.i 999
1363 }
1364 ;;
1365
1366 { .mfi
1367       ldfe P_1 = [table_ptr1],16              // Load P_1
1368       fmerge.s Abs_Arg = f0, Norm_Arg         // Get |x|
1369       mov GR_exp_m2tom33 = 0x2ffff - 33       // Form signexp of -2^-33
1370 }
1371 { .mfi
1372       ldfe d_1 = [table_ptr2],16              // Load d_1 for 2^24 <= |x| < 2^63
1373       nop.f 999
1374       mov GR_exp_2tom33 = 0xffff - 33         // Form signexp of 2^-33
1375 }
1376 ;;
1377
1378 { .mmi
1379       ldfe P_2 = [table_ptr1],16              // Load P_2
1380       ldfe d_2 = [table_ptr2],16              // Load d_2 for 2^24 <= |x| < 2^63
1381       cmp.ge p8,p0 = GR_exp_x, GR_exp_2_to_24 // Is |x| >= 2^24
1382 }
1383 ;;
1384
1385 // Use special scaling to right shift so N=Arg * 2/pi is in rightmost bits
1386 // Branch to Cases 3 or 4 if Arg <= -2**24 or Arg >= 2**24
1387 { .mfb
1388       ldfe   P_3 = [table_ptr1],16            // Load P_3
1389       fma.s1      N_fix = Norm_Arg, FR_inv_pi_2to63, FR_rshf_2to64
1390 (p8)  br.cond.spnt TANL_LARGER_ARG            // Branch if 2^24 <= |x| < 2^63
1391 }
1392 ;;
1393
1394 // Here if 0 < |x| < 2^24
1395 //     ARGUMENT REDUCTION CODE - CASE 1 and 2
1396 //
1397 { .mmf
1398       setf.exp TWO_TO_NEG33 = GR_exp_2tom33      // Form 2^-33
1399       setf.exp NEGTWO_TO_NEG33 = GR_exp_m2tom33  // Form -2^-33
1400       fmerge.s r = Norm_Arg,Norm_Arg          // Assume r=x, ok if |x| < pi/4
1401 }
1402 ;;
1403
1404 //
1405 // If |Arg| < pi/4,  set PR_8, else  pi/4 <=|Arg| < 2^24 - set PR_9.
1406 //
1407 //     Case 2: Convert integer N_fix back to normalized floating-point value.
1408 { .mfi
1409       getf.sig sig_r = Norm_Arg               // Get sig_r if 1/4 <= |x| < pi/4
1410       fcmp.lt.s1 p8,p9= Abs_Arg,PI_BY_4       // Test |x| < pi/4
1411       mov GR_exp_2tom2 = 0xffff - 2           // Form signexp of 2^-2
1412 }
1413 { .mfi
1414       ldfps TWO_TO_NEG2, NEGTWO_TO_NEG2 = [table_ptr2] // Load 2^-2, -2^-2
1415       fms.s1 N = N_fix, FR_2tom64, FR_rshf    // Use scaling to get N floated
1416       mov N_fix_gr = r0                       // Assume N=0, ok if |x| < pi/4
1417 }
1418 ;;
1419
1420 //
1421 //     Case 1: Is |r| < 2**(-2).
1422 //     Arg is the same as r in this case.
1423 //     r = Arg
1424 //     c = 0
1425 //
1426 //     Case 2: Place integer part of N in GP register.
1427 { .mfi
1428 (p9)  getf.sig N_fix_gr = N_fix
1429       fmerge.s c = f0, f0                     // Assume c=0, ok if |x| < pi/4
1430       cmp.lt p10, p0 = GR_exp_x, GR_exp_2tom2 // Test if |x| < 1/4
1431 }
1432 ;;
1433
1434 { .mfi
1435       setf.sig B_mask1 = bmask1               // Form mask to get 5 msb of r
1436       nop.f 999
1437       mov exp_r = GR_exp_x                    // Get exp_r if 1/4 <= |x| < pi/4
1438 }
1439 { .mbb
1440       setf.sig B_mask2 = bmask2               // Form mask to form B from r
1441 (p10) br.cond.spnt TANL_SMALL_R               // Branch if 0 < |x| < 1/4
1442 (p8)  br.cond.spnt TANL_NORMAL_R              // Branch if 1/4 <= |x| < pi/4
1443 }
1444 ;;
1445
1446 // Here if pi/4 <= |x| < 2^24
1447 //
1448 //     Case 1: PR_3 is only affected  when PR_1 is set.
1449 //
1450 //
1451 //     Case 2: w = N * P_2
1452 //     Case 2: s_val = -N * P_1  + Arg
1453 //
1454
1455 { .mfi
1456       nop.m 999
1457       fnma.s1 s_val = N, P_1, Norm_Arg
1458       nop.i 999
1459 }
1460 { .mfi
1461       nop.m 999
1462       fmpy.s1 w = N, P_2                     // w = N * P_2 for |s| >= 2^-33
1463       nop.i 999
1464 }
1465 ;;
1466
1467 //     Case 2_reduce: w = N * P_3 (change sign)
1468 { .mfi
1469       nop.m 999
1470       fmpy.s1 w2 = N, P_3                    // w = N * P_3 for |s| < 2^-33
1471       nop.i 999
1472 }
1473 ;;
1474
1475 //     Case 1_reduce: r = s + w (change sign)
1476 { .mfi
1477       nop.m 999
1478       fsub.s1 r = s_val, w                   // r = s_val - w for |s| >= 2^-33
1479       nop.i 999
1480 }
1481 ;;
1482
1483 //     Case 2_reduce: U_1 = N * P_2 + w
1484 { .mfi
1485       nop.m 999
1486       fma.s1  U_1 = N, P_2, w2              // U_1 = N * P_2 + w for |s| < 2^-33
1487       nop.i 999
1488 }
1489 ;;
1490
1491 //
1492 //     Decide between case_1 and case_2 reduce:
1493 //     Case 1_reduce:  |s| >= 2**(-33)
1494 //     Case 2_reduce:  |s| < 2**(-33)
1495 //
1496 { .mfi
1497       nop.m 999
1498       fcmp.lt.s1 p9, p8 = s_val, TWO_TO_NEG33
1499       nop.i 999
1500 }
1501 ;;
1502
1503 { .mfi
1504       nop.m 999
1505 (p9)  fcmp.gt.s1 p9, p8 = s_val, NEGTWO_TO_NEG33
1506       nop.i 999
1507 }
1508 ;;
1509
1510 //     Case 1_reduce: c = s - r
1511 { .mfi
1512       nop.m 999
1513       fsub.s1 c = s_val, r                     // c = s_val - r for |s| >= 2^-33
1514       nop.i 999
1515 }
1516 ;;
1517
1518 //     Case 2_reduce: r is complete here - continue to calculate c .
1519 //     r = s - U_1
1520 { .mfi
1521       nop.m 999
1522 (p9)  fsub.s1 r = s_val, U_1
1523       nop.i 999
1524 }
1525 { .mfi
1526       nop.m 999
1527 (p9)  fms.s1 U_2 = N, P_2, U_1
1528       nop.i 999
1529 }
1530 ;;
1531
1532 //
1533 //     Case 1_reduce: Is |r| < 2**(-2), if so set PR_10
1534 //     else set PR_13.
1535 //
1536
1537 { .mfi
1538       nop.m 999
1539       fand B = B_mask1, r
1540       nop.i 999
1541 }
1542 { .mfi
1543       nop.m 999
1544 (p8)  fcmp.lt.unc.s1 p10, p13 = r, TWO_TO_NEG2
1545       nop.i 999
1546 }
1547 ;;
1548
1549 { .mfi
1550 (p8)  getf.sig sig_r = r               // Get signif of r if |s| >= 2^-33
1551       nop.f 999
1552       nop.i 999
1553 }
1554 ;;
1555
1556 { .mfi
1557 (p8)  getf.exp exp_r = r               // Extract signexp of r if |s| >= 2^-33
1558 (p10) fcmp.gt.s1 p10, p13 = r, NEGTWO_TO_NEG2
1559       nop.i 999
1560 }
1561 ;;
1562
1563 //     Case 1_reduce: c is complete here.
1564 //     Case 1: Branch to SMALL_R or NORMAL_R.
1565 //     c = c + w (w has not been negated.)
1566 { .mfi
1567       nop.m 999
1568 (p8)  fsub.s1 c = c, w                         // c = c - w for |s| >= 2^-33
1569       nop.i 999
1570 }
1571 { .mbb
1572       nop.m 999
1573 (p10) br.cond.spnt TANL_SMALL_R     // Branch if pi/4 < |x| < 2^24 and |r|<1/4
1574 (p13) br.cond.sptk TANL_NORMAL_R_A  // Branch if pi/4 < |x| < 2^24 and |r|>=1/4
1575 }
1576 ;;
1577
1578
1579 // Here if pi/4 < |x| < 2^24 and |s| < 2^-33
1580 //
1581 //     Is i_1 = lsb of N_fix_gr even or odd?
1582 //     if i_1 == 0, set p11, else set p12.
1583 //
1584 { .mfi
1585       nop.m 999
1586       fsub.s1 s_val = s_val, r
1587       add N_fix_gr = N_fix_gr, cot_flag // N = N + 1 (for cotl)
1588 }
1589 { .mfi
1590       nop.m 999
1591 //
1592 //     Case 2_reduce:
1593 //     U_2 = N * P_2 - U_1
1594 //     Not needed until later.
1595 //
1596       fadd.s1 U_2 = U_2, w2
1597 //
1598 //     Case 2_reduce:
1599 //     s = s - r
1600 //     U_2 = U_2 + w
1601 //
1602       nop.i 999
1603 }
1604 ;;
1605
1606 //
1607 //     Case 2_reduce:
1608 //     c = c - U_2
1609 //     c is complete here
1610 //     Argument reduction ends here.
1611 //
1612 { .mfi
1613       nop.m 999
1614       fmpy.s1 rsq = r, r
1615       tbit.z p11, p12 = N_fix_gr, 0 ;;    // Set p11 if N even, p12 if odd
1616 }
1617
1618 { .mfi
1619       nop.m 999
1620 (p12) frcpa.s1 S_hi,p0 = f1, r
1621       nop.i 999
1622 }
1623 { .mfi
1624       nop.m 999
1625       fsub.s1 c = s_val, U_1
1626       nop.i 999
1627 }
1628 ;;
1629
1630 { .mmi
1631       add table_ptr1 = 160, table_base ;;  // Point to tanl_table_p1
1632       ldfe P1_1 = [table_ptr1],144
1633       nop.i 999 ;;
1634 }
1635 //
1636 //     Load P1_1 and point to Q1_1 .
1637 //
1638 { .mfi
1639       ldfe Q1_1 = [table_ptr1]
1640 //
1641 //     N even: rsq = r * Z
1642 //     N odd:  S_hi = frcpa(r)
1643 //
1644 (p12) fmerge.ns S_hi = S_hi, S_hi
1645       nop.i 999
1646 }
1647 { .mfi
1648       nop.m 999
1649 //
1650 //     Case 2_reduce:
1651 //     c = s - U_1
1652 //
1653 (p9)  fsub.s1 c = c, U_2
1654       nop.i 999 ;;
1655 }
1656 { .mfi
1657       nop.m 999
1658 (p12) fma.s1  poly1 = S_hi, r, f1
1659       nop.i 999 ;;
1660 }
1661 { .mfi
1662       nop.m 999
1663 //
1664 //     N odd:  Change sign of S_hi
1665 //
1666 (p11) fmpy.s1 rsq = rsq, P1_1
1667       nop.i 999 ;;
1668 }
1669 { .mfi
1670       nop.m 999
1671 (p12) fma.s1 S_hi = S_hi, poly1, S_hi
1672       nop.i 999 ;;
1673 }
1674 { .mfi
1675       nop.m 999
1676 //
1677 //     N even: rsq = rsq * P1_1
1678 //     N odd:  poly1 =  1.0 +  S_hi * r    16 bits partial  account for necessary
1679 //
1680 (p11) fma.s1 Poly = r, rsq, c
1681       nop.i 999 ;;
1682 }
1683 { .mfi
1684       nop.m 999
1685 //
1686 //     N even: Poly = c  + r * rsq
1687 //     N odd:  S_hi  = S_hi + S_hi*poly1  16 bits account for necessary
1688 //
1689 (p12) fma.s1 poly1 = S_hi, r, f1
1690 (p11) tbit.z.unc p14, p15 = cot_flag, 0 ;; // p14=1 for tanl; p15=1 for cotl
1691 }
1692 { .mfi
1693       nop.m 999
1694 //
1695 //     N even: Result = Poly + r
1696 //     N odd:  poly1  = 1.0 + S_hi * r        32 bits partial
1697 //
1698 (p14) fadd.s0 Result = r, Poly             // for tanl
1699       nop.i 999
1700 }
1701 { .mfi
1702       nop.m 999
1703 (p15) fms.s0 Result = r, mOne, Poly        // for cotl
1704       nop.i 999
1705 }
1706 ;;
1707
1708 { .mfi
1709       nop.m 999
1710 (p12) fma.s1  S_hi = S_hi, poly1, S_hi
1711       nop.i 999 ;;
1712 }
1713 { .mfi
1714       nop.m 999
1715 //
1716 //     N even: Result1 = Result + r
1717 //     N odd:   S_hi  = S_hi * poly1 + S_hi   32 bits
1718 //
1719 (p12) fma.s1 poly1 = S_hi, r, f1
1720       nop.i 999 ;;
1721 }
1722 { .mfi
1723       nop.m 999
1724 //
1725 //     N odd:  poly1  =  S_hi * r + 1.0       64 bits partial
1726 //
1727 (p12) fma.s1 S_hi = S_hi, poly1, S_hi
1728       nop.i 999 ;;
1729 }
1730 { .mfi
1731       nop.m 999
1732 //
1733 //     N odd:  poly1  =  S_hi * poly + 1.0    64 bits
1734 //
1735 (p12) fma.s1 poly1 = S_hi, r, f1
1736       nop.i 999 ;;
1737 }
1738 { .mfi
1739       nop.m 999
1740 //
1741 //     N odd:  poly1  =  S_hi * r + 1.0
1742 //
1743 (p12) fma.s1 poly1 = S_hi, c, poly1
1744       nop.i 999 ;;
1745 }
1746 { .mfi
1747       nop.m 999
1748 //
1749 //     N odd:  poly1  =  S_hi * c + poly1
1750 //
1751 (p12) fmpy.s1 S_lo = S_hi, poly1
1752       nop.i 999 ;;
1753 }
1754 { .mfi
1755       nop.m 999
1756 //
1757 //     N odd:  S_lo  =  S_hi *  poly1
1758 //
1759 (p12) fma.s1 S_lo = Q1_1, r, S_lo
1760 (p12) tbit.z.unc p14, p15 = cot_flag, 0 // p14=1 for tanl; p15=1 for cotl
1761 }
1762 { .mfi
1763       nop.m 999
1764 //
1765 //     N odd:  Result =  S_hi + S_lo
1766 //
1767       fmpy.s0 fp_tmp = fp_tmp, fp_tmp  // Dummy mult to set inexact
1768       nop.i 999 ;;
1769 }
1770 { .mfi
1771       nop.m 999
1772 //
1773 //     N odd:  S_lo  =  S_lo + Q1_1 * r
1774 //
1775 (p14) fadd.s0 Result = S_hi, S_lo          // for tanl
1776       nop.i 999
1777 }
1778 { .mfb
1779       nop.m 999
1780 (p15) fms.s0 Result = S_hi, mOne, S_lo     // for cotl
1781       br.ret.sptk b0 ;;          // Exit for pi/4 <= |x| < 2^24 and |s| < 2^-33
1782 }
1783
1784
1785 TANL_LARGER_ARG:
1786 // Here if 2^24 <= |x| < 2^63
1787 //
1788 // ARGUMENT REDUCTION CODE - CASE 3 and 4
1789 //
1790
1791 { .mmf
1792       mov GR_exp_2tom14 = 0xffff - 14          // Form signexp of 2^-14
1793       mov GR_exp_m2tom14 = 0x2ffff - 14        // Form signexp of -2^-14
1794       fmpy.s1 N_0 = Norm_Arg, Inv_P_0
1795 }
1796 ;;
1797
1798 { .mmi
1799       setf.exp TWO_TO_NEG14 = GR_exp_2tom14    // Form 2^-14
1800       setf.exp NEGTWO_TO_NEG14 = GR_exp_m2tom14// Form -2^-14
1801       nop.i 999
1802 }
1803 ;;
1804
1805
1806 //
1807 //    Adjust table_ptr1 to beginning of table.
1808 //    N_0 = Arg * Inv_P_0
1809 //
1810 { .mmi
1811       add table_ptr2 = 144, table_base ;;     // Point to 2^-2
1812       ldfps TWO_TO_NEG2, NEGTWO_TO_NEG2 = [table_ptr2]
1813       nop.i 999
1814 }
1815 ;;
1816
1817 //
1818 //    N_0_fix  = integer part of N_0 .
1819 //
1820 //
1821 //    Make N_0 the integer part.
1822 //
1823 { .mfi
1824       nop.m 999
1825       fcvt.fx.s1 N_0_fix = N_0
1826       nop.i 999 ;;
1827 }
1828 { .mfi
1829       setf.sig B_mask1 = bmask1               // Form mask to get 5 msb of r
1830       fcvt.xf N_0 = N_0_fix
1831       nop.i 999 ;;
1832 }
1833 { .mfi
1834       setf.sig B_mask2 = bmask2               // Form mask to form B from r
1835       fnma.s1 ArgPrime = N_0, P_0, Norm_Arg
1836       nop.i 999
1837 }
1838 { .mfi
1839       nop.m 999
1840       fmpy.s1 w = N_0, d_1
1841       nop.i 999 ;;
1842 }
1843 //
1844 //    ArgPrime = -N_0 * P_0 + Arg
1845 //    w  = N_0 * d_1
1846 //
1847 //
1848 //    N = ArgPrime * 2/pi
1849 //
1850 //      fcvt.fx.s1 N_fix = N
1851 // Use special scaling to right shift so N=Arg * 2/pi is in rightmost bits
1852 // Branch to Cases 3 or 4 if Arg <= -2**24 or Arg >= 2**24
1853 { .mfi
1854       nop.m 999
1855       fma.s1      N_fix = ArgPrime, FR_inv_pi_2to63, FR_rshf_2to64
1856
1857       nop.i 999 ;;
1858 }
1859 //     Convert integer N_fix back to normalized floating-point value.
1860 { .mfi
1861       nop.m 999
1862       fms.s1 N = N_fix, FR_2tom64, FR_rshf    // Use scaling to get N floated
1863       nop.i 999
1864 }
1865 ;;
1866
1867 //
1868 //    N is the integer part of the reduced-reduced argument.
1869 //    Put the integer in a GP register.
1870 //
1871 { .mfi
1872       getf.sig N_fix_gr = N_fix
1873       nop.f 999
1874       nop.i 999
1875 }
1876 ;;
1877
1878 //
1879 //    s_val = -N*P_1 + ArgPrime
1880 //    w = -N*P_2 + w
1881 //
1882 { .mfi
1883       nop.m 999
1884       fnma.s1 s_val = N, P_1, ArgPrime
1885       nop.i 999
1886 }
1887 { .mfi
1888       nop.m 999
1889       fnma.s1 w = N, P_2, w
1890       nop.i 999
1891 }
1892 ;;
1893
1894 //    Case 4: V_hi = N * P_2
1895 //    Case 4: U_hi = N_0 * d_1
1896 { .mfi
1897       nop.m 999
1898       fmpy.s1 V_hi = N, P_2               // V_hi = N * P_2 for |s| < 2^-14
1899       nop.i 999
1900 }
1901 { .mfi
1902       nop.m 999
1903       fmpy.s1 U_hi = N_0, d_1             // U_hi = N_0 * d_1 for |s| < 2^-14
1904       nop.i 999
1905 }
1906 ;;
1907
1908 //    Case 3: r = s_val + w (Z complete)
1909 //    Case 4: w = N * P_3
1910 { .mfi
1911       nop.m 999
1912       fadd.s1 r = s_val, w                // r = s_val + w for |s| >= 2^-14
1913       nop.i 999
1914 }
1915 { .mfi
1916       nop.m 999
1917       fmpy.s1 w2 = N, P_3                 // w = N * P_3 for |s| < 2^-14
1918       nop.i 999
1919 }
1920 ;;
1921
1922 //    Case 4: A =  U_hi + V_hi
1923 //    Note: Worry about switched sign of V_hi, so subtract instead of add.
1924 //    Case 4: V_lo = -N * P_2 - V_hi (U_hi is in place of V_hi in writeup)
1925 //    Note: the (-) is still missing for V_hi.
1926 { .mfi
1927       nop.m 999
1928       fsub.s1 A = U_hi, V_hi           // A = U_hi - V_hi for |s| < 2^-14
1929       nop.i 999
1930 }
1931 { .mfi
1932       nop.m 999
1933       fnma.s1 V_lo = N, P_2, V_hi      // V_lo = V_hi - N * P_2 for |s| < 2^-14
1934       nop.i 999
1935 }
1936 ;;
1937
1938 //    Decide between case 3 and 4:
1939 //    Case 3:  |s| >= 2**(-14)     Set p10
1940 //    Case 4:  |s| <  2**(-14)     Set p11
1941 //
1942 //    Case 4: U_lo = N_0 * d_1 - U_hi
1943 { .mfi
1944       nop.m 999
1945       fms.s1 U_lo = N_0, d_1, U_hi     // U_lo = N_0*d_1 - U_hi for |s| < 2^-14
1946       nop.i 999
1947 }
1948 { .mfi
1949       nop.m 999
1950       fcmp.lt.s1 p11, p10 = s_val, TWO_TO_NEG14
1951       nop.i 999
1952 }
1953 ;;
1954
1955 //    Case 4: We need abs of both U_hi and V_hi - dont
1956 //    worry about switched sign of V_hi.
1957 { .mfi
1958       nop.m 999
1959       fabs V_hiabs = V_hi              // |V_hi| for |s| < 2^-14
1960       nop.i 999
1961 }
1962 { .mfi
1963       nop.m 999
1964 (p11) fcmp.gt.s1 p11, p10 = s_val, NEGTWO_TO_NEG14
1965       nop.i 999
1966 }
1967 ;;
1968
1969 //    Case 3: c = s_val - r
1970 { .mfi
1971       nop.m 999
1972       fabs U_hiabs = U_hi              // |U_hi| for |s| < 2^-14
1973       nop.i 999
1974 }
1975 { .mfi
1976       nop.m 999
1977       fsub.s1 c = s_val, r             // c = s_val - r    for |s| >= 2^-14
1978       nop.i 999
1979 }
1980 ;;
1981
1982 // For Case 3, |s| >= 2^-14, determine if |r| < 1/4
1983 //
1984 //    Case 4: C_hi = s_val + A
1985 //
1986 { .mfi
1987       nop.m 999
1988 (p11) fadd.s1 C_hi = s_val, A              // C_hi = s_val + A for |s| < 2^-14
1989       nop.i 999
1990 }
1991 { .mfi
1992       nop.m 999
1993 (p10) fcmp.lt.unc.s1 p14, p15 = r, TWO_TO_NEG2
1994       nop.i 999
1995 }
1996 ;;
1997
1998 { .mfi
1999       getf.sig sig_r = r               // Get signif of r if |s| >= 2^-33
2000       fand B = B_mask1, r
2001       nop.i 999
2002 }
2003 ;;
2004
2005 //    Case 4: t = U_lo + V_lo
2006 { .mfi
2007       getf.exp exp_r = r               // Extract signexp of r if |s| >= 2^-33
2008 (p11) fadd.s1 t = U_lo, V_lo               // t = U_lo + V_lo for |s| < 2^-14
2009       nop.i 999
2010 }
2011 { .mfi
2012       nop.m 999
2013 (p14) fcmp.gt.s1 p14, p15 = r, NEGTWO_TO_NEG2
2014       nop.i 999
2015 }
2016 ;;
2017
2018 //    Case 3: c = (s - r) + w (c complete)
2019 { .mfi
2020       nop.m 999
2021 (p10) fadd.s1 c = c, w              // c = c + w for |s| >= 2^-14
2022       nop.i 999
2023 }
2024 { .mbb
2025       nop.m 999
2026 (p14) br.cond.spnt TANL_SMALL_R     // Branch if 2^24 <= |x| < 2^63 and |r|< 1/4
2027 (p15) br.cond.sptk TANL_NORMAL_R_A  // Branch if 2^24 <= |x| < 2^63 and |r|>=1/4
2028 }
2029 ;;
2030
2031
2032 // Here if 2^24 <= |x| < 2^63 and |s| < 2^-14  >>>>>>>  Case 4.
2033 //
2034 //    Case 4: Set P_12 if U_hiabs >= V_hiabs
2035 //    Case 4: w = w + N_0 * d_2
2036 //    Note: the (-) is now incorporated in w .
2037 { .mfi
2038       add table_ptr1 = 160, table_base           // Point to tanl_table_p1
2039       fcmp.ge.unc.s1 p12, p13 = U_hiabs, V_hiabs
2040       nop.i 999
2041 }
2042 { .mfi
2043       nop.m 999
2044       fms.s1 w2 = N_0, d_2, w2
2045       nop.i 999
2046 }
2047 ;;
2048
2049 //    Case 4: C_lo = s_val - C_hi
2050 { .mfi
2051       ldfe P1_1 = [table_ptr1], 16               // Load P1_1
2052       fsub.s1 C_lo = s_val, C_hi
2053       nop.i 999
2054 }
2055 ;;
2056
2057 //
2058 //    Case 4: a = U_hi - A
2059 //            a = V_hi - A (do an add to account for missing (-) on V_hi
2060 //
2061 { .mfi
2062       ldfe P1_2 = [table_ptr1], 128              // Load P1_2
2063 (p12) fsub.s1 a = U_hi, A
2064       nop.i 999
2065 }
2066 { .mfi
2067       nop.m 999
2068 (p13) fadd.s1 a = V_hi, A
2069       nop.i 999
2070 }
2071 ;;
2072
2073 //    Case 4: t = U_lo + V_lo  + w
2074 { .mfi
2075       ldfe Q1_1 = [table_ptr1], 16               // Load Q1_1
2076       fadd.s1 t = t, w2
2077       nop.i 999
2078 }
2079 ;;
2080
2081 //    Case 4: a = (U_hi - A)  + V_hi
2082 //            a = (V_hi - A)  + U_hi
2083 //    In each case account for negative missing form V_hi .
2084 //
2085 { .mfi
2086       ldfe Q1_2 = [table_ptr1], 16               // Load Q1_2
2087 (p12) fsub.s1 a = a, V_hi
2088       nop.i 999
2089 }
2090 { .mfi
2091       nop.m 999
2092 (p13) fsub.s1 a = U_hi, a
2093       nop.i 999
2094 }
2095 ;;
2096
2097 //
2098 //    Case 4: C_lo = (s_val - C_hi) + A
2099 //
2100 { .mfi
2101       nop.m 999
2102       fadd.s1 C_lo = C_lo, A
2103       nop.i 999 ;;
2104 }
2105 //
2106 //    Case 4: t = t + a
2107 //
2108 { .mfi
2109       nop.m 999
2110       fadd.s1 t = t, a
2111       nop.i 999
2112 }
2113 ;;
2114
2115 //    Case 4: C_lo = C_lo + t
2116 //    Case 4: r = C_hi + C_lo
2117 { .mfi
2118       nop.m 999
2119       fadd.s1 C_lo = C_lo, t
2120       nop.i 999
2121 }
2122 ;;
2123
2124 { .mfi
2125       nop.m 999
2126       fadd.s1 r = C_hi, C_lo
2127       nop.i 999
2128 }
2129 ;;
2130
2131 //
2132 //    Case 4: c = C_hi - r
2133 //
2134 { .mfi
2135       nop.m 999
2136       fsub.s1 c = C_hi, r
2137       nop.i 999
2138 }
2139 { .mfi
2140       nop.m 999
2141       fmpy.s1 rsq = r, r
2142       add N_fix_gr = N_fix_gr, cot_flag // N = N + 1 (for cotl)
2143 }
2144 ;;
2145
2146 //    Case 4: c = c + C_lo  finished.
2147 //
2148 //    Is i_1 = lsb of N_fix_gr even or odd?
2149 //    if i_1 == 0, set PR_11, else set PR_12.
2150 //
2151 { .mfi
2152       nop.m 999
2153       fadd.s1 c = c , C_lo
2154       tbit.z p11, p12 =  N_fix_gr, 0
2155 }
2156 ;;
2157
2158 // r and c have been computed.
2159 { .mfi
2160       nop.m 999
2161 (p12) frcpa.s1 S_hi, p0 = f1, r
2162       nop.i 999
2163 }
2164 { .mfi
2165       nop.m 999
2166 //
2167 //    N odd: Change sign of S_hi
2168 //
2169 (p11) fma.s1 Poly = rsq, P1_2, P1_1
2170       nop.i 999 ;;
2171 }
2172 { .mfi
2173       nop.m 999
2174 (p12) fma.s1 P = rsq, Q1_2, Q1_1
2175       nop.i 999
2176 }
2177 { .mfi
2178       nop.m 999
2179 //
2180 //    N odd:  Result  =  S_hi + S_lo      (User supplied rounding mode for C1)
2181 //
2182        fmpy.s0 fp_tmp = fp_tmp, fp_tmp  // Dummy mult to set inexact
2183       nop.i 999 ;;
2184 }
2185 { .mfi
2186       nop.m 999
2187 //
2188 //    N even: rsq = r * r
2189 //    N odd:  S_hi = frcpa(r)
2190 //
2191 (p12) fmerge.ns S_hi = S_hi, S_hi
2192       nop.i 999
2193 }
2194 { .mfi
2195       nop.m 999
2196 //
2197 //    N even: rsq = rsq * P1_2 + P1_1
2198 //    N odd:  poly1 =  1.0 +  S_hi * r    16 bits partial  account for necessary
2199 //
2200 (p11) fmpy.s1 Poly = rsq, Poly
2201       nop.i 999 ;;
2202 }
2203 { .mfi
2204       nop.m 999
2205 (p12) fma.s1 poly1 = S_hi, r,f1
2206 (p11) tbit.z.unc p14, p15 = cot_flag, 0 // p14=1 for tanl; p15=1 for cotl
2207 }
2208 { .mfi
2209       nop.m 999
2210 //
2211 //    N even: Poly =  Poly * rsq
2212 //    N odd:  S_hi  = S_hi + S_hi*poly1  16 bits account for necessary
2213 //
2214 (p11) fma.s1 Poly = r, Poly, c
2215       nop.i 999 ;;
2216 }
2217 { .mfi
2218       nop.m 999
2219 (p12) fma.s1 S_hi = S_hi, poly1, S_hi
2220       nop.i 999
2221 }
2222 { .mfi
2223       nop.m 999
2224 //
2225 //    N odd:   S_hi  = S_hi * poly1 + S_hi   32 bits
2226 //
2227 (p14) fadd.s0 Result = r, Poly          // for tanl
2228       nop.i 999 ;;
2229 }
2230
2231 .pred.rel "mutex",p15,p12
2232 { .mfi
2233       nop.m 999
2234 (p15) fms.s0 Result = r, mOne, Poly     // for cotl
2235       nop.i 999
2236 }
2237 { .mfi
2238       nop.m 999
2239 (p12) fma.s1 poly1 =  S_hi, r, f1
2240       nop.i 999 ;;
2241 }
2242 { .mfi
2243       nop.m 999
2244 //
2245 //    N even: Poly = Poly * r + c
2246 //    N odd:  poly1  = 1.0 + S_hi * r        32 bits partial
2247 //
2248 (p12) fma.s1 S_hi = S_hi, poly1, S_hi
2249       nop.i 999 ;;
2250 }
2251 { .mfi
2252       nop.m 999
2253 (p12) fma.s1 poly1 = S_hi, r, f1
2254       nop.i 999 ;;
2255 }
2256 { .mfi
2257       nop.m 999
2258 //
2259 //    N even: Result = Poly + r  (Rounding mode S0)
2260 //    N odd:  poly1  =  S_hi * r + 1.0       64 bits partial
2261 //
2262 (p12) fma.s1 S_hi = S_hi, poly1, S_hi
2263       nop.i 999 ;;
2264 }
2265 { .mfi
2266       nop.m 999
2267 //
2268 //    N odd:  poly1  =  S_hi * poly + S_hi    64 bits
2269 //
2270 (p12) fma.s1 poly1 = S_hi, r, f1
2271       nop.i 999 ;;
2272 }
2273 { .mfi
2274       nop.m 999
2275 //
2276 //    N odd:  poly1  =  S_hi * r + 1.0
2277 //
2278 (p12) fma.s1 poly1 = S_hi, c, poly1
2279       nop.i 999 ;;
2280 }
2281 { .mfi
2282       nop.m 999
2283 //
2284 //    N odd:  poly1  =  S_hi * c + poly1
2285 //
2286 (p12) fmpy.s1 S_lo = S_hi, poly1
2287       nop.i 999 ;;
2288 }
2289 { .mfi
2290       nop.m 999
2291 //
2292 //    N odd:  S_lo  =  S_hi *  poly1
2293 //
2294 (p12) fma.s1 S_lo = P, r, S_lo
2295 (p12) tbit.z.unc p14, p15 = cot_flag, 0 ;; // p14=1 for tanl; p15=1 for cotl
2296 }
2297
2298 { .mfi
2299       nop.m 999
2300 (p14) fadd.s0 Result = S_hi, S_lo           // for tanl
2301       nop.i 999
2302 }
2303 { .mfb
2304       nop.m 999
2305 //
2306 //    N odd:  S_lo  =  S_lo + r * P
2307 //
2308 (p15) fms.s0 Result = S_hi, mOne, S_lo      // for cotl
2309       br.ret.sptk b0 ;;      // Exit for 2^24 <= |x| < 2^63 and |s| < 2^-14
2310 }
2311
2312
2313 TANL_SMALL_R:
2314 // Here if |r| < 1/4
2315 // r and c have been computed.
2316 // *****************************************************************
2317 // *****************************************************************
2318 // *****************************************************************
2319 //    N odd:  S_hi = frcpa(r)
2320 //    Get [i_1] - lsb of N_fix_gr.  Set p11 if N even, p12 if N odd.
2321 //    N even: rsq = r * r
2322 { .mfi
2323       add table_ptr1 = 160, table_base    // Point to tanl_table_p1
2324       frcpa.s1 S_hi, p0 = f1, r           // S_hi for N odd
2325       add N_fix_gr = N_fix_gr, cot_flag   // N = N + 1 (for cotl)
2326 }
2327 { .mfi
2328       add table_ptr2 = 400, table_base    // Point to Q1_7
2329       fmpy.s1 rsq = r, r
2330       nop.i 999
2331 }
2332 ;;
2333
2334 { .mmi
2335       ldfe P1_1 = [table_ptr1], 16
2336 ;;
2337       ldfe P1_2 = [table_ptr1], 16
2338       tbit.z p11, p12 = N_fix_gr, 0
2339 }
2340 ;;
2341
2342
2343 { .mfi
2344       ldfe P1_3 = [table_ptr1], 96
2345       nop.f 999
2346       nop.i 999
2347 }
2348 ;;
2349
2350 { .mfi
2351 (p11) ldfe P1_9 = [table_ptr1], -16
2352 (p12) fmerge.ns S_hi = S_hi, S_hi
2353       nop.i 999
2354 }
2355 { .mfi
2356       nop.m 999
2357 (p11) fmpy.s1 r_to_the_8 = rsq, rsq
2358       nop.i 999
2359 }
2360 ;;
2361
2362 //
2363 //    N even: Poly2 = P1_7 + Poly2 * rsq
2364 //    N odd:  poly2 = Q1_5 + poly2 * rsq
2365 //
2366 { .mfi
2367 (p11) ldfe P1_8 = [table_ptr1], -16
2368 (p11) fadd.s1 CORR = rsq, f1
2369       nop.i 999
2370 }
2371 ;;
2372
2373 //
2374 //    N even: Poly1 = P1_2 + P1_3 * rsq
2375 //    N odd:  poly1 =  1.0 +  S_hi * r
2376 //    16 bits partial  account for necessary (-1)
2377 //
2378 { .mmi
2379 (p11) ldfe P1_7 = [table_ptr1], -16
2380 ;;
2381 (p11) ldfe P1_6 = [table_ptr1], -16
2382       nop.i 999
2383 }
2384 ;;
2385
2386 //
2387 //    N even: Poly1 = P1_1 + Poly1 * rsq
2388 //    N odd:  S_hi  =  S_hi + S_hi * poly1)     16 bits account for necessary
2389 //
2390 //
2391 //    N even: Poly2 = P1_5 + Poly2 * rsq
2392 //    N odd:  poly2 = Q1_3 + poly2 * rsq
2393 //
2394 { .mfi
2395 (p11) ldfe P1_5 = [table_ptr1], -16
2396 (p11) fmpy.s1 r_to_the_8 = r_to_the_8, r_to_the_8
2397       nop.i 999
2398 }
2399 { .mfi
2400       nop.m 999
2401 (p12) fma.s1 poly1 =  S_hi, r, f1
2402       nop.i 999
2403 }
2404 ;;
2405
2406 //
2407 //    N even: Poly1 =  Poly1 * rsq
2408 //    N odd:  poly1  = 1.0 + S_hi * r         32 bits partial
2409 //
2410
2411 //
2412 //    N even: CORR =  CORR * c
2413 //    N odd:  S_hi  =  S_hi * poly1 + S_hi    32 bits
2414 //
2415
2416 //
2417 //    N even: Poly2 = P1_6 + Poly2 * rsq
2418 //    N odd:  poly2 = Q1_4 + poly2 * rsq
2419 //
2420
2421 { .mmf
2422 (p11) ldfe P1_4 = [table_ptr1], -16
2423       nop.m 999
2424 (p11) fmpy.s1 CORR =  CORR, c
2425 }
2426 ;;
2427
2428 { .mfi
2429       nop.m 999
2430 (p11) fma.s1 Poly1 = P1_3, rsq, P1_2
2431       nop.i 999 ;;
2432 }
2433 { .mfi
2434 (p12) ldfe Q1_7 = [table_ptr2], -16
2435 (p12) fma.s1 S_hi = S_hi, poly1, S_hi
2436       nop.i 999 ;;
2437 }
2438 { .mfi
2439 (p12) ldfe Q1_6 = [table_ptr2], -16
2440 (p11) fma.s1 Poly2 = P1_9, rsq, P1_8
2441       nop.i 999 ;;
2442 }
2443 { .mmi
2444 (p12) ldfe Q1_5 = [table_ptr2], -16 ;;
2445 (p12) ldfe Q1_4 = [table_ptr2], -16
2446       nop.i 999 ;;
2447 }
2448 { .mfi
2449 (p12) ldfe Q1_3 = [table_ptr2], -16
2450 //
2451 //    N even: Poly2 = P1_8 + P1_9 * rsq
2452 //    N odd:  poly2 = Q1_6 + Q1_7 * rsq
2453 //
2454 (p11) fma.s1 Poly1 = Poly1, rsq, P1_1
2455       nop.i 999 ;;
2456 }
2457 { .mfi
2458 (p12) ldfe Q1_2 = [table_ptr2], -16
2459 (p12) fma.s1 poly1 = S_hi, r, f1
2460       nop.i 999 ;;
2461 }
2462 { .mfi
2463 (p12) ldfe Q1_1 = [table_ptr2], -16
2464 (p11) fma.s1 Poly2 = Poly2, rsq, P1_7
2465       nop.i 999 ;;
2466 }
2467 { .mfi
2468       nop.m 999
2469 //
2470 //    N even: CORR =  rsq + 1
2471 //    N even: r_to_the_8 =  rsq * rsq
2472 //
2473 (p11) fmpy.s1 Poly1 = Poly1, rsq
2474       nop.i 999 ;;
2475 }
2476 { .mfi
2477       nop.m 999
2478 (p12) fma.s1 S_hi = S_hi, poly1, S_hi
2479       nop.i 999
2480 }
2481 { .mfi
2482       nop.m 999
2483 (p12) fma.s1 poly2 = Q1_7, rsq, Q1_6
2484       nop.i 999 ;;
2485 }
2486 { .mfi
2487       nop.m 999
2488 (p11) fma.s1 Poly2 = Poly2, rsq, P1_6
2489       nop.i 999 ;;
2490 }
2491 { .mfi
2492       nop.m 999
2493 (p12) fma.s1 poly1 = S_hi, r, f1
2494       nop.i 999
2495 }
2496 { .mfi
2497       nop.m 999
2498 (p12) fma.s1 poly2 = poly2, rsq, Q1_5
2499       nop.i 999 ;;
2500 }
2501 { .mfi
2502       nop.m 999
2503 (p11) fma.s1 Poly2= Poly2, rsq, P1_5
2504       nop.i 999 ;;
2505 }
2506 { .mfi
2507       nop.m 999
2508 (p12) fma.s1 S_hi =  S_hi, poly1, S_hi
2509       nop.i 999
2510 }
2511 { .mfi
2512       nop.m 999
2513 (p12) fma.s1 poly2 = poly2, rsq, Q1_4
2514       nop.i 999 ;;
2515 }
2516 { .mfi
2517       nop.m 999
2518 //
2519 //    N even: r_to_the_8 = r_to_the_8 * r_to_the_8
2520 //    N odd:  poly1  =  S_hi * r + 1.0       64 bits partial
2521 //
2522 (p11) fma.s1 Poly2 = Poly2, rsq, P1_4
2523       nop.i 999 ;;
2524 }
2525 { .mfi
2526       nop.m 999
2527 //
2528 //    N even: Poly = CORR + Poly * r
2529 //    N odd:  P = Q1_1 + poly2 * rsq
2530 //
2531 (p12) fma.s1 poly1 = S_hi, r, f1
2532       nop.i 999
2533 }
2534 { .mfi
2535       nop.m 999
2536 (p12) fma.s1 poly2 = poly2, rsq, Q1_3
2537       nop.i 999 ;;
2538 }
2539 { .mfi
2540       nop.m 999
2541 //
2542 //    N even: Poly2 = P1_4 + Poly2 * rsq
2543 //    N odd:  poly2 = Q1_2 + poly2 * rsq
2544 //
2545 (p11) fma.s1 Poly = Poly2, r_to_the_8, Poly1
2546       nop.i 999 ;;
2547 }
2548 { .mfi
2549       nop.m 999
2550 (p12) fma.s1 poly1 = S_hi, c, poly1
2551       nop.i 999
2552 }
2553 { .mfi
2554       nop.m 999
2555 (p12) fma.s1 poly2 = poly2, rsq, Q1_2
2556       nop.i 999 ;;
2557 }
2558
2559 { .mfi
2560       nop.m 999
2561 //
2562 //    N even: Poly = Poly1 + Poly2 * r_to_the_8
2563 //    N odd:  S_hi =  S_hi * poly1 + S_hi    64 bits
2564 //
2565 (p11) fma.s1 Poly = Poly, r, CORR
2566       nop.i 999 ;;
2567 }
2568 { .mfi
2569       nop.m 999
2570 //
2571 //    N even: Result =  r + Poly  (User supplied rounding mode)
2572 //    N odd:  poly1  =  S_hi * c + poly1
2573 //
2574 (p12) fmpy.s1 S_lo = S_hi, poly1
2575 (p11) tbit.z.unc p14, p15 = cot_flag, 0 // p14=1 for tanl; p15=1 for cotl
2576 }
2577 { .mfi
2578       nop.m 999
2579 (p12) fma.s1 P = poly2, rsq, Q1_1
2580       nop.i 999 ;;
2581 }
2582 { .mfi
2583       nop.m 999
2584 //
2585 //    N odd:  poly1  =  S_hi * r + 1.0
2586 //
2587 //
2588 //    N odd:  S_lo  =  S_hi *  poly1
2589 //
2590 (p14) fadd.s0 Result = Poly, r          // for tanl
2591       nop.i 999
2592 }
2593 { .mfi
2594       nop.m 999
2595 (p15) fms.s0 Result = Poly, mOne, r     // for cotl
2596       nop.i 999 ;;
2597 }
2598
2599 { .mfi
2600       nop.m 999
2601 //
2602 //    N odd:  S_lo  =  Q1_1 * c + S_lo
2603 //
2604 (p12) fma.s1 S_lo = Q1_1, c, S_lo
2605       nop.i 999
2606 }
2607 { .mfi
2608       nop.m 999
2609       fmpy.s0 fp_tmp = fp_tmp, fp_tmp  // Dummy mult to set inexact
2610       nop.i 999 ;;
2611 }
2612 { .mfi
2613       nop.m 999
2614 //
2615 //    N odd:  Result =  S_lo + r * P
2616 //
2617 (p12) fma.s1 Result = P, r, S_lo
2618 (p12) tbit.z.unc p14, p15 = cot_flag, 0 ;; // p14=1 for tanl; p15=1 for cotl
2619 }
2620
2621 //
2622 //    N odd:  Result = Result + S_hi  (user supplied rounding mode)
2623 //
2624 { .mfi
2625       nop.m 999
2626 (p14) fadd.s0 Result = Result, S_hi         // for tanl
2627       nop.i 999
2628 }
2629 { .mfb
2630       nop.m 999
2631 (p15) fms.s0 Result = Result, mOne, S_hi    // for cotl
2632       br.ret.sptk b0 ;;              // Exit |r| < 1/4 path
2633 }
2634
2635
2636 TANL_NORMAL_R:
2637 // Here if 1/4 <= |x| < pi/4  or  if |x| >= 2^63 and |r| >= 1/4
2638 // *******************************************************************
2639 // *******************************************************************
2640 // *******************************************************************
2641 //
2642 //    r and c have been computed.
2643 //
2644 { .mfi
2645       nop.m 999
2646       fand B = B_mask1, r
2647       nop.i 999
2648 }
2649 ;;
2650
2651 TANL_NORMAL_R_A:
2652 // Enter here if pi/4 <= |x| < 2^63 and |r| >= 1/4
2653 //    Get the 5 bits or r for the lookup.   1.xxxxx ....
2654 { .mmi
2655       add table_ptr1 = 416, table_base     // Point to tanl_table_p2
2656       mov GR_exp_2tom65 = 0xffff - 65      // Scaling constant for B
2657       extr.u lookup = sig_r, 58, 5
2658 }
2659 ;;
2660
2661 { .mmi
2662       ldfe P2_1 = [table_ptr1], 16
2663       setf.exp TWO_TO_NEG65 = GR_exp_2tom65  // 2^-65 for scaling B if exp_r=-2
2664       add N_fix_gr = N_fix_gr, cot_flag      // N = N + 1 (for cotl)
2665 }
2666 ;;
2667
2668 .pred.rel "mutex",p11,p12
2669 //    B =  2^63 * 1.xxxxx 100...0
2670 { .mfi
2671       ldfe P2_2 = [table_ptr1], 16
2672       for B = B_mask2, B
2673       mov table_offset = 512               // Assume table offset is 512
2674 }
2675 ;;
2676
2677 { .mfi
2678       ldfe P2_3 = [table_ptr1], 16
2679       fmerge.s  Pos_r = f1, r
2680       tbit.nz p8,p9 = exp_r, 0
2681 }
2682 ;;
2683
2684 //    Is  B = 2** -2 or  B= 2** -1? If 2**-1, then
2685 //    we want an offset of 512 for table addressing.
2686 { .mii
2687       add table_ptr2 = 1296, table_base     // Point to tanl_table_cm2
2688 (p9)  shladd table_offset = lookup, 4, table_offset
2689 (p8)  shladd table_offset = lookup, 4, r0
2690 }
2691 ;;
2692
2693 { .mmi
2694       add table_ptr1 = table_ptr1, table_offset  // Point to T_hi
2695       add table_ptr2 = table_ptr2, table_offset  // Point to C_hi
2696       add table_ptr3 = 2128, table_base     // Point to tanl_table_scim2
2697 }
2698 ;;
2699
2700 { .mmi
2701       ldfd T_hi = [table_ptr1], 8                // Load T_hi
2702 ;;
2703       ldfd C_hi = [table_ptr2], 8                // Load C_hi
2704       add table_ptr3 = table_ptr3, table_offset  // Point to SC_inv
2705 }
2706 ;;
2707
2708 //
2709 //    x = |r| - B
2710 //
2711 //   Convert B so it has the same exponent as Pos_r before subtracting
2712 { .mfi
2713       ldfs T_lo = [table_ptr1]                   // Load T_lo
2714 (p9)  fnma.s1 x = B, FR_2tom64, Pos_r
2715       nop.i 999
2716 }
2717 { .mfi
2718       nop.m 999
2719 (p8)  fnma.s1 x = B, TWO_TO_NEG65, Pos_r
2720       nop.i 999
2721 }
2722 ;;
2723
2724 { .mfi
2725       ldfs C_lo = [table_ptr2]                   // Load C_lo
2726       nop.f 999
2727       nop.i 999
2728 }
2729 ;;
2730
2731 { .mfi
2732       ldfe SC_inv = [table_ptr3]                 // Load SC_inv
2733       fmerge.s  sgn_r = r, f1
2734       tbit.z p11, p12 = N_fix_gr, 0              // p11 if N even, p12 if odd
2735
2736 }
2737 ;;
2738
2739 //
2740 //    xsq = x * x
2741 //    N even: Tx = T_hi * x
2742 //
2743 //    N even: Tx1 = Tx + 1
2744 //    N odd:  Cx1 = 1 - Cx
2745 //
2746
2747 { .mfi
2748       nop.m 999
2749       fmpy.s1 xsq = x, x
2750       nop.i 999
2751 }
2752 { .mfi
2753       nop.m 999
2754 (p11) fmpy.s1 Tx = T_hi, x
2755       nop.i 999
2756 }
2757 ;;
2758
2759 //
2760 //    N odd: Cx = C_hi * x
2761 //
2762 { .mfi
2763       nop.m 999
2764 (p12) fmpy.s1 Cx = C_hi, x
2765       nop.i 999
2766 }
2767 ;;
2768 //
2769 //    N even and odd: P = P2_3 + P2_2 * xsq
2770 //
2771 { .mfi
2772       nop.m 999
2773       fma.s1 P = P2_3, xsq, P2_2
2774       nop.i 999
2775 }
2776 { .mfi
2777       nop.m 999
2778 (p11) fadd.s1 Tx1 = Tx, f1
2779       nop.i 999 ;;
2780 }
2781 { .mfi
2782       nop.m 999
2783 //
2784 //    N even: D = C_hi - tanx
2785 //    N odd: D = T_hi + tanx
2786 //
2787 (p11) fmpy.s1 CORR = SC_inv, T_hi
2788       nop.i 999
2789 }
2790 { .mfi
2791       nop.m 999
2792       fmpy.s1 Sx = SC_inv, x
2793       nop.i 999 ;;
2794 }
2795 { .mfi
2796       nop.m 999
2797 (p12) fmpy.s1 CORR = SC_inv, C_hi
2798       nop.i 999 ;;
2799 }
2800 { .mfi
2801       nop.m 999
2802 (p12) fsub.s1 V_hi = f1, Cx
2803       nop.i 999 ;;
2804 }
2805 { .mfi
2806       nop.m 999
2807       fma.s1 P = P, xsq, P2_1
2808       nop.i 999
2809 }
2810 { .mfi
2811       nop.m 999
2812 //
2813 //    N even and odd: P = P2_1 + P * xsq
2814 //
2815 (p11) fma.s1 V_hi = Tx, Tx1, f1
2816       nop.i 999 ;;
2817 }
2818 { .mfi
2819       nop.m 999
2820 //
2821 //    N even: Result  = sgn_r * tail + T_hi (user rounding mode for C1)
2822 //    N odd:  Result  = sgn_r * tail + C_hi (user rounding mode for C1)
2823 //
2824       fmpy.s0 fp_tmp = fp_tmp, fp_tmp  // Dummy mult to set inexact
2825       nop.i 999 ;;
2826 }
2827 { .mfi
2828       nop.m 999
2829       fmpy.s1 CORR = CORR, c
2830       nop.i 999 ;;
2831 }
2832 { .mfi
2833       nop.m 999
2834 (p12) fnma.s1 V_hi = Cx,V_hi,f1
2835       nop.i 999 ;;
2836 }
2837 { .mfi
2838       nop.m 999
2839 //
2840 //    N even: V_hi = Tx * Tx1 + 1
2841 //    N odd: Cx1 = 1 - Cx * Cx1
2842 //
2843       fmpy.s1 P = P, xsq
2844       nop.i 999
2845 }
2846 { .mfi
2847       nop.m 999
2848 //
2849 //    N even and odd: P = P * xsq
2850 //
2851 (p11) fmpy.s1 V_hi = V_hi, T_hi
2852       nop.i 999 ;;
2853 }
2854 { .mfi
2855       nop.m 999
2856 //
2857 //    N even and odd: tail = P * tail + V_lo
2858 //
2859 (p11) fmpy.s1 T_hi = sgn_r, T_hi
2860       nop.i 999 ;;
2861 }
2862 { .mfi
2863       nop.m 999
2864       fmpy.s1 CORR = CORR, sgn_r
2865       nop.i 999 ;;
2866 }
2867 { .mfi
2868       nop.m 999
2869 (p12) fmpy.s1 V_hi = V_hi,C_hi
2870       nop.i 999 ;;
2871 }
2872 { .mfi
2873       nop.m 999
2874 //
2875 //    N even: V_hi = T_hi * V_hi
2876 //    N odd: V_hi  = C_hi * V_hi
2877 //
2878       fma.s1 tanx = P, x, x
2879       nop.i 999
2880 }
2881 { .mfi
2882       nop.m 999
2883 (p12) fnmpy.s1 C_hi = sgn_r, C_hi
2884       nop.i 999 ;;
2885 }
2886 { .mfi
2887       nop.m 999
2888 //
2889 //    N even: V_lo = 1 - V_hi + C_hi
2890 //    N odd: V_lo = 1 - V_hi + T_hi
2891 //
2892 (p11) fadd.s1 CORR = CORR, T_lo
2893       nop.i 999
2894 }
2895 { .mfi
2896       nop.m 999
2897 (p12) fsub.s1 CORR = CORR, C_lo
2898       nop.i 999 ;;
2899 }
2900 { .mfi
2901       nop.m 999
2902 //
2903 //    N even and odd: tanx = x + x * P
2904 //    N even and odd: Sx = SC_inv * x
2905 //
2906 (p11) fsub.s1 D = C_hi, tanx
2907       nop.i 999
2908 }
2909 { .mfi
2910       nop.m 999
2911 (p12) fadd.s1 D = T_hi, tanx
2912       nop.i 999 ;;
2913 }
2914 { .mfi
2915       nop.m 999
2916 //
2917 //    N odd: CORR = SC_inv * C_hi
2918 //    N even: CORR = SC_inv * T_hi
2919 //
2920       fnma.s1 D = V_hi, D, f1
2921       nop.i 999 ;;
2922 }
2923 { .mfi
2924       nop.m 999
2925 //
2926 //    N even and odd: D = 1 - V_hi * D
2927 //    N even and odd: CORR = CORR * c
2928 //
2929       fma.s1 V_hi = V_hi, D, V_hi
2930       nop.i 999 ;;
2931 }
2932 { .mfi
2933       nop.m 999
2934 //
2935 //    N even and odd: V_hi = V_hi + V_hi * D
2936 //    N even and odd: CORR = sgn_r * CORR
2937 //
2938 (p11) fnma.s1 V_lo = V_hi, C_hi, f1
2939       nop.i 999
2940 }
2941 { .mfi
2942       nop.m 999
2943 (p12) fnma.s1 V_lo = V_hi, T_hi, f1
2944       nop.i 999 ;;
2945 }
2946 { .mfi
2947       nop.m 999
2948 //
2949 //    N even: CORR = COOR + T_lo
2950 //    N odd: CORR = CORR - C_lo
2951 //
2952 (p11) fma.s1 V_lo = tanx, V_hi, V_lo
2953       tbit.nz p15, p0 = cot_flag, 0       // p15=1 if we compute cotl
2954 }
2955 { .mfi
2956       nop.m 999
2957 (p12) fnma.s1 V_lo = tanx, V_hi, V_lo
2958       nop.i 999 ;;
2959 }
2960
2961 { .mfi
2962       nop.m 999
2963 (p15) fms.s1 T_hi = f0, f0, T_hi        // to correct result's sign for cotl
2964       nop.i 999
2965 }
2966 { .mfi
2967       nop.m 999
2968 (p15) fms.s1 C_hi = f0, f0, C_hi        // to correct result's sign for cotl
2969       nop.i 999
2970 };;
2971
2972 { .mfi
2973       nop.m 999
2974 (p15) fms.s1 sgn_r = f0, f0, sgn_r      // to correct result's sign for cotl
2975       nop.i 999
2976 };;
2977
2978 { .mfi
2979       nop.m 999
2980 //
2981 //    N even: V_lo = V_lo + V_hi * tanx
2982 //    N odd: V_lo = V_lo - V_hi * tanx
2983 //
2984 (p11) fnma.s1 V_lo = C_lo, V_hi, V_lo
2985       nop.i 999
2986 }
2987 { .mfi
2988       nop.m 999
2989 (p12) fnma.s1 V_lo = T_lo, V_hi, V_lo
2990       nop.i 999 ;;
2991 }
2992 { .mfi
2993       nop.m 999
2994 //
2995 //    N  even: V_lo = V_lo - V_hi * C_lo
2996 //    N  odd: V_lo = V_lo - V_hi * T_lo
2997 //
2998       fmpy.s1 V_lo = V_hi, V_lo
2999       nop.i 999 ;;
3000 }
3001 { .mfi
3002       nop.m 999
3003 //
3004 //    N even and odd: V_lo = V_lo * V_hi
3005 //
3006       fadd.s1 tail = V_hi, V_lo
3007       nop.i 999 ;;
3008 }
3009 { .mfi
3010       nop.m 999
3011 //
3012 //    N even and odd: tail = V_hi + V_lo
3013 //
3014       fma.s1 tail = tail, P, V_lo
3015       nop.i 999 ;;
3016 }
3017 { .mfi
3018       nop.m 999
3019 //
3020 //    N even: T_hi = sgn_r * T_hi
3021 //    N odd : C_hi = -sgn_r * C_hi
3022 //
3023       fma.s1 tail = tail, Sx, CORR
3024       nop.i 999 ;;
3025 }
3026 { .mfi
3027       nop.m 999
3028 //
3029 //    N even and odd: tail = Sx * tail + CORR
3030 //
3031       fma.s1 tail = V_hi, Sx, tail
3032       nop.i 999 ;;
3033 }
3034 { .mfi
3035       nop.m 999
3036 //
3037 //    N even an odd: tail = Sx * V_hi + tail
3038 //
3039 (p11) fma.s0 Result = sgn_r, tail, T_hi
3040       nop.i 999
3041 }
3042 { .mfb
3043       nop.m 999
3044 (p12) fma.s0 Result = sgn_r, tail, C_hi
3045       br.ret.sptk b0 ;;                 // Exit for 1/4 <= |r| < pi/4
3046 }
3047
3048 TANL_DENORMAL:
3049 // Here if x denormal
3050 { .mfb
3051       getf.exp GR_signexp_x = Norm_Arg          // Get sign and exponent of x
3052       nop.f 999
3053       br.cond.sptk TANL_COMMON                  // Return to common code
3054 }
3055 ;;
3056
3057
3058 TANL_SPECIAL:
3059 TANL_UNSUPPORTED:
3060 //
3061 //     Code for NaNs, Unsupporteds, Infs, or +/- zero ?
3062 //     Invalid raised for Infs and SNaNs.
3063 //
3064
3065 { .mfi
3066       nop.m 999
3067       fmerge.s  f10 = f8, f8            // Save input for error call
3068       tbit.nz p6, p7 = cot_flag, 0      // p6=1 if we compute cotl
3069 }
3070 ;;
3071
3072 { .mfi
3073       nop.m 999
3074 (p6)  fclass.m p6, p7 = f8, 0x7         // Test for zero (cotl only)
3075       nop.i 999
3076 }
3077 ;;
3078
3079 .pred.rel "mutex", p6, p7
3080 { .mfi
3081 (p6)  mov GR_Parameter_Tag = 225        // (cotl)
3082 (p6)  frcpa.s0  f8, p0 = f1, f8         // cotl(+-0) = +-Inf
3083       nop.i 999
3084 }
3085 { .mfb
3086       nop.m 999
3087 (p7)  fmpy.s0 f8 = f8, f0
3088 (p7)  br.ret.sptk b0
3089 }
3090 ;;
3091
3092 GLOBAL_IEEE754_END(tanl)
3093 libm_alias_ldouble_other (__tan, tan)
3094
3095
3096 LOCAL_LIBM_ENTRY(__libm_error_region)
3097 .prologue
3098
3099 // (1)
3100 { .mfi
3101       add           GR_Parameter_Y=-32,sp        // Parameter 2 value
3102       nop.f         0
3103 .save   ar.pfs,GR_SAVE_PFS
3104       mov           GR_SAVE_PFS=ar.pfs           // Save ar.pfs
3105 }
3106 { .mfi
3107 .fframe 64
3108       add sp=-64,sp                              // Create new stack
3109       nop.f 0
3110       mov GR_SAVE_GP=gp                          // Save gp
3111 };;
3112
3113 // (2)
3114 { .mmi
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
3119 };;
3120
3121 .body
3122 // (3)
3123 { .mib
3124       stfe [GR_Parameter_X] = f10                // STORE Parameter 1 on stack
3125       add   GR_Parameter_RESULT = 0,GR_Parameter_Y  // Parameter 3 address
3126       nop.b 0
3127 }
3128 { .mib
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
3132 };;
3133 { .mmi
3134       nop.m 0
3135       nop.m 0
3136       add   GR_Parameter_RESULT = 48,sp
3137 };;
3138
3139 // (4)
3140 { .mmi
3141       ldfe  f8 = [GR_Parameter_RESULT]           // Get return result off stack
3142 .restore sp
3143       add   sp = 64,sp                           // Restore stack pointer
3144       mov   b0 = GR_SAVE_B0                      // Restore return address
3145 };;
3146 { .mib
3147       mov   gp = GR_SAVE_GP                      // Restore gp
3148       mov   ar.pfs = GR_SAVE_PFS                 // Restore ar.pfs
3149       br.ret.sptk     b0                         // Return
3150 };;
3151
3152 LOCAL_LIBM_END(__libm_error_region)
3153
3154 .type   __libm_error_support#,@function
3155 .global __libm_error_support#
3156
3157
3158 // *******************************************************************
3159 // *******************************************************************
3160 // *******************************************************************
3161 //
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:
3165 //       On input:
3166 //         (Arg or x) is in f8
3167 //       On output:
3168 //         r is in f8
3169 //         c is in f9
3170 //         N is in r8
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
3173 //     function.
3174 //
3175 // *******************************************************************
3176 // *******************************************************************
3177 // *******************************************************************
3178
3179 LOCAL_LIBM_ENTRY(__libm_callout)
3180 TANL_ARG_TOO_LARGE:
3181 .prologue
3182 { .mfi
3183       add table_ptr2 = 144, table_base        // Point to 2^-2
3184       nop.f 999
3185 .save   ar.pfs,GR_SAVE_PFS
3186       mov  GR_SAVE_PFS=ar.pfs                 // Save ar.pfs
3187 }
3188 ;;
3189
3190 //     Load 2^-2, -2^-2
3191 { .mmi
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
3196 };;
3197
3198 .body
3199 //
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
3203 //
3204 { .mib
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#
3208 }
3209 ;;
3210
3211 //
3212 //     Is |r| < 2**(-2)
3213 //
3214 { .mfi
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
3218 }
3219 ;;
3220
3221 { .mfi
3222       getf.exp exp_r = r                     // Extract signexp of r
3223       nop.f 999
3224       mov    b0 = GR_SAVE_B0                 // Restore return address
3225 }
3226 ;;
3227
3228 //
3229 //     Get N_fix_gr
3230 //
3231 { .mfi
3232       mov   N_fix_gr = r8
3233 (p6)  fcmp.gt.unc.s1  p6, p0 = r, NEGTWO_TO_NEG2
3234       mov   ar.pfs = GR_SAVE_PFS             // Restore pfs
3235 }
3236 ;;
3237
3238 { .mbb
3239       nop.m 999
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
3242 }
3243 ;;
3244
3245 LOCAL_LIBM_END(__libm_callout)
3246
3247 .type __libm_pi_by_2_reduce#,@function
3248 .global __libm_pi_by_2_reduce#