sysdeps/ia64/fpu/e_pow.S

   1 .file "pow.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 // 2/02/00  Initial version
  43 // 2/03/00  Added p12 to definite over/under path. With odd power we did not
  44 //          maintain the sign of x in this path.
  45 // 4/04/00  Unwind support added
  46 // 4/19/00  pow(+-1,inf) now returns NaN
  47 //          pow(+-val, +-inf) returns 0 or inf, but now does not call error support
  48 //          Added s1 to fcvt.fx because invalid flag was incorrectly set.
  49 // 8/15/00  Bundle added after call to __libm_error_support to properly
  50 //          set [the previously overwritten] GR_Parameter_RESULT.
  51 // 9/07/00  Improved performance by eliminating bank conflicts and other stalls,
  52 //          and tweaking the critical path
  53 // 9/08/00  Per c99, pow(+-1,inf) now returns 1, and pow(+1,nan) returns 1
  54 // 9/28/00  Updated NaN**0 path
  55 // 1/20/01  Fixed denormal flag settings.
  56 // 2/12/01  Improved speed.
  57 //
  58 // API
  59 //==============================================================
  60 // double pow(double)
  61 // float  powf(float)
  62 //
  63 // Overview of operation
  64 //==============================================================
  65 //
  66 // Three steps...
  67 // 1. Log(x)
  68 // 2. y Log(x)
  69 // 3. exp(y log(x))
  70 //
  71 // This means we work with the absolute value of x and merge in the sign later.
  72 //      Log(x) = G + delta + r -rsq/2 + p
  73 // G,delta depend on the exponent of x and table entries. The table entries are
  74 // indexed by the exponent of x, called K.
  75 //
  76 // The G and delta come out of the reduction; r is the reduced x.
  77 //
  78 // B = frcpa(x)
  79 // xB-1 is small means that B is the approximate inverse of x.
  80 //
  81 //      Log(x) = Log( (1/B)(Bx) )
  82 //             = Log(1/B) + Log(Bx)
  83 //             = Log(1/B) + Log( 1 + (Bx-1))
  84 //
  85 //      x  = 2^K 1.x_1x_2.....x_52
  86 //      B= frcpa(x) = 2^-k Cm
  87 //      Log(1/B) = Log(1/(2^-K Cm))
  88 //      Log(1/B) = Log((2^K/ Cm))
  89 //      Log(1/B) = K Log(2) + Log(1/Cm)
  90 //
  91 //      Log(x)   = K Log(2) + Log(1/Cm) + Log( 1 + (Bx-1))
  92 //
  93 // If you take the significand of x, set the exponent to true 0, then Cm is
  94 // the frcpa. We tabulate the Log(1/Cm) values. There are 256 of them.
  95 // The frcpa table is indexed by 8 bits, the x_1 thru x_8.
  96 // m = x_1x_2...x_8 is an 8-bit index.
  97 //
  98 //      Log(1/Cm) = log(1/frcpa(1+m/256)) where m goes from 0 to 255.
  99 //
 100 // We tabluate as two doubles, T and t, where T +t is the value itself.
 101 //
 102 //      Log(x)   = (K Log(2)_hi + T) + (Log(2)_hi + t) + Log( 1 + (Bx-1))
 103 //      Log(x)   =  G + delta           + Log( 1 + (Bx-1))
 104 //
 105 // The Log( 1 + (Bx-1)) can be calculated as a series in r = Bx-1.
 106 //
 107 //      Log( 1 + (Bx-1)) = r - rsq/2 + p
 108 //
 109 // Then,
 110 //
 111 //      yLog(x) = yG + y delta + y(r-rsq/2) + yp
 112 //      yLog(x) = Z1 + e3      + Z2         + Z3 + (e2 + e3)
 113 //
 114 //
 115 //     exp(yLog(x)) = exp(Z1 + Z2 + Z3) exp(e1 + e2 + e3)
 116 //
 117 //
 118 //       exp(Z3) is another series.
 119 //       exp(e1 + e2 + e3) is approximated as f3 = 1 + (e1 + e2 + e3)
 120 //
 121 //       Z1 (128/log2) = number of log2/128 in Z1 is N1
 122 //       Z2 (128/log2) = number of log2/128 in Z2 is N2
 123 //
 124 //       s1 = Z1 - N1 log2/128
 125 //       s2 = Z2 - N2 log2/128
 126 //
 127 //       s = s1 + s2
 128 //       N = N1 + N2
 129 //
 130 //       exp(Z1 + Z2) = exp(Z)
 131 //       exp(Z)       = exp(s) exp(N log2/128)
 132 //
 133 //       exp(r)       = exp(Z - N log2/128)
 134 //
 135 //      r = s + d = (Z - N (log2/128)_hi) -N (log2/128)_lo
 136 //                =  Z - N (log2/128)
 137 //
 138 //      Z         = s+d +N (log2/128)
 139 //
 140 //      exp(Z)    = exp(s) (1+d) exp(N log2/128)
 141 //
 142 //      N = M 128 + n
 143 //
 144 //      N log2/128 = M log2 + n log2/128
 145 //
 146 //      n is 8 binary digits = n_7n_6...n_1
 147 //
 148 //      n log2/128 = n_7n_6n_5 16 log2/128 + n_4n_3n_2n_1 log2/128
 149 //      n log2/128 = n_7n_6n_5 log2/8 + n_4n_3n_2n_1 log2/128
 150 //      n log2/128 = I2 log2/8 + I1 log2/128
 151 //
 152 //      N log2/128 = M log2 + I2 log2/8 + I1 log2/128
 153 //
 154 //      exp(Z)    = exp(s) (1+d) exp(log(2^M) + log(2^I2/8) + log(2^I1/128))
 155 //      exp(Z)    = exp(s) (1+d1) (1+d2)(2^M) 2^I2/8 2^I1/128
 156 //      exp(Z)    = exp(s) f1 f2 (2^M) 2^I2/8 2^I1/128
 157 //
 158 // I1, I2 are table indices. Use a series for exp(s).
 159 // Then get exp(Z)
 160 //
 161 //     exp(yLog(x)) = exp(Z1 + Z2 + Z3) exp(e1 + e2 + e3)
 162 //     exp(yLog(x)) = exp(Z) exp(Z3) f3
 163 //     exp(yLog(x)) = exp(Z)f3 exp(Z3)
 164 //     exp(yLog(x)) = A exp(Z3)
 165 //
 166 // We actually calculate exp(Z3) -1.
 167 // Then,
 168 //     exp(yLog(x)) = A + A( exp(Z3)   -1)
 169 //
 170
 171 // Table Generation
 172 //==============================================================
 173
 174 // The log values
 175 // ==============
 176 // The operation (K*log2_hi) must be exact. K is the true exponent of x.
 177 // If we allow gradual underflow (denormals), K can be represented in 12 bits
 178 // (as a two's complement number). We assume 13 bits as an engineering precaution.
 179 //
 180 //           +------------+----------------+-+
 181 //           |  13 bits   | 50 bits        | |
 182 //           +------------+----------------+-+
 183 //           0            1                66
 184 //                        2                34
 185 //
 186 // So we want the lsb(log2_hi) to be 2^-50
 187 // We get log2 as a quad-extended (15-bit exponent, 128-bit significand)
 188 //
 189 //      0 fffe b17217f7d1cf79ab c9e3b39803f2f6af (4...)
 190 //
 191 // Consider numbering the bits left to right, starting at 0 thru 127.
 192 // Bit 0 is the 2^-1 bit; bit 49 is the 2^-50 bit.
 193 //
 194 //  ...79ab
 195 //     0111 1001 1010 1011
 196 //     44
 197 //     89
 198 //
 199 // So if we shift off the rightmost 14 bits, then (shift back only
 200 // the top half) we get
 201 //
 202 //      0 fffe b17217f7d1cf4000 e6af278ece600fcb dabc000000000000
 203 //
 204 // Put the right 64-bit signficand in an FR register, convert to double;
 205 // it is exact. Put the next 128 bits into a quad register and round to double.
 206 // The true exponent of the low part is -51.
 207 //
 208 // hi is 0 fffe b17217f7d1cf4000
 209 // lo is 0 ffcc e6af278ece601000
 210 //
 211 // Convert to double memory format and get
 212 //
 213 // hi is 0x3fe62e42fefa39e8
 214 // lo is 0x3cccd5e4f1d9cc02
 215 //
 216 // log2_hi + log2_lo is an accurate value for log2.
 217 //
 218 //
 219 // The T and t values
 220 // ==================
 221 // A similar method is used to generate the T and t values.
 222 //
 223 // K * log2_hi + T  must be exact.
 224 //
 225 // Smallest T,t
 226 // ----------
 227 // The smallest T,t is
 228 //       T                   t
 229 // data8 0x3f60040155d58800, 0x3c93bce0ce3ddd81  log(1/frcpa(1+0/256))=  +1.95503e-003
 230 //
 231 // The exponent is 0x3f6 (biased)  or -9 (true).
 232 // For the smallest T value, what we want is to clip the significand such that
 233 // when it is shifted right by 9, its lsb is in the bit for 2^-51. The 9 is the specific
 234 // for the first entry. In general, it is 0xffff - (biased 15-bit exponent).
 235
 236 // Independently, what we have calculated is the table value as a quad precision number.
 237 // Table entry 1 is
 238 // 0 fff6 80200aaeac44ef38 338f77605fdf8000
 239 //
 240 // We store this quad precision number in a data structure that is
 241 //    sign:           1
 242 //    exponent:      15
 243 //    signficand_hi: 64 (includes explicit bit)
 244 //    signficand_lo: 49
 245 // Because the explicit bit is included, the significand is 113 bits.
 246 //
 247 // Consider significand_hi for table entry 1.
 248 //
 249 //
 250 // +-+--- ... -------+--------------------+
 251 // | |
 252 // +-+--- ... -------+--------------------+
 253 // 0 1               4444444455555555556666
 254 //                   2345678901234567890123
 255 //
 256 // Labeled as above, bit 0 is 2^0, bit 1 is 2^-1, etc.
 257 // Bit 42 is 2^-42. If we shift to the right by 9, the bit in
 258 // bit 42 goes in 51.
 259 //
 260 // So what we want to do is shift bits 43 thru 63 into significand_lo.
 261 // This is shifting bit 42 into bit 63, taking care to retain the shifted-off bits.
 262 // Then shifting (just with signficaand_hi) back into bit 42.
 263 //
 264 // The shift_value is 63-42 = 21. In general, this is
 265 //      63 - (51 -(0xffff - 0xfff6))
 266 // For this example, it is
 267 //      63 - (51 - 9) = 63 - 42  = 21
 268 //
 269 // This means we are shifting 21 bits into significand_lo.  We must maintain more
 270 // that a 128-bit signficand not to lose bits. So before the shift we put the 128-bit
 271 // significand into a 256-bit signficand and then shift.
 272 // The 256-bit significand has four parts: hh, hl, lh, and ll.
 273 //
 274 // Start off with
 275 //      hh         hl         lh         ll
 276 //      <64>       <49><15_0> <64_0>     <64_0>
 277 //
 278 // After shift by 21 (then return for significand_hi),
 279 //      <43><21_0> <21><43>   <6><58_0>  <64_0>
 280 //
 281 // Take the hh part and convert to a double. There is no rounding here.
 282 // The conversion is exact. The true exponent of the high part is the same as the
 283 // true exponent of the input quad.
 284 //
 285 // We have some 64 plus significand bits for the low part. In this example, we have
 286 // 70 bits. We want to round this to a double. Put them in a quad and then do a quad fnorm.
 287 // For this example the true exponent of the low part is
 288 //      true_exponent_of_high - 43 = true_exponent_of_high - (64-21)
 289 // In general, this is
 290 //      true_exponent_of_high - (64 - shift_value)
 291 //
 292 //
 293 // Largest T,t
 294 // ----------
 295 // The largest T,t is
 296 // data8 0x3fe62643fecf9742, 0x3c9e3147684bd37d    log(1/frcpa(1+255/256))=  +6.92171e-001
 297 //
 298 // Table entry 256 is
 299 // 0 fffe b1321ff67cba178c 51da12f4df5a0000
 300 //
 301 // The shift value is
 302 //      63 - (51 -(0xffff - 0xfffe)) = 13
 303 //
 304 // The true exponent of the low part is
 305 //      true_exponent_of_high - (64 - shift_value)
 306 //      -1 - (64-13) = -52
 307 // Biased as a double, this is 0x3cb
 308 //
 309 //
 310 //
 311 // So then lsb(T) must be >= 2^-51
 312 // msb(Klog2_hi) <= 2^12
 313 //
 314 //              +--------+---------+
 315 //              |       51 bits    | <== largest T
 316 //              +--------+---------+
 317 //              | 9 bits | 42 bits | <== smallest T
 318 // +------------+----------------+-+
 319 // |  13 bits   | 50 bits        | |
 320 // +------------+----------------+-+
 321
 322
 323
 324 // Special Cases
 325 //==============================================================
 326
 327 //                                   double     float
 328 // overflow                          error 24   30
 329
 330 // underflow                         error 25   31
 331
 332 // X zero  Y zero
 333 //  +0     +0                 +1     error 26   32
 334 //  -0     +0                 +1     error 26   32
 335 //  +0     -0                 +1     error 26   32
 336 //  -0     -0                 +1     error 26   32
 337
 338 // X zero  Y negative
 339 //  +0     -odd integer       +inf   error 27   33  divide-by-zero
 340 //  -0     -odd integer       -inf   error 27   33  divide-by-zero
 341 //  +0     !-odd integer      +inf   error 27   33  divide-by-zero
 342 //  -0     !-odd integer      +inf   error 27   33  divide-by-zero
 343 //  +0     -inf               +inf   error 27   33  divide-by-zero
 344 //  -0     -inf               +inf   error 27   33  divide-by-zero
 345
 346 // X zero  Y positve
 347 //  +0     +odd integer       +0
 348 //  -0     +odd integer       -0
 349 //  +0     !+odd integer      +0
 350 //  -0     !+odd integer      +0
 351 //  +0     +inf               +0
 352 //  -0     +inf               +0
 353 //  +0     Y NaN              quiet Y               invalid if Y SNaN
 354 //  -0     Y NaN              quiet Y               invalid if Y SNaN
 355
 356 // X one
 357 //  -1     Y inf              +1
 358 //  -1     Y NaN              quiet Y               invalid if Y SNaN
 359 //  +1     Y NaN              +1                    invalid if Y SNaN
 360 //  +1     Y any else         +1
 361
 362 // X -     Y not integer      QNAN   error 28   34  invalid
 363
 364 // X NaN   Y 0                +1     error 29   35
 365 // X NaN   Y NaN              quiet X               invalid if X or Y SNaN
 366 // X NaN   Y any else         quiet X               invalid if X SNaN
 367 // X !+1   Y NaN              quiet Y               invalid if Y SNaN
 368
 369
 370 // X +inf  Y >0               +inf
 371 // X -inf  Y >0, !odd integer +inf
 372 // X -inf  Y >0, odd integer  -inf
 373
 374 // X +inf  Y <0               +0
 375 // X -inf  Y <0, !odd integer +0
 376 // X -inf  Y <0, odd integer  -0
 377
 378 // X +inf  Y =0               +1
 379 // X -inf  Y =0               +1
 380
 381 // |X|<1   Y +inf             +0
 382 // |X|<1   Y -inf             +inf
 383 // |X|>1   Y +inf             +inf
 384 // |X|>1   Y -inf             +0
 385
 386 // X any   Y =0               +1
 387
 388 #include "libm_support.h"
 389
 390 // Assembly macros
 391 //==============================================================
 392
 393 // integer registers used
 394
 395 pow_AD_Tt                 = r33
 396 pow_GR_FFF7               = r34
 397 pow_GR_exp_Y              = r34 // duplicate
 398 pow_GR_17ones             = r35
 399
 400 pow_AD_P                  = r36
 401 pow_AD_Q                  = r37
 402 pow_AD_tbl1               = r38
 403 pow_AD_tbl2               = r39
 404 pow_GR_exp_X              = r40
 405 pow_GR_true_exp_X         = r40 // duplicate
 406
 407 pow_GR_offset             = r41
 408 pow_GR_exp_Xm1            = r42
 409 pow_GR_sig_X              = r43
 410 pow_GR_signexp_X          = r44
 411
 412 pow_GR_signexp_Xm1        = r46
 413 pow_GR_int_W1             = r47
 414 pow_GR_int_W2             = r48
 415 pow_GR_int_N              = r49
 416 pow_GR_index1             = r50
 417
 418 pow_GR_index2             = r51
 419 pow_AD_T1                 = r52
 420 pow_AD_T2                 = r53
 421 pow_GR_gt_ln              = r53 // duplicate
 422 pow_int_GR_M              = r54
 423 pow_GR_10033              = r55
 424
 425 pow_GR_16ones             = r56
 426 pow_GR_sig_int_Y          = r57
 427 pow_GR_sign_Y_Gpr         = r58
 428 pow_GR_17ones_m1          = r59
 429 pow_GR_one                = r60
 430 pow_GR_sign_Y             = r60
 431
 432 pow_GR_signexp_Y_Gpr      = r61
 433 pow_GR_exp_Y_Gpr          = r62
 434 pow_GR_true_exp_Y_Gpr     = r63
 435 pow_GR_signexp_Y          = r64
 436
 437 GR_SAVE_B0                = r65
 438 GR_SAVE_GP                = r66
 439 GR_SAVE_PFS               = r67
 440
 441 GR_Parameter_X            = r68
 442 GR_Parameter_Y            = r69
 443 GR_Parameter_RESULT       = r70
 444 pow_GR_tag                = r71
 445
 446
 447 // floating point registers used
 448
 449 POW_B                     = f32
 450 POW_NORM_X                = f33
 451 POW_Xm1                   = f34
 452 POW_r1                    = f34
 453 POW_P4                    = f35
 454
 455 POW_P5                    = f36
 456 POW_NORM_Y                = f37
 457 POW_Q2                    = f38
 458 POW_Q3                    = f39
 459 POW_P2                    = f40
 460
 461 POW_P3                    = f41
 462 POW_P0                    = f42
 463 POW_log2_lo               = f43
 464 POW_r                     = f44
 465 POW_Q0_half               = f45
 466
 467 POW_Q1                    = f46
 468 POW_log2_hi               = f48
 469 POW_Q4                    = f49
 470 POW_P1                    = f50
 471
 472 POW_log2_by_128_hi        = f51
 473 POW_inv_log2_by_128       = f52
 474 POW_rsq                   = f53
 475 POW_Yrcub                 = f54
 476 POW_log2_by_128_lo        = f55
 477
 478 POW_v6                    = f56
 479 POW_v4                    = f58
 480 POW_v2                    = f59
 481 POW_T                     = f60
 482
 483 POW_Tt                    = f61
 484 POW_RSHF                  = f62
 485 POW_v21ps                 = f63
 486 POW_s4                    = f64
 487
 488 POW_U                     = f66
 489 POW_G                     = f67
 490 POW_delta                 = f68
 491 POW_v3                    = f69
 492 POW_V                     = f70
 493
 494 POW_p                     = f71
 495 POW_Z1                    = f72
 496 POW_e3                    = f73
 497 POW_e2                    = f74
 498 POW_Z2                    = f75
 499
 500 POW_e1                    = f76
 501 POW_W1                    = f77
 502 POW_UmZ2                  = f78
 503 POW_W2                    = f79
 504 POW_Z3                    = f80
 505
 506 POW_int_W1                = f81
 507 POW_e12                   = f82
 508 POW_int_W2                = f83
 509 POW_UmZ2pV                = f84
 510 POW_Z3sq                  = f85
 511
 512 POW_e123                  = f86
 513 POW_N1float               = f87
 514 POW_N2float               = f88
 515 POW_f3                    = f89
 516 POW_q                     = f90
 517
 518 POW_s1                    = f91
 519 POW_Nfloat                = f92
 520 POW_s2                    = f93
 521 POW_f2                    = f94
 522 POW_f1                    = f95
 523
 524 POW_T1                    = f96
 525 POW_T2                    = f97
 526 POW_2M                    = f98
 527 POW_s                     = f99
 528 POW_f12                   = f100
 529
 530 POW_ssq                   = f101
 531 POW_T1T2                  = f102
 532 POW_1ps                   = f103
 533 POW_A                     = f104
 534 POW_es                    = f105
 535
 536 POW_int_K                 = f107
 537 POW_K                     = f108
 538 POW_f123                  = f109
 539 POW_Gpr                   = f110
 540
 541 POW_Y_Gpr                 = f111
 542 POW_int_Y                 = f112
 543
 544 POW_float_int_Y           = f116
 545 POW_ftz_urm_f8            = f117
 546 POW_wre_urm_f8            = f118
 547 POW_abs_A                 = f119
 548 POW_gt_pln                = f120
 549
 550 POW_xsq                   = f121
 551
 552 POW_twoV                  = f122
 553 POW_Xp1                   = f123
 554
 555 // Data tables
 556 //==============================================================
 557
 558 #ifdef _LIBC
 559 .rodata
 560 #else
 561 .data
 562 #endif
 563
 564 .align 16
 565
 566 pow_table_P:
 567 ASM_TYPE_DIRECTIVE(pow_table_P,@object)
 568 data8 0x8000F7B249FF332D, 0x0000BFFC  // P_5
 569 data8 0xAAAAAAA9E7902C7F, 0x0000BFFC  // P_3
 570 data8 0x80000000000018E5, 0x0000BFFD  // P_1
 571 data8 0xb8aa3b295c17f0bc, 0x00004006  // inv_ln2_by_128
 572
 573
 574 data8 0x3FA5555555554A9E // Q_2
 575 data8 0x3F8111124F4DD9F9 // Q_3
 576 data8 0x3FE0000000000000 // Q_0
 577 data8 0x3FC5555555554733 // Q_1
 578 data8 0x3F56C16D9360FFA0 // Q_4
 579 data8 0x43e8000000000000 // Right shift constant for exp
 580 data8 0xc9e3b39803f2f6af, 0x00003fb7  // ln2_by_128_lo
 581 data8 0x0000000000000000 // pad to eliminate bank conflicts with pow_table_Q
 582 data8 0x0000000000000000 // pad to eliminate bank conflicts with pow_table_Q
 583 ASM_SIZE_DIRECTIVE(pow_table_P)
 584
 585 pow_table_Q:
 586 ASM_TYPE_DIRECTIVE(pow_table_Q,@object)
 587 data8 0x9249FE7F0DC423CF, 0x00003FFC  // P_4
 588 data8 0xCCCCCCCC4ED2BA7F, 0x00003FFC  // P_2
 589 data8 0xAAAAAAAAAAAAB505, 0x00003FFD  // P_0
 590 data8 0x3fe62e42fefa39e8, 0x3cccd5e4f1d9cc02 // log2 hi lo =  +6.93147e-001
 591 data8 0xb17217f7d1cf79ab, 0x00003ff7  // ln2_by_128_hi
 592 ASM_SIZE_DIRECTIVE(pow_table_Q)
 593
 594
 595 pow_Tt:
 596 ASM_TYPE_DIRECTIVE(pow_Tt,@object)
 597 data8 0x3f60040155d58800, 0x3c93bce0ce3ddd81 // log(1/frcpa(1+0/256))=  +1.95503e-003
 598 data8 0x3f78121214586a00, 0x3cb540e0a5cfc9bc // log(1/frcpa(1+1/256))=  +5.87661e-003
 599 data8 0x3f841929f9683200, 0x3cbdf1d57404da1f // log(1/frcpa(1+2/256))=  +9.81362e-003
 600 data8 0x3f8c317384c75f00, 0x3c69806208c04c22 // log(1/frcpa(1+3/256))=  +1.37662e-002
 601 data8 0x3f91a6b91ac73380, 0x3c7874daa716eb32 // log(1/frcpa(1+4/256))=  +1.72376e-002
 602 data8 0x3f95ba9a5d9ac000, 0x3cacbb84e08d78ac // log(1/frcpa(1+5/256))=  +2.12196e-002
 603 data8 0x3f99d2a807432580, 0x3cbcf80538b441e1 // log(1/frcpa(1+6/256))=  +2.52177e-002
 604 data8 0x3f9d6b2725979800, 0x3c6095e5c8f8f359 // log(1/frcpa(1+7/256))=  +2.87291e-002
 605 data8 0x3fa0c58fa19dfa80, 0x3cb4c5d4e9d0dda2 // log(1/frcpa(1+8/256))=  +3.27573e-002
 606 data8 0x3fa2954c78cbce00, 0x3caa932b860ab8d6 // log(1/frcpa(1+9/256))=  +3.62953e-002
 607 data8 0x3fa4a94d2da96c40, 0x3ca670452b76bbd5 // log(1/frcpa(1+10/256))=  +4.03542e-002
 608 data8 0x3fa67c94f2d4bb40, 0x3ca84104f9941798 // log(1/frcpa(1+11/256))=  +4.39192e-002
 609 data8 0x3fa85188b630f040, 0x3cb40a882cbf0153 // log(1/frcpa(1+12/256))=  +4.74971e-002
 610 data8 0x3faa6b8abe73af40, 0x3c988d46e25c9059 // log(1/frcpa(1+13/256))=  +5.16017e-002
 611 data8 0x3fac441e06f72a80, 0x3cae3e930a1a2a96 // log(1/frcpa(1+14/256))=  +5.52072e-002
 612 data8 0x3fae1e6713606d00, 0x3c8a796f6283b580 // log(1/frcpa(1+15/256))=  +5.88257e-002
 613 data8 0x3faffa6911ab9300, 0x3c5193070351e88a // log(1/frcpa(1+16/256))=  +6.24574e-002
 614 data8 0x3fb0ec139c5da600, 0x3c623f2a75eb992d // log(1/frcpa(1+17/256))=  +6.61022e-002
 615 data8 0x3fb1dbd2643d1900, 0x3ca649b2ef8927f0 // log(1/frcpa(1+18/256))=  +6.97605e-002
 616 data8 0x3fb2cc7284fe5f00, 0x3cbc5e86599513e2 // log(1/frcpa(1+19/256))=  +7.34321e-002
 617 data8 0x3fb3bdf5a7d1ee60, 0x3c90bd4bb69dada3 // log(1/frcpa(1+20/256))=  +7.71173e-002
 618 data8 0x3fb4b05d7aa012e0, 0x3c54e377c9b8a54f // log(1/frcpa(1+21/256))=  +8.08161e-002
 619 data8 0x3fb580db7ceb5700, 0x3c7fdb2f98354cde // log(1/frcpa(1+22/256))=  +8.39975e-002
 620 data8 0x3fb674f089365a60, 0x3cb9994c9d3301c1 // log(1/frcpa(1+23/256))=  +8.77219e-002
 621 data8 0x3fb769ef2c6b5680, 0x3caaec639db52a79 // log(1/frcpa(1+24/256))=  +9.14602e-002
 622 data8 0x3fb85fd927506a40, 0x3c9f9f99a3cf8e25 // log(1/frcpa(1+25/256))=  +9.52125e-002
 623 data8 0x3fb9335e5d594980, 0x3ca15c3abd47d99a // log(1/frcpa(1+26/256))=  +9.84401e-002
 624 data8 0x3fba2b0220c8e5e0, 0x3cb4ca639adf6fc3 // log(1/frcpa(1+27/256))=  +1.02219e-001
 625 data8 0x3fbb0004ac1a86a0, 0x3ca7cb81bf959a59 // log(1/frcpa(1+28/256))=  +1.05469e-001
 626 data8 0x3fbbf968769fca00, 0x3cb0c646c121418e // log(1/frcpa(1+29/256))=  +1.09274e-001
 627 data8 0x3fbccfedbfee13a0, 0x3ca0465fce24ab4b // log(1/frcpa(1+30/256))=  +1.12548e-001
 628 data8 0x3fbda727638446a0, 0x3c82803f4e2e6603 // log(1/frcpa(1+31/256))=  +1.15832e-001
 629 data8 0x3fbea3257fe10f60, 0x3cb986a3f2313d1a // log(1/frcpa(1+32/256))=  +1.19677e-001
 630 data8 0x3fbf7be9fedbfde0, 0x3c97d16a6a621cf4 // log(1/frcpa(1+33/256))=  +1.22985e-001
 631 data8 0x3fc02ab352ff25f0, 0x3c9cc6baad365600 // log(1/frcpa(1+34/256))=  +1.26303e-001
 632 data8 0x3fc097ce579d2040, 0x3cb9ba16d329440b // log(1/frcpa(1+35/256))=  +1.29633e-001
 633 data8 0x3fc1178e8227e470, 0x3cb7bc671683f8e6 // log(1/frcpa(1+36/256))=  +1.33531e-001
 634 data8 0x3fc185747dbecf30, 0x3c9d1116f66d2345 // log(1/frcpa(1+37/256))=  +1.36885e-001
 635 data8 0x3fc1f3b925f25d40, 0x3c8162c9ef939ac6 // log(1/frcpa(1+38/256))=  +1.40250e-001
 636 data8 0x3fc2625d1e6ddf50, 0x3caad3a1ec384fc3 // log(1/frcpa(1+39/256))=  +1.43627e-001
 637 data8 0x3fc2d1610c868130, 0x3cb3ad997036941b // log(1/frcpa(1+40/256))=  +1.47015e-001
 638 data8 0x3fc340c597411420, 0x3cbc2308262c7998 // log(1/frcpa(1+41/256))=  +1.50414e-001
 639 data8 0x3fc3b08b6757f2a0, 0x3cb2170d6cdf0526 // log(1/frcpa(1+42/256))=  +1.53825e-001
 640 data8 0x3fc40dfb08378000, 0x3c9bb453c4f7b685 // log(1/frcpa(1+43/256))=  +1.56677e-001
 641 data8 0x3fc47e74e8ca5f70, 0x3cb836a48fdfce9d // log(1/frcpa(1+44/256))=  +1.60109e-001
 642 data8 0x3fc4ef51f6466de0, 0x3ca07a43919aa64b // log(1/frcpa(1+45/256))=  +1.63553e-001
 643 data8 0x3fc56092e02ba510, 0x3ca85006899d97b0 // log(1/frcpa(1+46/256))=  +1.67010e-001
 644 data8 0x3fc5d23857cd74d0, 0x3ca30a5ba6e7abbe // log(1/frcpa(1+47/256))=  +1.70478e-001
 645 data8 0x3fc6313a37335d70, 0x3ca905586f0ac97e // log(1/frcpa(1+48/256))=  +1.73377e-001
 646 data8 0x3fc6a399dabbd380, 0x3c9b2c6657a96684 // log(1/frcpa(1+49/256))=  +1.76868e-001
 647 data8 0x3fc70337dd3ce410, 0x3cb50bc52f55cdd8 // log(1/frcpa(1+50/256))=  +1.79786e-001
 648 data8 0x3fc77654128f6120, 0x3cad2eb7c9a39efe // log(1/frcpa(1+51/256))=  +1.83299e-001
 649 data8 0x3fc7e9d82a0b0220, 0x3cba127e90393c01 // log(1/frcpa(1+52/256))=  +1.86824e-001
 650 data8 0x3fc84a6b759f5120, 0x3cbd7fd52079f706 // log(1/frcpa(1+53/256))=  +1.89771e-001
 651 data8 0x3fc8ab47d5f5a300, 0x3cbfae141751a3de // log(1/frcpa(1+54/256))=  +1.92727e-001
 652 data8 0x3fc91fe490965810, 0x3cb69cf30a1c319e // log(1/frcpa(1+55/256))=  +1.96286e-001
 653 data8 0x3fc981634011aa70, 0x3ca5bb3d208bc42a // log(1/frcpa(1+56/256))=  +1.99261e-001
 654 data8 0x3fc9f6c407089660, 0x3ca04d68658179a0 // log(1/frcpa(1+57/256))=  +2.02843e-001
 655 data8 0x3fca58e729348f40, 0x3c99f5411546c286 // log(1/frcpa(1+58/256))=  +2.05838e-001
 656 data8 0x3fcabb55c31693a0, 0x3cb9a5350eb327d5 // log(1/frcpa(1+59/256))=  +2.08842e-001
 657 data8 0x3fcb1e104919efd0, 0x3c18965fcce7c406 // log(1/frcpa(1+60/256))=  +2.11855e-001
 658 data8 0x3fcb94ee93e367c0, 0x3cb503716da45184 // log(1/frcpa(1+61/256))=  +2.15483e-001
 659 data8 0x3fcbf851c0675550, 0x3cbdf1b3f7ab5378 // log(1/frcpa(1+62/256))=  +2.18516e-001
 660 data8 0x3fcc5c0254bf23a0, 0x3ca7aab9ed0b1d7b // log(1/frcpa(1+63/256))=  +2.21558e-001
 661 data8 0x3fccc000c9db3c50, 0x3c92a7a2a850072a // log(1/frcpa(1+64/256))=  +2.24609e-001
 662 data8 0x3fcd244d99c85670, 0x3c9f6019120edf4c // log(1/frcpa(1+65/256))=  +2.27670e-001
 663 data8 0x3fcd88e93fb2f450, 0x3c6affb96815e081 // log(1/frcpa(1+66/256))=  +2.30741e-001
 664 data8 0x3fcdedd437eaef00, 0x3c72553595897976 // log(1/frcpa(1+67/256))=  +2.33820e-001
 665 data8 0x3fce530effe71010, 0x3c90913b020fa182 // log(1/frcpa(1+68/256))=  +2.36910e-001
 666 data8 0x3fceb89a1648b970, 0x3c837ba4045bfd25 // log(1/frcpa(1+69/256))=  +2.40009e-001
 667 data8 0x3fcf1e75fadf9bd0, 0x3cbcea6d13e0498d // log(1/frcpa(1+70/256))=  +2.43117e-001
 668 data8 0x3fcf84a32ead7c30, 0x3ca5e3a67b3c6d77 // log(1/frcpa(1+71/256))=  +2.46235e-001
 669 data8 0x3fcfeb2233ea07c0, 0x3cba0c6f0049c5a6 // log(1/frcpa(1+72/256))=  +2.49363e-001
 670 data8 0x3fd028f9c7035c18, 0x3cb0a30b06677ff6 // log(1/frcpa(1+73/256))=  +2.52501e-001
 671 data8 0x3fd05c8be0d96358, 0x3ca0f1c77ccb5865 // log(1/frcpa(1+74/256))=  +2.55649e-001
 672 data8 0x3fd085eb8f8ae790, 0x3cbd513f45fe7a97 // log(1/frcpa(1+75/256))=  +2.58174e-001
 673 data8 0x3fd0b9c8e32d1910, 0x3c927449047ca006 // log(1/frcpa(1+76/256))=  +2.61339e-001
 674 data8 0x3fd0edd060b78080, 0x3c89b52d8435f53e // log(1/frcpa(1+77/256))=  +2.64515e-001
 675 data8 0x3fd122024cf00638, 0x3cbdd976fabda4bd // log(1/frcpa(1+78/256))=  +2.67701e-001
 676 data8 0x3fd14be2927aecd0, 0x3cb02f90ad0bc471 // log(1/frcpa(1+79/256))=  +2.70257e-001
 677 data8 0x3fd180618ef18ad8, 0x3cbd003792c71a98 // log(1/frcpa(1+80/256))=  +2.73461e-001
 678 data8 0x3fd1b50bbe2fc638, 0x3ca9ae64c6403ead // log(1/frcpa(1+81/256))=  +2.76675e-001
 679 data8 0x3fd1df4cc7cf2428, 0x3cb43f0455f7e395 // log(1/frcpa(1+82/256))=  +2.79254e-001
 680 data8 0x3fd214456d0eb8d0, 0x3cb0fbd748d75d30 // log(1/frcpa(1+83/256))=  +2.82487e-001
 681 data8 0x3fd23ec5991eba48, 0x3c906edd746b77e2 // log(1/frcpa(1+84/256))=  +2.85081e-001
 682 data8 0x3fd2740d9f870af8, 0x3ca9802e6a00a670 // log(1/frcpa(1+85/256))=  +2.88333e-001
 683 data8 0x3fd29ecdabcdfa00, 0x3cacecef70890cfa // log(1/frcpa(1+86/256))=  +2.90943e-001
 684 data8 0x3fd2d46602adcce8, 0x3cb97911955f3521 // log(1/frcpa(1+87/256))=  +2.94214e-001
 685 data8 0x3fd2ff66b04ea9d0, 0x3cb12dabe191d1c9 // log(1/frcpa(1+88/256))=  +2.96838e-001
 686 data8 0x3fd335504b355a30, 0x3cbdf9139df924ec // log(1/frcpa(1+89/256))=  +3.00129e-001
 687 data8 0x3fd360925ec44f58, 0x3cb253e68977a1e3 // log(1/frcpa(1+90/256))=  +3.02769e-001
 688 data8 0x3fd38bf1c3337e70, 0x3cb3d283d2a2da21 // log(1/frcpa(1+91/256))=  +3.05417e-001
 689 data8 0x3fd3c25277333180, 0x3cadaa5b035eae27 // log(1/frcpa(1+92/256))=  +3.08735e-001
 690 data8 0x3fd3edf463c16838, 0x3cb983d680d3c108 // log(1/frcpa(1+93/256))=  +3.11399e-001
 691 data8 0x3fd419b423d5e8c0, 0x3cbc86dd921c139d // log(1/frcpa(1+94/256))=  +3.14069e-001
 692 data8 0x3fd44591e0539f48, 0x3c86a76d6dc2782e // log(1/frcpa(1+95/256))=  +3.16746e-001
 693 data8 0x3fd47c9175b6f0a8, 0x3cb59a2e013c6b5f // log(1/frcpa(1+96/256))=  +3.20103e-001
 694 data8 0x3fd4a8b341552b08, 0x3c93f1e86e468694 // log(1/frcpa(1+97/256))=  +3.22797e-001
 695 data8 0x3fd4d4f390890198, 0x3cbf5e4ea7c5105a // log(1/frcpa(1+98/256))=  +3.25498e-001
 696 data8 0x3fd501528da1f960, 0x3cbf58da53e9ad10 // log(1/frcpa(1+99/256))=  +3.28206e-001
 697 data8 0x3fd52dd06347d4f0, 0x3cb98a28cebf6eef // log(1/frcpa(1+100/256))=  +3.30921e-001
 698 data8 0x3fd55a6d3c7b8a88, 0x3c9c76b67c2d1fd4 // log(1/frcpa(1+101/256))=  +3.33644e-001
 699 data8 0x3fd5925d2b112a58, 0x3c9029616a4331b8 // log(1/frcpa(1+102/256))=  +3.37058e-001
 700 data8 0x3fd5bf406b543db0, 0x3c9fb8292ecfc820 // log(1/frcpa(1+103/256))=  +3.39798e-001
 701 data8 0x3fd5ec433d5c35a8, 0x3cb71a1229d17eec // log(1/frcpa(1+104/256))=  +3.42545e-001
 702 data8 0x3fd61965cdb02c18, 0x3cbba94fe1dbb8d2 // log(1/frcpa(1+105/256))=  +3.45300e-001
 703 data8 0x3fd646a84935b2a0, 0x3c9ee496d2c9ae57 // log(1/frcpa(1+106/256))=  +3.48063e-001
 704 data8 0x3fd6740add31de90, 0x3cb1da3a6c7a9dfd // log(1/frcpa(1+107/256))=  +3.50833e-001
 705 data8 0x3fd6a18db74a58c0, 0x3cb494c257add8dc // log(1/frcpa(1+108/256))=  +3.53610e-001
 706 data8 0x3fd6cf31058670e8, 0x3cb0b244a70a8da9 // log(1/frcpa(1+109/256))=  +3.56396e-001
 707 data8 0x3fd6f180e852f0b8, 0x3c9db7aefa866720 // log(1/frcpa(1+110/256))=  +3.58490e-001
 708 data8 0x3fd71f5d71b894e8, 0x3cbe91c4bf324957 // log(1/frcpa(1+111/256))=  +3.61289e-001
 709 data8 0x3fd74d5aefd66d58, 0x3cb06b3d9bfac023 // log(1/frcpa(1+112/256))=  +3.64096e-001
 710 data8 0x3fd77b79922bd378, 0x3cb727d8804491f4 // log(1/frcpa(1+113/256))=  +3.66911e-001
 711 data8 0x3fd7a9b9889f19e0, 0x3ca2ef22df5bc543 // log(1/frcpa(1+114/256))=  +3.69734e-001
 712 data8 0x3fd7d81b037eb6a0, 0x3cb8fd3ba07a7ece // log(1/frcpa(1+115/256))=  +3.72565e-001
 713 data8 0x3fd8069e33827230, 0x3c8bd1e25866e61a // log(1/frcpa(1+116/256))=  +3.75404e-001
 714 data8 0x3fd82996d3ef8bc8, 0x3ca5aab9f5928928 // log(1/frcpa(1+117/256))=  +3.77538e-001
 715 data8 0x3fd85855776dcbf8, 0x3ca56f33337789d6 // log(1/frcpa(1+118/256))=  +3.80391e-001
 716 data8 0x3fd8873658327cc8, 0x3cbb8ef0401db49d // log(1/frcpa(1+119/256))=  +3.83253e-001
 717 data8 0x3fd8aa75973ab8c8, 0x3cbb9961f509a680 // log(1/frcpa(1+120/256))=  +3.85404e-001
 718 data8 0x3fd8d992dc8824e0, 0x3cb220512a53732d // log(1/frcpa(1+121/256))=  +3.88280e-001
 719 data8 0x3fd908d2ea7d9510, 0x3c985f0e513bfb5c // log(1/frcpa(1+122/256))=  +3.91164e-001
 720 data8 0x3fd92c59e79c0e50, 0x3cb82e073fd30d63 // log(1/frcpa(1+123/256))=  +3.93332e-001
 721 data8 0x3fd95bd750ee3ed0, 0x3ca4aa7cdb6dd8a8 // log(1/frcpa(1+124/256))=  +3.96231e-001
 722 data8 0x3fd98b7811a3ee58, 0x3caa93a5b660893e // log(1/frcpa(1+125/256))=  +3.99138e-001
 723 data8 0x3fd9af47f33d4068, 0x3cac294b3b3190ba // log(1/frcpa(1+126/256))=  +4.01323e-001
 724 data8 0x3fd9df270c1914a0, 0x3cbe1a58fd0cd67e // log(1/frcpa(1+127/256))=  +4.04245e-001
 725 data8 0x3fda0325ed14fda0, 0x3cb1efa7950fb57e // log(1/frcpa(1+128/256))=  +4.06442e-001
 726 data8 0x3fda33440224fa78, 0x3c8915fe75e7d477 // log(1/frcpa(1+129/256))=  +4.09379e-001
 727 data8 0x3fda57725e80c380, 0x3ca72bd1062b1b7f // log(1/frcpa(1+130/256))=  +4.11587e-001
 728 data8 0x3fda87d0165dd198, 0x3c91f7845f58dbad // log(1/frcpa(1+131/256))=  +4.14539e-001
 729 data8 0x3fdaac2e6c03f890, 0x3cb6f237a911c509 // log(1/frcpa(1+132/256))=  +4.16759e-001
 730 data8 0x3fdadccc6fdf6a80, 0x3c90ddc4b7687169 // log(1/frcpa(1+133/256))=  +4.19726e-001
 731 data8 0x3fdb015b3eb1e790, 0x3c692dd7d90e1e8e // log(1/frcpa(1+134/256))=  +4.21958e-001
 732 data8 0x3fdb323a3a635948, 0x3c6f85655cbe14de // log(1/frcpa(1+135/256))=  +4.24941e-001
 733 data8 0x3fdb56fa04462908, 0x3c95252d841994de // log(1/frcpa(1+136/256))=  +4.27184e-001
 734 data8 0x3fdb881aa659bc90, 0x3caa53a745a3642f // log(1/frcpa(1+137/256))=  +4.30182e-001
 735 data8 0x3fdbad0bef3db160, 0x3cb32f2540dcc16a // log(1/frcpa(1+138/256))=  +4.32437e-001
 736 data8 0x3fdbd21297781c28, 0x3cbd8e891e106f1d // log(1/frcpa(1+139/256))=  +4.34697e-001
 737 data8 0x3fdc039236f08818, 0x3c809435af522ba7 // log(1/frcpa(1+140/256))=  +4.37718e-001
 738 data8 0x3fdc28cb1e4d32f8, 0x3cb3944752fbd81e // log(1/frcpa(1+141/256))=  +4.39990e-001
 739 data8 0x3fdc4e19b84723c0, 0x3c9a465260cd3fe5 // log(1/frcpa(1+142/256))=  +4.42267e-001
 740 data8 0x3fdc7ff9c74554c8, 0x3c92447d5b6ca369 // log(1/frcpa(1+143/256))=  +4.45311e-001
 741 data8 0x3fdca57b64e9db00, 0x3cb44344a8a00c82 // log(1/frcpa(1+144/256))=  +4.47600e-001
 742 data8 0x3fdccb130a5ceba8, 0x3cbefaddfb97b73f // log(1/frcpa(1+145/256))=  +4.49895e-001
 743 data8 0x3fdcf0c0d18f3268, 0x3cbd3e7bfee57898 // log(1/frcpa(1+146/256))=  +4.52194e-001
 744 data8 0x3fdd232075b5a200, 0x3c9222599987447c // log(1/frcpa(1+147/256))=  +4.55269e-001
 745 data8 0x3fdd490246defa68, 0x3cabafe9a767a80d // log(1/frcpa(1+148/256))=  +4.57581e-001
 746 data8 0x3fdd6efa918d25c8, 0x3cb58a2624e1c6fd // log(1/frcpa(1+149/256))=  +4.59899e-001
 747 data8 0x3fdd9509707ae528, 0x3cbdc3babce578e7 // log(1/frcpa(1+150/256))=  +4.62221e-001
 748 data8 0x3fddbb2efe92c550, 0x3cb0ac0943c434a4 // log(1/frcpa(1+151/256))=  +4.64550e-001
 749 data8 0x3fddee2f3445e4a8, 0x3cbba9d07ce820e8 // log(1/frcpa(1+152/256))=  +4.67663e-001
 750 data8 0x3fde148a1a2726c8, 0x3cb6537e3375b205 // log(1/frcpa(1+153/256))=  +4.70004e-001
 751 data8 0x3fde3afc0a49ff38, 0x3cbfed5518dbc20e // log(1/frcpa(1+154/256))=  +4.72350e-001
 752 data8 0x3fde6185206d5168, 0x3cb6572601f73d5c // log(1/frcpa(1+155/256))=  +4.74702e-001
 753 data8 0x3fde882578823d50, 0x3c9b24abd4584d1a // log(1/frcpa(1+156/256))=  +4.77060e-001
 754 data8 0x3fdeaedd2eac9908, 0x3cb0ceb5e4d2c8f7 // log(1/frcpa(1+157/256))=  +4.79423e-001
 755 data8 0x3fded5ac5f436be0, 0x3ca72f21f1f5238e // log(1/frcpa(1+158/256))=  +4.81792e-001
 756 data8 0x3fdefc9326d16ab8, 0x3c85081a1639a45c // log(1/frcpa(1+159/256))=  +4.84166e-001
 757 data8 0x3fdf2391a21575f8, 0x3cbf11015bdd297a // log(1/frcpa(1+160/256))=  +4.86546e-001
 758 data8 0x3fdf4aa7ee031928, 0x3cb3795bc052a2d1 // log(1/frcpa(1+161/256))=  +4.88932e-001
 759 data8 0x3fdf71d627c30bb0, 0x3c35c61f0f5a88f3 // log(1/frcpa(1+162/256))=  +4.91323e-001
 760 data8 0x3fdf991c6cb3b378, 0x3c97d99419be6028 // log(1/frcpa(1+163/256))=  +4.93720e-001
 761 data8 0x3fdfc07ada69a908, 0x3cbfe9341ded70b1 // log(1/frcpa(1+164/256))=  +4.96123e-001
 762 data8 0x3fdfe7f18eb03d38, 0x3cb85718a640c33f // log(1/frcpa(1+165/256))=  +4.98532e-001
 763 data8 0x3fe007c053c5002c, 0x3cb3addc9c065f09 // log(1/frcpa(1+166/256))=  +5.00946e-001
 764 data8 0x3fe01b942198a5a0, 0x3c9d5aa4c77da6ac // log(1/frcpa(1+167/256))=  +5.03367e-001
 765 data8 0x3fe02f74400c64e8, 0x3cb5a0ee4450ef52 // log(1/frcpa(1+168/256))=  +5.05793e-001
 766 data8 0x3fe04360be7603ac, 0x3c9dd00c35630fe0 // log(1/frcpa(1+169/256))=  +5.08225e-001
 767 data8 0x3fe05759ac47fe30, 0x3cbd063e1f0bd82c // log(1/frcpa(1+170/256))=  +5.10663e-001
 768 data8 0x3fe06b5f1911cf50, 0x3cae8da674af5289 // log(1/frcpa(1+171/256))=  +5.13107e-001
 769 data8 0x3fe078bf0533c568, 0x3c62241edf5fd1f7 // log(1/frcpa(1+172/256))=  +5.14740e-001
 770 data8 0x3fe08cd9687e7b0c, 0x3cb3007febcca227 // log(1/frcpa(1+173/256))=  +5.17194e-001
 771 data8 0x3fe0a10074cf9018, 0x3ca496e84603816b // log(1/frcpa(1+174/256))=  +5.19654e-001
 772 data8 0x3fe0b5343a234474, 0x3cb46098d14fc90a // log(1/frcpa(1+175/256))=  +5.22120e-001
 773 data8 0x3fe0c974c89431cc, 0x3cac0a7cdcbb86c6 // log(1/frcpa(1+176/256))=  +5.24592e-001
 774 data8 0x3fe0ddc2305b9884, 0x3cb2f753210410ff // log(1/frcpa(1+177/256))=  +5.27070e-001
 775 data8 0x3fe0eb524bafc918, 0x3c88affd6682229e // log(1/frcpa(1+178/256))=  +5.28726e-001
 776 data8 0x3fe0ffb54213a474, 0x3cadeefbab9af993 // log(1/frcpa(1+179/256))=  +5.31214e-001
 777 data8 0x3fe114253da97d9c, 0x3cbaf1c2b8bc160a // log(1/frcpa(1+180/256))=  +5.33709e-001
 778 data8 0x3fe128a24f1d9afc, 0x3cb9cf4df375e650 // log(1/frcpa(1+181/256))=  +5.36210e-001
 779 data8 0x3fe1365252bf0864, 0x3c985a621d4be111 // log(1/frcpa(1+182/256))=  +5.37881e-001
 780 data8 0x3fe14ae558b4a92c, 0x3ca104c4aa8977d1 // log(1/frcpa(1+183/256))=  +5.40393e-001
 781 data8 0x3fe15f85a19c7658, 0x3cbadf26e540f375 // log(1/frcpa(1+184/256))=  +5.42910e-001
 782 data8 0x3fe16d4d38c119f8, 0x3cb3aea11caec416 // log(1/frcpa(1+185/256))=  +5.44592e-001
 783 data8 0x3fe18203c20dd130, 0x3cba82d1211d1d6d // log(1/frcpa(1+186/256))=  +5.47121e-001
 784 data8 0x3fe196c7bc4b1f38, 0x3cb6267acc4f4f4a // log(1/frcpa(1+187/256))=  +5.49656e-001
 785 data8 0x3fe1a4a738b7a33c, 0x3c858930213c987d // log(1/frcpa(1+188/256))=  +5.51349e-001
 786 data8 0x3fe1b981c0c9653c, 0x3c9bc2a4a30f697b // log(1/frcpa(1+189/256))=  +5.53895e-001
 787 data8 0x3fe1ce69e8bb1068, 0x3cb7ae6199cf2a00 // log(1/frcpa(1+190/256))=  +5.56447e-001
 788 data8 0x3fe1dc619de06944, 0x3c6b50bb38388177 // log(1/frcpa(1+191/256))=  +5.58152e-001
 789 data8 0x3fe1f160a2ad0da0, 0x3cbd05b2778a5e1d // log(1/frcpa(1+192/256))=  +5.60715e-001
 790 data8 0x3fe2066d7740737c, 0x3cb32e828f9c6bd6 // log(1/frcpa(1+193/256))=  +5.63285e-001
 791 data8 0x3fe2147dba47a390, 0x3cbd579851b8b672 // log(1/frcpa(1+194/256))=  +5.65001e-001
 792 data8 0x3fe229a1bc5ebac0, 0x3cbb321be5237ce8 // log(1/frcpa(1+195/256))=  +5.67582e-001
 793 data8 0x3fe237c1841a502c, 0x3cb3b56e0915ea64 // log(1/frcpa(1+196/256))=  +5.69306e-001
 794 data8 0x3fe24cfce6f80d98, 0x3cb34a4d1a422919 // log(1/frcpa(1+197/256))=  +5.71898e-001
 795 data8 0x3fe25b2c55cd5760, 0x3cb237401ea5015e // log(1/frcpa(1+198/256))=  +5.73630e-001
 796 data8 0x3fe2707f4d5f7c40, 0x3c9d30f20acc8341 // log(1/frcpa(1+199/256))=  +5.76233e-001
 797 data8 0x3fe285e0842ca380, 0x3cbc4d866d5f21c0 // log(1/frcpa(1+200/256))=  +5.78842e-001
 798 data8 0x3fe294294708b770, 0x3cb85e14d5dc54fa // log(1/frcpa(1+201/256))=  +5.80586e-001
 799 data8 0x3fe2a9a2670aff0c, 0x3c7e6f8f468bbf91 // log(1/frcpa(1+202/256))=  +5.83207e-001
 800 data8 0x3fe2b7fb2c8d1cc0, 0x3c930ffcf63c8b65 // log(1/frcpa(1+203/256))=  +5.84959e-001
 801 data8 0x3fe2c65a6395f5f4, 0x3ca0afe20b53d2d2 // log(1/frcpa(1+204/256))=  +5.86713e-001
 802 data8 0x3fe2dbf557b0df40, 0x3cb646be1188fbc9 // log(1/frcpa(1+205/256))=  +5.89350e-001
 803 data8 0x3fe2ea64c3f97654, 0x3c96516fa8df33b2 // log(1/frcpa(1+206/256))=  +5.91113e-001
 804 data8 0x3fe3001823684d70, 0x3cb96d64e16d1360 // log(1/frcpa(1+207/256))=  +5.93762e-001
 805 data8 0x3fe30e97e9a8b5cc, 0x3c98ef96bc97cca0 // log(1/frcpa(1+208/256))=  +5.95531e-001
 806 data8 0x3fe32463ebdd34e8, 0x3caef1dc9a56c1bf // log(1/frcpa(1+209/256))=  +5.98192e-001
 807 data8 0x3fe332f4314ad794, 0x3caa4f0ac5d5fa11 // log(1/frcpa(1+210/256))=  +5.99970e-001
 808 data8 0x3fe348d90e7464cc, 0x3cbe7889f0516acd // log(1/frcpa(1+211/256))=  +6.02643e-001
 809 data8 0x3fe35779f8c43d6c, 0x3ca96bbab7245411 // log(1/frcpa(1+212/256))=  +6.04428e-001
 810 data8 0x3fe36621961a6a98, 0x3ca31f32262db9fb // log(1/frcpa(1+213/256))=  +6.06217e-001
 811 data8 0x3fe37c299f3c3668, 0x3cb15c72c107ee29 // log(1/frcpa(1+214/256))=  +6.08907e-001
 812 data8 0x3fe38ae2171976e4, 0x3cba42a2554b2dd4 // log(1/frcpa(1+215/256))=  +6.10704e-001
 813 data8 0x3fe399a157a603e4, 0x3cb99c62286d8919 // log(1/frcpa(1+216/256))=  +6.12504e-001
 814 data8 0x3fe3afccfe77b9d0, 0x3ca11048f96a43bd // log(1/frcpa(1+217/256))=  +6.15210e-001
 815 data8 0x3fe3be9d503533b4, 0x3ca4022f47588c3e // log(1/frcpa(1+218/256))=  +6.17018e-001
 816 data8 0x3fe3cd7480b4a8a0, 0x3cb4ba7afc2dc56a // log(1/frcpa(1+219/256))=  +6.18830e-001
 817 data8 0x3fe3e3c43918f76c, 0x3c859673d064b8ba // log(1/frcpa(1+220/256))=  +6.21554e-001
 818 data8 0x3fe3f2acb27ed6c4, 0x3cb55c6b452a16a8 // log(1/frcpa(1+221/256))=  +6.23373e-001
 819 data8 0x3fe4019c2125ca90, 0x3cb8c367879c5a31 // log(1/frcpa(1+222/256))=  +6.25197e-001
 820 data8 0x3fe4181061389720, 0x3cb2c17a79c5cc6c // log(1/frcpa(1+223/256))=  +6.27937e-001
 821 data8 0x3fe42711518df544, 0x3ca5f38d47012fc5 // log(1/frcpa(1+224/256))=  +6.29769e-001
 822 data8 0x3fe436194e12b6bc, 0x3cb9854d65a9b426 // log(1/frcpa(1+225/256))=  +6.31604e-001
 823 data8 0x3fe445285d68ea68, 0x3ca3ff9b3a81cd81 // log(1/frcpa(1+226/256))=  +6.33442e-001
 824 data8 0x3fe45bcc464c8938, 0x3cb0a2d8011a6c05 // log(1/frcpa(1+227/256))=  +6.36206e-001
 825 data8 0x3fe46aed21f117fc, 0x3c8a2be41f8e9f3d // log(1/frcpa(1+228/256))=  +6.38053e-001
 826 data8 0x3fe47a1527e8a2d0, 0x3cba4a83594fab09 // log(1/frcpa(1+229/256))=  +6.39903e-001
 827 data8 0x3fe489445efffcc8, 0x3cbf306a23dcbcde // log(1/frcpa(1+230/256))=  +6.41756e-001
 828 data8 0x3fe4a018bcb69834, 0x3ca46c9285029fd1 // log(1/frcpa(1+231/256))=  +6.44543e-001
 829 data8 0x3fe4af5a0c9d65d4, 0x3cbbc1db897580e3 // log(1/frcpa(1+232/256))=  +6.46405e-001
 830 data8 0x3fe4bea2a5bdbe84, 0x3cb84d880d7ef775 // log(1/frcpa(1+233/256))=  +6.48271e-001
 831 data8 0x3fe4cdf28f10ac44, 0x3cb3ec4b7893ce1f // log(1/frcpa(1+234/256))=  +6.50140e-001
 832 data8 0x3fe4dd49cf994058, 0x3c897224d59d3408 // log(1/frcpa(1+235/256))=  +6.52013e-001
 833 data8 0x3fe4eca86e64a680, 0x3cbccf620f24f0cd // log(1/frcpa(1+236/256))=  +6.53889e-001
 834 data8 0x3fe503c43cd8eb68, 0x3c3f872c65971084 // log(1/frcpa(1+237/256))=  +6.56710e-001
 835 data8 0x3fe513356667fc54, 0x3cb9ca64cc3d52c8 // log(1/frcpa(1+238/256))=  +6.58595e-001
 836 data8 0x3fe522ae0738a3d4, 0x3cbe708164c75968 // log(1/frcpa(1+239/256))=  +6.60483e-001
 837 data8 0x3fe5322e26867854, 0x3cb9988ba4aea615 // log(1/frcpa(1+240/256))=  +6.62376e-001
 838 data8 0x3fe541b5cb979808, 0x3ca1662e3a6b95f5 // log(1/frcpa(1+241/256))=  +6.64271e-001
 839 data8 0x3fe55144fdbcbd60, 0x3cb3acd4ca45c1e0 // log(1/frcpa(1+242/256))=  +6.66171e-001
 840 data8 0x3fe560dbc45153c4, 0x3cb4988947959fed // log(1/frcpa(1+243/256))=  +6.68074e-001
 841 data8 0x3fe5707a26bb8c64, 0x3cb3017fe6607ba9 // log(1/frcpa(1+244/256))=  +6.69980e-001
 842 data8 0x3fe587f60ed5b8fc, 0x3cbe7a3266366ed4 // log(1/frcpa(1+245/256))=  +6.72847e-001
 843 data8 0x3fe597a7977c8f30, 0x3ca1e12b9959a90e // log(1/frcpa(1+246/256))=  +6.74763e-001
 844 data8 0x3fe5a760d634bb88, 0x3cb7c365e53d9602 // log(1/frcpa(1+247/256))=  +6.76682e-001
 845 data8 0x3fe5b721d295f10c, 0x3cb716c2551ccbf0 // log(1/frcpa(1+248/256))=  +6.78605e-001
 846 data8 0x3fe5c6ea94431ef8, 0x3ca02b2ed0e28261 // log(1/frcpa(1+249/256))=  +6.80532e-001
 847 data8 0x3fe5d6bb22ea86f4, 0x3caf43a8bbb2f974 // log(1/frcpa(1+250/256))=  +6.82462e-001
 848 data8 0x3fe5e6938645d38c, 0x3cbcedc98821b333 // log(1/frcpa(1+251/256))=  +6.84397e-001
 849 data8 0x3fe5f673c61a2ed0, 0x3caa385eef5f2789 // log(1/frcpa(1+252/256))=  +6.86335e-001
 850 data8 0x3fe6065bea385924, 0x3cb11624f165c5b4 // log(1/frcpa(1+253/256))=  +6.88276e-001
 851 data8 0x3fe6164bfa7cc068, 0x3cbad884f87073fa // log(1/frcpa(1+254/256))=  +6.90222e-001
 852 data8 0x3fe62643fecf9740, 0x3cb78c51da12f4df // log(1/frcpa(1+255/256))=  +6.92171e-001
 853 ASM_SIZE_DIRECTIVE(pow_Tt)
 854
 855
 856 // Table 1 is 2^(index_1/128) where
 857 // index_1 goes from 0 to 15
 858 pow_tbl1:
 859 ASM_TYPE_DIRECTIVE(pow_tbl1,@object)
 860 data8 0x8000000000000000 , 0x00003FFF
 861 data8 0x80B1ED4FD999AB6C , 0x00003FFF
 862 data8 0x8164D1F3BC030773 , 0x00003FFF
 863 data8 0x8218AF4373FC25EC , 0x00003FFF
 864 data8 0x82CD8698AC2BA1D7 , 0x00003FFF
 865 data8 0x8383594EEFB6EE37 , 0x00003FFF
 866 data8 0x843A28C3ACDE4046 , 0x00003FFF
 867 data8 0x84F1F656379C1A29 , 0x00003FFF
 868 data8 0x85AAC367CC487B15 , 0x00003FFF
 869 data8 0x8664915B923FBA04 , 0x00003FFF
 870 data8 0x871F61969E8D1010 , 0x00003FFF
 871 data8 0x87DB357FF698D792 , 0x00003FFF
 872 data8 0x88980E8092DA8527 , 0x00003FFF
 873 data8 0x8955EE03618E5FDD , 0x00003FFF
 874 data8 0x8A14D575496EFD9A , 0x00003FFF
 875 data8 0x8AD4C6452C728924 , 0x00003FFF
 876 ASM_SIZE_DIRECTIVE(pow_tbl1)
 877
 878
 879 // Table 2 is 2^(index_1/8) where
 880 // index_2 goes from 0 to 7
 881 pow_tbl2:
 882 ASM_TYPE_DIRECTIVE(pow_tbl2,@object)
 883 data8 0x8000000000000000 , 0x00003FFF
 884 data8 0x8B95C1E3EA8BD6E7 , 0x00003FFF
 885 data8 0x9837F0518DB8A96F , 0x00003FFF
 886 data8 0xA5FED6A9B15138EA , 0x00003FFF
 887 data8 0xB504F333F9DE6484 , 0x00003FFF
 888 data8 0xC5672A115506DADD , 0x00003FFF
 889 data8 0xD744FCCAD69D6AF4 , 0x00003FFF
 890 data8 0xEAC0C6E7DD24392F , 0x00003FFF
 891 ASM_SIZE_DIRECTIVE(pow_tbl2)
 892
 893 .global pow
 894
 895 .section .text
 896 .proc  pow
 897 .align 32
 898
 899 pow:
 900
 901 { .mfi
 902           alloc         r32=ar.pfs,1,35,4,0
 903           fms.s1 POW_Xm1 = f8,f1,f1   // Will be used for r1 if x>0
 904           mov           pow_GR_17ones  = 0x1FFFF
 905 }
 906 { .mfi
 907 (p0)      addl          pow_AD_P   = @ltoff(pow_table_P), gp
 908           fma.s1 POW_Xp1 = f8,f1,f1   // Will be used for r1 if x<0
 909           nop.i 999
 910 ;;
 911 }
 912
 913
 914 // Get exponent of x.  Will be used to calculate K.
 915 { .mfi
 916           getf.exp      pow_GR_signexp_X    = f8
 917           frcpa.s1      POW_B, p6   = f1,f8
 918           nop.i 999
 919 }
 920 { .mfi
 921           ld8 pow_AD_P = [pow_AD_P]
 922           fma.s1        POW_NORM_X     = f8,f1,f0
 923           mov           pow_GR_FFF7    = 0xFFF7
 924 }
 925 ;;
 926
 927
 928
 929 // Get significand of x.  Will be used to get index to fetch T, Tt.
 930 // p13 = TRUE ==> X is unorm
 931 // DOUBLE 0x10033  exponent limit at which y is an integer
 932 // SINGLE 0x10016
 933 { .mfi
 934           getf.sig      pow_GR_sig_X        = f8
 935           fclass.m  p13,p0          = f8, 0x0b  // Test for x unorm
 936           addl pow_GR_10033                 = 0x10033, r0
 937 }
 938 { .mfi
 939           mov           pow_GR_16ones   = 0xFFFF
 940           fma.s1        POW_NORM_Y     = f9,f1,f0
 941           nop.i 999
 942 }
 943 ;;
 944
 945
 946 // p14 = TRUE ==> X is ZERO
 947 { .mfi
 948           adds          pow_AD_Tt       = pow_Tt - pow_table_P,  pow_AD_P
 949           fclass.m  p14,p15          = f8, 0x07
 950           and           pow_GR_exp_X        = pow_GR_signexp_X, pow_GR_17ones
 951 }
 952 { .mfi
 953           adds          pow_AD_Q       = pow_table_Q - pow_table_P,  pow_AD_P
 954           nop.f 999
 955           nop.i 999
 956 }
 957 ;;
 958
 959 { .mfi
 960           ldfe          POW_P5         = [pow_AD_P], 16
 961           fcmp.lt.s1 p8,p9 = f8, f0    // Test for x<0
 962           shl           pow_GR_offset       = pow_GR_sig_X, 1
 963 }
 964 { .mib
 965           ldfe          POW_P4         = [pow_AD_Q], 16
 966           sub       pow_GR_true_exp_X       = pow_GR_exp_X, pow_GR_16ones
 967 (p13)     br.cond.spnt L(POW_X_DENORM)
 968 }
 969 ;;
 970
 971
 972 // Continue normal and denormal paths here
 973 L(POW_COMMON):
 974 // p11 = TRUE ==> Y is a NAN
 975 { .mfi
 976           ldfe          POW_P3         = [pow_AD_P], 16
 977           fclass.m.unc  p11,p0         = f9, 0xc3
 978           shr.u     pow_GR_offset           = pow_GR_offset,56
 979 }
 980 { .mfi
 981           ldfe          POW_P2         = [pow_AD_Q], 16
 982           nop.f 999
 983           nop.i 999
 984 }
 985 ;;
 986
 987
 988
 989 // Compute xsq to decide later if |x|=1
 990 // p11 = TRUE ==> Y is a NaN
 991 { .mfi
 992           setf.sig POW_int_K                = pow_GR_true_exp_X
 993 (p15)     fms.s1        POW_r          = POW_B, POW_NORM_X,f1
 994           shladd pow_AD_Tt = pow_GR_offset, 4, pow_AD_Tt
 995 }
 996 { .mfi
 997           nop.m 999
 998 (p8)      fnma.s1        POW_Xm1       = POW_Xp1,f1,f0
 999           nop.i 999
1000 }
1001 ;;
1002
1003
1004
1005 // p12 = TRUE ==> X is ZERO and Y is ZERO
1006 { .mfi
1007           ldfe          POW_P1         = [pow_AD_P], 16
1008 (p14)     fclass.m.unc  p12,p0              = f9, 0x07
1009           nop.i 999
1010 }
1011 { .mfb
1012           ldfe          POW_P0         = [pow_AD_Q], 16
1013           fma.s1        POW_xsq = POW_NORM_X, POW_NORM_X, f0
1014 (p11)     br.cond.spnt   L(POW_Y_NAN)
1015 }
1016 ;;
1017
1018
1019 .pred.rel "mutex",p8,p9
1020 // Get exponent of |x|-1 to use in comparison to 2^-8
1021 { .mmf
1022 (p8)      getf.exp      pow_GR_signexp_Xm1  = POW_Xp1
1023 (p9)      getf.exp      pow_GR_signexp_Xm1  = POW_Xm1
1024           fcvt.fx.s1   POW_int_Y            = POW_NORM_Y
1025 }
1026 ;;
1027
1028
1029 // p11 = TRUE ==> X is a NAN
1030 { .mfi
1031           ldfpd         POW_log2_hi, POW_log2_lo  = [pow_AD_Q], 16
1032           fclass.m.unc  p11,p0              = f8, 0xc3
1033           nop.i 999
1034 }
1035 { .mib
1036           ldfpd  POW_T, POW_Tt              = [pow_AD_Tt], 16
1037           nop.i 999
1038 (p12)     br.cond.spnt L(POW_X_0_Y_0)
1039 }
1040 ;;
1041
1042
1043 // p14 = TRUE ==> X is zero
1044 //    p15 = TRUE ==> X is zero AND Y is negative
1045 //    p10 = TRUE ==> X is zero AND Y is >= zero
1046 { .mfi
1047           ldfe          POW_inv_log2_by_128 = [pow_AD_P], 16
1048 (p14)     fcmp.lt.unc.s1 p15, p10           = f9,f0
1049           nop.i 999
1050 }
1051 { .mfi
1052           nop.m 999
1053           nop.f 999
1054           and       pow_GR_exp_Xm1          = pow_GR_signexp_Xm1, pow_GR_17ones
1055 }
1056 ;;
1057
1058
1059 // Determine if we will use the |x| near 1 path (p6) or normal path (p7)
1060 // p12 = TRUE ==> X is a NAN and Y is a zero
1061 // p13 = TRUE ==> X is a NAN and Y is anything else
1062 { .mfi
1063           getf.exp  pow_GR_signexp_Y        = POW_NORM_Y
1064 (p11)     fclass.m.unc  p12,p13             = f9, 0x07
1065           cmp.lt.unc p6,p7                  = pow_GR_exp_Xm1, pow_GR_FFF7
1066 }
1067 { .mfi
1068           ldfpd  POW_Q2, POW_Q3             = [pow_AD_P], 16
1069           fma.s1        POW_rsq             = POW_r, POW_r,f0
1070           nop.i 999
1071 ;;
1072 }
1073
1074 // If on the x near 1 path, assign r1 to r and r1*r1 to rsq
1075 { .mfi
1076           ldfpd  POW_Q0_half, POW_Q1             = [pow_AD_P], 16
1077 (p6)      fma.s1    POW_r                 = POW_r1, f1, f0
1078           nop.i 999
1079 }
1080 { .mfi
1081           nop.m 999
1082 (p6)      fma.s1    POW_rsq                 = POW_r1, POW_r1, f0
1083           nop.i 999
1084 ;;
1085 }
1086
1087
1088 { .mfi
1089           ldfpd   POW_Q4, POW_RSHF          = [pow_AD_P], 16
1090 (p7)      fma.s1 POW_v6                     = POW_r,  POW_P5, POW_P4
1091           and pow_GR_exp_Y                   = pow_GR_signexp_Y, pow_GR_17ones
1092 }
1093 { .mfb
1094           nop.m 999
1095 (p6)      fma.s1 POW_v6                     = POW_r1, POW_P5, POW_P4
1096 (p12)     br.cond.spnt L(POW_X_NAN_Y_0)
1097 }
1098 ;;
1099
1100
1101 { .mfi
1102           nop.m 999
1103 (p7)      fma.s1 POW_v4                     = POW_P3, POW_r,  POW_P2
1104           andcm pow_GR_sign_Y               = pow_GR_signexp_Y, pow_GR_17ones
1105 }
1106 { .mfb
1107           nop.m 999
1108 (p6)      fma.s1 POW_v4                     = POW_P3, POW_r1, POW_P2
1109 (p12)     br.cond.spnt L(POW_X_NAN_Y_0)
1110 }
1111 ;;
1112
1113 { .mfi
1114           nop.m 999
1115           fcvt.xf POW_K                     = POW_int_K
1116           nop.i 999
1117 }
1118 { .mfb
1119           nop.m 999
1120 (p13)     fma.d f8                           = f8,f1,f0
1121 (p13)     br.ret.spnt  b0    // Exit if x nan, y anything but zero
1122 }
1123 ;;
1124
1125 // p10 = TRUE ==> X is zero  AND Y is positive
1126 //  p8  = TRUE ==> X is zero  AND Y is outside integer range (treat as even int)
1127 //                   return +0
1128 //  p9  = TRUE ==> X is zero  AND Y is within integer range (may not be integer)
1129 { .mfi
1130 (p10)     cmp.gt.unc p8,p9                  =  pow_GR_exp_Y, pow_GR_10033
1131 (p6)      fmerge.s POW_delta                 = f0,f0
1132           nop.i 999
1133 }
1134 { .mfi
1135           nop.m 999
1136 (p6)      fma.s1 POW_G                       = f0,f0,f0
1137           nop.i 999
1138 }
1139 ;;
1140
1141 { .mfi
1142           getf.sig pow_GR_sig_int_Y         = POW_int_Y
1143           fnma.s1 POW_twoV                   = POW_NORM_Y, POW_rsq,f0
1144           nop.i 999
1145 }
1146 { .mfi
1147           nop.m 999
1148           fma.s1 POW_U                      = POW_NORM_Y,POW_r,f0
1149           nop.i 999
1150 }
1151 ;;
1152
1153 { .mfi
1154           ldfe      POW_log2_by_128_lo      = [pow_AD_P], 16
1155 (p6)      fma.s1 POW_v2                     = POW_P1, POW_r1, POW_P0
1156           nop.i 999
1157 }
1158 { .mfi
1159           ldfe          POW_log2_by_128_hi  = [pow_AD_Q], 16
1160 (p7)      fma.s1 POW_v2                     = POW_P1, POW_r,  POW_P0
1161           nop.i 999
1162 }
1163 ;;
1164
1165
1166 { .mfi
1167           nop.m 999
1168           fcvt.xf   POW_float_int_Y               = POW_int_Y
1169           nop.i 999
1170 }
1171 { .mfi
1172           nop.m 999
1173           fma.s1 POW_v3                     = POW_v6, POW_rsq,  POW_v4
1174           adds          pow_AD_tbl1       = pow_tbl1 - pow_Tt,  pow_AD_Q
1175 }
1176 ;;
1177
1178 { .mfi
1179           nop.m 999
1180 (p7)      fma.s1 POW_delta                  = POW_K, POW_log2_lo, POW_Tt
1181           nop.i 999
1182 }
1183 { .mfi
1184           nop.m 999
1185 (p7)      fma.s1 POW_G                      = POW_K, POW_log2_hi, POW_T
1186           adds pow_AD_tbl2                  = pow_tbl2 - pow_tbl1,  pow_AD_tbl1
1187 }
1188 ;;
1189
1190
1191 { .mfi
1192           nop.m 999
1193           fms.s1 POW_e2                     = POW_NORM_Y, POW_r, POW_U
1194           nop.i 999
1195 }
1196 { .mfi
1197           nop.m 999
1198           fma.s1 POW_Z2                     = POW_twoV, POW_Q0_half, POW_U
1199           nop.i 999
1200 }
1201 ;;
1202
1203 // p11 = TRUE ==> X is NEGATIVE
1204 // p8  = TRUE ==> X is zero  AND Y is outside intger range (treat as even int)
1205 //                return +0
1206 { .mfi
1207           nop.m 999
1208           fclass.m.unc  p11,p0              = f8, 0x1a
1209           nop.i 999
1210 }
1211 { .mfb
1212           nop.m 999
1213 (p8)      fma.d f8                          = f0,f0,f0
1214 (p8)      br.ret.spnt b0
1215 }
1216 ;;
1217
1218 { .mfi
1219           nop.m 999
1220           fma.s1 POW_Yrcub                 = POW_rsq, POW_U, f0
1221           nop.i 999
1222 }
1223 { .mfi
1224           nop.m 999
1225           fma.s1 POW_p                      = POW_rsq, POW_v3, POW_v2
1226           nop.i 999
1227 }
1228 ;;
1229
1230
1231 // p11 = TRUE ==> X is NEGATIVE
1232 //    p12 = TRUE ==> X is NEGATIVE  AND  Y  already int
1233 //    p13 = TRUE ==> X is NEGATIVE  AND  Y possible int
1234 { .mfi
1235           nop.m 999
1236           fma.s1 POW_Z1                     = POW_NORM_Y, POW_G, f0
1237 (p11)     cmp.ge.unc  p12,p13                = pow_GR_exp_Y, pow_GR_10033
1238 }
1239 { .mfi
1240           nop.m 999
1241           fma.s1 POW_e3                     = POW_NORM_Y, POW_delta, f0
1242           nop.i 999
1243 }
1244 ;;
1245
1246 // p9  = TRUE ==> X is zero  AND Y is within integer range (may not be integer)
1247 //    p6 = TRUE ==>  X is zero  AND  Y is an integer (may be even or odd)
1248 //    p7 = TRUE ==>  X is zero  AND  Y is NOT an integer, return +0
1249 { .mfi
1250           nop.m 999
1251 (p9)      fcmp.eq.unc.s1 p6,p7             = POW_float_int_Y,  POW_NORM_Y
1252           nop.i 999
1253 }
1254 { .mfi
1255           nop.m 999
1256           fma.s1 POW_Gpr                    = POW_G, f1, POW_r
1257           nop.i 999
1258 }
1259 ;;
1260
1261 // By adding RSHF (1.1000...*2^63) we put integer part in rightmost significand
1262 { .mfi
1263           nop.m 999
1264           fma.s1 POW_W2  = POW_Z2, POW_inv_log2_by_128, POW_RSHF
1265           nop.i 999
1266 }
1267 { .mfi
1268           nop.m 999
1269           fms.s1 POW_UmZ2                   = POW_U, f1, POW_Z2
1270           nop.i 999
1271 }
1272 ;;
1273
1274
1275 // If x=0 and y>0, test y and flag denormal
1276 // p6  = TRUE ==>  X is zero  AND  Y is an integer (may be even or odd)
1277 //    p8 = TRUE ==>  X is zero  AND  Y is an odd  integer
1278 //    p9 = TRUE ==>  X is zero  AND  Y is an even integer
1279 { .mfi
1280           nop.m 999
1281 (p10)     fcmp.eq.s0 p15,p0 = f9,f0
1282 (p6)      tbit.nz.unc  p8,p9                = pow_GR_sig_int_Y,0
1283 }
1284 { .mfi
1285           nop.m 999
1286           fma.s1 POW_Z3                      = POW_p, POW_Yrcub, f0
1287           nop.i 999
1288 }
1289 ;;
1290
1291 // By adding RSHF (1.1000...*2^63) we put integer part in rightmost significand
1292 { .mfi
1293           nop.m 999
1294           fms.s1 POW_e1                     = POW_NORM_Y, POW_G, POW_Z1
1295           nop.i 999
1296 }
1297 { .mfi
1298           nop.m 999
1299           fma.s1 POW_W1  = POW_Z1, POW_inv_log2_by_128, POW_RSHF
1300           nop.i 999
1301 }
1302 ;;
1303
1304 { .mfi
1305           nop.m 999
1306 (p7)      fma.d f8  = f0,f0,f0  // Result +0 if x zero and y not integer
1307           nop.i 999
1308 }
1309 { .mfb
1310           nop.m 999
1311           fma.s1 POW_Y_Gpr                  = POW_NORM_Y, POW_Gpr, f0
1312 (p8)      br.ret.spnt b0        // Exit if x zero and y odd integer
1313 }
1314 ;;
1315
1316 // By subtracting RSHF we get rounded integer POW_N2float
1317 // p15 = TRUE ==> X_0_Y_NEG
1318 { .mfi
1319           nop.m 999
1320           fms.s1 POW_N2float  = POW_W2, f1, POW_RSHF
1321           nop.i 999
1322 }
1323 { .mfb
1324           nop.m 999
1325           fma.s1 POW_UmZ2pV                 = POW_twoV,POW_Q0_half,POW_UmZ2
1326 (p15)     br.cond.spnt L(POW_X_0_Y_NEG)
1327 }
1328 ;;
1329
1330
1331
1332 { .mfi
1333           nop.m 999
1334           fma.s1 POW_Z3sq                   = POW_Z3, POW_Z3, f0
1335           nop.i 999
1336 }
1337 { .mfb
1338           nop.m 999
1339           fma.s1 POW_v4                     = POW_Z3, POW_Q3, POW_Q2
1340 (p7)      br.ret.spnt b0     // Exit if x zero and y not an integer
1341 }
1342 ;;
1343
1344
1345
1346 // Extract rounded integer from rightmost significand of POW_W2
1347 // By subtracting RSHF we get rounded integer POW_N1float
1348 { .mfi
1349           getf.sig pow_GR_int_W2             = POW_W2
1350           fms.s1 POW_N1float  = POW_W1, f1, POW_RSHF
1351           nop.i 999
1352 }
1353 { .mfi
1354           nop.m 999
1355           fma.s1 POW_v2                     = POW_Z3, POW_Q1, POW_Q0_half
1356           nop.i 999
1357 }
1358 ;;
1359
1360
1361
1362
1363 // p13 = TRUE ==> X is NEGATIVE  AND  Y possible int
1364 //     p10 = TRUE ==> X is NEG and Y is an int
1365 //     p12 = TRUE ==> X is NEG and Y is not an int
1366 { .mfi
1367           nop.m 999
1368 (p13)     fcmp.eq.unc.s1 p10,p12             = POW_float_int_Y,  POW_NORM_Y
1369           nop.i 999
1370 }
1371 { .mfb
1372           nop.m 999
1373 (p9)      fma.d f8  = f0,f0,f0   // Result +0 if x zero and y even integer
1374 (p9)      br.ret.spnt b0    // Exit if x zero and y even integer
1375 }
1376 ;;
1377
1378
1379 { .mfi
1380           nop.m 999
1381           fnma.s1 POW_s2 = POW_N2float, POW_log2_by_128_hi, POW_Z2
1382           nop.i 999
1383 }
1384 { .mfi
1385           nop.m 999
1386           fma.s1 POW_e2                     = POW_e2,f1,POW_UmZ2pV
1387           nop.i 999
1388 }
1389 ;;
1390
1391 // Extract rounded integer from rightmost significand of POW_W1
1392 // Test if x inf
1393 { .mfi
1394           getf.sig pow_GR_int_W1             = POW_W1
1395           fclass.m.unc p15,p0 = POW_NORM_X,  0x23
1396           nop.i 999
1397 }
1398 { .mfb
1399           nop.m 999
1400           fnma.s1 POW_f2  = POW_N2float, POW_log2_by_128_lo, f1
1401 (p12)     br.cond.spnt L(POW_X_NEG_Y_NONINT)  // Branch if x neg, y not integer
1402 }
1403 ;;
1404
1405 // p12 = TRUE ==> X is NEGATIVE  AND Y is an odd integer
1406 { .mfi
1407           getf.exp pow_GR_signexp_Y_Gpr       = POW_Y_Gpr
1408           fma.s1 POW_v3                     = POW_Z3sq, POW_Q4, POW_v4
1409 (p10)     tbit.nz.unc  p12,p0                = pow_GR_sig_int_Y,0
1410 }
1411 ;;
1412
1413
1414 { .mfi
1415           add pow_GR_int_N                   = pow_GR_int_W1, pow_GR_int_W2
1416           fnma.s1 POW_f1  = POW_N1float, POW_log2_by_128_lo, f1
1417           nop.i 999
1418 }
1419 { .mfb
1420           nop.m 999
1421           fnma.s1 POW_s1  = POW_N1float, POW_log2_by_128_hi, POW_Z1
1422 (p15)     br.cond.spnt L(POW_X_INF)
1423 }
1424 ;;
1425
1426
1427 // Test x and y and flag denormal
1428 { .mfi
1429           and pow_GR_index1                  = 0x0f, pow_GR_int_N
1430           fcmp.eq.s0 p15,p0 = f8,f9
1431           shr r2                             = pow_GR_int_N, 7
1432 }
1433 { .mfi
1434           and pow_GR_exp_Y_Gpr               = pow_GR_signexp_Y_Gpr, pow_GR_17ones
1435           nop.f 999
1436           and pow_GR_index2                  = 0x70, pow_GR_int_N
1437 }
1438 ;;
1439
1440
1441
1442 { .mfi
1443           shladd pow_AD_T1                   = pow_GR_index1, 4, pow_AD_tbl1
1444           fcmp.eq.s1 p7,p0 = POW_NORM_Y, f1  // Test for y=1.0
1445           sub pow_GR_true_exp_Y_Gpr          = pow_GR_exp_Y_Gpr, pow_GR_16ones
1446 }
1447 { .mfi
1448           addl pow_int_GR_M                  = 0xFFFF, r2
1449           fma.s1  POW_e12                     = POW_e1,f1,POW_e2
1450           add pow_AD_T2                      = pow_AD_tbl2, pow_GR_index2
1451 }
1452 ;;
1453
1454
1455 { .mmi
1456           ldfe POW_T1                        = [pow_AD_T1],16
1457           setf.exp POW_2M                    = pow_int_GR_M
1458           andcm pow_GR_sign_Y_Gpr            = pow_GR_signexp_Y_Gpr, pow_GR_17ones
1459 }
1460 ;;
1461
1462
1463 { .mfb
1464           ldfe POW_T2                        = [pow_AD_T2],16
1465           fma.s1 POW_q                       = POW_Z3sq, POW_v3, POW_v2
1466 (p7)      br.ret.spnt b0        // Early exit if y=1.0, result is x
1467 }
1468 ;;
1469
1470
1471 // double: p8 TRUE ==> |Y(G + r)| >= 10
1472 // single: p8 TRUE ==> |Y(G + r)| >= 7
1473
1474 // double
1475 //     -2^10  -2^9             2^9   2^10
1476 // -----+-----+----+ ... +-----+-----+-----
1477 //  p8  |             p9             |  p8
1478 //      |     |       p10      |     |
1479 // single
1480 //     -2^7   -2^6             2^6   2^7
1481 // -----+-----+----+ ... +-----+-----+-----
1482 //  p8  |             p9             |  p8
1483 //      |     |       p10      |     |
1484
1485
1486 { .mfi
1487 (p0)      cmp.le.unc p8,p9                   = 10, pow_GR_true_exp_Y_Gpr
1488           fma.s1 POW_s                       = POW_s1, f1, POW_s2
1489           nop.i 999
1490 }
1491 { .mfi
1492           nop.m 999
1493           fma.s1 POW_f12                     = POW_f1, POW_f2,f0
1494           nop.i 999
1495 }
1496 ;;
1497
1498
1499 { .mfi
1500           nop.f 999
1501 (p9)      cmp.le.unc p0,p10                  = 9, pow_GR_true_exp_Y_Gpr
1502 }
1503 ;;
1504
1505
1506
1507 { .mfb
1508           nop.m 999
1509           fma.s1 POW_e123                    = POW_e12, f1, POW_e3
1510 (p8)      br.cond.spnt L(POW_OVER_UNDER_X_NOT_INF)
1511 }
1512 ;;
1513
1514
1515 { .mmf
1516           fma.s1 POW_q                       = POW_Z3sq, POW_q, POW_Z3
1517 }
1518 ;;
1519
1520
1521 { .mfi
1522           nop.m 999
1523           fma.s1 POW_ssq                     = POW_s, POW_s, f0
1524           nop.i 999
1525 }
1526 { .mfi
1527           nop.m 999
1528           fma.s1 POW_v4                      = POW_s, POW_Q3, POW_Q2
1529           nop.i 999
1530 }
1531 ;;
1532
1533 { .mfi
1534           nop.m 999
1535           fma.s1 POW_v2                      = POW_s, POW_Q1, POW_Q0_half
1536           nop.i 999
1537 }
1538 { .mfi
1539           nop.m 999
1540           fma.s1 POW_1ps                     = f1,f1,POW_s
1541           nop.i 999
1542 }
1543 ;;
1544
1545 { .mfi
1546           nop.m 999
1547           fma.s1 POW_f3                      = POW_e123,f1,f1
1548           nop.i 999
1549 }
1550 ;;
1551
1552 { .mfi
1553           nop.m 999
1554           fma.s1 POW_T1T2                    = POW_T1, POW_T2, f0
1555           nop.i 999
1556 }
1557 ;;
1558
1559 { .mfi
1560           nop.m 999
1561           fma.s1 POW_v3                     = POW_ssq, POW_Q4, POW_v4
1562           nop.i 999
1563 }
1564 ;;
1565
1566 { .mfi
1567           nop.m 999
1568           fma.s1 POW_v21ps                  = POW_ssq, POW_v2, POW_1ps
1569           nop.i 999
1570 }
1571 { .mfi
1572           nop.m 999
1573           fma.s1 POW_s4                     = POW_ssq, POW_ssq, f0
1574           nop.i 999
1575 }
1576 ;;
1577
1578 { .mfi
1579           nop.m 999
1580           fma.s1 POW_f123                    = POW_f12, POW_f3, f0
1581           nop.i 999
1582 }
1583 ;;
1584
1585 { .mfi
1586           nop.m 999
1587           fma.s1 POW_A                      =  POW_2M, POW_T1T2, f0
1588           nop.i 999
1589 }
1590 ;;
1591
1592
1593
1594 { .mfi
1595           nop.m 999
1596 (p12)     fmerge.s POW_f123 = f8,POW_f123  // if x neg, y odd int
1597           nop.i 999
1598 }
1599 { .mfi
1600           nop.m 999
1601 //          fma.s1 POW_es                     = POW_ssq,  POW_v3, POW_v2
1602           nop.i 999
1603 }
1604 ;;
1605
1606 { .mfi
1607           nop.m 999
1608           fma.s1 POW_es                     = POW_s4,  POW_v3, POW_v21ps
1609           nop.i 999
1610 }
1611 ;;
1612
1613
1614 { .mfi
1615           nop.m 999
1616           fma.s1 POW_A                      = POW_A, POW_f123, f0
1617           nop.i 999
1618 }
1619 { .mfi
1620           nop.m 999
1621 //          fma.s1 POW_es                     = POW_es, POW_ssq, POW_1ps
1622           nop.i 999
1623 }
1624 ;;
1625
1626
1627 { .mfi
1628           nop.m 999
1629           fma.s1 POW_A                      = POW_A, POW_es,f0
1630           nop.i 999
1631 }
1632 ;;
1633
1634
1635
1636 { .mfb
1637           nop.m 999
1638 (p10)     fma.d f8                          = POW_A, POW_q, POW_A
1639 (p10)     br.ret.sptk     b0
1640 }
1641 ;;
1642
1643
1644
1645
1646
1647 // POSSIBLE_OVER_UNDER
1648 // p6 = TRUE ==> Y negative
1649
1650 { .mfi
1651         nop.m 999
1652         fmerge.s POW_abs_A                = f0, POW_A
1653         cmp.eq.unc p0,p6                  = pow_GR_sign_Y, r0
1654 }
1655 ;;
1656
1657 { .mib
1658         nop.m 999
1659         nop.i 999
1660 (p6)    br.cond.spnt L(POW_POSSIBLE_UNDER)
1661 }
1662 ;;
1663
1664 // POSSIBLE_OVER
1665 // We got an answer.
1666 // overflow is a possibility, not a certainty
1667
1668
1669 // We define an overflow when the answer with
1670 //    WRE set
1671 //    user-defined rounding mode
1672
1673 // double
1674 // Largest double is 7FE (biased double)
1675 //                   7FE - 3FF + FFFF = 103FE
1676 // Create + largest_double_plus_ulp
1677 // Create - largest_double_plus_ulp
1678 // Calculate answer with WRE set.
1679
1680 // single
1681 // Largest single is FE (biased double)
1682 //                   FE - 7F + FFFF = 1007E
1683 // Create + largest_single_plus_ulp
1684 // Create - largest_single_plus_ulp
1685 // Calculate answer with WRE set.
1686
1687 // Cases when answer is ldn+1  are as follows:
1688 //  ldn                   ldn+1
1689 // --+----------|----------+------------
1690 //              |
1691 //    +inf          +inf      -inf
1692 //                  RN         RN
1693 //                             RZ
1694
1695
1696 // Put in s2 (td set, wre set)
1697 { .mfi
1698         mov           pow_GR_gt_ln                 = 0x103ff
1699         fsetc.s2 0x7F,0x42
1700         nop.i 999
1701 }
1702 ;;
1703
1704
1705 { .mfi
1706         setf.exp POW_gt_pln                        = pow_GR_gt_ln
1707         fma.d.s2 POW_wre_urm_f8                    = POW_abs_A, POW_q, POW_abs_A
1708         nop.i 999 ;;
1709 }
1710
1711 // Return s2 to default
1712 { .mfi
1713         nop.m 999
1714         fsetc.s2 0x7F,0x40
1715         nop.i 999
1716 }
1717 ;;
1718
1719
1720 // p7 = TRUE ==> yes, we have an overflow
1721 { .mfi
1722         nop.m 999
1723         fcmp.ge.unc.s1 p7, p0                    =  POW_wre_urm_f8, POW_gt_pln
1724         nop.i 999
1725 }
1726 ;;
1727
1728
1729
1730 { .mfb
1731 (p7)   mov pow_GR_tag                            = 24
1732        fma.d f8                                  = POW_A, POW_q, POW_A
1733 (p7)   br.cond.spnt __libm_error_region
1734 }
1735 { .mfb
1736        nop.m 999
1737        nop.f 999
1738 (p0)   br.ret.sptk     b0
1739 }
1740 ;;
1741
1742
1743 L(POW_POSSIBLE_UNDER):
1744 // We got an answer. input was < -2^9 but > -2^10 (double)
1745 // We got an answer. input was < -2^6 but > -2^7  (float)
1746 // underflow is a possibility, not a certainty
1747
1748 // We define an underflow when the answer with
1749 //    ftz set
1750 // is zero (tiny numbers become zero)
1751 // Notice (from below) that if we have an unlimited exponent range,
1752 // then there is an extra machine number E between the largest denormal and
1753 // the smallest normal.
1754 // So if with unbounded exponent we round to E or below, then we are
1755 // tiny and underflow has occurred.
1756 // But notice that you can be in a situation where we are tiny, namely
1757 // rounded to E, but when the exponent is bounded we round to smallest
1758 // normal. So the answer can be the smallest normal with underflow.
1759 //                           E
1760 // -----+--------------------+--------------------+-----
1761 //      |                    |                    |
1762 //   1.1...10 2^-3fff    1.1...11 2^-3fff    1.0...00 2^-3ffe
1763 //   0.1...11 2^-3ffe                                   (biased, 1)
1764 //    largest dn                               smallest normal
1765
1766
1767 // Put in s2 (td set, ftz set)
1768 { .mfi
1769         nop.m 999
1770         fsetc.s2 0x7F,0x41
1771         nop.i 999
1772 }
1773 ;;
1774
1775
1776
1777 { .mfi
1778         nop.m 999
1779         fma.d.s2 POW_ftz_urm_f8                    = POW_A, POW_q, POW_A
1780         nop.i 999
1781 }
1782 ;;
1783
1784
1785 // Return s2 to default
1786 { .mfi
1787         nop.m 999
1788         fsetc.s2 0x7F,0x40
1789         nop.i 999
1790 }
1791 ;;
1792
1793
1794 // p7 = TRUE ==> yes, we have an underflow
1795 { .mfi
1796         nop.m 999
1797         fcmp.eq.unc.s1 p7, p0                     =  POW_ftz_urm_f8, f0
1798         nop.i 999
1799 }
1800 ;;
1801
1802
1803
1804
1805 { .mfb
1806 (p7)    mov pow_GR_tag                           = 25
1807         fma.d f8                                 = POW_A, POW_q, POW_A
1808 (p7)    br.cond.spnt __libm_error_region
1809 }
1810 ;;
1811
1812
1813 { .mfb
1814         nop.m 999
1815         nop.f 999
1816         br.ret.sptk     b0
1817 }
1818 ;;
1819
1820
1821 L(POW_X_DENORM):
1822 // Here if x unorm. Use the NORM_X for getf instructions, and the back
1823 // to normal path
1824 { .mfi
1825         getf.exp      pow_GR_signexp_X    = POW_NORM_X
1826         nop.f 999
1827         nop.i 999
1828 }
1829 ;;
1830
1831 { .mfi
1832         getf.sig      pow_GR_sig_X        = POW_NORM_X
1833         nop.f 999
1834         nop.i 999
1835 }
1836 ;;
1837
1838 { .mfi
1839         and           pow_GR_exp_X        = pow_GR_signexp_X, pow_GR_17ones
1840         nop.f 999
1841 }
1842 ;;
1843
1844 { .mib
1845         sub       pow_GR_true_exp_X       = pow_GR_exp_X, pow_GR_16ones
1846         shl           pow_GR_offset       = pow_GR_sig_X, 1
1847         br.cond.sptk    L(POW_COMMON)
1848 }
1849 ;;
1850
1851
1852 L(POW_X_0_Y_0):
1853 // When X is +-0 and Y is +-0, IEEE returns 1.0
1854 // We call error support with this value
1855
1856 { .mfb
1857          mov pow_GR_tag                     = 26
1858          fma.d f8                           = f1,f1,f0
1859          br.cond.sptk __libm_error_region
1860 }
1861 ;;
1862
1863
1864
1865
1866 L(POW_X_INF):
1867 // When X is +-inf and Y is +-, IEEE returns
1868
1869 // overflow
1870 // X +inf  Y +inf             +inf
1871 // X -inf  Y +inf             +inf
1872
1873 // X +inf  Y >0               +inf
1874 // X -inf  Y >0, !odd integer +inf     <== (-inf)^0.5 = +inf !!
1875 // X -inf  Y >0,  odd integer  -inf
1876
1877 // underflow
1878 // X +inf  Y -inf             +0
1879 // X -inf  Y -inf             +0
1880
1881 // X +inf  Y <0               +0
1882 // X -inf  Y <0, !odd integer +0
1883 // X -inf  Y <0, odd integer  -0
1884
1885 // X + inf Y=+0                +1
1886 // X + inf Y=-0                +1
1887 // X - inf Y=+0                +1
1888 // X - inf Y=-0                +1
1889
1890 // p13 == Y negative
1891 // p14 == Y positive
1892
1893 // p6 == Y is a floating point number outside the integer.
1894 //       Hence it is an integer and is even.
1895 //       p13 == (Y negative)
1896 //          return +inf
1897 //       p14 == (Y positive)
1898 //          return +0
1899
1900
1901
1902 // p7 == Y is a floating point number within the integer range.
1903 //      p9  == (int_Y = NORM_Y), Y is an integer, which may be odd or even.
1904 //           p11 odd
1905 //              p13 == (Y negative)
1906 //                 return (sign_of_x)inf
1907 //              p14 == (Y positive)
1908 //                 return (sign_of_x)0
1909 //           pxx even
1910 //              p13 == (Y negative)
1911 //                 return +inf
1912 //              p14 == (Y positive)
1913 //                 return +0
1914
1915 //      pxx == Y is not an integer
1916 //           p13 == (Y negative)
1917 //                 return +inf
1918 //           p14 == (Y positive)
1919 //                 return +0
1920 //
1921
1922 // If x=inf, test y and flag denormal
1923 { .mfi
1924           nop.m 999
1925           fcmp.eq.s0 p10,p11 = f9,f0
1926           nop.i 999
1927 }
1928 ;;
1929
1930 { .mfi
1931           nop.m 999
1932           fcmp.lt p13,p14                    = POW_NORM_Y,f0
1933           cmp.gt.unc  p6,p7                  = pow_GR_exp_Y, pow_GR_10033
1934 }
1935 { .mfi
1936           nop.m 999
1937           fclass.m p12,p0                    = f9, 0x23
1938           nop.i 999
1939 }
1940 ;;
1941
1942
1943 { .mfi
1944           nop.m 999
1945           fclass.m p15,p0                    = f9, 0x07 //@zero
1946           nop.i 999
1947 }
1948 ;;
1949
1950 { .mfb
1951           nop.m 999
1952 (p15)     fmerge.s f8 = f1,f1
1953 (p15)     br.ret.spnt b0
1954 }
1955 ;;
1956
1957
1958 { .mfi
1959 (p13)     mov pow_GR_tag                     = 25
1960 (p14)     frcpa.s1 f8,p10                       = f1,f0
1961           nop.i 999
1962 }
1963 { .mfb
1964 (p14)     mov pow_GR_tag                     = 24
1965 (p13)     fma.s1 f8                          = f0,f0,f0
1966 (p12)     br.ret.spnt b0
1967 }
1968 ;;
1969
1970
1971
1972 { .mfb
1973           nop.m 999
1974 (p7)      fcmp.eq.unc.s1 p9,p0              = POW_float_int_Y,  POW_NORM_Y
1975           nop.b 999
1976 }
1977 ;;
1978
1979 { .mfi
1980           nop.m 999
1981           nop.f 999
1982 (p9)      tbit.nz.unc p11,p0                 = pow_GR_sig_int_Y,0
1983 }
1984 ;;
1985
1986 { .mfb
1987           nop.m 999
1988 (p11)     fmerge.s f8 = POW_NORM_X,f8
1989           br.ret.sptk b0
1990 }
1991 ;;
1992
1993
1994
1995 L(POW_X_0_Y_NEG):
1996 // When X is +-0 and Y is negative, IEEE returns
1997 // X     Y           answer
1998 // +0    -odd int    +inf
1999 // -0    -odd int    -inf
2000
2001 // +0    !-odd int   +inf
2002 // -0    !-odd int   +inf
2003
2004
2005 // p6 == Y is a floating point number outside the integer.
2006 //       Hence it is an integer and is even.
2007 //       return +inf
2008
2009 // p7 == Y is a floating point number within the integer range.
2010 //      p9  == (int_Y = NORM_Y), Y is an integer, which may be odd or even.
2011 //           p11 odd
2012 //              return (sign_of_x)inf
2013 //           p12 even
2014 //              return +inf
2015 //      p10 == Y is not an integer
2016 //         return +inf
2017 //
2018 //
2019
2020 { .mfi
2021           nop.m 999
2022           nop.f 999
2023           cmp.gt.unc  p6,p7                  = pow_GR_exp_Y, pow_GR_10033
2024 }
2025 ;;
2026
2027
2028 { .mfi
2029           mov pow_GR_tag                     = 27
2030 (p7)      fcmp.eq.unc.s1 p9,p10              = POW_float_int_Y,  POW_NORM_Y
2031           nop.i 999
2032 }
2033 ;;
2034
2035
2036 { .mfb
2037           nop.m 999
2038 (p6)      frcpa.s0 f8,p13                       = f1, f0
2039 (p6)      br.cond.sptk __libm_error_region
2040 }
2041 ;;
2042
2043 { .mfb
2044           nop.m 999
2045 (p10)     frcpa.s0 f8,p13                       = f1, f0
2046 (p10)     br.cond.sptk __libm_error_region
2047 }
2048 ;;
2049
2050
2051
2052 { .mib
2053           nop.m 999
2054 (p9)      tbit.nz.unc p11,p12                = pow_GR_sig_int_Y,0
2055           nop.b 999
2056 }
2057 ;;
2058
2059
2060
2061 { .mfi
2062           nop.m 999
2063 (p12)     frcpa.s0 f8,p13                      = f1,f0
2064           nop.i 999
2065 }
2066 ;;
2067
2068 { .mfb
2069           nop.m 999
2070 (p11)     frcpa f8,p13                      = f1,f8
2071           br.cond.sptk __libm_error_region
2072 }
2073 ;;
2074
2075
2076
2077
2078 L(POW_X_NEG_Y_NONINT):
2079 // When X is negative and Y is a non-integer, IEEE
2080 // returns a qnan indefinite.
2081 // We call error support with this value
2082
2083 { .mfb
2084          mov pow_GR_tag                     = 28
2085          frcpa f8,p6                        = f0,f0
2086          br.cond.sptk __libm_error_region
2087 }
2088 ;;
2089
2090
2091
2092
2093 L(POW_X_NAN_Y_0):
2094 // When X is a NAN and Y is zero, IEEE returns 1.
2095 // We call error support with this value.
2096
2097 { .mfi
2098          nop.m 0
2099          fma.d.s0 f10 = f8,f1,f0
2100          nop.i 0
2101 }
2102 { .mfb
2103          mov pow_GR_tag                     = 29
2104          fma.d.s0 f8 = f0,f0,f1
2105          br.cond.sptk __libm_error_region
2106 }
2107 ;;
2108
2109
2110 L(POW_OVER_UNDER_X_NOT_INF):
2111
2112 // p8 is TRUE for overflow
2113 // p9 is TRUE for underflow
2114
2115 // if y is infinity, we should not over/underflow
2116
2117
2118 { .mfi
2119           nop.m 999
2120           fcmp.eq.unc.s1     p14, p13        = POW_xsq,f1
2121           cmp.eq.unc p8,p9                   = pow_GR_sign_Y_Gpr, r0
2122 }
2123 ;;
2124
2125 { .mfi
2126           nop.m 999
2127 (p14)     fclass.m.unc       p15, p0         = f9, 0x23
2128           nop.i 999
2129 }
2130 { .mfi
2131           nop.m 999
2132 (p13)     fclass.m.unc       p11,p0         = f9, 0x23
2133           nop.i 999
2134 }
2135 ;;
2136
2137 // p15 = TRUE if |x|=1, y=inf, return +1
2138 { .mfb
2139           nop.m 999
2140 (p15)     fma.d              f8              = f1,f1,f0
2141 (p15)     br.ret.spnt b0
2142 }
2143 ;;
2144
2145 .pred.rel "mutex",p8,p9
2146 {  .mfb
2147 (p8)      setf.exp           f8              = pow_GR_17ones
2148 (p9)      fmerge.s           f8              = f0,f0
2149 (p11)     br.ret.sptk b0
2150 }
2151
2152 { .mfb
2153           nop.m 999
2154           nop.f 999
2155           br.cond.sptk L(POW_OVER_UNDER_ERROR)
2156 }
2157 ;;
2158
2159 L(POW_Y_NAN):
2160
2161 // Is x = +1 then result is +1, else result is quiet Y
2162 { .mfi
2163        nop.m 999
2164        fcmp.eq.s1         p10,p9               = POW_NORM_X, f1
2165        nop.i 999
2166 }
2167 ;;
2168
2169 { .mfi
2170        nop.m 999
2171 (p10)  fcmp.eq.s0 p6,p0 = f9,f1   // Set invalid, even if x=+1
2172        nop.i 999
2173 }
2174 ;;
2175
2176 { .mfi
2177        nop.m 999
2178 (p10)  fma.d f8 = f1,f1,f0
2179        nop.i 999
2180 }
2181 { .mfb
2182        nop.m 999
2183 (p9)   fma.d f8 = f9,f8,f0
2184        br.ret.sptk b0
2185 }
2186 ;;
2187
2188
2189 L(POW_OVER_UNDER_ERROR):
2190
2191 { .mfi
2192           nop.m 999
2193           fmerge.s f10                      = POW_NORM_X,POW_NORM_X
2194           nop.i 999
2195 }
2196 { .mfi
2197           sub   pow_GR_17ones_m1            = pow_GR_17ones, r0, 1
2198           nop.f 999
2199           mov pow_GR_one                    = 0x1
2200 }
2201 ;;
2202
2203 // overflow
2204 { .mmb
2205 (p8)     mov pow_GR_tag                     = 24
2206 (p8)     setf.exp f11                       = pow_GR_17ones_m1
2207          nop.b 999
2208 }
2209 ;;
2210
2211
2212 // underflow
2213 { .mmi
2214 (p9)    mov pow_GR_tag                     = 25
2215 (p9)    setf.exp f11                       = pow_GR_one
2216         nop.i 999
2217 }
2218 ;;
2219
2220
2221 // p12 x is negative and y is an odd integer
2222
2223
2224 { .mfi
2225         nop.m 999
2226         fma.d f8                               = f11, f11, f0
2227         nop.i 999
2228 }
2229 ;;
2230
2231 { .mfi
2232         nop.m 999
2233 (p12)   fmerge.ns f8                           = f8, f8
2234         nop.i 999
2235 }
2236 ;;
2237
2238
2239 .endp pow
2240 ASM_SIZE_DIRECTIVE(pow)
2241
2242
2243 // Stack operations when calling error support.
2244 //       (1)               (2)                          (3) (call)              (4)
2245 //   sp   -> +          psp -> +                     psp -> +                   sp -> +
2246 //           |                 |                            |                         |
2247 //           |                 | <- GR_Y               R3 ->| <- GR_RESULT            | -> f8
2248 //           |                 |                            |                         |
2249 //           | <-GR_Y      Y2->|                       Y2 ->| <- GR_Y                 |
2250 //           |                 |                            |                         |
2251 //           |                 | <- GR_X               X1 ->|                         |
2252 //           |                 |                            |                         |
2253 //  sp-64 -> +          sp ->  +                     sp ->  +                         +
2254 //    save ar.pfs          save b0                                               restore gp
2255 //    save gp                                                                    restore ar.pfs
2256
2257
2258
2259 .proc __libm_error_region
2260 __libm_error_region:
2261
2262 // Answer is inf for overflow and 0 for underflow.
2263 .prologue
2264 // (1)
2265 { .mfi
2266         add   GR_Parameter_Y=-32,sp             // Parameter 2 value
2267         nop.f 0
2268 .save   ar.pfs,GR_SAVE_PFS
2269         mov  GR_SAVE_PFS=ar.pfs                 // Save ar.pfs
2270 }
2271 { .mfi
2272 .fframe 64
2273         add sp=-64,sp                          // Create new stack
2274         nop.f 0
2275         mov GR_SAVE_GP=gp                      // Save gp
2276 };;
2277
2278
2279 // (2)
2280 { .mmi
2281         stfd [GR_Parameter_Y] = POW_NORM_Y,16 // STORE Parameter 2 on stack
2282         add GR_Parameter_X = 16,sp            // Parameter 1 address
2283 .save   b0, GR_SAVE_B0
2284         mov GR_SAVE_B0=b0                     // Save b0
2285 };;
2286
2287 .body
2288 // (3)
2289 { .mib
2290         stfd [GR_Parameter_X] = POW_NORM_X              // STORE Parameter 1 on stack
2291         add   GR_Parameter_RESULT = 0,GR_Parameter_Y    // Parameter 3 address
2292         nop.b 0
2293 }
2294 { .mib
2295         stfd [GR_Parameter_Y] = f8                      // STORE Parameter 3 on stack
2296         add   GR_Parameter_Y = -16,GR_Parameter_Y
2297         br.call.sptk b0=__libm_error_support#           // Call error handling function
2298 };;
2299 { .mmi
2300         nop.m 0
2301         nop.m 0
2302         add   GR_Parameter_RESULT = 48,sp
2303 };;
2304
2305 // (4)
2306 { .mmi
2307         ldfd  f8 = [GR_Parameter_RESULT]       // Get return result off stack
2308 .restore sp
2309         add   sp = 64,sp                       // Restore stack pointer
2310         mov   b0 = GR_SAVE_B0                  // Restore return address
2311 };;
2312 { .mib
2313         mov   gp = GR_SAVE_GP                  // Restore gp
2314         mov   ar.pfs = GR_SAVE_PFS             // Restore ar.pfs
2315         br.ret.sptk     b0                     // Return
2316 };;
2317
2318 .endp __libm_error_region
2319 ASM_SIZE_DIRECTIVE(__libm_error_region)
2320
2321 .type   __libm_error_support#,@function
2322 .global __libm_error_support#