sysdeps/ia64/fpu/e_cosh.S

   1 .file "cosh.s"
   2
   3
   4 // Copyright (c) 2000 - 2005, 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 // 04/04/00 Unwind support added
  44 // 08/15/00 Bundle added after call to __libm_error_support to properly
  45 //          set [the previously overwritten] GR_Parameter_RESULT.
  46 // 05/07/01 Reworked to improve speed of all paths
  47 // 05/20/02 Cleaned up namespace and sf0 syntax
  48 // 11/15/02 Improved speed with new algorithm
  49 // 03/31/05 Reformatted delimiters between data tables
  50
  51 // API
  52 //==============================================================
  53 // double cosh(double)
  54
  55 // Overview of operation
  56 //==============================================================
  57 // Case 1:  0 < |x| < 0.25
  58 //  Evaluate cosh(x) by a 12th order polynomial
  59 //  Care is take for the order of multiplication; and A2 is not exactly 1/4!,
  60 //  A3 is not exactly 1/6!, etc.
  61 //  cosh(x) = 1 + (A1*x^2 + A2*x^4 + A3*x^6 + A4*x^8 + A5*x^10 + A6*x^12)
  62 //
  63 // Case 2:  0.25 < |x| < 710.47586
  64 //  Algorithm is based on the identity cosh(x) = ( exp(x) + exp(-x) ) / 2.
  65 //  The algorithm for exp is described as below.  There are a number of
  66 //  economies from evaluating both exp(x) and exp(-x).  Although we
  67 //  are evaluating both quantities, only where the quantities diverge do we
  68 //  duplicate the computations.  The basic algorithm for exp(x) is described
  69 //  below.
  70 //
  71 // Take the input x. w is "how many log2/128 in x?"
  72 //  w = x * 128/log2
  73 //  n = int(w)
  74 //  x = n log2/128 + r + delta
  75
  76 //  n = 128M + index_1 + 2^4 index_2
  77 //  x = M log2 + (log2/128) index_1 + (log2/8) index_2 + r + delta
  78
  79 //  exp(x) = 2^M  2^(index_1/128)  2^(index_2/8) exp(r) exp(delta)
  80 //       Construct 2^M
  81 //       Get 2^(index_1/128) from table_1;
  82 //       Get 2^(index_2/8)   from table_2;
  83 //       Calculate exp(r) by 5th order polynomial
  84 //          r = x - n (log2/128)_high
  85 //          delta = - n (log2/128)_low
  86 //       Calculate exp(delta) as 1 + delta
  87
  88
  89 // Special values
  90 //==============================================================
  91 // cosh(+0)    = 1.0
  92 // cosh(-0)    = 1.0
  93
  94 // cosh(+qnan) = +qnan
  95 // cosh(-qnan) = -qnan
  96 // cosh(+snan) = +qnan
  97 // cosh(-snan) = -qnan
  98
  99 // cosh(-inf)  = +inf
 100 // cosh(+inf)  = +inf
 101
 102 // Overflow and Underflow
 103 //=======================
 104 // cosh(x) = largest double normal when
 105 //     x = 710.47586 = 0x408633ce8fb9f87d
 106 //
 107 // There is no underflow.
 108
 109 // Registers used
 110 //==============================================================
 111 // Floating Point registers used:
 112 // f8, input, output
 113 // f6 -> f15,  f32 -> f61
 114
 115 // General registers used:
 116 // r14 -> r40
 117
 118 // Predicate registers used:
 119 // p6 -> p15
 120
 121 // Assembly macros
 122 //==============================================================
 123
 124 rRshf                 = r14
 125 rN_neg                = r14
 126 rAD_TB1               = r15
 127 rAD_TB2               = r16
 128 rAD_P                 = r17
 129 rN                    = r18
 130 rIndex_1              = r19
 131 rIndex_2_16           = r20
 132 rM                    = r21
 133 rBiased_M             = r21
 134 rSig_inv_ln2          = r22
 135 rIndex_1_neg          = r22
 136 rExp_bias             = r23
 137 rExp_bias_minus_1     = r23
 138 rExp_mask             = r24
 139 rTmp                  = r24
 140 rGt_ln                = r24
 141 rIndex_2_16_neg       = r24
 142 rM_neg                = r25
 143 rBiased_M_neg         = r25
 144 rRshf_2to56           = r26
 145 rAD_T1_neg            = r26
 146 rExp_2tom56           = r28
 147 rAD_T2_neg            = r28
 148 rAD_T1                = r29
 149 rAD_T2                = r30
 150 rSignexp_x            = r31
 151 rExp_x                = r31
 152
 153 GR_SAVE_B0            = r33
 154 GR_SAVE_PFS           = r34
 155 GR_SAVE_GP            = r35
 156 GR_SAVE_SP            = r36
 157
 158 GR_Parameter_X        = r37
 159 GR_Parameter_Y        = r38
 160 GR_Parameter_RESULT   = r39
 161 GR_Parameter_TAG      = r40
 162
 163
 164 FR_X                  = f10
 165 FR_Y                  = f1
 166 FR_RESULT             = f8
 167
 168 fRSHF_2TO56           = f6
 169 fINV_LN2_2TO63        = f7
 170 fW_2TO56_RSH          = f9
 171 f2TOM56               = f11
 172 fP5                   = f12
 173 fP4                   = f13
 174 fP3                   = f14
 175 fP2                   = f15
 176
 177 fLn2_by_128_hi        = f33
 178 fLn2_by_128_lo        = f34
 179
 180 fRSHF                 = f35
 181 fNfloat               = f36
 182 fNormX                = f37
 183 fR                    = f38
 184 fF                    = f39
 185
 186 fRsq                  = f40
 187 f2M                   = f41
 188 fS1                   = f42
 189 fT1                   = f42
 190 fS2                   = f43
 191 fT2                   = f43
 192 fS                    = f43
 193 fWre_urm_f8           = f44
 194 fAbsX                 = f44
 195
 196 fMIN_DBL_OFLOW_ARG    = f45
 197 fMAX_DBL_NORM_ARG     = f46
 198 fXsq                  = f47
 199 fX4                   = f48
 200 fGt_pln               = f49
 201 fTmp                  = f49
 202
 203 fP54                  = f50
 204 fP5432                = f50
 205 fP32                  = f51
 206 fP                    = f52
 207 fP54_neg              = f53
 208 fP5432_neg            = f53
 209 fP32_neg              = f54
 210 fP_neg                = f55
 211 fF_neg                = f56
 212
 213 f2M_neg               = f57
 214 fS1_neg               = f58
 215 fT1_neg               = f58
 216 fS2_neg               = f59
 217 fT2_neg               = f59
 218 fS_neg                = f59
 219 fExp                  = f60
 220 fExp_neg              = f61
 221
 222 fA6                   = f50
 223 fA65                  = f50
 224 fA6543                = f50
 225 fA654321              = f50
 226 fA5                   = f51
 227 fA4                   = f52
 228 fA43                  = f52
 229 fA3                   = f53
 230 fA2                   = f54
 231 fA21                  = f54
 232 fA1                   = f55
 233
 234 // Data tables
 235 //==============================================================
 236
 237 RODATA
 238 .align 16
 239
 240 // ************* DO NOT CHANGE ORDER OF THESE TABLES ********************
 241
 242 // double-extended 1/ln(2)
 243 // 3fff b8aa 3b29 5c17 f0bb be87fed0691d3e88
 244 // 3fff b8aa 3b29 5c17 f0bc
 245 // For speed the significand will be loaded directly with a movl and setf.sig
 246 //   and the exponent will be bias+63 instead of bias+0.  Thus subsequent
 247 //   computations need to scale appropriately.
 248 // The constant 128/ln(2) is needed for the computation of w.  This is also
 249 //   obtained by scaling the computations.
 250 //
 251 // Two shifting constants are loaded directly with movl and setf.d.
 252 //   1. fRSHF_2TO56 = 1.1000..00 * 2^(63-7)
 253 //        This constant is added to x*1/ln2 to shift the integer part of
 254 //        x*128/ln2 into the rightmost bits of the significand.
 255 //        The result of this fma is fW_2TO56_RSH.
 256 //   2. fRSHF       = 1.1000..00 * 2^(63)
 257 //        This constant is subtracted from fW_2TO56_RSH * 2^(-56) to give
 258 //        the integer part of w, n, as a floating-point number.
 259 //        The result of this fms is fNfloat.
 260
 261
 262 LOCAL_OBJECT_START(exp_table_1)
 263 data8 0x408633ce8fb9f87e // smallest dbl overflow arg
 264 data8 0x408633ce8fb9f87d // largest dbl arg to give normal dbl result
 265 data8 0xb17217f7d1cf79ab , 0x00003ff7 // ln2/128 hi
 266 data8 0xc9e3b39803f2f6af , 0x00003fb7 // ln2/128 lo
 267 //
 268 // Table 1 is 2^(index_1/128) where
 269 // index_1 goes from 0 to 15
 270 //
 271 data8 0x8000000000000000 , 0x00003FFF
 272 data8 0x80B1ED4FD999AB6C , 0x00003FFF
 273 data8 0x8164D1F3BC030773 , 0x00003FFF
 274 data8 0x8218AF4373FC25EC , 0x00003FFF
 275 data8 0x82CD8698AC2BA1D7 , 0x00003FFF
 276 data8 0x8383594EEFB6EE37 , 0x00003FFF
 277 data8 0x843A28C3ACDE4046 , 0x00003FFF
 278 data8 0x84F1F656379C1A29 , 0x00003FFF
 279 data8 0x85AAC367CC487B15 , 0x00003FFF
 280 data8 0x8664915B923FBA04 , 0x00003FFF
 281 data8 0x871F61969E8D1010 , 0x00003FFF
 282 data8 0x87DB357FF698D792 , 0x00003FFF
 283 data8 0x88980E8092DA8527 , 0x00003FFF
 284 data8 0x8955EE03618E5FDD , 0x00003FFF
 285 data8 0x8A14D575496EFD9A , 0x00003FFF
 286 data8 0x8AD4C6452C728924 , 0x00003FFF
 287 LOCAL_OBJECT_END(exp_table_1)
 288
 289 // Table 2 is 2^(index_1/8) where
 290 // index_2 goes from 0 to 7
 291 LOCAL_OBJECT_START(exp_table_2)
 292 data8 0x8000000000000000 , 0x00003FFF
 293 data8 0x8B95C1E3EA8BD6E7 , 0x00003FFF
 294 data8 0x9837F0518DB8A96F , 0x00003FFF
 295 data8 0xA5FED6A9B15138EA , 0x00003FFF
 296 data8 0xB504F333F9DE6484 , 0x00003FFF
 297 data8 0xC5672A115506DADD , 0x00003FFF
 298 data8 0xD744FCCAD69D6AF4 , 0x00003FFF
 299 data8 0xEAC0C6E7DD24392F , 0x00003FFF
 300 LOCAL_OBJECT_END(exp_table_2)
 301
 302 LOCAL_OBJECT_START(exp_p_table)
 303 data8 0x3f8111116da21757 //P5
 304 data8 0x3fa55555d787761c //P4
 305 data8 0x3fc5555555555414 //P3
 306 data8 0x3fdffffffffffd6a //P2
 307 LOCAL_OBJECT_END(exp_p_table)
 308
 309 LOCAL_OBJECT_START(cosh_p_table)
 310 data8 0x8FA02AC65BCBD5BC, 0x00003FE2  // A6
 311 data8 0xD00D00D1021D7370, 0x00003FEF  // A4
 312 data8 0xAAAAAAAAAAAAAB80, 0x00003FFA  // A2
 313 data8 0x93F27740C0C2F1CC, 0x00003FE9  // A5
 314 data8 0xB60B60B60B4FE884, 0x00003FF5  // A3
 315 data8 0x8000000000000000, 0x00003FFE  // A1
 316 LOCAL_OBJECT_END(cosh_p_table)
 317
 318
 319 .section .text
 320 GLOBAL_IEEE754_ENTRY(cosh)
 321
 322 { .mlx
 323       getf.exp        rSignexp_x = f8  // Must recompute if x unorm
 324       movl            rSig_inv_ln2 = 0xb8aa3b295c17f0bc  // significand of 1/ln2
 325 }
 326 { .mlx
 327       addl            rAD_TB1    = @ltoff(exp_table_1), gp
 328       movl            rRshf_2to56 = 0x4768000000000000   // 1.10000 2^(63+56)
 329 }
 330 ;;
 331
 332 { .mfi
 333       ld8             rAD_TB1    = [rAD_TB1]
 334       fclass.m        p6,p0 = f8,0x0b  // Test for x=unorm
 335       mov             rExp_mask = 0x1ffff
 336 }
 337 { .mfi
 338       mov             rExp_bias = 0xffff
 339       fnorm.s1        fNormX   = f8
 340       mov             rExp_2tom56 = 0xffff-56
 341 }
 342 ;;
 343
 344 // Form two constants we need
 345 //  1/ln2 * 2^63  to compute  w = x * 1/ln2 * 128
 346 //  1.1000..000 * 2^(63+63-7) to right shift int(w) into the significand
 347
 348 { .mfi
 349       setf.sig        fINV_LN2_2TO63 = rSig_inv_ln2 // form 1/ln2 * 2^63
 350       fclass.m        p8,p0 = f8,0x07  // Test for x=0
 351       nop.i 999
 352 }
 353 { .mlx
 354       setf.d          fRSHF_2TO56 = rRshf_2to56 // Form const 1.100 * 2^(63+56)
 355       movl            rRshf = 0x43e8000000000000 // 1.10000 2^63 for right shift
 356 }
 357 ;;
 358
 359 { .mfi
 360       ldfpd           fMIN_DBL_OFLOW_ARG, fMAX_DBL_NORM_ARG = [rAD_TB1],16
 361       fclass.m        p10,p0 = f8,0x1e3  // Test for x=inf, nan, NaT
 362       nop.i           0
 363 }
 364 { .mfb
 365       setf.exp        f2TOM56 = rExp_2tom56 // form 2^-56 for scaling Nfloat
 366       nop.f           0
 367 (p6)  br.cond.spnt    COSH_UNORM            // Branch if x=unorm
 368 }
 369 ;;
 370
 371 COSH_COMMON:
 372 { .mfi
 373       ldfe            fLn2_by_128_hi  = [rAD_TB1],16
 374       nop.f           0
 375       nop.i           0
 376 }
 377 { .mfb
 378       setf.d          fRSHF = rRshf // Form right shift const 1.100 * 2^63
 379 (p8)  fma.d.s0        f8 = f1,f1,f0           // quick exit for x=0
 380 (p8)  br.ret.spnt     b0
 381 }
 382 ;;
 383
 384 { .mfi
 385       ldfe            fLn2_by_128_lo  = [rAD_TB1],16
 386       nop.f           0
 387       nop.i           0
 388 }
 389 { .mfb
 390       and             rExp_x = rExp_mask, rSignexp_x // Biased exponent of x
 391 (p10) fma.d.s0        f8 = f8,f8,f0  // Result if x=inf, nan, NaT
 392 (p10) br.ret.spnt     b0               // quick exit for x=inf, nan, NaT
 393 }
 394 ;;
 395
 396 // After that last load rAD_TB1 points to the beginning of table 1
 397 { .mfi
 398       nop.m           0
 399       fcmp.eq.s0      p6,p0 = f8, f0       // Dummy to set D
 400       sub             rExp_x = rExp_x, rExp_bias // True exponent of x
 401 }
 402 ;;
 403
 404 { .mfi
 405       nop.m           0
 406       fmerge.s        fAbsX = f0, fNormX   // Form |x|
 407       nop.i           0
 408 }
 409 { .mfb
 410       cmp.gt          p7, p0 = -2, rExp_x      // Test |x| < 2^(-2)
 411       fma.s1          fXsq = fNormX, fNormX, f0  // x*x for small path
 412 (p7)  br.cond.spnt    COSH_SMALL               // Branch if 0 < |x| < 2^-2
 413 }
 414 ;;
 415
 416 // W = X * Inv_log2_by_128
 417 // By adding 1.10...0*2^63 we shift and get round_int(W) in significand.
 418 // We actually add 1.10...0*2^56 to X * Inv_log2 to do the same thing.
 419
 420 { .mfi
 421       add             rAD_P = 0x180, rAD_TB1
 422       fma.s1          fW_2TO56_RSH  = fNormX, fINV_LN2_2TO63, fRSHF_2TO56
 423       add             rAD_TB2 = 0x100, rAD_TB1
 424 }
 425 ;;
 426
 427 // Divide arguments into the following categories:
 428 //  Certain Safe                - 0.25 <= |x| <= MAX_DBL_NORM_ARG
 429 //  Possible Overflow       p14 - MAX_DBL_NORM_ARG < |x| < MIN_DBL_OFLOW_ARG
 430 //  Certain Overflow        p15 - MIN_DBL_OFLOW_ARG <= |x| < +inf
 431 //
 432 // If the input is really a double arg, then there will never be
 433 // "Possible Overflow" arguments.
 434 //
 435
 436 { .mfi
 437       ldfpd           fP5, fP4  = [rAD_P] ,16
 438       fcmp.ge.s1      p15,p14 = fAbsX,fMIN_DBL_OFLOW_ARG
 439       nop.i           0
 440 }
 441 ;;
 442
 443 // Nfloat = round_int(W)
 444 // The signficand of fW_2TO56_RSH contains the rounded integer part of W,
 445 // as a twos complement number in the lower bits (that is, it may be negative).
 446 // That twos complement number (called N) is put into rN.
 447
 448 // Since fW_2TO56_RSH is scaled by 2^56, it must be multiplied by 2^-56
 449 // before the shift constant 1.10000 * 2^63 is subtracted to yield fNfloat.
 450 // Thus, fNfloat contains the floating point version of N
 451
 452 { .mfi
 453       ldfpd           fP3, fP2  = [rAD_P]
 454 (p14) fcmp.gt.unc.s1  p14,p0 = fAbsX,fMAX_DBL_NORM_ARG
 455       nop.i           0
 456 }
 457 { .mfb
 458       nop.m           0
 459       fms.s1          fNfloat = fW_2TO56_RSH, f2TOM56, fRSHF
 460 (p15) br.cond.spnt    COSH_CERTAIN_OVERFLOW
 461 }
 462 ;;
 463
 464 { .mfi
 465       getf.sig        rN        = fW_2TO56_RSH
 466       nop.f           0
 467       mov             rExp_bias_minus_1 = 0xfffe
 468 }
 469 ;;
 470
 471 // rIndex_1 has index_1
 472 // rIndex_2_16 has index_2 * 16
 473 // rBiased_M has M
 474
 475 // rM has true M
 476 // r = x - Nfloat * ln2_by_128_hi
 477 // f = 1 - Nfloat * ln2_by_128_lo
 478 { .mfi
 479       and             rIndex_1 = 0x0f, rN
 480       fnma.s1         fR   = fNfloat, fLn2_by_128_hi, fNormX
 481       shr             rM = rN,  0x7
 482 }
 483 { .mfi
 484       and             rIndex_2_16 = 0x70, rN
 485       fnma.s1         fF   = fNfloat, fLn2_by_128_lo, f1
 486       sub             rN_neg = r0, rN
 487 }
 488 ;;
 489
 490 { .mmi
 491       and             rIndex_1_neg = 0x0f, rN_neg
 492       add             rBiased_M = rExp_bias_minus_1, rM
 493       shr             rM_neg = rN_neg,  0x7
 494 }
 495 { .mmi
 496       and             rIndex_2_16_neg = 0x70, rN_neg
 497       add             rAD_T2 = rAD_TB2, rIndex_2_16
 498       shladd          rAD_T1 = rIndex_1, 4, rAD_TB1
 499 }
 500 ;;
 501
 502 // rAD_T1 has address of T1
 503 // rAD_T2 has address if T2
 504
 505 { .mmi
 506       setf.exp        f2M = rBiased_M
 507       ldfe            fT2  = [rAD_T2]
 508       nop.i           0
 509 }
 510 { .mmi
 511       add             rBiased_M_neg = rExp_bias_minus_1, rM_neg
 512       add             rAD_T2_neg = rAD_TB2, rIndex_2_16_neg
 513       shladd          rAD_T1_neg = rIndex_1_neg, 4, rAD_TB1
 514 }
 515 ;;
 516
 517 // Create Scale = 2^M
 518 // Load T1 and T2
 519 { .mmi
 520       ldfe            fT1  = [rAD_T1]
 521       nop.m           0
 522       nop.i           0
 523 }
 524 { .mmf
 525       setf.exp        f2M_neg = rBiased_M_neg
 526       ldfe            fT2_neg  = [rAD_T2_neg]
 527       fma.s1          fF_neg   = fNfloat, fLn2_by_128_lo, f1
 528 }
 529 ;;
 530
 531 { .mfi
 532       nop.m           0
 533       fma.s1          fRsq = fR, fR, f0
 534       nop.i           0
 535 }
 536 { .mfi
 537       ldfe            fT1_neg  = [rAD_T1_neg]
 538       fma.s1          fP54 = fR, fP5, fP4
 539       nop.i           0
 540 }
 541 ;;
 542
 543 { .mfi
 544       nop.m           0
 545       fma.s1          fP32 = fR, fP3, fP2
 546       nop.i           0
 547 }
 548 { .mfi
 549       nop.m           0
 550       fnma.s1         fP54_neg = fR, fP5, fP4
 551       nop.i           0
 552 }
 553 ;;
 554
 555 { .mfi
 556       nop.m           0
 557       fnma.s1         fP32_neg = fR, fP3, fP2
 558       nop.i           0
 559 }
 560 ;;
 561
 562 { .mfi
 563       nop.m           0
 564       fma.s1          fP5432  = fRsq, fP54, fP32
 565       nop.i           0
 566 }
 567 { .mfi
 568       nop.m           0
 569       fma.s1          fS2  = fF,fT2,f0
 570       nop.i           0
 571 }
 572 ;;
 573
 574 { .mfi
 575       nop.m           0
 576       fma.s1          fS1  = f2M,fT1,f0
 577       nop.i           0
 578 }
 579 { .mfi
 580       nop.m           0
 581       fma.s1          fP5432_neg  = fRsq, fP54_neg, fP32_neg
 582       nop.i           0
 583 }
 584 ;;
 585
 586 { .mfi
 587       nop.m           0
 588       fma.s1          fS1_neg  = f2M_neg,fT1_neg,f0
 589       nop.i           0
 590 }
 591 { .mfi
 592       nop.m           0
 593       fma.s1          fS2_neg  = fF_neg,fT2_neg,f0
 594       nop.i           0
 595 }
 596 ;;
 597
 598 { .mfi
 599       nop.m           0
 600       fma.s1          fP     = fRsq, fP5432, fR
 601       nop.i           0
 602 }
 603 { .mfi
 604       nop.m           0
 605       fma.s1          fS   = fS1,fS2,f0
 606       nop.i           0
 607 }
 608 ;;
 609
 610 { .mfi
 611       nop.m           0
 612       fms.s1          fP_neg     = fRsq, fP5432_neg, fR
 613       nop.i           0
 614 }
 615 { .mfi
 616       nop.m           0
 617       fma.s1          fS_neg   = fS1_neg,fS2_neg,f0
 618       nop.i           0
 619 }
 620 ;;
 621
 622 { .mfb
 623       nop.m           0
 624       fmpy.s0         fTmp = fLn2_by_128_lo, fLn2_by_128_lo // Force inexact
 625 (p14) br.cond.spnt    COSH_POSSIBLE_OVERFLOW
 626 }
 627 ;;
 628
 629 { .mfi
 630       nop.m           0
 631       fma.s1          fExp = fS, fP, fS
 632       nop.i           0
 633 }
 634 { .mfi
 635       nop.m           0
 636       fma.s1          fExp_neg = fS_neg, fP_neg, fS_neg
 637       nop.i           0
 638 }
 639 ;;
 640
 641 { .mfb
 642       nop.m           0
 643       fma.d.s0        f8 = fExp, f1, fExp_neg
 644       br.ret.sptk     b0                  // Normal path exit
 645 }
 646 ;;
 647
 648 // Here if 0 < |x| < 0.25
 649 COSH_SMALL:
 650 { .mmf
 651       add             rAD_T1 = 0x1a0, rAD_TB1
 652       add             rAD_T2 = 0x1d0, rAD_TB1
 653 }
 654 ;;
 655
 656 { .mmf
 657       ldfe            fA6 = [rAD_T1],16
 658       ldfe            fA5 = [rAD_T2],16
 659       nop.f           0
 660 }
 661 ;;
 662
 663 { .mmi
 664       ldfe            fA4 = [rAD_T1],16
 665       ldfe            fA3 = [rAD_T2],16
 666       nop.i           0
 667 }
 668 ;;
 669
 670 { .mmi
 671       ldfe            fA2 = [rAD_T1],16
 672       ldfe            fA1 = [rAD_T2],16
 673       nop.i           0
 674 }
 675 ;;
 676
 677 { .mfi
 678       nop.m           0
 679       fma.s1          fX4 = fXsq, fXsq, f0
 680       nop.i           0
 681 }
 682 ;;
 683
 684 { .mfi
 685       nop.m           0
 686       fma.s1          fA65 = fXsq, fA6, fA5
 687       nop.i           0
 688 }
 689 { .mfi
 690       nop.m           0
 691       fma.s1          fA43 = fXsq, fA4, fA3
 692       nop.i           0
 693 }
 694 ;;
 695
 696 { .mfi
 697       nop.m           0
 698       fma.s1          fA21 = fXsq, fA2, fA1
 699       nop.i           0
 700 }
 701 ;;
 702
 703 { .mfi
 704       nop.m           0
 705       fma.s1          fA6543 = fX4, fA65, fA43
 706       nop.i           0
 707 }
 708 ;;
 709
 710 { .mfi
 711       nop.m           0
 712       fma.s1          fA654321 = fX4, fA6543, fA21
 713       nop.i           0
 714 }
 715 ;;
 716
 717 // Dummy multiply to generate inexact
 718 { .mfi
 719       nop.m           0
 720       fmpy.s0         fTmp = fA6, fA6
 721       nop.i           0
 722 }
 723 { .mfb
 724       nop.m           0
 725       fma.d.s0        f8 = fA654321, fXsq, f1
 726       br.ret.sptk     b0                // Exit if 0 < |x| < 0.25
 727 }
 728 ;;
 729
 730
 731 COSH_POSSIBLE_OVERFLOW:
 732
 733 // Here if fMAX_DBL_NORM_ARG < |x| < fMIN_DBL_OFLOW_ARG
 734 // This cannot happen if input is a double, only if input higher precision.
 735 // Overflow is a possibility, not a certainty.
 736
 737 // Recompute result using status field 2 with user's rounding mode,
 738 // and wre set.  If result is larger than largest double, then we have
 739 // overflow
 740
 741 { .mfi
 742       mov             rGt_ln  = 0x103ff // Exponent for largest dbl + 1 ulp
 743       fsetc.s2        0x7F,0x42         // Get user's round mode, set wre
 744       nop.i           0
 745 }
 746 ;;
 747
 748 { .mfi
 749       setf.exp        fGt_pln = rGt_ln  // Create largest double + 1 ulp
 750       fma.d.s2        fWre_urm_f8 = fS, fP, fS    // Result with wre set
 751       nop.i           0
 752 }
 753 ;;
 754
 755 { .mfi
 756       nop.m           0
 757       fsetc.s2        0x7F,0x40                   // Turn off wre in sf2
 758       nop.i           0
 759 }
 760 ;;
 761
 762 { .mfi
 763       nop.m           0
 764       fcmp.ge.s1      p6, p0 =  fWre_urm_f8, fGt_pln // Test for overflow
 765       nop.i           0
 766 }
 767 ;;
 768
 769 { .mfb
 770       nop.m           0
 771       nop.f           0
 772 (p6)  br.cond.spnt    COSH_CERTAIN_OVERFLOW // Branch if overflow
 773 }
 774 ;;
 775
 776 { .mfb
 777       nop.m           0
 778       fma.d.s0        f8 = fS, fP, fS
 779       br.ret.sptk     b0                     // Exit if really no overflow
 780 }
 781 ;;
 782
 783 COSH_CERTAIN_OVERFLOW:
 784 { .mmi
 785       sub             rTmp = rExp_mask, r0, 1
 786 ;;
 787       setf.exp        fTmp = rTmp
 788       nop.i           0
 789 }
 790 ;;
 791
 792 { .mfi
 793       alloc           r32=ar.pfs,1,4,4,0
 794       fmerge.s        FR_X = f8,f8
 795       nop.i           0
 796 }
 797 { .mfb
 798       mov             GR_Parameter_TAG = 64
 799       fma.d.s0        FR_RESULT = fTmp, fTmp, f0    // Set I,O and +INF result
 800       br.cond.sptk    __libm_error_region
 801 }
 802 ;;
 803
 804 // Here if x unorm
 805 COSH_UNORM:
 806 { .mfb
 807       getf.exp        rSignexp_x = fNormX    // Must recompute if x unorm
 808       fcmp.eq.s0      p6, p0 = f8, f0        // Set D flag
 809       br.cond.sptk    COSH_COMMON
 810 }
 811 ;;
 812
 813 GLOBAL_IEEE754_END(cosh)
 814
 815
 816 LOCAL_LIBM_ENTRY(__libm_error_region)
 817 .prologue
 818 { .mfi
 819         add   GR_Parameter_Y=-32,sp             // Parameter 2 value
 820         nop.f 0
 821 .save   ar.pfs,GR_SAVE_PFS
 822         mov  GR_SAVE_PFS=ar.pfs                 // Save ar.pfs
 823 }
 824 { .mfi
 825 .fframe 64
 826         add sp=-64,sp                           // Create new stack
 827         nop.f 0
 828         mov GR_SAVE_GP=gp                       // Save gp
 829 };;
 830 { .mmi
 831         stfd [GR_Parameter_Y] = FR_Y,16         // STORE Parameter 2 on stack
 832         add GR_Parameter_X = 16,sp              // Parameter 1 address
 833 .save   b0, GR_SAVE_B0
 834         mov GR_SAVE_B0=b0                       // Save b0
 835 };;
 836 .body
 837 { .mib
 838         stfd [GR_Parameter_X] = FR_X            // STORE Parameter 1 on stack
 839         add   GR_Parameter_RESULT = 0,GR_Parameter_Y  // Parameter 3 address
 840         nop.b 0
 841 }
 842 { .mib
 843         stfd [GR_Parameter_Y] = FR_RESULT       // STORE Parameter 3 on stack
 844         add   GR_Parameter_Y = -16,GR_Parameter_Y
 845         br.call.sptk b0=__libm_error_support#   // Call error handling function
 846 };;
 847 { .mmi
 848         add   GR_Parameter_RESULT = 48,sp
 849         nop.m 0
 850         nop.i 0
 851 };;
 852 { .mmi
 853         ldfd  f8 = [GR_Parameter_RESULT]       // Get return result off stack
 854 .restore sp
 855         add   sp = 64,sp                       // Restore stack pointer
 856         mov   b0 = GR_SAVE_B0                  // Restore return address
 857 };;
 858 { .mib
 859         mov   gp = GR_SAVE_GP                  // Restore gp
 860         mov   ar.pfs = GR_SAVE_PFS             // Restore ar.pfs
 861         br.ret.sptk     b0                     // Return
 862 };;
 863
 864 LOCAL_LIBM_END(__libm_error_region)
 865 .type   __libm_error_support#,@function
 866 .global __libm_error_support#