sysdeps/ia64/fpu/e_coshl.S

   1 .file "coshl.s"
   2
   3
   4 // Copyright (c) 2000 - 2002, 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 // 01/23/01 Set inexact flag for large args.
  47 // 05/07/01 Reworked to improve speed of all paths
  48 // 05/20/02 Cleaned up namespace and sf0 syntax
  49 // 12/06/02 Improved performance
  50 //
  51 // API
  52 //==============================================================
  53 // long double = coshl(long double)
  54 // input  floating point f8
  55 // output floating point f8
  56 //
  57 // Registers used
  58 //==============================================================
  59 // general registers:
  60 // r14 -> r40
  61 // predicate registers used:
  62 // p6 -> p11
  63 // floating-point registers used:
  64 // f9 -> f15; f32 -> f90;
  65 // f8 has input, then output
  66 //
  67 // Overview of operation
  68 //==============================================================
  69 // There are seven paths
  70 // 1. 0 < |x| < 0.25          COSH_BY_POLY
  71 // 2. 0.25 <=|x| < 32         COSH_BY_TBL
  72 // 3. 32 <= |x| < 11357.21655 COSH_BY_EXP (merged path with COSH_BY_TBL)
  73 // 4. |x| >= 11357.21655      COSH_HUGE
  74 // 5. x=0                     Done with early exit
  75 // 6. x=inf,nan               Done with early exit
  76 // 7. x=denormal              COSH_DENORM
  77 //
  78 // For double extended we get overflow for x >= 400c b174 ddc0 31ae c0ea
  79 //                                           >= 11357.21655
  80 //
  81 //
  82 // 1. COSH_BY_POLY   0 < |x| < 0.25
  83 // ===============
  84 // Evaluate cosh(x) by a 12th order polynomial
  85 // Care is take for the order of multiplication; and P2 is not exactly 1/4!,
  86 // P3 is not exactly 1/6!, etc.
  87 // cosh(x) = 1 + (P1*x^2 + P2*x^4 + P3*x^6 + P4*x^8 + P5*x^10 + P6*x^12)
  88 //
  89 // 2. COSH_BY_TBL   0.25 <= |x| < 32.0
  90 // =============
  91 // cosh(x) = cosh(B+R)
  92 //         = cosh(B)cosh(R) + sinh(B)sinh(R)
  93 //
  94 // ax = |x| = M*log2/64 + R
  95 // B = M*log2/64
  96 // M = 64*N + j
  97 //   We will calculate M and get N as (M-j)/64
  98 //   The division is a shift.
  99 // exp(B)  = exp(N*log2 + j*log2/64)
 100 //         = 2^N * 2^(j*log2/64)
 101 // cosh(B) = 1/2(e^B + e^-B)
 102 //         = 1/2(2^N * 2^(j*log2/64) + 2^-N * 2^(-j*log2/64))
 103 // cosh(B) = (2^(N-1) * 2^(j*log2/64) + 2^(-N-1) * 2^(-j*log2/64))
 104 // sinh(B) = (2^(N-1) * 2^(j*log2/64) - 2^(-N-1) * 2^(-j*log2/64))
 105 // 2^(j*log2/64) is stored as Tjhi + Tjlo , j= -32,....,32
 106 // Tjhi is double-extended (80-bit) and Tjlo is single(32-bit)
 107 //
 108 // R = ax - M*log2/64
 109 // R = ax - M*log2_by_64_hi - M*log2_by_64_lo
 110 // exp(R) = 1 + R +R^2(1/2! + R(1/3! + R(1/4! + ... + R(1/n!)...)
 111 //        = 1 + p_odd + p_even
 112 //        where the p_even uses the A coefficients and the p_even uses
 113 //        the B coefficients
 114 //
 115 // So sinh(R) = 1 + p_odd + p_even -(1 -p_odd -p_even)/2 = p_odd
 116 //    cosh(R) = 1 + p_even
 117 //    cosh(B) = C_hi + C_lo
 118 //    sinh(B) = S_hi
 119 // cosh(x) = cosh(B)cosh(R) + sinh(B)sinh(R)
 120 //
 121 // 3. COSH_BY_EXP   32.0 <= |x| < 11357.21655  ( 400c b174 ddc0 31ae c0ea )
 122 // ==============
 123 // Can approximate result by exp(x)/2 in this region.
 124 // Y_hi = Tjhi
 125 // Y_lo = Tjhi * (p_odd + p_even) + Tjlo
 126 // cosh(x) = Y_hi + Y_lo
 127 //
 128 // 4. COSH_HUGE     |x| >= 11357.21655  ( 400c b174 ddc0 31ae c0ea )
 129 // ============
 130 // Set error tag and call error support
 131 //
 132 //
 133 // Assembly macros
 134 //==============================================================
 135 r_ad5                 = r14
 136 r_rshf_2to57          = r15
 137 r_exp_denorm          = r15
 138 r_ad_mJ_lo            = r15
 139 r_ad_J_lo             = r16
 140 r_2Nm1                = r17
 141 r_2mNm1               = r18
 142 r_exp_x               = r18
 143 r_ad_J_hi             = r19
 144 r_ad2o                = r19
 145 r_ad_mJ_hi            = r20
 146 r_mj                  = r21
 147 r_ad2e                = r22
 148 r_ad3                 = r23
 149 r_ad1                 = r24
 150 r_Mmj                 = r24
 151 r_rshf                = r25
 152 r_M                   = r25
 153 r_N                   = r25
 154 r_jshf                = r26
 155 r_exp_2tom57          = r26
 156 r_j                   = r26
 157 r_exp_mask            = r27
 158 r_signexp_x           = r28
 159 r_signexp_0_5         = r28
 160 r_exp_0_25            = r29
 161 r_sig_inv_ln2         = r30
 162 r_exp_32              = r30
 163 r_exp_huge            = r30
 164 r_ad4                 = r31
 165
 166 GR_SAVE_PFS           = r34
 167 GR_SAVE_B0            = r35
 168 GR_SAVE_GP            = r36
 169
 170 GR_Parameter_X        = r37
 171 GR_Parameter_Y        = r38
 172 GR_Parameter_RESULT   = r39
 173 GR_Parameter_TAG      = r40
 174
 175
 176 f_ABS_X               = f9
 177 f_X2                  = f10
 178 f_X4                  = f11
 179 f_tmp                 = f14
 180 f_RSHF                = f15
 181
 182 f_Inv_log2by64        = f32
 183 f_log2by64_lo         = f33
 184 f_log2by64_hi         = f34
 185 f_A1                  = f35
 186
 187 f_A2                  = f36
 188 f_A3                  = f37
 189 f_Rcub                = f38
 190 f_M_temp              = f39
 191 f_R_temp              = f40
 192
 193 f_Rsq                 = f41
 194 f_R                   = f42
 195 f_M                   = f43
 196 f_B1                  = f44
 197 f_B2                  = f45
 198
 199 f_B3                  = f46
 200 f_peven_temp1         = f47
 201 f_peven_temp2         = f48
 202 f_peven               = f49
 203 f_podd_temp1          = f50
 204
 205 f_podd_temp2          = f51
 206 f_podd                = f52
 207 f_poly65              = f53
 208 f_poly6543            = f53
 209 f_poly6to1            = f53
 210 f_poly43              = f54
 211 f_poly21              = f55
 212
 213 f_X3                  = f56
 214 f_INV_LN2_2TO63       = f57
 215 f_RSHF_2TO57          = f58
 216 f_2TOM57              = f59
 217 f_smlst_oflow_input   = f60
 218
 219 f_pre_result          = f61
 220 f_huge                = f62
 221 f_spos                = f63
 222 f_sneg                = f64
 223 f_Tjhi                = f65
 224
 225 f_Tjlo                = f66
 226 f_Tmjhi               = f67
 227 f_Tmjlo               = f68
 228 f_S_hi                = f69
 229 f_SC_hi_temp          = f70
 230
 231 f_C_lo_temp1          = f71
 232 f_C_lo_temp2          = f72
 233 f_C_lo_temp3          = f73
 234 f_C_lo_temp4          = f73
 235 f_C_lo                = f74
 236 f_C_hi                = f75
 237
 238 f_Y_hi                = f77
 239 f_Y_lo_temp           = f78
 240 f_Y_lo                = f79
 241 f_NORM_X              = f80
 242
 243 f_P1                  = f81
 244 f_P2                  = f82
 245 f_P3                  = f83
 246 f_P4                  = f84
 247 f_P5                  = f85
 248
 249 f_P6                  = f86
 250 f_Tjhi_spos           = f87
 251 f_Tjlo_spos           = f88
 252 f_huge                = f89
 253 f_signed_hi_lo        = f90
 254
 255
 256 // Data tables
 257 //==============================================================
 258
 259 // DO NOT CHANGE ORDER OF THESE TABLES
 260 RODATA
 261
 262 .align 16
 263 LOCAL_OBJECT_START(cosh_arg_reduction)
 264 //   data8 0xB8AA3B295C17F0BC, 0x00004005  // 64/log2 -- signif loaded with setf
 265    data8 0xB17217F7D1000000, 0x00003FF8  // log2/64 high part
 266    data8 0xCF79ABC9E3B39804, 0x00003FD0  // log2/64 low part
 267    data8 0xb174ddc031aec0ea, 0x0000400c  // Smallest x to overflow (11357.21655)
 268 LOCAL_OBJECT_END(cosh_arg_reduction)
 269
 270 LOCAL_OBJECT_START(cosh_p_table)
 271    data8 0x8FA02AC65BCBD5BC, 0x00003FE2  // P6
 272    data8 0xD00D00D1021D7370, 0x00003FEF  // P4
 273    data8 0xAAAAAAAAAAAAAB80, 0x00003FFA  // P2
 274    data8 0x93F27740C0C2F1CC, 0x00003FE9  // P5
 275    data8 0xB60B60B60B4FE884, 0x00003FF5  // P3
 276    data8 0x8000000000000000, 0x00003FFE  // P1
 277 LOCAL_OBJECT_END(cosh_p_table)
 278
 279 LOCAL_OBJECT_START(cosh_ab_table)
 280    data8 0xAAAAAAAAAAAAAAAC, 0x00003FFC  // A1
 281    data8 0x88888888884ECDD5, 0x00003FF8  // A2
 282    data8 0xD00D0C6DCC26A86B, 0x00003FF2  // A3
 283    data8 0x8000000000000002, 0x00003FFE  // B1
 284    data8 0xAAAAAAAAAA402C77, 0x00003FFA  // B2
 285    data8 0xB60B6CC96BDB144D, 0x00003FF5  // B3
 286 LOCAL_OBJECT_END(cosh_ab_table)
 287
 288 LOCAL_OBJECT_START(cosh_j_hi_table)
 289    data8 0xB504F333F9DE6484, 0x00003FFE
 290    data8 0xB6FD91E328D17791, 0x00003FFE
 291    data8 0xB8FBAF4762FB9EE9, 0x00003FFE
 292    data8 0xBAFF5AB2133E45FB, 0x00003FFE
 293    data8 0xBD08A39F580C36BF, 0x00003FFE
 294    data8 0xBF1799B67A731083, 0x00003FFE
 295    data8 0xC12C4CCA66709456, 0x00003FFE
 296    data8 0xC346CCDA24976407, 0x00003FFE
 297    data8 0xC5672A115506DADD, 0x00003FFE
 298    data8 0xC78D74C8ABB9B15D, 0x00003FFE
 299    data8 0xC9B9BD866E2F27A3, 0x00003FFE
 300    data8 0xCBEC14FEF2727C5D, 0x00003FFE
 301    data8 0xCE248C151F8480E4, 0x00003FFE
 302    data8 0xD06333DAEF2B2595, 0x00003FFE
 303    data8 0xD2A81D91F12AE45A, 0x00003FFE
 304    data8 0xD4F35AABCFEDFA1F, 0x00003FFE
 305    data8 0xD744FCCAD69D6AF4, 0x00003FFE
 306    data8 0xD99D15C278AFD7B6, 0x00003FFE
 307    data8 0xDBFBB797DAF23755, 0x00003FFE
 308    data8 0xDE60F4825E0E9124, 0x00003FFE
 309    data8 0xE0CCDEEC2A94E111, 0x00003FFE
 310    data8 0xE33F8972BE8A5A51, 0x00003FFE
 311    data8 0xE5B906E77C8348A8, 0x00003FFE
 312    data8 0xE8396A503C4BDC68, 0x00003FFE
 313    data8 0xEAC0C6E7DD24392F, 0x00003FFE
 314    data8 0xED4F301ED9942B84, 0x00003FFE
 315    data8 0xEFE4B99BDCDAF5CB, 0x00003FFE
 316    data8 0xF281773C59FFB13A, 0x00003FFE
 317    data8 0xF5257D152486CC2C, 0x00003FFE
 318    data8 0xF7D0DF730AD13BB9, 0x00003FFE
 319    data8 0xFA83B2DB722A033A, 0x00003FFE
 320    data8 0xFD3E0C0CF486C175, 0x00003FFE
 321    data8 0x8000000000000000, 0x00003FFF // Center of table
 322    data8 0x8164D1F3BC030773, 0x00003FFF
 323    data8 0x82CD8698AC2BA1D7, 0x00003FFF
 324    data8 0x843A28C3ACDE4046, 0x00003FFF
 325    data8 0x85AAC367CC487B15, 0x00003FFF
 326    data8 0x871F61969E8D1010, 0x00003FFF
 327    data8 0x88980E8092DA8527, 0x00003FFF
 328    data8 0x8A14D575496EFD9A, 0x00003FFF
 329    data8 0x8B95C1E3EA8BD6E7, 0x00003FFF
 330    data8 0x8D1ADF5B7E5BA9E6, 0x00003FFF
 331    data8 0x8EA4398B45CD53C0, 0x00003FFF
 332    data8 0x9031DC431466B1DC, 0x00003FFF
 333    data8 0x91C3D373AB11C336, 0x00003FFF
 334    data8 0x935A2B2F13E6E92C, 0x00003FFF
 335    data8 0x94F4EFA8FEF70961, 0x00003FFF
 336    data8 0x96942D3720185A00, 0x00003FFF
 337    data8 0x9837F0518DB8A96F, 0x00003FFF
 338    data8 0x99E0459320B7FA65, 0x00003FFF
 339    data8 0x9B8D39B9D54E5539, 0x00003FFF
 340    data8 0x9D3ED9A72CFFB751, 0x00003FFF
 341    data8 0x9EF5326091A111AE, 0x00003FFF
 342    data8 0xA0B0510FB9714FC2, 0x00003FFF
 343    data8 0xA27043030C496819, 0x00003FFF
 344    data8 0xA43515AE09E6809E, 0x00003FFF
 345    data8 0xA5FED6A9B15138EA, 0x00003FFF
 346    data8 0xA7CD93B4E965356A, 0x00003FFF
 347    data8 0xA9A15AB4EA7C0EF8, 0x00003FFF
 348    data8 0xAB7A39B5A93ED337, 0x00003FFF
 349    data8 0xAD583EEA42A14AC6, 0x00003FFF
 350    data8 0xAF3B78AD690A4375, 0x00003FFF
 351    data8 0xB123F581D2AC2590, 0x00003FFF
 352    data8 0xB311C412A9112489, 0x00003FFF
 353    data8 0xB504F333F9DE6484, 0x00003FFF
 354 LOCAL_OBJECT_END(cosh_j_hi_table)
 355
 356 LOCAL_OBJECT_START(cosh_j_lo_table)
 357    data4 0x1EB2FB13
 358    data4 0x1CE2CBE2
 359    data4 0x1DDC3CBC
 360    data4 0x1EE9AA34
 361    data4 0x9EAEFDC1
 362    data4 0x9DBF517B
 363    data4 0x1EF88AFB
 364    data4 0x1E03B216
 365    data4 0x1E78AB43
 366    data4 0x9E7B1747
 367    data4 0x9EFE3C0E
 368    data4 0x9D36F837
 369    data4 0x9DEE53E4
 370    data4 0x9E24AE8E
 371    data4 0x1D912473
 372    data4 0x1EB243BE
 373    data4 0x1E669A2F
 374    data4 0x9BBC610A
 375    data4 0x1E761035
 376    data4 0x9E0BE175
 377    data4 0x1CCB12A1
 378    data4 0x1D1BFE90
 379    data4 0x1DF2F47A
 380    data4 0x1EF22F22
 381    data4 0x9E3F4A29
 382    data4 0x1EC01A5B
 383    data4 0x1E8CAC3A
 384    data4 0x9DBB3FAB
 385    data4 0x1EF73A19
 386    data4 0x9BB795B5
 387    data4 0x1EF84B76
 388    data4 0x9EF5818B
 389    data4 0x00000000 // Center of table
 390    data4 0x1F77CACA
 391    data4 0x1EF8A91D
 392    data4 0x1E57C976
 393    data4 0x9EE8DA92
 394    data4 0x1EE85C9F
 395    data4 0x1F3BF1AF
 396    data4 0x1D80CA1E
 397    data4 0x9D0373AF
 398    data4 0x9F167097
 399    data4 0x1EB70051
 400    data4 0x1F6EB029
 401    data4 0x1DFD6D8E
 402    data4 0x9EB319B0
 403    data4 0x1EBA2BEB
 404    data4 0x1F11D537
 405    data4 0x1F0D5A46
 406    data4 0x9E5E7BCA
 407    data4 0x9F3AAFD1
 408    data4 0x9E86DACC
 409    data4 0x9F3EDDC2
 410    data4 0x1E496E3D
 411    data4 0x9F490BF6
 412    data4 0x1DD1DB48
 413    data4 0x1E65EBFB
 414    data4 0x9F427496
 415    data4 0x1F283C4A
 416    data4 0x1F4B0047
 417    data4 0x1F130152
 418    data4 0x9E8367C0
 419    data4 0x9F705F90
 420    data4 0x1EFB3C53
 421    data4 0x1F32FB13
 422 LOCAL_OBJECT_END(cosh_j_lo_table)
 423
 424
 425 .section .text
 426 GLOBAL_IEEE754_ENTRY(coshl)
 427
 428 { .mlx
 429       getf.exp        r_signexp_x = f8   // Get signexp of x, must redo if unorm
 430       movl            r_sig_inv_ln2 = 0xb8aa3b295c17f0bc // significand of 1/ln2
 431 }
 432 { .mlx
 433       addl            r_ad1 = @ltoff(cosh_arg_reduction), gp
 434       movl            r_rshf_2to57 = 0x4778000000000000 // 1.10000 2^(63+57)
 435 }
 436 ;;
 437
 438 { .mfi
 439       ld8             r_ad1 = [r_ad1]
 440       fmerge.s        f_ABS_X    = f0,f8
 441       mov             r_exp_0_25 = 0x0fffd    // Form exponent for 0.25
 442 }
 443 { .mfi
 444       nop.m           0
 445       fnorm.s1        f_NORM_X = f8
 446       mov             r_exp_2tom57 = 0xffff-57
 447 }
 448 ;;
 449
 450 { .mfi
 451       setf.d          f_RSHF_2TO57 = r_rshf_2to57 // Form const 1.100 * 2^120
 452       fclass.m        p10,p0 = f8, 0x0b           // Test for denorm
 453       mov             r_exp_mask = 0x1ffff
 454 }
 455 { .mlx
 456       setf.sig        f_INV_LN2_2TO63 = r_sig_inv_ln2 // Form 1/ln2 * 2^63
 457       movl            r_rshf = 0x43e8000000000000 // 1.1000 2^63 for right shift
 458 }
 459 ;;
 460
 461 { .mfi
 462       nop.m           0
 463       fclass.m        p7,p0 = f8, 0x07  // Test if x=0
 464       nop.i           0
 465 }
 466 { .mfi
 467       setf.exp        f_2TOM57 = r_exp_2tom57 // Form 2^-57 for scaling
 468       nop.f           0
 469       add             r_ad3 = 0x90, r_ad1  // Point to ab_table
 470 }
 471 ;;
 472
 473 { .mfi
 474       setf.d          f_RSHF = r_rshf     // Form right shift const 1.100 * 2^63
 475       fclass.m        p6,p0 = f8, 0xe3     // Test if x nan, inf
 476       add             r_ad4 = 0x2f0, r_ad1 // Point to j_hi_table midpoint
 477 }
 478 { .mib
 479       add             r_ad2e = 0x20, r_ad1 // Point to p_table
 480       nop.i           0
 481 (p10) br.cond.spnt    COSH_DENORM          // Branch if x denorm
 482 }
 483 ;;
 484
 485 // Common path -- return here from COSH_DENORM if x is unnorm
 486 COSH_COMMON:
 487 { .mfi
 488       ldfe            f_smlst_oflow_input = [r_ad2e],16
 489 (p7)  fma.s0          f8 = f1, f1, f0      // Result = 1.0 if x=0
 490       add             r_ad5 = 0x580, r_ad1 // Point to j_lo_table midpoint
 491 }
 492 { .mib
 493       ldfe            f_log2by64_hi  = [r_ad1],16
 494       and             r_exp_x = r_exp_mask, r_signexp_x
 495 (p7)  br.ret.spnt     b0                  // Exit if x=0
 496 }
 497 ;;
 498
 499 // Get the A coefficients for COSH_BY_TBL
 500 { .mfi
 501       ldfe            f_A1 = [r_ad3],16
 502       fcmp.lt.s1      p8,p9 = f8,f0           // Test for x<0
 503       cmp.lt          p7,p0 = r_exp_x, r_exp_0_25  // Test x < 0.25
 504 }
 505 { .mfb
 506       add             r_ad2o = 0x30, r_ad2e  // Point to p_table odd coeffs
 507 (p6)  fma.s0          f8 = f8,f8,f0          // Result for x nan, inf
 508 (p6)  br.ret.spnt     b0                     // Exit for x nan, inf
 509 }
 510 ;;
 511
 512 // Calculate X2 = ax*ax for COSH_BY_POLY
 513 { .mfi
 514       ldfe            f_log2by64_lo  = [r_ad1],16
 515       nop.f           0
 516       nop.i           0
 517 }
 518 { .mfb
 519       ldfe            f_A2 = [r_ad3],16
 520       fma.s1          f_X2 = f_NORM_X, f_NORM_X, f0
 521 (p7)  br.cond.spnt    COSH_BY_POLY
 522 }
 523 ;;
 524
 525 // Here if |x| >= 0.25
 526 COSH_BY_TBL:
 527 // ******************************************************
 528 // STEP 1 (TBL and EXP) - Argument reduction
 529 // ******************************************************
 530 // Get the following constants.
 531 // Inv_log2by64
 532 // log2by64_hi
 533 // log2by64_lo
 534
 535
 536 // We want 2^(N-1) and 2^(-N-1). So bias N-1 and -N-1 and
 537 // put them in an exponent.
 538 // f_spos = 2^(N-1) and f_sneg = 2^(-N-1)
 539 // 0xffff + (N-1)  = 0xffff +N -1
 540 // 0xffff - (N +1) = 0xffff -N -1
 541
 542
 543 // Calculate M and keep it as integer and floating point.
 544 // M = round-to-integer(x*Inv_log2by64)
 545 // f_M = M = truncate(ax/(log2/64))
 546 // Put the integer representation of M in r_M
 547 //    and the floating point representation of M in f_M
 548
 549 // Get the remaining A,B coefficients
 550 { .mmi
 551       ldfe            f_A3 = [r_ad3],16
 552       nop.m           0
 553       nop.i           0
 554 }
 555 ;;
 556
 557 // Use constant (1.100*2^(63-6)) to get rounded M into rightmost significand
 558 // |x| * 64 * 1/ln2 * 2^(63-6) + 1.1000 * 2^(63+(63-6))
 559 { .mfi
 560       nop.m           0
 561       fma.s1          f_M_temp = f_ABS_X, f_INV_LN2_2TO63, f_RSHF_2TO57
 562       mov             r_signexp_0_5 = 0x0fffe // signexp of +0.5
 563 }
 564 ;;
 565
 566 // Test for |x| >= overflow limit
 567 { .mfi
 568       ldfe            f_B1 = [r_ad3],16
 569       fcmp.ge.s1      p6,p0 = f_ABS_X, f_smlst_oflow_input
 570       nop.i           0
 571 }
 572 ;;
 573
 574 { .mfi
 575       ldfe            f_B2 = [r_ad3],16
 576       nop.f           0
 577       mov             r_exp_32 = 0x10004
 578 }
 579 ;;
 580
 581 // Subtract RSHF constant to get rounded M as a floating point value
 582 // M_temp * 2^(63-6) - 2^63
 583 { .mfb
 584       ldfe            f_B3 = [r_ad3],16
 585       fms.s1          f_M = f_M_temp, f_2TOM57, f_RSHF
 586 (p6)  br.cond.spnt    COSH_HUGE  // Branch if result will overflow
 587 }
 588 ;;
 589
 590 { .mfi
 591       getf.sig        r_M = f_M_temp
 592       nop.f           0
 593       cmp.ge          p7,p6 = r_exp_x, r_exp_32 // Test if x >= 32
 594 }
 595 ;;
 596
 597 // Calculate j. j is the signed extension of the six lsb of M. It
 598 // has a range of -32 thru 31.
 599
 600 // Calculate R
 601 // ax - M*log2by64_hi
 602 // R = (ax - M*log2by64_hi) - M*log2by64_lo
 603
 604 { .mfi
 605       nop.m           0
 606       fnma.s1         f_R_temp = f_M, f_log2by64_hi, f_ABS_X
 607       and             r_j = 0x3f, r_M
 608 }
 609 ;;
 610
 611 { .mii
 612       nop.m           0
 613       shl             r_jshf = r_j, 0x2 // Shift j so can sign extend it
 614 ;;
 615       sxt1            r_jshf = r_jshf
 616 }
 617 ;;
 618
 619 { .mii
 620       nop.m           0
 621       shr             r_j = r_jshf, 0x2    // Now j has range -32 to 31
 622       nop.i           0
 623 }
 624 ;;
 625
 626 { .mmi
 627       shladd          r_ad_J_hi = r_j, 4, r_ad4 // pointer to Tjhi
 628       sub             r_Mmj = r_M, r_j          // M-j
 629       sub             r_mj = r0, r_j            // Form -j
 630 }
 631 ;;
 632
 633 // The TBL and EXP branches are merged and predicated
 634 // If TBL, p6 true, 0.25 <= |x| < 32
 635 // If EXP, p7 true, 32 <= |x| < overflow_limit
 636 //
 637 // N = (M-j)/64
 638 { .mfi
 639       ldfe            f_Tjhi = [r_ad_J_hi]
 640       fnma.s1         f_R = f_M, f_log2by64_lo, f_R_temp
 641       shr             r_N = r_Mmj, 0x6            // N = (M-j)/64
 642 }
 643 { .mfi
 644       shladd          r_ad_mJ_hi = r_mj, 4, r_ad4 // pointer to Tmjhi
 645       nop.f           0
 646       shladd          r_ad_mJ_lo = r_mj, 2, r_ad5 // pointer to Tmjlo
 647 }
 648 ;;
 649
 650 { .mfi
 651       sub             r_2mNm1 = r_signexp_0_5, r_N // signexp 2^(-N-1)
 652       nop.f           0
 653       shladd          r_ad_J_lo = r_j, 2, r_ad5   // pointer to Tjlo
 654 }
 655 { .mfi
 656       ldfe            f_Tmjhi = [r_ad_mJ_hi]
 657       nop.f           0
 658       add             r_2Nm1 = r_signexp_0_5, r_N // signexp 2^(N-1)
 659 }
 660 ;;
 661
 662 { .mmf
 663       ldfs            f_Tmjlo = [r_ad_mJ_lo]
 664       setf.exp        f_sneg = r_2mNm1            // Form 2^(-N-1)
 665       nop.f           0
 666 }
 667 ;;
 668
 669 { .mmf
 670       ldfs            f_Tjlo  = [r_ad_J_lo]
 671       setf.exp        f_spos = r_2Nm1             // Form 2^(N-1)
 672       nop.f           0
 673 }
 674 ;;
 675
 676 // ******************************************************
 677 // STEP 2 (TBL and EXP)
 678 // ******************************************************
 679 // Calculate Rsquared and Rcubed in preparation for p_even and p_odd
 680
 681 { .mmf
 682       nop.m           0
 683       nop.m           0
 684       fma.s1          f_Rsq  = f_R, f_R, f0
 685 }
 686 ;;
 687
 688
 689 // Calculate p_even
 690 // B_2 + Rsq *B_3
 691 // B_1 + Rsq * (B_2 + Rsq *B_3)
 692 // p_even = Rsq * (B_1 + Rsq * (B_2 + Rsq *B_3))
 693 { .mfi
 694       nop.m           0
 695       fma.s1          f_peven_temp1 = f_Rsq, f_B3, f_B2
 696       nop.i           0
 697 }
 698 // Calculate p_odd
 699 // A_2 + Rsq *A_3
 700 // A_1 + Rsq * (A_2 + Rsq *A_3)
 701 // podd = R + Rcub * (A_1 + Rsq * (A_2 + Rsq *A_3))
 702 { .mfi
 703       nop.m           0
 704       fma.s1          f_podd_temp1 = f_Rsq, f_A3, f_A2
 705       nop.i           0
 706 }
 707 ;;
 708
 709 { .mfi
 710       nop.m           0
 711       fma.s1          f_Rcub = f_Rsq, f_R, f0
 712       nop.i           0
 713 }
 714 ;;
 715
 716 //
 717 // If TBL,
 718 // Calculate S_hi and S_lo, and C_hi
 719 // SC_hi_temp = sneg * Tmjhi
 720 // S_hi = spos * Tjhi - SC_hi_temp
 721 // S_hi = spos * Tjhi - (sneg * Tmjhi)
 722 // C_hi = spos * Tjhi + SC_hi_temp
 723 // C_hi = spos * Tjhi + (sneg * Tmjhi)
 724
 725 { .mfi
 726       nop.m           0
 727 (p6)  fma.s1          f_SC_hi_temp = f_sneg, f_Tmjhi, f0
 728       nop.i           0
 729 }
 730 ;;
 731
 732 // If TBL,
 733 // C_lo_temp3 = sneg * Tmjlo
 734 // C_lo_temp4 = spos * Tjlo + C_lo_temp3
 735 // C_lo_temp4 = spos * Tjlo + (sneg * Tmjlo)
 736 { .mfi
 737       nop.m           0
 738 (p6)  fma.s1          f_C_lo_temp3 =  f_sneg, f_Tmjlo, f0
 739       nop.i           0
 740 }
 741 ;;
 742
 743 { .mfi
 744       nop.m           0
 745       fma.s1          f_peven_temp2 = f_Rsq, f_peven_temp1, f_B1
 746       nop.i           0
 747 }
 748 { .mfi
 749       nop.m           0
 750       fma.s1          f_podd_temp2 = f_Rsq, f_podd_temp1, f_A1
 751       nop.i           0
 752 }
 753 ;;
 754
 755 // If EXP,
 756 // Compute 2^(N-1) * Tjhi and 2^(N-1) * Tjlo
 757 { .mfi
 758       nop.m           0
 759 (p7)  fma.s1          f_Tjhi_spos = f_Tjhi, f_spos, f0
 760       nop.i           0
 761 }
 762 { .mfi
 763       nop.m           0
 764 (p7)  fma.s1          f_Tjlo_spos = f_Tjlo, f_spos, f0
 765       nop.i           0
 766 }
 767 ;;
 768
 769 { .mfi
 770       nop.m           0
 771 (p6)  fma.s1          f_C_hi = f_spos, f_Tjhi, f_SC_hi_temp
 772       nop.i           0
 773 }
 774 ;;
 775
 776 { .mfi
 777       nop.m           0
 778 (p6)  fms.s1          f_S_hi = f_spos, f_Tjhi, f_SC_hi_temp
 779       nop.i           0
 780 }
 781 { .mfi
 782       nop.m           0
 783 (p6)  fma.s1          f_C_lo_temp4 = f_spos, f_Tjlo, f_C_lo_temp3
 784       nop.i           0
 785 }
 786 ;;
 787
 788 { .mfi
 789       nop.m           0
 790       fma.s1          f_peven = f_Rsq, f_peven_temp2, f0
 791       nop.i           0
 792 }
 793 { .mfi
 794       nop.m           0
 795       fma.s1          f_podd = f_podd_temp2, f_Rcub, f_R
 796       nop.i           0
 797 }
 798 ;;
 799
 800 // If TBL,
 801 // C_lo_temp1 =  spos * Tjhi - C_hi
 802 // C_lo_temp2 =  sneg * Tmjlo + C_lo_temp1
 803 // C_lo_temp2 =  sneg * Tmjlo + (spos * Tjhi - C_hi)
 804
 805 { .mfi
 806       nop.m           0
 807 (p6)  fms.s1          f_C_lo_temp1 =  f_spos, f_Tjhi,  f_C_hi
 808       nop.i           0
 809 }
 810 ;;
 811
 812 { .mfi
 813       nop.m           0
 814 (p6)  fma.s1          f_C_lo_temp2 = f_sneg, f_Tmjhi, f_C_lo_temp1
 815       nop.i           0
 816 }
 817 ;;
 818
 819 // If EXP,
 820 // Y_hi = 2^(N-1) * Tjhi
 821 // Y_lo = 2^(N-1) * Tjhi * (p_odd + p_even) + 2^(N-1) * Tjlo
 822 { .mfi
 823       nop.m           0
 824 (p7)  fma.s1          f_Y_lo_temp =  f_peven, f1, f_podd
 825       nop.i           0
 826 }
 827 ;;
 828
 829 // If TBL,
 830 // C_lo = C_lo_temp4 + C_lo_temp2
 831 { .mfi
 832       nop.m           0
 833 (p6)  fma.s1          f_C_lo = f_C_lo_temp4, f1, f_C_lo_temp2
 834       nop.i           0
 835 }
 836 ;;
 837
 838 // If TBL,
 839 // Y_hi = C_hi
 840 // Y_lo = S_hi*p_odd + (C_hi*p_even + C_lo)
 841 { .mfi
 842       nop.m           0
 843 (p6)  fma.s1          f_Y_lo_temp = f_C_hi, f_peven, f_C_lo
 844       nop.i           0
 845 }
 846 ;;
 847
 848 { .mfi
 849       nop.m           0
 850 (p7)  fma.s1          f_Y_lo = f_Tjhi_spos, f_Y_lo_temp, f_Tjlo_spos
 851       nop.i           0
 852 }
 853 ;;
 854
 855 // Dummy multiply to generate inexact
 856 { .mfi
 857       nop.m           0
 858       fmpy.s0         f_tmp = f_B2, f_B2
 859       nop.i           0
 860 }
 861 { .mfi
 862       nop.m           0
 863 (p6)  fma.s1          f_Y_lo = f_S_hi, f_podd, f_Y_lo_temp
 864       nop.i           0
 865 }
 866 ;;
 867
 868 // f8 = answer = Y_hi + Y_lo
 869 { .mfi
 870       nop.m           0
 871 (p7)  fma.s0          f8 = f_Y_lo,  f1, f_Tjhi_spos
 872       nop.i           0
 873 }
 874 ;;
 875
 876 // f8 = answer = Y_hi + Y_lo
 877 { .mfb
 878       nop.m           0
 879 (p6)  fma.s0          f8 = f_Y_lo, f1, f_C_hi
 880       br.ret.sptk     b0      // Exit for COSH_BY_TBL and COSH_BY_EXP
 881 }
 882 ;;
 883
 884
 885 // Here if 0 < |x| < 0.25
 886 COSH_BY_POLY:
 887 { .mmf
 888       ldfe            f_P6 = [r_ad2e],16
 889       ldfe            f_P5 = [r_ad2o],16
 890       nop.f           0
 891 }
 892 ;;
 893
 894 { .mmi
 895       ldfe            f_P4 = [r_ad2e],16
 896       ldfe            f_P3 = [r_ad2o],16
 897       nop.i           0
 898 }
 899 ;;
 900
 901 { .mmi
 902       ldfe            f_P2 = [r_ad2e],16
 903       ldfe            f_P1 = [r_ad2o],16
 904       nop.i           0
 905 }
 906 ;;
 907
 908 { .mfi
 909       nop.m           0
 910       fma.s1          f_X3 = f_NORM_X, f_X2, f0
 911       nop.i           0
 912 }
 913 { .mfi
 914       nop.m           0
 915       fma.s1          f_X4 = f_X2, f_X2, f0
 916       nop.i           0
 917 }
 918 ;;
 919
 920 { .mfi
 921       nop.m           0
 922       fma.s1          f_poly65 = f_X2, f_P6, f_P5
 923       nop.i           0
 924 }
 925 { .mfi
 926       nop.m           0
 927       fma.s1          f_poly43 = f_X2, f_P4, f_P3
 928       nop.i           0
 929 }
 930 ;;
 931
 932 { .mfi
 933       nop.m           0
 934       fma.s1          f_poly21 = f_X2, f_P2, f_P1
 935       nop.i           0
 936 }
 937 ;;
 938
 939 { .mfi
 940       nop.m           0
 941       fma.s1          f_poly6543 = f_X4, f_poly65, f_poly43
 942       nop.i           0
 943 }
 944 ;;
 945
 946 { .mfi
 947       nop.m           0
 948       fma.s1          f_poly6to1 = f_X4, f_poly6543, f_poly21
 949       nop.i           0
 950 }
 951 ;;
 952
 953 // Dummy multiply to generate inexact
 954 { .mfi
 955       nop.m           0
 956       fmpy.s0         f_tmp = f_P6, f_P6
 957       nop.i           0
 958 }
 959 { .mfb
 960       nop.m           0
 961       fma.s0          f8 = f_poly6to1, f_X2, f1
 962       br.ret.sptk     b0                // Exit COSH_BY_POLY
 963 }
 964 ;;
 965
 966
 967 // Here if x denorm or unorm
 968 COSH_DENORM:
 969 // Determine if x really a denorm and not a unorm
 970 { .mmf
 971       getf.exp        r_signexp_x = f_NORM_X
 972       mov             r_exp_denorm = 0x0c001   // Real denorms have exp < this
 973       fmerge.s        f_ABS_X = f0, f_NORM_X
 974 }
 975 ;;
 976
 977 { .mfi
 978       nop.m           0
 979       fcmp.eq.s0      p10,p0 = f8, f0  // Set denorm flag
 980       nop.i           0
 981 }
 982 ;;
 983
 984 // Set p8 if really a denorm
 985 { .mmi
 986       and             r_exp_x = r_exp_mask, r_signexp_x
 987 ;;
 988       cmp.lt          p8,p9 = r_exp_x, r_exp_denorm
 989       nop.i           0
 990 }
 991 ;;
 992
 993 // Identify denormal operands.
 994 { .mfb
 995       nop.m           0
 996 (p8)  fma.s0          f8 =  f8,f8,f1 // If x denorm, result=1+x^2
 997 (p9)  br.cond.sptk    COSH_COMMON    // Return to main path if x unorm
 998 }
 999 ;;
