sysdeps/ia64/fpu/e_asin.S

   1 .file "asin.s"
   2
   3
   4 // Copyright (c) 2000 - 2003 Intel Corporation
   5 // All rights reserved.
   6 //
   7 // Contributed 2000 by the Intel Numerics Group, Intel Corporation
   8 //
   9 // Redistribution and use in source and binary forms, with or without
  10 // modification, are permitted provided that the following conditions are
  11 // met:
  12 //
  13 // * Redistributions of source code must retain the above copyright
  14 // notice, this list of conditions and the following disclaimer.
  15 //
  16 // * Redistributions in binary form must reproduce the above copyright
  17 // notice, this list of conditions and the following disclaimer in the
  18 // documentation and/or other materials provided with the distribution.
  19 //
  20 // * The name of Intel Corporation may not be used to endorse or promote
  21 // products derived from this software without specific prior written
  22 // permission.
  23
  24 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
  25 // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
  26 // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
  27 // A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL INTEL OR ITS
  28 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
  29 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
  30 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
  31 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY
  32 // OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY OR TORT (INCLUDING
  33 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
  34 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  35 //
  36 // Intel Corporation is the author of this code, and requests that all
  37 // problem reports or change requests be submitted to it directly at
  38 // http://www.intel.com/software/products/opensource/libraries/num.htm.
  39
  40 // History
  41 //==============================================================
  42 // 02/02/00 Initial version
  43 // 08/17/00 New and much faster algorithm.
  44 // 08/31/00 Avoided bank conflicts on loads, shortened |x|=1 path,
  45 //          fixed mfb split issue stalls.
  46 // 12/19/00 Fixed small arg cases to force inexact, or inexact and underflow.
  47 // 08/02/02 New and much faster algorithm II
  48 // 02/06/03 Reordered header: .section, .global, .proc, .align
  49
  50 // Description
  51 //=========================================
  52 // The asin function computes the principal value of the arc sine of x.
  53 // asin(0) returns 0, asin(1) returns pi/2, asin(-1) returns -pi/2.
  54 // A doman error occurs for arguments not in the range [-1,+1].
  55 //
  56 // The asin function returns the arc sine in the range [-pi/2, +pi/2] radians.
  57 //
  58 // There are 8 paths:
  59 // 1. x = +/-0.0
  60 //    Return asin(x) = +/-0.0
  61 //
  62 // 2. 0.0 < |x| < 0.625
  63 //    Return asin(x) = x + x^3 *PolA(x^2)
  64 //    where PolA(x^2) = A3 + A5*x^2 + A7*x^4 +...+ A35*x^32
  65 //
  66 // 3. 0.625 <=|x| < 1.0
  67 //    Return asin(x) = sign(x) * ( Pi/2 - sqrt(R) * PolB(R))
  68 //    Where R = 1 - |x|,
  69 //          PolB(R) = B0 + B1*R + B2*R^2 +...+B12*R^12
  70 //
  71 //    sqrt(R) is approximated using the following sequence:
  72 //        y0 = (1 + eps)/sqrt(R) - initial approximation by frsqrta,
  73 //             |eps| < 2^(-8)
  74 //        Then 3 iterations are used to refine the result:
  75 //        H0 = 0.5*y0
  76 //        S0 = R*y0
  77 //
  78 //        d0 = 0.5 - H0*S0
  79 //        H1 = H0 + d0*H0
  80 //        S1 = S0 + d0*S0
  81 //
  82 //        d1 = 0.5 - H1*S1
  83 //        H2 = H1 + d0*H1
  84 //        S2 = S1 + d0*S1
  85 //
  86 //        d2 = 0.5 - H2*S2
  87 //        S3 = S3 + d2*S3
  88 //
  89 //        S3 approximates sqrt(R) with enough accuracy for this algorithm
  90 //
  91 //    So, the result should be reconstracted as follows:
  92 //    asin(x) = sign(x) * (Pi/2 - S3*PolB(R))
  93 //
  94 //    But for optimization perposes the reconstruction step is slightly
  95 //    changed:
  96 //    asin(x) = sign(x)*(Pi/2 - PolB(R)*S2) + sign(x)*d2*S2*PolB(R)
  97 //
  98 // 4. |x| = 1.0
  99 //    Return asin(x) = sign(x)*Pi/2
 100 //
 101 // 5. 1.0 < |x| <= +INF
 102 //    A doman error occurs for arguments not in the range [-1,+1]
 103 //
 104 // 6. x = [S,Q]NaN
 105 //    Return asin(x) = QNaN
 106 //
 107 // 7. x is denormal
 108 //    Return asin(x) = x + x^3,
 109 //
 110 // 8. x is unnormal
 111 //    Normalize input in f8 and return to the very beginning of the function
 112 //
 113 // Registers used
 114 //==============================================================
 115 // Floating Point registers used:
 116 // f8, input, output
 117 // f6, f7, f9 -> f15, f32 -> f63
 118
 119 // General registers used:
 120 // r3, r21 -> r31, r32 -> r38
 121
 122 // Predicate registers used:
 123 // p0, p6 -> p14
 124
 125 //
 126 // Assembly macros
 127 //=========================================
 128 // integer registers used
 129 // scratch
 130 rTblAddr                      = r3
 131
 132 rPiBy2Ptr                     = r21
 133 rTmpPtr3                      = r22
 134 rDenoBound                    = r23
 135 rOne                          = r24
 136 rAbsXBits                     = r25
 137 rHalf                         = r26
 138 r0625                         = r27
 139 rSign                         = r28
 140 rXBits                        = r29
 141 rTmpPtr2                      = r30
 142 rTmpPtr1                      = r31
 143
 144 // stacked
 145 GR_SAVE_PFS                   = r32
 146 GR_SAVE_B0                    = r33
 147 GR_SAVE_GP                    = r34
 148 GR_Parameter_X                = r35
 149 GR_Parameter_Y                = r36
 150 GR_Parameter_RESULT           = r37
 151 GR_Parameter_TAG              = r38
 152
 153 // floating point registers used
 154 FR_X                          = f10
 155 FR_Y                          = f1
 156 FR_RESULT                     = f8
 157
 158
 159 // scratch
 160 fXSqr                         = f6
 161 fXCube                        = f7
 162 fXQuadr                       = f9
 163 f1pX                          = f10
 164 f1mX                          = f11
 165 f1pXRcp                       = f12
 166 f1mXRcp                       = f13
 167 fH                            = f14
 168 fS                            = f15
 169 // stacked
 170 fA3                           = f32
 171 fB1                           = f32
 172 fA5                           = f33
 173 fB2                           = f33
 174 fA7                           = f34
 175 fPiBy2                        = f34
 176 fA9                           = f35
 177 fA11                          = f36
 178 fB10                          = f35
 179 fB11                          = f36
 180 fA13                          = f37
 181 fA15                          = f38
 182 fB4                           = f37
 183 fB5                           = f38
 184 fA17                          = f39
 185 fA19                          = f40
 186 fB6                           = f39
 187 fB7                           = f40
 188 fA21                          = f41
 189 fA23                          = f42
 190 fB3                           = f41
 191 fB8                           = f42
 192 fA25                          = f43
 193 fA27                          = f44
 194 fB9                           = f43
 195 fB12                          = f44
 196 fA29                          = f45
 197 fA31                          = f46
 198 fA33                          = f47
 199 fA35                          = f48
 200 fBaseP                        = f49
 201 fB0                           = f50
 202 fSignedS                      = f51
 203 fD                            = f52
 204 fHalf                         = f53
 205 fR                            = f54
 206 fCloseTo1Pol                  = f55
 207 fSignX                        = f56
 208 fDenoBound                    = f57
 209 fNormX                        = f58
 210 fX8                           = f59
 211 fRSqr                         = f60
 212 fRQuadr                       = f61
 213 fR8                           = f62
 214 fX16                          = f63
 215 // Data tables
 216 //==============================================================
 217 RODATA
 218 .align 16
 219 LOCAL_OBJECT_START(asin_base_range_table)
 220 // Ai: Polynomial coefficients for the asin(x), |x| < .625000
 221 // Bi: Polynomial coefficients for the asin(x), |x| > .625000
 222 data8 0xBFDAAB56C01AE468 //A29
 223 data8 0x3FE1C470B76A5B2B //A31
 224 data8 0xBFDC5FF82A0C4205 //A33
 225 data8 0x3FC71FD88BFE93F0 //A35
 226 data8 0xB504F333F9DE6487, 0x00003FFF //B0
 227 data8 0xAAAAAAAAAAAAFC18, 0x00003FFC //A3
 228 data8 0x3F9F1C71BC4A7823 //A9
 229 data8 0x3F96E8BBAAB216B2 //A11
 230 data8 0x3F91C4CA1F9F8A98 //A13
 231 data8 0x3F8C9DDCEDEBE7A6 //A15
 232 data8 0x3F877784442B1516 //A17
 233 data8 0x3F859C0491802BA2 //A19
 234 data8 0x9999999998C88B8F, 0x00003FFB //A5
 235 data8 0x3F6BD7A9A660BF5E //A21
 236 data8 0x3F9FC1659340419D //A23
 237 data8 0xB6DB6DB798149BDF, 0x00003FFA //A7
 238 data8 0xBFB3EF18964D3ED3 //A25
 239 data8 0x3FCD285315542CF2 //A27
 240 data8 0xF15BEEEFF7D2966A, 0x00003FFB //B1
 241 data8 0x3EF0DDA376D10FB3 //B10
 242 data8 0xBEB83CAFE05EBAC9 //B11
 243 data8 0x3F65FFB67B513644 //B4
 244 data8 0x3F5032FBB86A4501 //B5
 245 data8 0x3F392162276C7CBA //B6
 246 data8 0x3F2435949FD98BDF //B7
 247 data8 0xD93923D7FA08341C, 0x00003FF9 //B2
 248 data8 0x3F802995B6D90BDB //B3
 249 data8 0x3F10DF86B341A63F //B8
 250 data8 0xC90FDAA22168C235, 0x00003FFF // Pi/2
 251 data8 0x3EFA3EBD6B0ECB9D //B9
 252 data8 0x3EDE18BA080E9098 //B12
 253 LOCAL_OBJECT_END(asin_base_range_table)
 254
 255
 256 .section .text
 257 GLOBAL_LIBM_ENTRY(asin)
 258 asin_unnormal_back:
 259 { .mfi
 260       getf.d             rXBits = f8 // grab bits of input value
 261       // set p12 = 1 if x is a NaN, denormal, or zero
 262       fclass.m           p12, p0 = f8, 0xcf
 263       adds               rSign = 1, r0
 264 }
 265 { .mfi
 266       addl               rTblAddr = @ltoff(asin_base_range_table),gp
 267       // 1 - x = 1 - |x| for positive x
 268       fms.s1             f1mX = f1, f1, f8
 269       addl               rHalf = 0xFFFE, r0 // exponent of 1/2
 270 }
 271 ;;
 272 { .mfi
 273       addl               r0625 = 0x3FE4, r0 // high 16 bits of 0.625
 274       // set p8 = 1 if x < 0
 275       fcmp.lt.s1         p8, p9 = f8, f0
 276       shl                rSign = rSign, 63 // sign bit
 277 }
 278 { .mfi
 279       // point to the beginning of the table
 280       ld8                rTblAddr = [rTblAddr]
 281       // 1 + x = 1 - |x| for negative x
 282       fma.s1             f1pX = f1, f1, f8
 283       adds               rOne = 0x3FF, r0
 284 }
 285 ;;
 286 { .mfi
 287       andcm              rAbsXBits = rXBits, rSign // bits of |x|
 288       fmerge.s           fSignX = f8, f1 // signum(x)
 289       shl                r0625 = r0625, 48 // bits of DP representation of 0.625
 290 }
 291 { .mfb
 292       setf.exp           fHalf = rHalf // load A2 to FP reg
 293       fma.s1             fXSqr = f8, f8, f0 // x^2
 294       // branch on special path if x is a NaN, denormal, or zero
 295 (p12) br.cond.spnt       asin_special
 296 }
 297 ;;
 298 { .mfi
 299       adds               rPiBy2Ptr = 272, rTblAddr
 300       nop.f              0
 301       shl                rOne = rOne, 52 // bits of 1.0
 302 }
 303 { .mfi
 304       adds               rTmpPtr1 = 16, rTblAddr
 305       nop.f              0
 306       // set p6 = 1 if |x| < 0.625
 307       cmp.lt             p6, p7 = rAbsXBits, r0625
 308 }
 309 ;;
 310 { .mfi
 311       ldfpd              fA29, fA31 = [rTblAddr] // A29, fA31
 312       // 1 - x = 1 - |x| for positive x
 313 (p9)  fms.s1             fR = f1, f1, f8
 314       // point to coefficient of "near 1" polynomial
 315 (p7)  adds               rTmpPtr2 = 176, rTblAddr
 316 }
 317 { .mfi
 318       ldfpd              fA33, fA35 = [rTmpPtr1], 16 // A33, fA35
 319       // 1 + x = 1 - |x| for negative x
 320 (p8)  fma.s1             fR = f1, f1, f8
 321 (p6)  adds               rTmpPtr2 = 48, rTblAddr
 322 }
 323 ;;
 324 { .mfi
 325       ldfe               fB0 = [rTmpPtr1], 16 // B0
 326       nop.f              0
 327       nop.i              0
 328 }
 329 { .mib
 330       adds               rTmpPtr3 = 16, rTmpPtr2
 331       // set p10 = 1 if |x| = 1.0
 332       cmp.eq             p10, p0 = rAbsXBits, rOne
 333       // branch on special path for |x| = 1.0
 334 (p10) br.cond.spnt       asin_abs_1
 335 }
 336 ;;
 337 { .mfi
 338       ldfe               fA3 = [rTmpPtr2], 48 // A3 or B1
 339       nop.f              0
 340       adds               rTmpPtr1 = 64, rTmpPtr3
 341 }
 342 { .mib
 343       ldfpd              fA9, fA11 = [rTmpPtr3], 16 // A9, A11 or B10, B11
 344       // set p11 = 1 if |x| > 1.0
 345       cmp.gt             p11, p0 = rAbsXBits, rOne
 346       // branch on special path for |x| > 1.0
 347 (p11) br.cond.spnt       asin_abs_gt_1
 348 }
 349 ;;
 350 { .mfi
 351       ldfpd              fA17, fA19 = [rTmpPtr2], 16 // A17, A19 or B6, B7
 352       // initial approximation of 1 / sqrt(1 - x)
 353       frsqrta.s1         f1mXRcp, p0 = f1mX
 354       nop.i              0
 355 }
 356 { .mfi
 357       ldfpd              fA13, fA15 = [rTmpPtr3] // A13, A15 or B4, B5
 358       fma.s1             fXCube = fXSqr, f8, f0 // x^3
 359       nop.i              0
 360 }
 361 ;;
 362 { .mfi
 363       ldfe               fA5 = [rTmpPtr2], 48 // A5 or B2
 364       // initial approximation of 1 / sqrt(1 + x)
 365       frsqrta.s1         f1pXRcp, p0 = f1pX
 366       nop.i              0
 367 }
 368 { .mfi
 369       ldfpd              fA21, fA23 = [rTmpPtr1], 16 // A21, A23 or B3, B8
 370       fma.s1             fXQuadr = fXSqr, fXSqr, f0 // x^4
 371       nop.i              0
 372 }
 373 ;;
 374 { .mfi
 375       ldfe               fA7 = [rTmpPtr1] // A7 or Pi/2
 376       fma.s1             fRSqr = fR, fR, f0 // R^2
 377       nop.i              0
 378 }
 379 { .mfb
 380       ldfpd              fA25, fA27 = [rTmpPtr2] // A25, A27 or B9, B12
 381       nop.f              0
 382 (p6)  br.cond.spnt       asin_base_range;
 383 }
 384 ;;
 385
 386 { .mfi
 387       nop.m              0
 388 (p9)  fma.s1             fH = fHalf, f1mXRcp, f0 // H0 for x > 0
 389       nop.i              0
 390 }
 391 { .mfi
 392       nop.m              0
 393 (p9)  fma.s1             fS = f1mX, f1mXRcp, f0  // S0 for x > 0
 394       nop.i              0
 395 }
 396 ;;
 397 { .mfi
 398       nop.m              0
 399 (p8)  fma.s1             fH = fHalf, f1pXRcp, f0 // H0 for x < 0
 400       nop.i              0
 401 }
 402 { .mfi
 403       nop.m              0
 404 (p8)  fma.s1             fS = f1pX, f1pXRcp, f0  // S0 for x > 0
 405       nop.i              0
 406 }
 407 ;;
 408 { .mfi
 409       nop.m              0
 410       fma.s1             fRQuadr = fRSqr, fRSqr, f0 // R^4
 411       nop.i              0
 412 }
 413 ;;
 414 { .mfi
 415       nop.m              0
 416       fma.s1             fB11 = fB11, fR, fB10
 417       nop.i              0
 418 }
 419 { .mfi
 420       nop.m              0
 421       fma.s1             fB1 = fB1, fR, fB0
 422       nop.i              0
 423 }
 424 ;;
 425 { .mfi
 426       nop.m              0
 427       fma.s1             fB5 = fB5, fR, fB4
 428       nop.i              0
 429 }
 430 { .mfi
 431       nop.m              0
 432       fma.s1             fB7 = fB7, fR, fB6
 433       nop.i              0
 434 }
 435 ;;
 436 { .mfi
 437       nop.m              0
 438       fma.s1             fB3 = fB3, fR, fB2
 439       nop.i              0
 440 }
 441 ;;
 442 { .mfi
 443       nop.m              0
 444       fnma.s1            fD = fH, fS, fHalf // d0 = 1/2 - H0*S0
 445       nop.i              0
 446 }
 447 ;;
 448 { .mfi
 449       nop.m              0
 450       fma.s1             fR8 = fRQuadr, fRQuadr, f0 // R^4
 451       nop.i              0
 452 }
 453 { .mfi
 454       nop.m              0
 455       fma.s1             fB9 = fB9, fR, fB8
 456       nop.i              0
 457 }
 458 ;;
 459 {.mfi
 460       nop.m              0
 461       fma.s1             fB12 = fB12, fRSqr, fB11
 462       nop.i              0
 463 }
 464 {.mfi
 465       nop.m              0
 466       fma.s1             fB7 = fB7, fRSqr, fB5
 467       nop.i              0
 468 }
 469 ;;
 470 {.mfi
 471       nop.m              0
 472       fma.s1             fB3 = fB3, fRSqr, fB1
 473       nop.i              0
 474 }
 475 ;;
 476 { .mfi
 477       nop.m              0
 478       fma.s1             fH = fH, fD, fH // H1 = H0 + H0*d0
 479       nop.i              0
 480 }
 481 { .mfi
 482       nop.m              0
 483       fma.s1             fS = fS, fD, fS // S1 = S0 + S0*d0
 484       nop.i              0
 485 }
 486 ;;
 487 {.mfi
 488       nop.m              0
 489       fma.s1             fPiBy2 = fPiBy2, fSignX, f0 // signum(x)*Pi/2
 490       nop.i              0
 491 }
 492 ;;
 493 { .mfi
 494       nop.m              0
 495       fma.s1             fB12 = fB12, fRSqr, fB9
 496       nop.i              0
 497 }
 498 { .mfi
 499       nop.m              0
 500       fma.s1             fB7 = fB7, fRQuadr, fB3
 501       nop.i              0
 502 }
 503 ;;
 504 {.mfi
 505       nop.m              0
 506       fnma.s1            fD = fH, fS, fHalf // d1 = 1/2 - H1*S1
 507       nop.i              0
 508 }
 509 { .mfi
 510       nop.m              0
 511       fnma.s1            fSignedS = fSignX, fS, f0 // -signum(x)*S1
 512       nop.i              0
 513 }
 514 ;;
 515 { .mfi
 516       nop.m              0
 517       fma.s1             fCloseTo1Pol = fB12, fR8, fB7
 518       nop.i              0
 519 }
 520 ;;
 521 { .mfi
 522       nop.m              0
 523       fma.s1             fH = fH, fD, fH // H2 = H1 + H1*d1
 524       nop.i              0
 525 }
 526 { .mfi
 527       nop.m              0
 528       fma.s1             fS = fS, fD, fS // S2 = S1 + S1*d1
 529       nop.i              0
 530 }
 531 ;;
 532 { .mfi
 533       nop.m              0
 534       // -signum(x)* S2 = -signum(x)*(S1 + S1*d1)
 535       fma.s1             fSignedS = fSignedS, fD, fSignedS
 536       nop.i              0
 537 }
 538 ;;
 539 {.mfi
 540       nop.m              0
 541       fnma.s1            fD = fH, fS, fHalf // d2 = 1/2 - H2*S2
 542       nop.i              0
 543 }
 544 ;;
 545 { .mfi
 546       nop.m              0
 547       // signum(x)*(Pi/2 - PolB*S2)
 548       fma.s1             fPiBy2 = fSignedS, fCloseTo1Pol, fPiBy2
 549       nop.i              0
 550 }
 551 { .mfi
 552       nop.m              0
 553       // -signum(x)*PolB * S2
 554       fma.s1             fCloseTo1Pol = fSignedS, fCloseTo1Pol, f0
 555       nop.i              0
 556 }
 557 ;;
 558 { .mfb
 559       nop.m              0
 560       // final result for 0.625 <= |x| < 1
 561       fma.d.s0           f8 = fCloseTo1Pol, fD, fPiBy2
 562       // exit here for  0.625 <= |x| < 1
 563       br.ret.sptk        b0
 564 }
 565 ;;
 566
 567
 568 // here if |x| < 0.625
 569 .align 32
 570 asin_base_range:
 571 { .mfi
 572       nop.m              0
 573       fma.s1             fA33 = fA33, fXSqr, fA31
 574       nop.i              0
 575 }
 576 { .mfi
 577       nop.m              0
 578       fma.s1             fA15 = fA15, fXSqr, fA13
 579       nop.i              0
 580 }
 581 ;;
 582 { .mfi
 583       nop.m              0
 584       fma.s1             fA29 = fA29, fXSqr, fA27
 585       nop.i              0
 586 }
 587 { .mfi
 588       nop.m              0
 589       fma.s1             fA25 = fA25, fXSqr, fA23
 590       nop.i              0
 591 }
 592 ;;
 593 { .mfi
 594       nop.m              0
 595       fma.s1             fA21 = fA21, fXSqr, fA19
 596       nop.i              0
 597 }
 598 { .mfi
 599       nop.m              0
 600       fma.s1             fA9 = fA9, fXSqr, fA7
 601       nop.i              0
 602 }
 603 ;;
 604 { .mfi
 605       nop.m              0
 606       fma.s1             fA5 = fA5, fXSqr, fA3
 607       nop.i              0
 608 }
 609 ;;
 610 { .mfi
 611       nop.m              0
 612       fma.s1             fA35 = fA35, fXQuadr, fA33
 613       nop.i              0
 614 }
 615 { .mfi
 616       nop.m              0
 617       fma.s1             fA17 = fA17, fXQuadr, fA15
 618       nop.i              0
 619 }
 620 ;;
 621 { .mfi
 622       nop.m              0
 623       fma.s1             fX8 = fXQuadr, fXQuadr, f0 // x^8
 624       nop.i              0
 625 }
 626 { .mfi
 627       nop.m              0
 628       fma.s1             fA25 = fA25, fXQuadr, fA21
 629       nop.i              0
 630 }
 631 ;;
 632 { .mfi
 633       nop.m              0
 634       fma.s1             fA9 = fA9, fXQuadr, fA5
 635       nop.i              0
 636 }
 637 ;;
 638 { .mfi
 639       nop.m              0
 640       fma.s1             fA35 = fA35, fXQuadr, fA29
 641       nop.i              0
 642 }
 643 { .mfi
 644       nop.m              0
 645       fma.s1             fA17 = fA17, fXSqr, fA11
 646       nop.i              0
 647 }
 648 ;;
 649 { .mfi
 650       nop.m              0
 651       fma.s1             fX16 = fX8, fX8, f0 // x^16
 652       nop.i              0
 653 }
 654 ;;
 655 { .mfi
 656       nop.m              0
 657       fma.s1             fA35 = fA35, fX8, fA25
 658       nop.i              0
 659 }
 660 { .mfi
 661       nop.m              0
 662       fma.s1             fA17 = fA17, fX8, fA9
 663       nop.i              0
 664 }
 665 ;;
 666 { .mfi
 667       nop.m              0
 668       fma.s1             fBaseP = fA35, fX16, fA17
 669       nop.i              0
 670 }
 671 ;;
 672 { .mfb
 673       nop.m              0
 674       // final result for |x| < 0.625
 675       fma.d.s0           f8 = fBaseP, fXCube, f8
 676       // exit here for |x| < 0.625 path
 677       br.ret.sptk        b0
 678 }
 679 ;;
 680
 681 // here if |x| = 1
 682 // asin(x) = sign(x) * Pi/2
 683 .align 32
 684 asin_abs_1:
 685 { .mfi
 686       ldfe               fPiBy2 = [rPiBy2Ptr] // Pi/2
 687       nop.f              0
 688       nop.i              0
 689 }
 690 ;;
 691 {.mfb
 692       nop.m              0
 693       // result for |x| = 1.0
 694       fma.d.s0           f8 = fPiBy2, fSignX, f0
 695       // exit here for |x| = 1.0
 696       br.ret.sptk        b0
 697 }
 698 ;;
 699
 700 // here if x is a NaN, denormal, or zero
 701 .align 32
 702 asin_special:
 703 { .mfi
 704       nop.m              0
 705       // set p12 = 1 if x is a NaN
 706       fclass.m           p12, p0 = f8, 0xc3
 707       nop.i              0
 708 }
 709 { .mlx
 710       nop.m              0
 711       // smallest positive DP normalized number
 712       movl               rDenoBound = 0x0010000000000000
 713 }
 714 ;;
 715 { .mfi
 716       nop.m              0
 717       // set p13 = 1 if x = 0.0
 718       fclass.m           p13, p0 = f8, 0x07
 719       nop.i              0
 720 }
 721 { .mfi
 722       nop.m              0
 723       fnorm.s1           fNormX = f8
 724       nop.i              0
 725 }
 726 ;;
 727 { .mfb
 728       // load smallest normal to FP reg
 729       setf.d             fDenoBound = rDenoBound
 730       // answer if x is a NaN
 731 (p12) fma.d.s0           f8 = f8,f1,f0
 732       // exit here if x is a NaN
 733 (p12) br.ret.spnt        b0
 734 }
 735 ;;
 736 { .mfb
 737       nop.m              0
 738       nop.f              0
 739       // exit here if x = 0.0
 740 (p13) br.ret.spnt        b0
 741 }
 742 ;;
 743 // if we still here then x is denormal or unnormal
 744 { .mfi
 745       nop.m              0
 746       // absolute value of normalized x
 747       fmerge.s           fNormX = f1, fNormX
 748       nop.i              0
 749 }
 750 ;;
 751 { .mfi
 752       nop.m              0
 753       // set p14 = 1 if normalized x is greater than or
 754       // equal to the smallest denormalized value
 755       // So, if p14 is set to 1 it means that we deal with
 756       // unnormal rather than with "true" denormal
 757       fcmp.ge.s1         p14, p0 = fNormX, fDenoBound
 758       nop.i              0
 759 }
 760 ;;
 761 { .mfi
 762       nop.m              0
 763 (p14) fcmp.eq.s0         p6, p0 = f8, f0      // Set D flag if x unnormal
 764       nop.i              0
 765 }
 766 { .mfb
 767       nop.m              0
 768       // normalize unnormal input
 769 (p14) fnorm.s1           f8 = f8
 770       // return to the main path
 771 (p14) br.cond.sptk       asin_unnormal_back
 772 }
 773 ;;
 774 // if we still here it means that input is "true" denormal
 775 { .mfb
 776       nop.m              0
 777       // final result if x is denormal
 778       fma.d.s0           f8 = f8, fXSqr, f8
 779       // exit here if x is denormal
 780       br.ret.sptk        b0
 781 }
 782 ;;
 783
 784 // here if |x| > 1.0
 785 // error handler should be called
 786 .align 32
 787 asin_abs_gt_1:
 788 { .mfi
 789       alloc              r32 = ar.pfs, 0, 3, 4, 0 // get some registers
 790       fmerge.s           FR_X = f8,f8
 791       nop.i              0
 792 }
 793 { .mfb
 794       mov                GR_Parameter_TAG = 61 // error code
 795       frcpa.s0           FR_RESULT, p0 = f0,f0
 796       // call error handler routine
 797       br.cond.sptk       __libm_error_region
 798 }
 799 ;;
 800 GLOBAL_LIBM_END(asin)
 801
 802
 803
 804 LOCAL_LIBM_ENTRY(__libm_error_region)
 805 .prologue
 806 { .mfi
 807         add   GR_Parameter_Y=-32,sp             // Parameter 2 value
 808         nop.f 0
 809 .save   ar.pfs,GR_SAVE_PFS
 810         mov  GR_SAVE_PFS=ar.pfs                 // Save ar.pfs
 811 }
 812 { .mfi
 813 .fframe 64
 814         add sp=-64,sp                           // Create new stack
 815         nop.f 0
 816         mov GR_SAVE_GP=gp                       // Save gp
 817 };;
 818 { .mmi
 819         stfd [GR_Parameter_Y] = FR_Y,16         // STORE Parameter 2 on stack
 820         add GR_Parameter_X = 16,sp              // Parameter 1 address
 821 .save   b0, GR_SAVE_B0
 822         mov GR_SAVE_B0=b0                       // Save b0
 823 };;
 824 .body
 825 { .mib
 826         stfd [GR_Parameter_X] = FR_X                  // STORE Parameter 1 on stack
 827         add   GR_Parameter_RESULT = 0,GR_Parameter_Y  // Parameter 3 address
 828         nop.b 0
 829 }
 830 { .mib
 831         stfd [GR_Parameter_Y] = FR_RESULT             // STORE Parameter 3 on stack
 832         add   GR_Parameter_Y = -16,GR_Parameter_Y
 833         br.call.sptk b0=__libm_error_support#         // Call error handling function
 834 };;
 835 { .mmi
 836         add   GR_Parameter_RESULT = 48,sp
 837         nop.m 0
 838         nop.i 0
 839 };;
 840 { .mmi
 841         ldfd  f8 = [GR_Parameter_RESULT]       // Get return result off stack
 842 .restore sp
 843         add   sp = 64,sp                       // Restore stack pointer
 844         mov   b0 = GR_SAVE_B0                  // Restore return address
 845 };;
 846 { .mib
 847         mov   gp = GR_SAVE_GP                  // Restore gp
 848         mov   ar.pfs = GR_SAVE_PFS             // Restore ar.pfs
 849         br.ret.sptk     b0                     // Return
 850 };;
 851
 852 LOCAL_LIBM_END(__libm_error_region)
 853 .type   __libm_error_support#,@function
 854 .global __libm_error_support#