sysdeps/ia64/fpu/s_expm1l.S

   1 .file "exp_m1l.s"
   2
   3 // Copyright (C) 2000, 2001, Intel Corporation
   4 // All rights reserved.
   5 //
   6 // Contributed 2/2/2000 by John Harrison, Ted Kubaska, Bob Norin, Shane Story,
   7 // and Ping Tak Peter Tang of the Computational Software Lab, 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://developer.intel.com/opensource.
  39 //
  40 // History
  41 //==============================================================
  42 // 4/04/00  Unwind support added
  43 // 8/15/00  Bundle added after call to __libm_error_support to properly
  44 //          set [the previously overwritten] GR_Parameter_RESULT.
  45 //
  46 // *********************************************************************
  47 //
  48 // Function:   Combined expl(x) and expm1l(x), where
  49 //                        x
  50 //             expl(x) = e , for double-extended precision x values
  51 //                          x
  52 //             expm1l(x) = e  - 1  for double-extended precision x values
  53 //
  54 // *********************************************************************
  55 //
  56 // Resources Used:
  57 //
  58 //    Floating-Point Registers: f8  (Input and Return Value)
  59 //                              f9,f32-f61, f99-f102
  60 //
  61 //    General Purpose Registers:
  62 //      r32-r61
  63 //      r62-r65 (Used to pass arguments to error handling routine)
  64 //
  65 //    Predicate Registers:      p6-p15
  66 //
  67 // *********************************************************************
  68 //
  69 // IEEE Special Conditions:
  70 //
  71 //    Denormal  fault raised on denormal inputs
  72 //    Overflow exceptions raised when appropriate for exp and expm1
  73 //    Underflow exceptions raised when appropriate for exp and expm1
  74 //    (Error Handling Routine called for overflow and Underflow)
  75 //    Inexact raised when appropriate by algorithm
  76 //
  77 //    expl(inf) = inf
  78 //    expl(-inf) = +0
  79 //    expl(SNaN) = QNaN
  80 //    expl(QNaN) = QNaN
  81 //    expl(0) = 1
  82 //    expl(EM_special Values) = QNaN
  83 //    expl(inf) = inf
  84 //    expm1l(-inf) = -1
  85 //    expm1l(SNaN) = QNaN
  86 //    expm1l(QNaN) = QNaN
  87 //    expm1l(0) = 0
  88 //    expm1l(EM_special Values) = QNaN
  89 //
  90 // *********************************************************************
  91 //
  92 // Implementation and Algorithm Notes:
  93 //
  94 //  ker_exp_64( in_FR  : X,
  95 //            in_GR  : Flag,
  96 //            in_GR  : Expo_Range
  97 //            out_FR : Y_hi,
  98 //            out_FR : Y_lo,
  99 //            out_FR : scale,
 100 //            out_PR : Safe )
 101 //
 102 // On input, X is in register format and
 103 // Flag  = 0 for exp,
 104 // Flag  = 1 for expm1,
 105 //
 106 // On output, provided X and X_cor are real numbers, then
 107 //
 108 //   scale*(Y_hi + Y_lo)  approximates  expl(X)       if Flag is 0
 109 //   scale*(Y_hi + Y_lo)  approximates  expl(X)-1     if Flag is 1
 110 //
 111 // The accuracy is sufficient for a highly accurate 64 sig.
 112 // bit implementation.  Safe is set if there is no danger of
 113 // overflow/underflow when the result is composed from scale,
 114 // Y_hi and Y_lo. Thus, we can have a fast return if Safe is set.
 115 // Otherwise, one must prepare to handle the possible exception
 116 // appropriately.  Note that SAFE not set (false) does not mean
 117 // that overflow/underflow will occur; only the setting of SAFE
 118 // guarantees the opposite.
 119 //
 120 // **** High Level Overview ****
 121 //
 122 // The method consists of three cases.
 123 //
 124 // If           |X| < Tiny      use case exp_tiny;
 125 // else if      |X| < 2^(-6)    use case exp_small;
 126 // else         use case exp_regular;
 127 //
 128 // Case exp_tiny:
 129 //
 130 //   1 + X     can be used to approximate expl(X) or expl(X+X_cor);
 131 //   X + X^2/2 can be used to approximate expl(X) - 1
 132 //
 133 // Case exp_small:
 134 //
 135 //   Here, expl(X), expl(X+X_cor), and expl(X) - 1 can all be
 136 //   appproximated by a relatively simple polynomial.
 137 //
 138 //   This polynomial resembles the truncated Taylor series
 139 //
 140 //      expl(w) = 1 + w + w^2/2! + w^3/3! + ... + w^n/n!
 141 //
 142 // Case exp_regular:
 143 //
 144 //   Here we use a table lookup method. The basic idea is that in
 145 //   order to compute expl(X), we accurately decompose X into
 146 //
 147 //   X = N * log(2)/(2^12)  + r,        |r| <= log(2)/2^13.
 148 //
 149 //   Hence
 150 //
 151 //   expl(X) = 2^( N / 2^12 ) * expl(r).
 152 //
 153 //   The value 2^( N / 2^12 ) is obtained by simple combinations
 154 //   of values calculated beforehand and stored in table; expl(r)
 155 //   is approximated by a short polynomial because |r| is small.
 156 //
 157 //   We elaborate this method in 4 steps.
 158 //
 159 //   Step 1: Reduction
 160 //
 161 //   The value 2^12/log(2) is stored as a double-extended number
 162 //   L_Inv.
 163 //
 164 //   N := round_to_nearest_integer( X * L_Inv )
 165 //
 166 //   The value log(2)/2^12 is stored as two numbers L_hi and L_lo so
 167 //   that r can be computed accurately via
 168 //
 169 //   r := (X - N*L_hi) - N*L_lo
 170 //
 171 //   We pick L_hi such that N*L_hi is representable in 64 sig. bits
 172 //   and thus the FMA   X - N*L_hi   is error free. So r is the
 173 //   1 rounding error from an exact reduction with respect to
 174 //
 175 //   L_hi + L_lo.
 176 //
 177 //   In particular, L_hi has 30 significant bit and can be stored
 178 //   as a double-precision number; L_lo has 64 significant bits and
 179 //   stored as a double-extended number.
 180 //
 181 //   In the case Flag = 2, we further modify r by
 182 //
 183 //   r := r + X_cor.
 184 //
 185 //   Step 2: Approximation
 186 //
 187 //   expl(r) - 1 is approximated by a short polynomial of the form
 188 //
 189 //   r + A_1 r^2 + A_2 r^3 + A_3 r^4 .
 190 //
 191 //   Step 3: Composition from Table Values
 192 //
 193 //   The value 2^( N / 2^12 ) can be composed from a couple of tables
 194 //   of precalculated values. First, express N as three integers
 195 //   K, M_1, and M_2 as
 196 //
 197 //     N  =  K * 2^12  + M_1 * 2^6 + M_2
 198 //
 199 //   Where 0 <= M_1, M_2 < 2^6; and K can be positive or negative.
 200 //   When N is represented in 2's complement, M_2 is simply the 6
 201 //   lsb's, M_1 is the next 6, and K is simply N shifted right
 202 //   arithmetically (sign extended) by 12 bits.
 203 //
 204 //   Now, 2^( N / 2^12 ) is simply
 205 //
 206 //      2^K * 2^( M_1 / 2^6 ) * 2^( M_2 / 2^12 )
 207 //
 208 //   Clearly, 2^K needs no tabulation. The other two values are less
 209 //   trivial because if we store each accurately to more than working
 210 //   precision, than its product is too expensive to calculate. We
 211 //   use the following method.
 212 //
 213 //   Define two mathematical values, delta_1 and delta_2, implicitly
 214 //   such that
 215 //
 216 //     T_1 = expl( [M_1 log(2)/2^6]  -  delta_1 )
 217 //     T_2 = expl( [M_2 log(2)/2^12] -  delta_2 )
 218 //
 219 //   are representable as 24 significant bits. To illustrate the idea,
 220 //   we show how we define delta_1:
 221 //
 222 //     T_1     := round_to_24_bits( expl( M_1 log(2)/2^6 ) )
 223 //     delta_1  = (M_1 log(2)/2^6) - log( T_1 )
 224 //
 225 //   The last equality means mathematical equality. We then tabulate
 226 //
 227 //     W_1 := expl(delta_1) - 1
 228 //     W_2 := expl(delta_2) - 1
 229 //
 230 //   Both in double precision.
 231 //
 232 //   From the tabulated values T_1, T_2, W_1, W_2, we compose the values
 233 //   T and W via
 234 //
 235 //     T := T_1 * T_2                   ...exactly
 236 //     W := W_1 + (1 + W_1)*W_2
 237 //
 238 //   W approximates expl( delta ) - 1  where delta = delta_1 + delta_2.
 239 //   The mathematical product of T and (W+1) is an accurate representation
 240 //   of 2^(M_1/2^6) * 2^(M_2/2^12).
 241 //
 242 //   Step 4. Reconstruction
 243 //
 244 //   Finally, we can reconstruct expl(X), expl(X) - 1.
 245 //   Because
 246 //
 247 //      X = K * log(2) + (M_1*log(2)/2^6  - delta_1)
 248 //                     + (M_2*log(2)/2^12 - delta_2)
 249 //                     + delta_1 + delta_2 + r          ...accurately
 250 //   We have
 251 //
 252 //      expl(X) ~=~ 2^K * ( T + T*[expl(delta_1+delta_2+r) - 1] )
 253 //             ~=~ 2^K * ( T + T*[expl(delta + r) - 1]         )
 254 //             ~=~ 2^K * ( T + T*[(expl(delta)-1)
 255 //                               + expl(delta)*(expl(r)-1)]   )
 256 //             ~=~ 2^K * ( T + T*( W + (1+W)*poly(r) ) )
 257 //             ~=~ 2^K * ( Y_hi  +  Y_lo )
 258 //
 259 //   where Y_hi = T  and Y_lo = T*(W + (1+W)*poly(r))
 260 //
 261 //   For expl(X)-1, we have
 262 //
 263 //      expl(X)-1 ~=~ 2^K * ( Y_hi + Y_lo ) - 1
 264 //               ~=~ 2^K * ( Y_hi + Y_lo - 2^(-K) )
 265 //
 266 //   and we combine Y_hi + Y_lo - 2^(-N)  into the form of two
 267 //   numbers  Y_hi + Y_lo carefully.
 268 //
 269 //   **** Algorithm Details ****
 270 //
 271 //   A careful algorithm must be used to realize the mathematical ideas
 272 //   accurately. We describe each of the three cases. We assume SAFE
 273 //   is preset to be TRUE.
 274 //
 275 //   Case exp_tiny:
 276 //
 277 //   The important points are to ensure an accurate result under
 278 //   different rounding directions and a correct setting of the SAFE
 279 //   flag.
 280 //
 281 //   If Flag is 1, then
 282 //      SAFE  := False  ...possibility of underflow
 283 //      Scale := 1.0
 284 //      Y_hi  := X
 285 //      Y_lo  := 2^(-17000)
 286 //   Else
 287 //      Scale := 1.0
 288 //      Y_hi  := 1.0
 289 //      Y_lo  := X      ...for different rounding modes
 290 //   Endif
 291 //
 292 //   Case exp_small:
 293 //
 294 //   Here we compute a simple polynomial. To exploit parallelism, we split
 295 //   the polynomial into several portions.
 296 //
 297 //   Let r = X
 298 //
 299 //   If Flag is not 1   ...i.e. expl( argument )
 300 //
 301 //      rsq := r * r;
 302 //      r4  := rsq*rsq
 303 //      poly_lo := P_3 + r*(P_4 + r*(P_5 + r*P_6))
 304 //      poly_hi := r + rsq*(P_1 + r*P_2)
 305 //      Y_lo    := poly_hi + r4 * poly_lo
 306 //      set lsb(Y_lo) to 1
 307 //      Y_hi    := 1.0
 308 //      Scale   := 1.0
 309 //
 310 //   Else                       ...i.e. expl( argument ) - 1
 311 //
 312 //      rsq := r * r
 313 //      r4  := rsq * rsq
 314 //      r6  := rsq * r4
 315 //      poly_lo := r6*(Q_5 + r*(Q_6 + r*Q_7))
 316 //      poly_hi := Q_1 + r*(Q_2 + r*(Q_3 + r*Q_4))
 317 //      Y_lo    := rsq*poly_hi +  poly_lo
 318 //      set lsb(Y_lo) to 1
 319 //      Y_hi    := X
 320 //      Scale   := 1.0
 321 //
 322 //   Endif
 323 //
 324 //  Case exp_regular:
 325 //
 326 //  The previous description contain enough information except the
 327 //  computation of poly and the final Y_hi and Y_lo in the case for
 328 //  expl(X)-1.
 329 //
 330 //  The computation of poly for Step 2:
 331 //
 332 //   rsq := r*r
 333 //   poly := r + rsq*(A_1 + r*(A_2 + r*A_3))
 334 //
 335 //  For the case expl(X) - 1, we need to incorporate 2^(-K) into
 336 //  Y_hi and Y_lo at the end of Step 4.
 337 //
 338 //   If K > 10 then
 339 //      Y_lo := Y_lo - 2^(-K)
 340 //   Else
 341 //      If K < -10 then
 342 //       Y_lo := Y_hi + Y_lo
 343 //       Y_hi := -2^(-K)
 344 //      Else
 345 //       Y_hi := Y_hi - 2^(-K)
 346 //      End If
 347 //   End If
 348 //
 349
 350 #include "libm_support.h"
 351
 352 #ifdef _LIBC
 353 .rodata
 354 #else
 355 .data
 356 #endif
 357
 358 .align 64
 359 Constants_exp_64_Arg:
 360 ASM_TYPE_DIRECTIVE(Constants_exp_64_Arg,@object)
 361 data4 0x5C17F0BC,0xB8AA3B29,0x0000400B,0x00000000
 362 data4 0x00000000,0xB17217F4,0x00003FF2,0x00000000
 363 data4 0xF278ECE6,0xF473DE6A,0x00003FD4,0x00000000
 364 // /* Inv_L, L_hi, L_lo */
 365 ASM_SIZE_DIRECTIVE(Constants_exp_64_Arg)
 366
 367 .align 64
 368 Constants_exp_64_Exponents:
 369 ASM_TYPE_DIRECTIVE(Constants_exp_64_Exponents,@object)
 370 data4 0x0000007E,0x00000000,0xFFFFFF83,0xFFFFFFFF
 371 data4 0x000003FE,0x00000000,0xFFFFFC03,0xFFFFFFFF
 372 data4 0x00003FFE,0x00000000,0xFFFFC003,0xFFFFFFFF
 373 data4 0x00003FFE,0x00000000,0xFFFFC003,0xFFFFFFFF
 374 data4 0xFFFFFFE2,0xFFFFFFFF,0xFFFFFFC4,0xFFFFFFFF
 375 data4 0xFFFFFFBA,0xFFFFFFFF,0xFFFFFFBA,0xFFFFFFFF
 376 ASM_SIZE_DIRECTIVE(Constants_exp_64_Exponents)
 377
 378 .align 64
 379 Constants_exp_64_A:
 380 ASM_TYPE_DIRECTIVE(Constants_exp_64_A,@object)
 381 data4 0xB1B736A0,0xAAAAAAAB,0x00003FFA,0x00000000
 382 data4 0x90CD6327,0xAAAAAAAB,0x00003FFC,0x00000000
 383 data4 0xFFFFFFFF,0xFFFFFFFF,0x00003FFD,0x00000000
 384 // /* Reversed */
 385 ASM_SIZE_DIRECTIVE(Constants_exp_64_A)
 386
 387 .align 64
 388 Constants_exp_64_P:
 389 ASM_TYPE_DIRECTIVE(Constants_exp_64_P,@object)
 390 data4 0x43914A8A,0xD00D6C81,0x00003FF2,0x00000000
 391 data4 0x30304B30,0xB60BC4AC,0x00003FF5,0x00000000
 392 data4 0x7474C518,0x88888888,0x00003FF8,0x00000000
 393 data4 0x8DAE729D,0xAAAAAAAA,0x00003FFA,0x00000000
 394 data4 0xAAAAAF61,0xAAAAAAAA,0x00003FFC,0x00000000
 395 data4 0x000004C7,0x80000000,0x00003FFE,0x00000000
 396 // /* Reversed */
 397 ASM_SIZE_DIRECTIVE(Constants_exp_64_P)
 398
 399 .align 64
 400 Constants_exp_64_Q:
 401 ASM_TYPE_DIRECTIVE(Constants_exp_64_Q,@object)
 402 data4 0xA49EF6CA,0xD00D56F7,0x00003FEF,0x00000000
 403 data4 0x1C63493D,0xD00D59AB,0x00003FF2,0x00000000
 404 data4 0xFB50CDD2,0xB60B60B5,0x00003FF5,0x00000000
 405 data4 0x7BA68DC8,0x88888888,0x00003FF8,0x00000000
 406 data4 0xAAAAAC8D,0xAAAAAAAA,0x00003FFA,0x00000000
 407 data4 0xAAAAACCA,0xAAAAAAAA,0x00003FFC,0x00000000
 408 data4 0x00000000,0x80000000,0x00003FFE,0x00000000
 409 // /* Reversed */
 410 ASM_SIZE_DIRECTIVE(Constants_exp_64_Q)
 411
 412 .align 64
 413 Constants_exp_64_T1:
 414 ASM_TYPE_DIRECTIVE(Constants_exp_64_T1,@object)
 415 data4 0x3F800000,0x3F8164D2,0x3F82CD87,0x3F843A29
 416 data4 0x3F85AAC3,0x3F871F62,0x3F88980F,0x3F8A14D5
 417 data4 0x3F8B95C2,0x3F8D1ADF,0x3F8EA43A,0x3F9031DC
 418 data4 0x3F91C3D3,0x3F935A2B,0x3F94F4F0,0x3F96942D
 419 data4 0x3F9837F0,0x3F99E046,0x3F9B8D3A,0x3F9D3EDA
 420 data4 0x3F9EF532,0x3FA0B051,0x3FA27043,0x3FA43516
 421 data4 0x3FA5FED7,0x3FA7CD94,0x3FA9A15B,0x3FAB7A3A
 422 data4 0x3FAD583F,0x3FAF3B79,0x3FB123F6,0x3FB311C4
 423 data4 0x3FB504F3,0x3FB6FD92,0x3FB8FBAF,0x3FBAFF5B
 424 data4 0x3FBD08A4,0x3FBF179A,0x3FC12C4D,0x3FC346CD
 425 data4 0x3FC5672A,0x3FC78D75,0x3FC9B9BE,0x3FCBEC15
 426 data4 0x3FCE248C,0x3FD06334,0x3FD2A81E,0x3FD4F35B
 427 data4 0x3FD744FD,0x3FD99D16,0x3FDBFBB8,0x3FDE60F5
 428 data4 0x3FE0CCDF,0x3FE33F89,0x3FE5B907,0x3FE8396A
 429 data4 0x3FEAC0C7,0x3FED4F30,0x3FEFE4BA,0x3FF28177
 430 data4 0x3FF5257D,0x3FF7D0DF,0x3FFA83B3,0x3FFD3E0C
 431 ASM_SIZE_DIRECTIVE(Constants_exp_64_T1)
 432
 433 .align 64
 434 Constants_exp_64_T2:
 435 ASM_TYPE_DIRECTIVE(Constants_exp_64_T2,@object)
 436 data4 0x3F800000,0x3F80058C,0x3F800B18,0x3F8010A4
 437 data4 0x3F801630,0x3F801BBD,0x3F80214A,0x3F8026D7
 438 data4 0x3F802C64,0x3F8031F2,0x3F803780,0x3F803D0E
 439 data4 0x3F80429C,0x3F80482B,0x3F804DB9,0x3F805349
 440 data4 0x3F8058D8,0x3F805E67,0x3F8063F7,0x3F806987
 441 data4 0x3F806F17,0x3F8074A8,0x3F807A39,0x3F807FCA
 442 data4 0x3F80855B,0x3F808AEC,0x3F80907E,0x3F809610
 443 data4 0x3F809BA2,0x3F80A135,0x3F80A6C7,0x3F80AC5A
 444 data4 0x3F80B1ED,0x3F80B781,0x3F80BD14,0x3F80C2A8
 445 data4 0x3F80C83C,0x3F80CDD1,0x3F80D365,0x3F80D8FA
 446 data4 0x3F80DE8F,0x3F80E425,0x3F80E9BA,0x3F80EF50
 447 data4 0x3F80F4E6,0x3F80FA7C,0x3F810013,0x3F8105AA
 448 data4 0x3F810B41,0x3F8110D8,0x3F81166F,0x3F811C07
 449 data4 0x3F81219F,0x3F812737,0x3F812CD0,0x3F813269
 450 data4 0x3F813802,0x3F813D9B,0x3F814334,0x3F8148CE
 451 data4 0x3F814E68,0x3F815402,0x3F81599C,0x3F815F37
 452 ASM_SIZE_DIRECTIVE(Constants_exp_64_T2)
 453
 454 .align 64
 455 Constants_exp_64_W1:
 456 ASM_TYPE_DIRECTIVE(Constants_exp_64_W1,@object)
 457 data4 0x00000000,0x00000000,0x171EC4B4,0xBE384454
 458 data4 0x4AA72766,0xBE694741,0xD42518F8,0xBE5D32B6
 459 data4 0x3A319149,0x3E68D96D,0x62415F36,0xBE68F4DA
 460 data4 0xC9C86A3B,0xBE6DDA2F,0xF49228FE,0x3E6B2E50
 461 data4 0x1188B886,0xBE49C0C2,0x1A4C2F1F,0x3E64BFC2
 462 data4 0x2CB98B54,0xBE6A2FBB,0x9A55D329,0x3E5DC5DE
 463 data4 0x39A7AACE,0x3E696490,0x5C66DBA5,0x3E54728B
 464 data4 0xBA1C7D7D,0xBE62B0DB,0x09F1AF5F,0x3E576E04
 465 data4 0x1A0DD6A1,0x3E612500,0x795FBDEF,0xBE66A419
 466 data4 0xE1BD41FC,0xBE5CDE8C,0xEA54964F,0xBE621376
 467 data4 0x476E76EE,0x3E6370BE,0x3427EB92,0x3E390D1A
 468 data4 0x2BF82BF8,0x3E1336DE,0xD0F7BD9E,0xBE5FF1CB
 469 data4 0x0CEB09DD,0xBE60A355,0x0980F30D,0xBE5CA37E
 470 data4 0x4C082D25,0xBE5C541B,0x3B467D29,0xBE5BBECA
 471 data4 0xB9D946C5,0xBE400D8A,0x07ED374A,0xBE5E2A08
 472 data4 0x365C8B0A,0xBE66CB28,0xD3403BCA,0x3E3AAD5B
 473 data4 0xC7EA21E0,0x3E526055,0xE72880D6,0xBE442C75
 474 data4 0x85222A43,0x3E58B2BB,0x522C42BF,0xBE5AAB79
 475 data4 0x469DC2BC,0xBE605CB4,0xA48C40DC,0xBE589FA7
 476 data4 0x1AA42614,0xBE51C214,0xC37293F4,0xBE48D087
 477 data4 0xA2D673E0,0x3E367A1C,0x114F7A38,0xBE51BEBB
 478 data4 0x661A4B48,0xBE6348E5,0x1D3B9962,0xBDF52643
 479 data4 0x35A78A53,0x3E3A3B5E,0x1CECD788,0xBE46C46C
 480 data4 0x7857D689,0xBE60B7EC,0xD14F1AD7,0xBE594D3D
 481 data4 0x4C9A8F60,0xBE4F9C30,0x02DFF9D2,0xBE521873
 482 data4 0x55E6D68F,0xBE5E4C88,0x667F3DC4,0xBE62140F
 483 data4 0x3BF88747,0xBE36961B,0xC96EC6AA,0x3E602861
 484 data4 0xD57FD718,0xBE3B5151,0xFC4A627B,0x3E561CD0
 485 data4 0xCA913FEA,0xBE3A5217,0x9A5D193A,0x3E40A3CC
 486 data4 0x10A9C312,0xBE5AB713,0xC5F57719,0x3E4FDADB
 487 data4 0xDBDF59D5,0x3E361428,0x61B4180D,0x3E5DB5DB
 488 data4 0x7408D856,0xBE42AD5F,0x31B2B707,0x3E2A3148
 489 ASM_SIZE_DIRECTIVE(Constants_exp_64_W1)
 490
 491 .align 64
 492 Constants_exp_64_W2:
 493 ASM_TYPE_DIRECTIVE(Constants_exp_64_W2,@object)
 494 data4 0x00000000,0x00000000,0x37A3D7A2,0xBE641F25
 495 data4 0xAD028C40,0xBE68DD57,0xF212B1B6,0xBE5C77D8
 496 data4 0x1BA5B070,0x3E57878F,0x2ECAE6FE,0xBE55A36A
 497 data4 0x569DFA3B,0xBE620608,0xA6D300A3,0xBE53B50E
 498 data4 0x223F8F2C,0x3E5B5EF2,0xD6DE0DF4,0xBE56A0D9
 499 data4 0xEAE28F51,0xBE64EEF3,0x367EA80B,0xBE5E5AE2
 500 data4 0x5FCBC02D,0x3E47CB1A,0x9BDAFEB7,0xBE656BA0
 501 data4 0x805AFEE7,0x3E6E70C6,0xA3415EBA,0xBE6E0509
 502 data4 0x49BFF529,0xBE56856B,0x00508651,0x3E66DD33
 503 data4 0xC114BC13,0x3E51165F,0xC453290F,0x3E53333D
 504 data4 0x05539FDA,0x3E6A072B,0x7C0A7696,0xBE47CD87
 505 data4 0xEB05C6D9,0xBE668BF4,0x6AE86C93,0xBE67C3E3
 506 data4 0xD0B3E84B,0xBE533904,0x556B53CE,0x3E63E8D9
 507 data4 0x63A98DC8,0x3E212C89,0x032A7A22,0xBE33138F
 508 data4 0xBC584008,0x3E530FA9,0xCCB93C97,0xBE6ADF82
 509 data4 0x8370EA39,0x3E5F9113,0xFB6A05D8,0x3E5443A4
 510 data4 0x181FEE7A,0x3E63DACD,0xF0F67DEC,0xBE62B29D
 511 data4 0x3DDE6307,0x3E65C483,0xD40A24C1,0x3E5BF030
 512 data4 0x14E437BE,0x3E658B8F,0xED98B6C7,0xBE631C29
 513 data4 0x04CF7C71,0x3E6335D2,0xE954A79D,0x3E529EED
 514 data4 0xF64A2FB8,0x3E5D9257,0x854ED06C,0xBE6BED1B
 515 data4 0xD71405CB,0x3E5096F6,0xACB9FDF5,0xBE3D4893
 516 data4 0x01B68349,0xBDFEB158,0xC6A463B9,0x3E628D35
 517 data4 0xADE45917,0xBE559725,0x042FC476,0xBE68C29C
 518 data4 0x01E511FA,0xBE67593B,0x398801ED,0xBE4A4313
 519 data4 0xDA7C3300,0x3E699571,0x08062A9E,0x3E5349BE
 520 data4 0x755BB28E,0x3E5229C4,0x77A1F80D,0x3E67E426
 521 data4 0x6B69C352,0xBE52B33F,0x084DA57F,0xBE6B3550
 522 data4 0xD1D09A20,0xBE6DB03F,0x2161B2C1,0xBE60CBC4
 523 data4 0x78A2B771,0x3E56ED9C,0x9D0FA795,0xBE508E31
 524 data4 0xFD1A54E9,0xBE59482A,0xB07FD23E,0xBE2A17CE
 525 data4 0x17365712,0x3E68BF5C,0xB3785569,0x3E3956F9
 526 ASM_SIZE_DIRECTIVE(Constants_exp_64_W2)
 527
 528 GR_SAVE_PFS         = r59
 529 GR_SAVE_B0          = r60
 530 GR_SAVE_GP          = r61
 531 GR_Parameter_X      = r62
 532 GR_Parameter_Y      = r63
 533 GR_Parameter_RESULT = r64
 534 GR_Parameter_TAG    = r65
 535
 536 FR_X                = f9
 537 FR_Y                = f9
 538 FR_RESULT           = f99
 539
 540 .section .text
 541 .proc expm1l#
 542 .global expm1l#
 543 .align 64
 544 expm1l:
 545 #ifdef _LIBC
 546 .global __expm1l#
 547 __expm1l:
 548 #endif
 549 { .mii
 550 alloc r32 = ar.pfs,0,30,4,0
 551 (p0)  add r33 = 1, r0
 552 (p0)  cmp.eq.unc  p7, p0 =  r0, r0
 553 }
 554 { .mbb
 555         nop.m 999
 556 (p0)  br.cond.sptk exp_continue
 557         nop.b 999 ;;
 558 }
 559
 560 //
 561 //    Set p7 true for expm1
 562 //    Set Flag = r33 = 1 for expm1
 563 //
 564
 565 .endp expm1l
 566 ASM_SIZE_DIRECTIVE(expm1l)
 567
 568 #ifdef _LIBC
 569 libm_hidden_def (__expm1l)
 570 #endif
 571
 572 .section .text
 573 .proc expl#
 574 .global expl#
 575 .align 64
 576 expl:
 577 #ifdef _LIBC
 578 .global __ieee754_expl#
 579 __ieee754_expl:
 580 #endif
 581 { .mii
 582 alloc r32 = ar.pfs,0,30,4,0
 583 (p0)  add r33 = r0, r0
 584 (p0)  cmp.eq.unc  p0, p7 =  r0, r0 ;;
 585 }
 586 exp_continue:
 587 { .mfi
 588 (p0)  add r32 = 2,r0
 589 (p0)  fnorm.s1 f9 = f8
 590       nop.i 0
 591 }
 592 { .mfi
 593 (p0)  nop.m 0
 594 //
 595 //    Set p7 false for exp
 596 //    Set Flag = r33 = 0 for exp
 597 //
 598 (p0)  fclass.m.unc p6, p8 =  f8, 0x1E7
 599       nop.i 0;;
 600 }
 601 { .mfi
 602         nop.m 999
 603 (p0)  fclass.nm.unc p9, p0 =  f8, 0x1FF
 604       nop.i 0
 605 }
 606 { .mfi
 607         nop.m 999
 608 (p0)  mov f36 = f1
 609         nop.i 999 ;;
 610 }
 611 { .mfb
 612         nop.m 999
 613 //
 614 //    Identify NatVals, NaNs, Infs, and Zeros.
 615 //    Identify EM unsupporteds.
 616 //    Save special input registers
 617 (p0)  mov f32 = f0
 618 //
 619 //    Create FR_X_cor      = 0.0
 620 //           GR_Flag       = 0
 621 //           GR_Expo_Range = 2 (r32) for double-extended precision
 622 //           FR_Scale      = 1.0
 623 //
 624 (p6)  br.cond.spnt EXPL_64_SPECIAL ;;
 625 }
 626 { .mib
 627         nop.m 999
 628         nop.i 999
 629 (p9)  br.cond.spnt EXPL_64_UNSUPPORTED ;;
 630 }
 631 { .mfi
 632 (p0)  cmp.ne.unc p12, p13 = 0x01, r33
 633 //
 634 //    Branch out for special input values
 635 //
 636 (p0)  fcmp.lt.unc.s0 p9,p0 =  f8, f0
 637 (p0)  cmp.eq.unc  p15, p0 =  r0, r0
 638 }
 639 { .mmi
 640         nop.m 999
 641 //
 642 //    Raise possible denormal operand exception
 643 //    Normalize x
 644 //
 645 //    This function computes expl( x  + x_cor)
 646 //    Input  FR 1: FR_X
 647 //    Input  FR 2: FR_X_cor
 648 //    Input  GR 1: GR_Flag
 649 //    Input  GR 2: GR_Expo_Range
 650 //    Output FR 3: FR_Y_hi
 651 //    Output FR 4: FR_Y_lo
 652 //    Output FR 5: FR_Scale
 653 //    Output PR 1: PR_Safe
 654 (p0)  addl r34 = @ltoff(Constants_exp_64_Arg#),gp
 655 (p0)  addl r40 = @ltoff(Constants_exp_64_W1#),gp
 656 };;
 657 //
 658 //    Prepare to load constants
 659 //    Set Safe = True
 660 //
 661
 662 { .mmi
 663       ld8  r34 = [r34]
 664       ld8  r40 = [r40]
 665 (p0)  addl r41 = @ltoff(Constants_exp_64_W2#),gp
 666 };;
 667
 668 { .mmi
 669 (p0)  ldfe f37 = [r34],16
 670 (p0)  ld8 r41 = [r41] ;;
 671 }
 672
 673 //
 674 //    N = fcvt.fx(float_N)
 675 //    Set p14 if -6 > expo_X
 676 //
 677 //
 678 //    Bias = 0x0FFFF
 679 //    expo_X = expo_X and Mask
 680 //
 681
 682 { .mmi
 683 (p0)  ldfe f40 = [r34],16
 684       nop.m 999
 685 //
 686 //    Load L_lo
 687 //    Set p10 if 14 < expo_X
 688 //
 689 (p0)  addl r50 = @ltoff(Constants_exp_64_T1#),gp
 690 }
 691 { .mmi
 692         nop.m 999
 693         nop.m 999
 694 (p0)  addl r51 = @ltoff(Constants_exp_64_T2#),gp ;;
 695 }
 696 //
 697 //    Load W2_ptr
 698 //    Branch to SMALL is expo_X < -6
 699 //
 700
 701 {.mmi
 702 (p0)  ld8 r50 = [r50]
 703 (p0)  ld8 r51 = [r51]
 704 };;
 705
 706 { .mlx
 707 (p0)  ldfe f41 = [r34],16
 708 //
 709 //    float_N = X * L_Inv
 710 //    expo_X = exponent of X
 711 //    Mask = 0x1FFFF
 712 //
 713 (p0)  movl r58 = 0x0FFFF
 714 }
 715 { .mlx
 716         nop.m 999
 717 (p0)  movl r39 = 0x1FFFF ;;
 718 }
 719 { .mmi
 720 (p0)  getf.exp r37 = f9
 721         nop.m 999
 722 (p0)  addl r34 = @ltoff(Constants_exp_64_Exponents#),gp ;;
 723 }
 724 { .mii
 725 (p0)  ld8 r34 = [r34]
 726       nop.i 999
 727 (p0)  and  r37 = r37, r39 ;;
 728 }
 729 { .mmi
 730 (p0)  sub r37 = r37, r58 ;;
 731 (p0)  cmp.gt.unc  p14, p0 =  -6, r37
 732 (p0)  cmp.lt.unc  p10, p0 =  14, r37 ;;
 733 }
 734 { .mfi
 735 (p0)  nop.m 0
 736 //
 737 //    Load L_inv
 738 //    Set p12 true for Flag = 0 (exp)
 739 //    Set p13 true for Flag = 1 (expm1)
 740 //
 741 (p0)  fmpy.s1 f38 = f9, f37
 742         nop.i 999 ;;
 743 }
 744 { .mfb
 745         nop.m 999
 746 //
 747 //    Load L_hi
 748 //    expo_X = expo_X - Bias
 749 //    get W1_ptr
 750 //
 751 (p0)  fcvt.fx.s1 f39 = f38
 752 (p14) br.cond.spnt EXPL_SMALL ;;
 753 }
 754 { .mib
 755         nop.m 999
 756         nop.i 999
 757 (p10) br.cond.spnt EXPL_HUGE ;;
 758 }
 759 { .mmi
 760 (p0)  shladd r34 = r32,4,r34
 761       nop.m 999
 762 (p0)  addl r35 = @ltoff(Constants_exp_64_A#),gp ;;
 763 }
 764 //
 765 //    Load T_1,T_2
 766 //
 767 { .mmi
 768    nop.m 999
 769    ld8   r35 =[r35]
 770    nop.i 99
 771 };;
 772 { .mmb
 773 (p0)  ldfe f51 = [r35],16
 774 (p0)  ld8 r45 = [r34],8
 775         nop.b 999 ;;
 776 }
 777 //
 778 //    Set Safe = True  if k >= big_expo_neg
 779 //    Set Safe = False if k < big_expo_neg
 780 //
 781 { .mmb
 782 (p0)  ldfe f49 = [r35],16
 783 (p0)  ld8 r48 = [r34],0
 784         nop.b 999 ;;
 785 }
 786 { .mfi
 787         nop.m 999
 788 //
 789 //    Branch to HUGE is expo_X > 14
 790 //
 791 (p0)  fcvt.xf f38 = f39
 792         nop.i 999 ;;
 793 }
 794 { .mfi
 795 (p0)  getf.sig r52 = f39
 796         nop.f 999
 797         nop.i 999 ;;
 798 }
 799 { .mii
 800         nop.m 999
 801 (p0)  extr.u r43 = r52, 6, 6 ;;
 802 //
 803 //    r = r - float_N * L_lo
 804 //    K = extr(N_fix,12,52)
 805 //
 806 (p0)  shladd r40 = r43,3,r40 ;;
 807 }
 808 { .mfi
 809 (p0)  shladd r50 = r43,2,r50
 810 (p0)  fnma.s1 f42 = f40, f38, f9
 811 //
 812 //    float_N = float(N)
 813 //    N_fix = signficand N
 814 //
 815 (p0)  extr.u r42 = r52, 0, 6
 816 }
 817 { .mmi
 818 (p0)  ldfd  f43 = [r40],0 ;;
 819 (p0)  shladd r41 = r42,3,r41
 820 (p0)  shladd r51 = r42,2,r51
 821 }
 822 //
 823 //    W_1_p1 = 1 + W_1
 824 //
 825 { .mmi
 826 (p0)  ldfs  f44 = [r50],0 ;;
 827 (p0)  ldfd  f45 = [r41],0
 828 //
 829 //    M_2 = extr(N_fix,0,6)
 830 //    M_1 = extr(N_fix,6,6)
 831 //    r = X - float_N * L_hi
 832 //
 833 (p0)  extr r44 = r52, 12, 52
 834 }
 835 { .mmi
 836 (p0)  ldfs  f46 = [r51],0 ;;
 837 (p0)  sub r46 = r58, r44
 838 (p0)  cmp.gt.unc  p8, p15 =  r44, r45
 839 }
 840 //
 841 //    W = W_1 + W_1_p1*W_2
 842 //    Load  A_2
 843 //    Bias_m_K = Bias - K
 844 //
 845 { .mii
 846 (p0)  ldfe f40 = [r35],16
 847 //
 848 //    load A_1
 849 //    poly = A_2 + r*A_3
 850 //    rsq = r * r
 851 //    neg_2_mK = exponent of Bias_m_k
 852 //
 853 (p0)  add r47 = r58, r44 ;;
 854 //
 855 //    Set Safe = True  if k <= big_expo_pos
 856 //    Set Safe = False  if k >  big_expo_pos
 857 //    Load A_3
 858 //
 859 (p15) cmp.lt p8,p15 = r44,r48 ;;
 860 }
 861 { .mmf
 862 (p0)  setf.exp f61 = r46
 863 //
 864 //    Bias_p + K = Bias + K
 865 //    T = T_1 * T_2
 866 //
 867 (p0)  setf.exp f36 = r47
 868 (p0)  fnma.s1 f42 = f41, f38, f42 ;;
 869 }
 870 { .mfi
 871         nop.m 999
 872 //
 873 //    Load W_1,W_2
 874 //    Load big_exp_pos, load big_exp_neg
 875 //
 876 (p0)  fadd.s1 f47 = f43, f1
 877         nop.i 999 ;;
 878 }
 879 { .mfi
 880         nop.m 999
 881 (p0)  fma.s1 f52 = f42, f51, f49
 882         nop.i 999
 883 }
 884 { .mfi
 885         nop.m 999
 886 (p0)  fmpy.s1 f48 = f42, f42
 887         nop.i 999 ;;
 888 }
 889 { .mfi
 890         nop.m 999
 891 (p0)  fmpy.s1 f53 = f44, f46
 892         nop.i 999 ;;
 893 }
 894 { .mfi
 895         nop.m 999
 896 (p0)  fma.s1 f54 = f45, f47, f43
 897         nop.i 999
 898 }
 899 { .mfi
 900         nop.m 999
 901 (p0)  fneg f61 =  f61
 902         nop.i 999 ;;
 903 }
 904 { .mfi
 905         nop.m 999
 906 (p0)  fma.s1 f52 = f42, f52, f40
 907         nop.i 999 ;;
 908 }
 909 { .mfi
 910         nop.m 999
 911 (p0)  fadd.s1 f55 = f54, f1
 912         nop.i 999
 913 }
 914 { .mfi
 915         nop.m 999
 916 //
 917 //    W + Wp1 * poly
 918 //
 919 (p0)  mov f34 = f53
 920         nop.i 999 ;;
 921 }
 922 { .mfi
 923         nop.m 999
 924 //
 925 //    A_1 + r * poly
 926 //    Scale = setf_expl(Bias_p_k)
 927 //
 928 (p0)  fma.s1 f52 = f48, f52, f42
 929         nop.i 999 ;;
 930 }
 931 { .mfi
 932         nop.m 999
 933 //
 934 //    poly = r + rsq(A_1 + r*poly)
 935 //    Wp1 = 1 + W
 936 //    neg_2_mK = -neg_2_mK
 937 //
 938 (p0)  fma.s1 f35 = f55, f52, f54
 939         nop.i 999 ;;
 940 }
 941 { .mfb
 942         nop.m 999
 943 (p0)  fmpy.s1 f35 = f35, f53
 944 //
 945 //    Y_hi = T
 946 //    Y_lo = T * (W + Wp1*poly)
 947 //
 948 (p12) br.cond.sptk EXPL_MAIN ;;
 949 }
 950 //
 951 //    Branch if expl(x)
 952 //    Continue for expl(x-1)
 953 //
 954 { .mii
 955 (p0)  cmp.lt.unc  p12, p13 =  10, r44
 956         nop.i 999 ;;
 957 //
 958 //    Set p12 if 10 < K, Else p13
 959 //
 960 (p13) cmp.gt.unc  p13, p14 =  -10, r44 ;;
 961 }
 962 //
 963 //    K > 10:  Y_lo = Y_lo + neg_2_mK
 964 //    K <=10:  Set p13 if -10 > K, Else set p14
 965 //
 966 { .mfi
 967 (p13) cmp.eq  p15, p0 =  r0, r0
 968 (p14) fadd.s1 f34 = f61, f34
 969         nop.i 999 ;;
 970 }
 971 { .mfi
 972         nop.m 999
 973 (p12) fadd.s1 f35 = f35, f61
 974         nop.i 999 ;;
 975 }
 976 { .mfi
 977         nop.m 999
 978 (p13) fadd.s1 f35 = f35, f34
 979         nop.i 999
 980 }
 981 { .mfb
 982         nop.m 999
 983 //
 984 //    K <= 10 and K < -10, Set Safe = True
 985 //    K <= 10 and K < 10,   Y_lo = Y_hi + Y_lo
 986 //    K <= 10 and K > =-10, Y_hi = Y_hi + neg_2_mk
 987 //
 988 (p13) mov f34 = f61
 989 (p0)  br.cond.sptk EXPL_MAIN ;;
 990 }
 991 EXPL_SMALL:
 992 { .mmi
 993       nop.m 999
 994 (p0)  addl r34 = @ltoff(Constants_exp_64_Exponents#),gp
 995 (p12) addl r35 = @ltoff(Constants_exp_64_P#),gp ;;
 996 }
 997 .pred.rel "mutex",p12,p13
 998 { .mmi
 999 (p12) ld8  r35=[r35]
1000 nop.m 999
1001 (p13) addl r35 = @ltoff(Constants_exp_64_Q#),gp
1002 };;
1003 { .mmi
1004 (p13) ld8  r35=[r35]
1005 (p0) ld8  r34=[r34]
1006 nop.i 999
1007 };;
1008 { .mfi
1009 (p0)  add r34 = 0x48,r34
1010 //
1011 //    Return
1012 //    K <= 10 and K < 10,   Y_hi = neg_2_mk
1013 //
1014 //    /*******************************************************/
1015 //    /*********** Branch EXPL_SMALL  ************************/
1016 //    /*******************************************************/
1017 (p0)  mov f42 = f9
1018         nop.i 999 ;;
1019 }
1020 //
1021 //    Flag = 0
1022 //    r4 = rsq * rsq
1023 //
1024 { .mfi
1025 (p0)  ld8 r49 =[r34],0
1026         nop.f 999
1027         nop.i 999 ;;
1028 }
1029 { .mii
1030         nop.m 999
1031         nop.i 999 ;;
1032 //
1033 //    Flag = 1
1034 //
1035 (p0)  cmp.lt.unc  p14, p0 =  r37, r49 ;;
1036 }
1037 { .mfi
1038         nop.m 999
1039 //
1040 //    r = X
1041 //
1042 (p0)  fmpy.s1 f48 = f42, f42
1043         nop.i 999 ;;
1044 }
1045 { .mfb
1046         nop.m 999
1047 //
1048 //    rsq = r * r
1049 //
1050 (p0)  fmpy.s1 f50 = f48, f48
1051 //
1052 //    Is input very small?
1053 //
1054 (p14) br.cond.spnt EXPL_VERY_SMALL ;;
1055 }
1056 //
1057 //    Flag_not1: Y_hi = 1.0
1058 //    Flag is 1: r6 = rsq * r4
1059 //
1060 { .mfi
1061 (p12) ldfe f52 = [r35],16
1062 (p12) mov f34 = f1
1063 (p0)  add r53 = 0x1,r0 ;;
1064 }
1065 { .mfi
1066 (p13) ldfe f51 = [r35],16
1067 //
1068 //    Flag_not_1: Y_lo = poly_hi + r4 * poly_lo
1069 //
1070 (p13) mov f34 = f9
1071         nop.i 999 ;;
1072 }
1073 { .mmf
1074 (p12) ldfe f53 = [r35],16
1075 //
1076 //    For Flag_not_1, Y_hi = X
1077 //    Scale = 1
1078 //    Create 0x000...01
1079 //
1080 (p0)  setf.sig f37 = r53
1081 (p0)  mov f36 = f1 ;;
1082 }
1083 { .mmi
1084 (p13) ldfe f52 = [r35],16 ;;
1085 (p12) ldfe f54 = [r35],16
1086         nop.i 999 ;;
1087 }
1088 { .mfi
1089 (p13) ldfe f53 = [r35],16
1090 (p13) fmpy.s1 f58 = f48, f50
1091         nop.i 999 ;;
1092 }
1093 //
1094 //    Flag_not1: poly_lo = P_5 + r*P_6
1095 //    Flag_1: poly_lo = Q_6 + r*Q_7
1096 //
1097 { .mmi
1098 (p13) ldfe f54 = [r35],16 ;;
1099 (p12) ldfe f55 = [r35],16
1100         nop.i 999 ;;
1101 }
1102 { .mmi
1103 (p12) ldfe f56 = [r35],16 ;;
1104 (p13) ldfe f55 = [r35],16
1105         nop.i 999 ;;
1106 }
1107 { .mmi
1108 (p12) ldfe f57 = [r35],0 ;;
1109 (p13) ldfe f56 = [r35],16
1110         nop.i 999 ;;
1111 }
1112 { .mfi
1113 (p13) ldfe f57 = [r35],0
1114         nop.f 999
1115         nop.i 999 ;;
1116 }
1117 { .mfi
1118         nop.m 999
1119 //
1120 //    For  Flag_not_1, load p5,p6,p1,p2
1121 //    Else load p5,p6,p1,p2
1122 //
1123 (p12) fma.s1 f60 = f52, f42, f53
1124         nop.i 999 ;;
1125 }
1126 { .mfi
1127         nop.m 999
1128 (p13) fma.s1 f60 = f51, f42, f52
1129         nop.i 999 ;;
1130 }
1131 { .mfi
1132         nop.m 999
1133 (p12) fma.s1 f60 = f60, f42, f54
1134         nop.i 999 ;;
1135 }
1136 { .mfi
1137         nop.m 999
1138 (p12) fma.s1 f59 = f56, f42, f57
1139         nop.i 999 ;;
1140 }
1141 { .mfi
1142         nop.m 999
1143 (p13) fma.s1 f60 = f42, f60, f53
1144         nop.i 999 ;;
1145 }
1146 { .mfi
1147         nop.m 999
1148 (p12) fma.s1 f59 = f59, f48, f42
1149         nop.i 999 ;;
1150 }
1151 { .mfi
1152         nop.m 999
1153 //
1154 //    Flag_1: poly_lo = Q_5 + r*(Q_6 + r*Q_7)
1155 //    Flag_not1: poly_lo = P_4 + r*(P_5 + r*P_6)
1156 //    Flag_not1: poly_hi = (P_1 + r*P_2)
1157 //
1158 (p13) fmpy.s1 f60 = f60, f58
1159         nop.i 999 ;;
1160 }
1161 { .mfi
1162         nop.m 999
1163 (p12) fma.s1 f60 = f60, f42, f55
1164         nop.i 999 ;;
1165 }
1166 { .mfi
1167         nop.m 999
1168 //
1169 //    Flag_1: poly_lo = r6 *(Q_5 + ....)
1170 //    Flag_not1: poly_hi =  r + rsq *(P_1 + r*P_2)
1171 //
1172 (p12) fma.s1 f35 = f60, f50, f59
1173         nop.i 999
1174 }
1175 { .mfi
1176         nop.m 999
1177 (p13) fma.s1 f59 = f54, f42, f55
1178         nop.i 999 ;;
1179 }
1180 { .mfi
1181         nop.m 999
1182 //
1183 //    Flag_not1: Y_lo = rsq* poly_hi + poly_lo
1184 //    Flag_1: poly_lo = rsq* poly_hi + poly_lo
1185 //
1186 (p13) fma.s1 f59 = f59, f42, f56
1187         nop.i 999 ;;
1188 }
1189 { .mfi
1190         nop.m 999
1191 //
1192 //    Flag_not_1: (P_1 + r*P_2)
1193 //
1194 (p13) fma.s1 f59 = f59, f42, f57
1195         nop.i 999 ;;
1196 }
1197 { .mfi
1198         nop.m 999
1199 //
1200 //    Flag_not_1: poly_hi = r + rsq * (P_1 + r*P_2)
1201 //
1202 (p13) fma.s1 f35 = f59, f48, f60
1203         nop.i 999 ;;
1204 }
1205 { .mfi
1206         nop.m 999
1207 //
1208 //    Create 0.000...01
1209 //
1210 (p0)  for f37 = f35, f37
1211         nop.i 999 ;;
1212 }
1213 { .mfb
1214         nop.m 999
1215 //
1216 //    Set lsb of Y_lo to 1
1217 //
1218 (p0)  fmerge.se f35 = f35,f37
1219 (p0)  br.cond.sptk EXPL_MAIN ;;
1220 }
1221 EXPL_VERY_SMALL:
1222 { .mmi
1223         nop.m 999
1224         nop.m 999
1225 (p13) addl r34 = @ltoff(Constants_exp_64_Exponents#),gp
1226 }
1227 { .mfi
1228         nop.m 999
1229 (p12) mov f35 = f9
1230         nop.i 999 ;;
1231 }
1232 { .mfb
1233 (p13) ld8 r34 = [r34]
1234 (p12) mov f34 = f1
1235 (p12) br.cond.sptk EXPL_MAIN ;;
1236 }
1237 { .mlx
1238 (p13) add  r34 = 8,r34
1239 (p13) movl r39 = 0x0FFFE ;;
1240 }
1241 //
1242 //    Load big_exp_neg
1243 //    Create 1/2's exponent
1244 //
1245 { .mii
1246 (p13) setf.exp f56 = r39
1247 (p13) shladd r34 = r32,4,r34 ;;
1248         nop.i 999
1249 }
1250 //
1251 //    Negative exponents are stored after positive
1252 //
1253 { .mfi
1254 (p13) ld8 r45 = [r34],0
1255 //
1256 //    Y_hi = x
1257 //    Scale = 1
1258 //
1259 (p13) fmpy.s1 f35 = f9, f9
1260         nop.i 999 ;;
1261 }
1262 { .mfi
1263         nop.m 999
1264 //
1265 //    Reset Safe if necessary
1266 //    Create 1/2
1267 //
1268 (p13) mov f34 = f9
1269         nop.i 999 ;;
1270 }
1271 { .mfi
1272 (p13) cmp.lt.unc  p0, p15 =  r37, r45
1273 (p13) mov f36 = f1
1274         nop.i 999 ;;
1275 }
1276 { .mfb
1277         nop.m 999
1278 //
1279 //    Y_lo = x * x
1280 //
1281 (p13) fmpy.s1 f35 = f35, f56
1282 //
1283 //    Y_lo = x*x/2
1284 //
1285 (p13) br.cond.sptk EXPL_MAIN ;;
1286 }
1287 EXPL_HUGE:
1288 { .mfi
1289         nop.m 999
1290 (p0)  fcmp.gt.unc.s1 p14, p0 =  f9, f0
1291         nop.i 999
1292 }
1293 { .mlx
1294         nop.m 999
1295 (p0)  movl r39 = 0x15DC0 ;;
1296 }
1297 { .mfi
1298 (p14) setf.exp f34 = r39
1299 (p14) mov f35 = f1
1300 (p14) cmp.eq  p0, p15 =  r0, r0 ;;
1301 }
1302 { .mfb
1303         nop.m 999
1304 (p14) mov f36 = f34
1305 //
1306 //    If x > 0, Set Safe = False
1307 //    If x > 0, Y_hi = 2**(24,000)
1308 //    If x > 0, Y_lo = 1.0
1309 //    If x > 0, Scale = 2**(24,000)
1310 //
1311 (p14) br.cond.sptk EXPL_MAIN ;;
1312 }
1313 { .mlx
1314         nop.m 999
1315 (p12) movl r39 = 0xA240
1316 }
1317 { .mlx
1318         nop.m 999
1319 (p12) movl r38 = 0xA1DC ;;
1320 }
1321 { .mmb
1322 (p13) cmp.eq  p15, p14 =  r0, r0
1323 (p12) setf.exp f34 = r39
1324         nop.b 999 ;;
1325 }
1326 { .mlx
1327 (p12) setf.exp f35 = r38
1328 (p13) movl r39 = 0xFF9C
1329 }
1330 { .mfi
1331         nop.m 999
1332 (p13) fsub.s1 f34 = f0, f1
1333         nop.i 999 ;;
1334 }
1335 { .mfi
1336         nop.m 999
1337 (p12) mov f36 = f34
1338 (p12) cmp.eq  p0, p15 =  r0, r0 ;;
1339 }
1340 { .mfi
1341 (p13) setf.exp f35 = r39
1342 (p13) mov f36 = f1
1343         nop.i 999 ;;
1344 }
1345 EXPL_MAIN:
1346 { .mfi
1347 (p0)  cmp.ne.unc p12, p0 = 0x01, r33
1348 (p0)  fmpy.s1 f101 = f36, f35
1349         nop.i 999 ;;
1350 }
1351 { .mfb
1352         nop.m 999
1353 (p0)  fma.s0 f99 = f34, f36, f101
1354 (p15) br.cond.sptk EXPL_64_RETURN ;;
1355 }
1356 { .mfi
1357         nop.m 999
1358 (p0)  fsetc.s3 0x7F,0x01
1359         nop.i 999
1360 }
1361 { .mlx
1362         nop.m 999
1363 (p0)  movl r50 = 0x00000000013FFF ;;
1364 }
1365 //
1366 //    S0 user supplied status
1367 //    S2 user supplied status + WRE + TD  (Overflows)
1368 //    S3 user supplied status + RZ + TD   (Underflows)
1369 //
1370 //
1371 //    If (Safe) is true, then
1372 //        Compute result using user supplied status field.
1373 //        No overflow or underflow here, but perhaps inexact.
1374 //        Return
1375 //    Else
1376 //       Determine if overflow or underflow  was raised.
1377 //       Fetch +/- overflow threshold for IEEE single, double,
1378 //       double extended
1379 //
1380 { .mfi
1381 (p0)  setf.exp f60 = r50
1382 (p0)  fma.s3 f102 = f34, f36, f101
1383         nop.i 999
1384 }
1385 { .mfi
1386         nop.m 999
1387 (p0)  fsetc.s3 0x7F,0x40
1388         nop.i 999 ;;
1389 }
1390 { .mfi
1391         nop.m 999
1392 //
1393 //    For Safe, no need to check for over/under.
1394 //    For expm1, handle errors like exp.
1395 //
1396 (p0)  fsetc.s2 0x7F,0x42
1397         nop.i 999;;
1398 }
1399 { .mfi
1400         nop.m 999
1401 (p0)  fma.s2 f100 = f34, f36, f101
1402         nop.i 999 ;;
1403 }
1404 { .mfi
1405         nop.m 999
1406 (p0)  fsetc.s2 0x7F,0x40
1407         nop.i 999 ;;
1408 }
1409 { .mfi
1410         nop.m 999
1411 (p7)  fclass.m.unc   p12, p0 =  f102, 0x00F
1412         nop.i 999
1413 }
1414 { .mfi
1415         nop.m 999
1416 (p0)  fclass.m.unc   p11, p0 =  f102, 0x00F
1417         nop.i 999 ;;
1418 }
1419 { .mfi
1420         nop.m 999
1421 (p7)  fcmp.ge.unc.s1 p10, p0 =  f100, f60
1422         nop.i 999
1423 }
1424 { .mfi
1425         nop.m 999
1426 //
1427 //    Create largest double exponent + 1.
1428 //    Create smallest double exponent - 1.
1429 //
1430 (p0)  fcmp.ge.unc.s1 p8, p0 =  f100, f60
1431         nop.i 999 ;;
1432 }
1433 //
1434 //    fcmp:   resultS2 >= + overflow threshold  -> set (a) if true
1435 //    fcmp:   resultS2 <= - overflow threshold  -> set (b) if true
1436 //    fclass: resultS3 is denorm/unorm/0        -> set (d) if true
1437 //
1438 { .mib
1439 (p10) mov   GR_Parameter_TAG = 39
1440         nop.i 999
1441 (p10) br.cond.sptk __libm_error_region ;;
1442 }
1443 { .mib
1444 (p8)  mov   GR_Parameter_TAG = 12
1445         nop.i 999
1446 (p8)  br.cond.sptk __libm_error_region ;;
1447 }
1448 //
1449 //    Report that exp overflowed
1450 //
1451 { .mib
1452 (p12) mov   GR_Parameter_TAG = 40
1453         nop.i 999
1454 (p12) br.cond.sptk __libm_error_region ;;
1455 }
1456 { .mib
1457 (p11) mov   GR_Parameter_TAG = 13
1458         nop.i 999
1459 (p11) br.cond.sptk __libm_error_region ;;
1460 }
1461 { .mib
1462         nop.m 999
1463         nop.i 999
1464 //
1465 //    Report that exp underflowed
1466 //
1467 (p0)  br.cond.sptk EXPL_64_RETURN ;;
1468 }
1469 EXPL_64_SPECIAL:
1470 { .mfi
1471         nop.m 999
1472 (p0)  fclass.m.unc p6,  p0 =  f8, 0x0c3
1473         nop.i 999
1474 }
1475 { .mfi
1476         nop.m 999
1477 (p0)  fclass.m.unc p13, p8 =  f8, 0x007
1478         nop.i 999 ;;
1479 }
1480 { .mfi
1481         nop.m 999
1482 (p7)  fclass.m.unc p14, p0 =  f8, 0x007
1483         nop.i 999
1484 }
1485 { .mfi
1486         nop.m 999
1487 (p0)  fclass.m.unc p12, p9 =  f8, 0x021
1488         nop.i 999 ;;
1489 }
1490 { .mfi
1491         nop.m 999
1492 (p0)  fclass.m.unc p11, p0 =  f8, 0x022
1493         nop.i 999
1494 }
1495 { .mfi
1496         nop.m 999
1497 (p7)  fclass.m.unc p10, p0 =  f8, 0x022
1498         nop.i 999 ;;
1499 }
1500 { .mfi
1501         nop.m 999
1502 //
1503 //    Identify +/- 0, Inf, or -Inf
1504 //    Generate the right kind of NaN.
1505 //
1506 (p13) fadd.s0 f99 = f0, f1
1507         nop.i 999 ;;
1508 }
1509 { .mfi
1510         nop.m 999
1511 (p14) mov f99 = f8
1512         nop.i 999 ;;
1513 }
1514 { .mfb
1515         nop.m 999
1516 (p6)  fadd.s0 f99 = f8, f1
1517 //
1518 //    expl(+/-0) = 1
1519 //    expm1l(+/-0) = +/-0
1520 //    No exceptions raised
1521 //
1522 (p6)  br.cond.sptk EXPL_64_RETURN ;;
1523 }
1524 { .mib
1525         nop.m 999
1526         nop.i 999
1527 (p14) br.cond.sptk EXPL_64_RETURN ;;
1528 }
1529 { .mfi
1530         nop.m 999
1531 (p11) mov f99 = f0
1532         nop.i 999 ;;
1533 }
1534 { .mfb
1535         nop.m 999
1536 (p10) fsub.s1 f99 = f0, f1
1537 //
1538 //    expl(-Inf) = 0
1539 //    expm1l(-Inf) = -1
1540 //    No exceptions raised.
1541 //
1542 (p10) br.cond.sptk EXPL_64_RETURN ;;
1543 }
1544 { .mfb
1545         nop.m 999
1546 (p12) fmpy.s1 f99 = f8, f1
1547 //
1548 //    expl(+Inf) = Inf
1549 //    No exceptions raised.
1550 //
1551 (p0)  br.cond.sptk EXPL_64_RETURN ;;
1552 }
1553 EXPL_64_UNSUPPORTED:
1554 { .mfb
1555         nop.m 999
1556 (p0)  fmpy.s0 f99 = f8, f0
1557 (p0)  br.cond.sptk EXPL_64_RETURN ;;
1558 }
1559 EXPL_64_RETURN:
1560 { .mfb
1561       nop.m 999
1562 (p0)  mov   f8     = f99
1563 (p0)  br.ret.sptk   b0
1564 }
1565 .endp
1566 ASM_SIZE_DIRECTIVE(expl)
1567
1568 .proc __libm_error_region
1569 __libm_error_region:
1570 .prologue
1571 { .mfi
1572         add   GR_Parameter_Y=-32,sp             // Parameter 2 value
1573         nop.f 0
1574 .save   ar.pfs,GR_SAVE_PFS
1575         mov  GR_SAVE_PFS=ar.pfs                 // Save ar.pfs
1576 }
1577 { .mfi
1578 .fframe 64
1579         add sp=-64,sp                           // Create new stack
1580         nop.f 0
1581         mov GR_SAVE_GP=gp                       // Save gp
1582 };;
1583 { .mmi
1584         stfe [GR_Parameter_Y] = FR_Y,16         // Save Parameter 2 on stack
1585         add GR_Parameter_X = 16,sp              // Parameter 1 address
1586 .save   b0, GR_SAVE_B0
1587         mov GR_SAVE_B0=b0                       // Save b0
1588 };;
1589 .body
1590 { .mib
1591         stfe [GR_Parameter_X] = FR_X            // Store Parameter 1 on stack
1592         add   GR_Parameter_RESULT = 0,GR_Parameter_Y
1593         nop.b 0                                 // Parameter 3 address
1594 }
1595 { .mib
1596         stfe [GR_Parameter_Y] = FR_RESULT      // Store Parameter 3 on stack
1597         add   GR_Parameter_Y = -16,GR_Parameter_Y
1598         br.call.sptk b0=__libm_error_support#  // Call error handling function
1599 };;
1600 { .mmi
1601         nop.m 0
1602         nop.m 0
1603         add   GR_Parameter_RESULT = 48,sp
1604 };;
1605 { .mmi
1606         ldfe  f8 = [GR_Parameter_RESULT]       // Get return result off stack
1607 .restore sp
1608         add   sp = 64,sp                       // Restore stack pointer
1609         mov   b0 = GR_SAVE_B0                  // Restore return address
1610 };;
1611 { .mib
1612         mov   gp = GR_SAVE_GP                  // Restore gp
1613         mov   ar.pfs = GR_SAVE_PFS             // Restore ar.pfs
1614         br.ret.sptk     b0                     // Return
1615 };;
1616 .endp __libm_error_region
1617 ASM_SIZE_DIRECTIVE(__libm_error_region)
1618
1619 .type   __libm_error_support#,@function
1620 .global __libm_error_support#