1000
1001 { .mfb
1002       nop.m           0
1003       nop.f           0
1004       br.ret.sptk     b0            // Exit if x denorm
1005 }
1006 ;;
1007
1008
1009 // Here if |x| >= overflow limit
1010 COSH_HUGE:
1011 // for COSH_HUGE, put 24000 in exponent; take sign from input
1012 { .mmi
1013       mov             r_exp_huge = 0x15dbf
1014 ;;
1015       setf.exp        f_huge  = r_exp_huge
1016       nop.i           0
1017 }
1018 ;;
1019
1020 { .mfi
1021       alloc           r32 = ar.pfs,0,5,4,0
1022       fma.s1          f_signed_hi_lo = f_huge, f1, f1
1023       nop.i           0
1024 }
1025 ;;
1026
1027 { .mfi
1028       nop.m           0
1029       fma.s0          f_pre_result = f_signed_hi_lo, f_huge, f0
1030       mov             GR_Parameter_TAG = 63
1031 }
1032 ;;
1033
1034 GLOBAL_IEEE754_END(coshl)
1035
1036
1037 LOCAL_LIBM_ENTRY(__libm_error_region)
1038 .prologue
1039
1040 { .mfi
1041         add   GR_Parameter_Y=-32,sp              // Parameter 2 value
1042         nop.f 0
1043 .save   ar.pfs,GR_SAVE_PFS
1044         mov  GR_SAVE_PFS=ar.pfs                  // Save ar.pfs
1045 }
1046 { .mfi
1047 .fframe 64
1048         add sp=-64,sp                            // Create new stack
1049         nop.f 0
1050         mov GR_SAVE_GP=gp                        // Save gp
1051 };;
1052
1053 { .mmi
1054         stfe [GR_Parameter_Y] = f0,16            // STORE Parameter 2 on stack
1055         add GR_Parameter_X = 16,sp               // Parameter 1 address
1056 .save   b0, GR_SAVE_B0
1057         mov GR_SAVE_B0=b0                        // Save b0
1058 };;
1059
1060 .body
1061 { .mib
1062         stfe [GR_Parameter_X] = f8               // STORE Parameter 1 on stack
1063         add   GR_Parameter_RESULT = 0,GR_Parameter_Y   // Parameter 3 address
1064         nop.b 0
1065 }
1066 { .mib
1067         stfe [GR_Parameter_Y] = f_pre_result     // STORE Parameter 3 on stack
1068         add   GR_Parameter_Y = -16,GR_Parameter_Y
1069         br.call.sptk b0=__libm_error_support#    // Call error handling function
1070 };;
1071
1072 { .mmi
1073         add   GR_Parameter_RESULT = 48,sp
1074         nop.m 0
1075         nop.i 0
1076 };;
1077
1078 { .mmi
1079         ldfe  f8 = [GR_Parameter_RESULT]         // Get return result off stack
1080 .restore sp
1081         add   sp = 64,sp                         // Restore stack pointer
1082         mov   b0 = GR_SAVE_B0                    // Restore return address
1083 };;
1084
1085 { .mib
1086         mov   gp = GR_SAVE_GP                    // Restore gp
1087         mov   ar.pfs = GR_SAVE_PFS               // Restore ar.pfs
1088         br.ret.sptk     b0                       // Return
1089 };;
1090
1091 LOCAL_LIBM_END(__libm_error_region)
1092
1093
1094 .type   __libm_error_support#,@function
1095 .global __libm_error_support#