src/compiler/glsl/float64.glsl

   1 /*
   2  * The implementations contained in this file are heavily based on the
   3  * implementations found in the Berkeley SoftFloat library. As such, they are
   4  * licensed under the same 3-clause BSD license:
   5  *
   6  * License for Berkeley SoftFloat Release 3e
   7  *
   8  * John R. Hauser
   9  * 2018 January 20
  10  *
  11  * The following applies to the whole of SoftFloat Release 3e as well as to
  12  * each source file individually.
  13  *
  14  * Copyright 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018 The Regents of the
  15  * University of California.  All rights reserved.
  16  *
  17  * Redistribution and use in source and binary forms, with or without
  18  * modification, are permitted provided that the following conditions are met:
  19  *
  20  *  1. Redistributions of source code must retain the above copyright notice,
  21  *     this list of conditions, and the following disclaimer.
  22  *
  23  *  2. Redistributions in binary form must reproduce the above copyright
  24  *     notice, this list of conditions, and the following disclaimer in the
  25  *     documentation and/or other materials provided with the distribution.
  26  *
  27  *  3. Neither the name of the University nor the names of its contributors
  28  *     may be used to endorse or promote products derived from this software
  29  *     without specific prior written permission.
  30  *
  31  * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS "AS IS", AND ANY
  32  * EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
  33  * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, ARE
  34  * DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE FOR ANY
  35  * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
  36  * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
  37  * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
  38  * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
  39  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
  40  * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  41 */
  42
  43 #version 430
  44 #extension GL_ARB_gpu_shader_int64 : enable
  45 #extension GL_ARB_shader_bit_encoding : enable
  46 #extension GL_EXT_shader_integer_mix : enable
  47 #extension GL_MESA_shader_integer_functions : enable
  48
  49 #pragma warning(off)
  50
  51 /* Software IEEE floating-point rounding mode.
  52  * GLSL spec section "4.7.1 Range and Precision":
  53  * The rounding mode cannot be set and is undefined.
  54  * But here, we are able to define the rounding mode at the compilation time.
  55  */
  56 #define FLOAT_ROUND_NEAREST_EVEN    0
  57 #define FLOAT_ROUND_TO_ZERO         1
  58 #define FLOAT_ROUND_DOWN            2
  59 #define FLOAT_ROUND_UP              3
  60 #define FLOAT_ROUNDING_MODE         FLOAT_ROUND_NEAREST_EVEN
  61
  62 /* Absolute value of a Float64 :
  63  * Clear the sign bit
  64  */
  65 uint64_t
  66 __fabs64(uint64_t __a)
  67 {
  68    uvec2 a = unpackUint2x32(__a);
  69    a.y &= 0x7FFFFFFFu;
  70    return packUint2x32(a);
  71 }
  72
  73 /* Returns 1 if the double-precision floating-point value `a' is a NaN;
  74  * otherwise returns 0.
  75  */
  76 bool
  77 __is_nan(uint64_t __a)
  78 {
  79    uvec2 a = unpackUint2x32(__a);
  80    return (0xFFE00000u <= (a.y<<1)) &&
  81       ((a.x != 0u) || ((a.y & 0x000FFFFFu) != 0u));
  82 }
  83
  84 /* Negate value of a Float64 :
  85  * Toggle the sign bit
  86  */
  87 uint64_t
  88 __fneg64(uint64_t __a)
  89 {
  90    uvec2 a = unpackUint2x32(__a);
  91    a.y ^= (1u << 31);
  92    return packUint2x32(a);
  93 }
  94
  95 uint64_t
  96 __fsign64(uint64_t __a)
  97 {
  98    uvec2 a = unpackUint2x32(__a);
  99    uvec2 retval;
 100    retval.x = 0u;
 101    retval.y = mix((a.y & 0x80000000u) | 0x3FF00000u, 0u, (a.y << 1 | a.x) == 0u);
 102    return packUint2x32(retval);
 103 }
 104
 105 /* Returns the fraction bits of the double-precision floating-point value `a'.*/
 106 uint
 107 __extractFloat64FracLo(uint64_t a)
 108 {
 109    return unpackUint2x32(a).x;
 110 }
 111
 112 uint
 113 __extractFloat64FracHi(uint64_t a)
 114 {
 115    return unpackUint2x32(a).y & 0x000FFFFFu;
 116 }
 117
 118 /* Returns the exponent bits of the double-precision floating-point value `a'.*/
 119 int
 120 __extractFloat64Exp(uint64_t __a)
 121 {
 122    uvec2 a = unpackUint2x32(__a);
 123    return int((a.y>>20) & 0x7FFu);
 124 }
 125
 126 bool
 127 __feq64_nonnan(uint64_t __a, uint64_t __b)
 128 {
 129    uvec2 a = unpackUint2x32(__a);
 130    uvec2 b = unpackUint2x32(__b);
 131    return (a.x == b.x) &&
 132           ((a.y == b.y) || ((a.x == 0u) && (((a.y | b.y)<<1) == 0u)));
 133 }
 134
 135 /* Returns true if the double-precision floating-point value `a' is equal to the
 136  * corresponding value `b', and false otherwise.  The comparison is performed
 137  * according to the IEEE Standard for Floating-Point Arithmetic.
 138  */
 139 bool
 140 __feq64(uint64_t a, uint64_t b)
 141 {
 142    if (__is_nan(a) || __is_nan(b))
 143       return false;
 144
 145    return __feq64_nonnan(a, b);
 146 }
 147
 148 /* Returns true if the double-precision floating-point value `a' is not equal
 149  * to the corresponding value `b', and false otherwise.  The comparison is
 150  * performed according to the IEEE Standard for Floating-Point Arithmetic.
 151  */
 152 bool
 153 __fne64(uint64_t a, uint64_t b)
 154 {
 155    if (__is_nan(a) || __is_nan(b))
 156       return true;
 157
 158    return !__feq64_nonnan(a, b);
 159 }
 160
 161 /* Returns the sign bit of the double-precision floating-point value `a'.*/
 162 uint
 163 __extractFloat64Sign(uint64_t a)
 164 {
 165    return unpackUint2x32(a).y & 0x80000000u;
 166 }
 167
 168 /* Returns true if the 64-bit value formed by concatenating `a0' and `a1' is less
 169  * than the 64-bit value formed by concatenating `b0' and `b1'.  Otherwise,
 170  * returns false.
 171  */
 172 bool
 173 lt64(uint a0, uint a1, uint b0, uint b1)
 174 {
 175    return (a0 < b0) || ((a0 == b0) && (a1 < b1));
 176 }
 177
 178 bool
 179 __flt64_nonnan(uint64_t __a, uint64_t __b)
 180 {
 181    uvec2 a = unpackUint2x32(__a);
 182    uvec2 b = unpackUint2x32(__b);
 183    uint aSign = __extractFloat64Sign(__a);
 184    uint bSign = __extractFloat64Sign(__b);
 185    if (aSign != bSign)
 186       return (aSign != 0u) && ((((a.y | b.y)<<1) | a.x | b.x) != 0u);
 187
 188    return mix(lt64(a.y, a.x, b.y, b.x), lt64(b.y, b.x, a.y, a.x), aSign != 0u);
 189 }
 190
 191 /* Returns true if the double-precision floating-point value `a' is less than
 192  * the corresponding value `b', and false otherwise.  The comparison is performed
 193  * according to the IEEE Standard for Floating-Point Arithmetic.
 194  */
 195 bool
 196 __flt64(uint64_t a, uint64_t b)
 197 {
 198    if (__is_nan(a) || __is_nan(b))
 199       return false;
 200
 201    return __flt64_nonnan(a, b);
 202 }
 203
 204 /* Returns true if the double-precision floating-point value `a' is greater
 205  * than or equal to * the corresponding value `b', and false otherwise.  The
 206  * comparison is performed * according to the IEEE Standard for Floating-Point
 207  * Arithmetic.
 208  */
 209 bool
 210 __fge64(uint64_t a, uint64_t b)
 211 {
 212    if (__is_nan(a) || __is_nan(b))
 213       return false;
 214
 215    return !__flt64_nonnan(a, b);
 216 }
 217
 218 uint64_t
 219 __fsat64(uint64_t __a)
 220 {
 221    if (__flt64(__a, 0ul))
 222       return 0ul;
 223
 224    if (__fge64(__a, 0x3FF0000000000000ul /* 1.0 */))
 225       return 0x3FF0000000000000ul;
 226
 227    return __a;
 228 }
 229
 230 /* Adds the 64-bit value formed by concatenating `a0' and `a1' to the 64-bit
 231  * value formed by concatenating `b0' and `b1'.  Addition is modulo 2^64, so
 232  * any carry out is lost.  The result is broken into two 32-bit pieces which
 233  * are stored at the locations pointed to by `z0Ptr' and `z1Ptr'.
 234  */
 235 void
 236 __add64(uint a0, uint a1, uint b0, uint b1,
 237         out uint z0Ptr,
 238         out uint z1Ptr)
 239 {
 240    uint z1 = a1 + b1;
 241    z1Ptr = z1;
 242    z0Ptr = a0 + b0 + uint(z1 < a1);
 243 }
 244
 245
 246 /* Subtracts the 64-bit value formed by concatenating `b0' and `b1' from the
 247  * 64-bit value formed by concatenating `a0' and `a1'.  Subtraction is modulo
 248  * 2^64, so any borrow out (carry out) is lost.  The result is broken into two
 249  * 32-bit pieces which are stored at the locations pointed to by `z0Ptr' and
 250  * `z1Ptr'.
 251  */
 252 void
 253 __sub64(uint a0, uint a1, uint b0, uint b1,
 254         out uint z0Ptr,
 255         out uint z1Ptr)
 256 {
 257    z1Ptr = a1 - b1;
 258    z0Ptr = a0 - b0 - uint(a1 < b1);
 259 }
 260
 261 /* Shifts the 64-bit value formed by concatenating `a0' and `a1' right by the
 262  * number of bits given in `count'.  If any nonzero bits are shifted off, they
 263  * are "jammed" into the least significant bit of the result by setting the
 264  * least significant bit to 1.  The value of `count' can be arbitrarily large;
 265  * in particular, if `count' is greater than 64, the result will be either 0
 266  * or 1, depending on whether the concatenation of `a0' and `a1' is zero or
 267  * nonzero.  The result is broken into two 32-bit pieces which are stored at
 268  * the locations pointed to by `z0Ptr' and `z1Ptr'.
 269  */
 270 void
 271 __shift64RightJamming(uint a0,
 272                       uint a1,
 273                       int count,
 274                       out uint z0Ptr,
 275                       out uint z1Ptr)
 276 {
 277    uint z0;
 278    uint z1;
 279    int negCount = (-count) & 31;
 280
 281    z0 = mix(0u, a0, count == 0);
 282    z0 = mix(z0, (a0 >> count), count < 32);
 283
 284    z1 = uint((a0 | a1) != 0u); /* count >= 64 */
 285    uint z1_lt64 = (a0>>(count & 31)) | uint(((a0<<negCount) | a1) != 0u);
 286    z1 = mix(z1, z1_lt64, count < 64);
 287    z1 = mix(z1, (a0 | uint(a1 != 0u)), count == 32);
 288    uint z1_lt32 = (a0<<negCount) | (a1>>count) | uint ((a1<<negCount) != 0u);
 289    z1 = mix(z1, z1_lt32, count < 32);
 290    z1 = mix(z1, a1, count == 0);
 291    z1Ptr = z1;
 292    z0Ptr = z0;
 293 }
 294
 295 /* Shifts the 96-bit value formed by concatenating `a0', `a1', and `a2' right
 296  * by 32 _plus_ the number of bits given in `count'.  The shifted result is
 297  * at most 64 nonzero bits; these are broken into two 32-bit pieces which are
 298  * stored at the locations pointed to by `z0Ptr' and `z1Ptr'.  The bits shifted
 299  * off form a third 32-bit result as follows:  The _last_ bit shifted off is
 300  * the most-significant bit of the extra result, and the other 31 bits of the
 301  * extra result are all zero if and only if _all_but_the_last_ bits shifted off
 302  * were all zero.  This extra result is stored in the location pointed to by
 303  * `z2Ptr'.  The value of `count' can be arbitrarily large.
 304  *     (This routine makes more sense if `a0', `a1', and `a2' are considered
 305  * to form a fixed-point value with binary point between `a1' and `a2'.  This
 306  * fixed-point value is shifted right by the number of bits given in `count',
 307  * and the integer part of the result is returned at the locations pointed to
 308  * by `z0Ptr' and `z1Ptr'.  The fractional part of the result may be slightly
 309  * corrupted as described above, and is returned at the location pointed to by
 310  * `z2Ptr'.)
 311  */
 312 void
 313 __shift64ExtraRightJamming(uint a0, uint a1, uint a2,
 314                            int count,
 315                            out uint z0Ptr,
 316                            out uint z1Ptr,
 317                            out uint z2Ptr)
 318 {
 319    uint z0 = 0u;
 320    uint z1;
 321    uint z2;
 322    int negCount = (-count) & 31;
 323
 324    z2 = mix(uint(a0 != 0u), a0, count == 64);
 325    z2 = mix(z2, a0 << negCount, count < 64);
 326    z2 = mix(z2, a1 << negCount, count < 32);
 327
 328    z1 = mix(0u, (a0 >> (count & 31)), count < 64);
 329    z1 = mix(z1, (a0<<negCount) | (a1>>count), count < 32);
 330
 331    a2 = mix(a2 | a1, a2, count < 32);
 332    z0 = mix(z0, a0 >> count, count < 32);
 333    z2 |= uint(a2 != 0u);
 334
 335    z0 = mix(z0, 0u, (count == 32));
 336    z1 = mix(z1, a0, (count == 32));
 337    z2 = mix(z2, a1, (count == 32));
 338    z0 = mix(z0, a0, (count == 0));
 339    z1 = mix(z1, a1, (count == 0));
 340    z2 = mix(z2, a2, (count == 0));
 341    z2Ptr = z2;
 342    z1Ptr = z1;
 343    z0Ptr = z0;
 344 }
 345
 346 /* Shifts the 64-bit value formed by concatenating `a0' and `a1' left by the
 347  * number of bits given in `count'.  Any bits shifted off are lost.  The value
 348  * of `count' must be less than 32.  The result is broken into two 32-bit
 349  * pieces which are stored at the locations pointed to by `z0Ptr' and `z1Ptr'.
 350  */
 351 void
 352 __shortShift64Left(uint a0, uint a1,
 353                    int count,
 354                    out uint z0Ptr,
 355                    out uint z1Ptr)
 356 {
 357    z1Ptr = a1<<count;
 358    z0Ptr = mix((a0 << count | (a1 >> ((-count) & 31))), a0, count == 0);
 359 }
 360
 361 /* Packs the sign `zSign', the exponent `zExp', and the significand formed by
 362  * the concatenation of `zFrac0' and `zFrac1' into a double-precision floating-
 363  * point value, returning the result.  After being shifted into the proper
 364  * positions, the three fields `zSign', `zExp', and `zFrac0' are simply added
 365  * together to form the most significant 32 bits of the result.  This means
 366  * that any integer portion of `zFrac0' will be added into the exponent.  Since
 367  * a properly normalized significand will have an integer portion equal to 1,
 368  * the `zExp' input should be 1 less than the desired result exponent whenever
 369  * `zFrac0' and `zFrac1' concatenated form a complete, normalized significand.
 370  */
 371 uint64_t
 372 __packFloat64(uint zSign, int zExp, uint zFrac0, uint zFrac1)
 373 {
 374    uvec2 z;
 375
 376    z.y = zSign + (uint(zExp) << 20) + zFrac0;
 377    z.x = zFrac1;
 378    return packUint2x32(z);
 379 }
 380
 381 /* Takes an abstract floating-point value having sign `zSign', exponent `zExp',
 382  * and extended significand formed by the concatenation of `zFrac0', `zFrac1',
 383  * and `zFrac2', and returns the proper double-precision floating-point value
 384  * corresponding to the abstract input.  Ordinarily, the abstract value is
 385  * simply rounded and packed into the double-precision format, with the inexact
 386  * exception raised if the abstract input cannot be represented exactly.
 387  * However, if the abstract value is too large, the overflow and inexact
 388  * exceptions are raised and an infinity or maximal finite value is returned.
 389  * If the abstract value is too small, the input value is rounded to a
 390  * subnormal number, and the underflow and inexact exceptions are raised if the
 391  * abstract input cannot be represented exactly as a subnormal double-precision
 392  * floating-point number.
 393  *     The input significand must be normalized or smaller.  If the input
 394  * significand is not normalized, `zExp' must be 0; in that case, the result
 395  * returned is a subnormal number, and it must not require rounding.  In the
 396  * usual case that the input significand is normalized, `zExp' must be 1 less
 397  * than the "true" floating-point exponent.  The handling of underflow and
 398  * overflow follows the IEEE Standard for Floating-Point Arithmetic.
 399  */
 400 uint64_t
 401 __roundAndPackFloat64(uint zSign,
 402                       int zExp,
 403                       uint zFrac0,
 404                       uint zFrac1,
 405                       uint zFrac2)
 406 {
 407    bool roundNearestEven;
 408    bool increment;
 409
 410    roundNearestEven = FLOAT_ROUNDING_MODE == FLOAT_ROUND_NEAREST_EVEN;
 411    increment = int(zFrac2) < 0;
 412    if (!roundNearestEven) {
 413       if (FLOAT_ROUNDING_MODE == FLOAT_ROUND_TO_ZERO) {
 414          increment = false;
 415       } else {
 416          if (zSign != 0u) {
 417             increment = (FLOAT_ROUNDING_MODE == FLOAT_ROUND_DOWN) &&
 418                (zFrac2 != 0u);
 419          } else {
 420             increment = (FLOAT_ROUNDING_MODE == FLOAT_ROUND_UP) &&
 421                (zFrac2 != 0u);
 422          }
 423       }
 424    }
 425    if (0x7FD <= zExp) {
 426       if ((0x7FD < zExp) ||
 427          ((zExp == 0x7FD) &&
 428             (0x001FFFFFu == zFrac0 && 0xFFFFFFFFu == zFrac1) &&
 429                increment)) {
 430          if ((FLOAT_ROUNDING_MODE == FLOAT_ROUND_TO_ZERO) ||
 431             ((zSign != 0u) && (FLOAT_ROUNDING_MODE == FLOAT_ROUND_UP)) ||
 432                ((zSign == 0u) && (FLOAT_ROUNDING_MODE == FLOAT_ROUND_DOWN))) {
 433             return __packFloat64(zSign, 0x7FE, 0x000FFFFFu, 0xFFFFFFFFu);
 434          }
 435          return __packFloat64(zSign, 0x7FF, 0u, 0u);
 436       }
 437       if (zExp < 0) {
 438          __shift64ExtraRightJamming(
 439             zFrac0, zFrac1, zFrac2, -zExp, zFrac0, zFrac1, zFrac2);
 440          zExp = 0;
 441          if (roundNearestEven) {
 442             increment = zFrac2 < 0u;
 443          } else {
 444             if (zSign != 0u) {
 445                increment = (FLOAT_ROUNDING_MODE == FLOAT_ROUND_DOWN) &&
 446                   (zFrac2 != 0u);
 447             } else {
 448                increment = (FLOAT_ROUNDING_MODE == FLOAT_ROUND_UP) &&
 449                   (zFrac2 != 0u);
 450             }
 451          }
 452       }
 453    }
 454    if (increment) {
 455       __add64(zFrac0, zFrac1, 0u, 1u, zFrac0, zFrac1);
 456       zFrac1 &= ~((zFrac2 + uint(zFrac2 == 0u)) & uint(roundNearestEven));
 457    } else {
 458       zExp = mix(zExp, 0, (zFrac0 | zFrac1) == 0u);
 459    }
 460    return __packFloat64(zSign, zExp, zFrac0, zFrac1);
 461 }
 462
 463 uint64_t
 464 __roundAndPackUInt64(uint zSign, uint zFrac0, uint zFrac1, uint zFrac2)
 465 {
 466    bool roundNearestEven;
 467    bool increment;
 468    uint64_t default_nan = 0xFFFFFFFFFFFFFFFFUL;
 469
 470    roundNearestEven = FLOAT_ROUNDING_MODE == FLOAT_ROUND_NEAREST_EVEN;
 471
 472    if (zFrac2 >= 0x80000000u)
 473       increment = false;
 474
 475    if (!roundNearestEven) {
 476       if (zSign != 0u) {
 477          if ((FLOAT_ROUNDING_MODE == FLOAT_ROUND_DOWN) && (zFrac2 != 0u)) {
 478             increment = false;
 479          }
 480       } else {
 481          increment = (FLOAT_ROUNDING_MODE == FLOAT_ROUND_UP) &&
 482             (zFrac2 != 0u);
 483       }
 484    }
 485
 486    if (increment) {
 487       __add64(zFrac0, zFrac1, 0u, 1u, zFrac0, zFrac1);
 488       if ((zFrac0 | zFrac1) != 0u)
 489          zFrac1 &= ~(1u) + uint(zFrac2 == 0u) & uint(roundNearestEven);
 490    }
 491    return mix(packUint2x32(uvec2(zFrac1, zFrac0)), default_nan,
 492               (zSign != 0u && (zFrac0 | zFrac1) != 0u));
 493 }
 494
 495 int64_t
 496 __roundAndPackInt64(uint zSign, uint zFrac0, uint zFrac1, uint zFrac2)
 497 {
 498    bool roundNearestEven;
 499    bool increment;
 500    int64_t default_NegNaN = -0x7FFFFFFFFFFFFFFEL;
 501    int64_t default_PosNaN = 0xFFFFFFFFFFFFFFFFL;
 502
 503    roundNearestEven = FLOAT_ROUNDING_MODE == FLOAT_ROUND_NEAREST_EVEN;
 504
 505    if (zFrac2 >= 0x80000000u)
 506       increment = false;
 507
 508    if (!roundNearestEven) {
 509       if (zSign != 0u) {
 510          increment = ((FLOAT_ROUNDING_MODE == FLOAT_ROUND_DOWN) &&
 511             (zFrac2 != 0u));
 512       } else {
 513          increment = (FLOAT_ROUNDING_MODE == FLOAT_ROUND_UP) &&
 514             (zFrac2 != 0u);
 515       }
 516    }
 517
 518    if (increment) {
 519       __add64(zFrac0, zFrac1, 0u, 1u, zFrac0, zFrac1);
 520       if ((zFrac0 | zFrac1) != 0u)
 521          zFrac1 &= ~(1u) + uint(zFrac2 == 0u) & uint(roundNearestEven);
 522    }
 523
 524    int64_t absZ = mix(int64_t(packUint2x32(uvec2(zFrac1, zFrac0))),
 525                       -int64_t(packUint2x32(uvec2(zFrac1, zFrac0))),
 526                       zSign != 0u);
 527    int64_t nan = mix(default_PosNaN, default_NegNaN, zSign != 0u);
 528    return mix(absZ, nan, ((zSign != 0u) != (absZ < 0)) && bool(absZ));
 529 }
 530
 531 /* Returns the number of leading 0 bits before the most-significant 1 bit of
 532  * `a'.  If `a' is zero, 32 is returned.
 533  */
 534 int
 535 __countLeadingZeros32(uint a)
 536 {
 537    return 31 - findMSB(a);
 538 }
 539
 540 /* Takes an abstract floating-point value having sign `zSign', exponent `zExp',
 541  * and significand formed by the concatenation of `zSig0' and `zSig1', and
 542  * returns the proper double-precision floating-point value corresponding
 543  * to the abstract input.  This routine is just like `__roundAndPackFloat64'
 544  * except that the input significand has fewer bits and does not have to be
 545  * normalized.  In all cases, `zExp' must be 1 less than the "true" floating-
 546  * point exponent.
 547  */
 548 uint64_t
 549 __normalizeRoundAndPackFloat64(uint zSign,
 550                                int zExp,
 551                                uint zFrac0,
 552                                uint zFrac1)
 553 {
 554    int shiftCount;
 555    uint zFrac2;
 556
 557    if (zFrac0 == 0u) {
 558       zExp -= 32;
 559       zFrac0 = zFrac1;
 560       zFrac1 = 0u;
 561    }
 562
 563    shiftCount = __countLeadingZeros32(zFrac0) - 11;
 564    if (0 <= shiftCount) {
 565       zFrac2 = 0u;
 566       __shortShift64Left(zFrac0, zFrac1, shiftCount, zFrac0, zFrac1);
 567    } else {
 568       __shift64ExtraRightJamming(
 569          zFrac0, zFrac1, 0u, -shiftCount, zFrac0, zFrac1, zFrac2);
 570    }
 571    zExp -= shiftCount;
 572    return __roundAndPackFloat64(zSign, zExp, zFrac0, zFrac1, zFrac2);
 573 }
 574
 575 /* Takes two double-precision floating-point values `a' and `b', one of which
 576  * is a NaN, and returns the appropriate NaN result.
 577  */
 578 uint64_t
 579 __propagateFloat64NaN(uint64_t __a, uint64_t __b)
 580 {
 581    bool aIsNaN = __is_nan(__a);
 582    bool bIsNaN = __is_nan(__b);
 583    uvec2 a = unpackUint2x32(__a);
 584    uvec2 b = unpackUint2x32(__b);
 585    a.y |= 0x00080000u;
 586    b.y |= 0x00080000u;
 587
 588    return packUint2x32(mix(b, mix(a, b, bvec2(bIsNaN, bIsNaN)), bvec2(aIsNaN, aIsNaN)));
 589 }
 590
 591 /* Returns the result of adding the double-precision floating-point values
 592  * `a' and `b'.  The operation is performed according to the IEEE Standard for
 593  * Floating-Point Arithmetic.
 594  */
 595 uint64_t
 596 __fadd64(uint64_t a, uint64_t b)
 597 {
 598    uint aSign = __extractFloat64Sign(a);
 599    uint bSign = __extractFloat64Sign(b);
 600    uint aFracLo = __extractFloat64FracLo(a);
 601    uint aFracHi = __extractFloat64FracHi(a);
 602    uint bFracLo = __extractFloat64FracLo(b);
 603    uint bFracHi = __extractFloat64FracHi(b);
 604    int aExp = __extractFloat64Exp(a);
 605    int bExp = __extractFloat64Exp(b);
 606    uint zFrac0 = 0u;
 607    uint zFrac1 = 0u;
 608    int expDiff = aExp - bExp;
 609    if (aSign == bSign) {
 610       uint zFrac2 = 0u;
 611       int zExp;
 612       bool orig_exp_diff_is_zero = (expDiff == 0);
 613
 614       if (orig_exp_diff_is_zero) {
 615          if (aExp == 0x7FF) {
 616             bool propagate = (aFracHi | aFracLo | bFracHi | bFracLo) != 0u;
 617             return mix(a, __propagateFloat64NaN(a, b), propagate);
 618          }
 619          __add64(aFracHi, aFracLo, bFracHi, bFracLo, zFrac0, zFrac1);
 620          if (aExp == 0)
 621             return __packFloat64(aSign, 0, zFrac0, zFrac1);
 622          zFrac2 = 0u;
 623          zFrac0 |= 0x00200000u;
 624          zExp = aExp;
 625          __shift64ExtraRightJamming(
 626             zFrac0, zFrac1, zFrac2, 1, zFrac0, zFrac1, zFrac2);
 627       } else if (0 < expDiff) {
 628          if (aExp == 0x7FF) {
 629             bool propagate = (aFracHi | aFracLo) != 0u;
 630             return mix(a, __propagateFloat64NaN(a, b), propagate);
 631          }
 632
 633          expDiff = mix(expDiff, expDiff - 1, bExp == 0);
 634          bFracHi = mix(bFracHi | 0x00100000u, bFracHi, bExp == 0);
 635          __shift64ExtraRightJamming(
 636             bFracHi, bFracLo, 0u, expDiff, bFracHi, bFracLo, zFrac2);
 637          zExp = aExp;
 638       } else if (expDiff < 0) {
 639          if (bExp == 0x7FF) {
 640             bool propagate = (bFracHi | bFracLo) != 0u;
 641             return mix(__packFloat64(aSign, 0x7ff, 0u, 0u), __propagateFloat64NaN(a, b), propagate);
 642          }
 643          expDiff = mix(expDiff, expDiff + 1, aExp == 0);
 644          aFracHi = mix(aFracHi | 0x00100000u, aFracHi, aExp == 0);
 645          __shift64ExtraRightJamming(
 646             aFracHi, aFracLo, 0u, - expDiff, aFracHi, aFracLo, zFrac2);
 647          zExp = bExp;
 648       }
 649       if (!orig_exp_diff_is_zero) {
 650          aFracHi |= 0x00100000u;
 651          __add64(aFracHi, aFracLo, bFracHi, bFracLo, zFrac0, zFrac1);
 652          --zExp;
 653          if (!(zFrac0 < 0x00200000u)) {
 654             __shift64ExtraRightJamming(zFrac0, zFrac1, zFrac2, 1, zFrac0, zFrac1, zFrac2);
 655             ++zExp;
 656          }
 657       }
 658       return __roundAndPackFloat64(aSign, zExp, zFrac0, zFrac1, zFrac2);
 659
 660    } else {
 661       int zExp;
 662
 663       __shortShift64Left(aFracHi, aFracLo, 10, aFracHi, aFracLo);
 664       __shortShift64Left(bFracHi, bFracLo, 10, bFracHi, bFracLo);
 665       if (0 < expDiff) {
 666          if (aExp == 0x7FF) {
 667             bool propagate = (aFracHi | aFracLo) != 0u;
 668             return mix(a, __propagateFloat64NaN(a, b), propagate);
 669          }
 670          expDiff = mix(expDiff, expDiff - 1, bExp == 0);
 671          bFracHi = mix(bFracHi | 0x40000000u, bFracHi, bExp == 0);
 672          __shift64RightJamming(bFracHi, bFracLo, expDiff, bFracHi, bFracLo);
 673          aFracHi |= 0x40000000u;
 674          __sub64(aFracHi, aFracLo, bFracHi, bFracLo, zFrac0, zFrac1);
 675          zExp = aExp;
 676          --zExp;
 677          return __normalizeRoundAndPackFloat64(aSign, zExp - 10, zFrac0, zFrac1);
 678       }
 679       if (expDiff < 0) {
 680          if (bExp == 0x7FF) {
 681             bool propagate = (bFracHi | bFracLo) != 0u;
 682             return mix(__packFloat64(aSign ^ 0x80000000u, 0x7ff, 0u, 0u), __propagateFloat64NaN(a, b), propagate);
 683          }
 684          expDiff = mix(expDiff, expDiff + 1, aExp == 0);
 685          aFracHi = mix(aFracHi | 0x40000000u, aFracHi, aExp == 0);
 686          __shift64RightJamming(aFracHi, aFracLo, - expDiff, aFracHi, aFracLo);
 687          bFracHi |= 0x40000000u;
 688          __sub64(bFracHi, bFracLo, aFracHi, aFracLo, zFrac0, zFrac1);
 689          zExp = bExp;
 690          aSign ^= 0x80000000u;
 691          --zExp;
 692          return __normalizeRoundAndPackFloat64(aSign, zExp - 10, zFrac0, zFrac1);
 693       }
 694       if (aExp == 0x7FF) {
 695           bool propagate = (aFracHi | aFracLo | bFracHi | bFracLo) != 0u;
 696          return mix(0xFFFFFFFFFFFFFFFFUL, __propagateFloat64NaN(a, b), propagate);
 697       }
 698       bExp = mix(bExp, 1, aExp == 0);
 699       aExp = mix(aExp, 1, aExp == 0);
 700       bool zexp_normal = false;
 701       bool blta = true;
 702       if (bFracHi < aFracHi) {
 703          __sub64(aFracHi, aFracLo, bFracHi, bFracLo, zFrac0, zFrac1);
 704          zexp_normal = true;
 705       }
 706       else if (aFracHi < bFracHi) {
 707          __sub64(bFracHi, bFracLo, aFracHi, aFracLo, zFrac0, zFrac1);
 708          blta = false;
 709          zexp_normal = true;
 710       }
 711       else if (bFracLo < aFracLo) {
 712          __sub64(aFracHi, aFracLo, bFracHi, bFracLo, zFrac0, zFrac1);
 713          zexp_normal = true;
 714       }
 715       else if (aFracLo < bFracLo) {
 716          __sub64(bFracHi, bFracLo, aFracHi, aFracLo, zFrac0, zFrac1);
 717           blta = false;
 718           zexp_normal = true;
 719       }
 720       zExp = mix(bExp, aExp, blta);
 721       aSign = mix(aSign ^ 0x80000000u, aSign, blta);
 722       uint64_t retval_0 = __packFloat64(uint(FLOAT_ROUNDING_MODE == FLOAT_ROUND_DOWN) << 31, 0, 0u, 0u);
 723       uint64_t retval_1 = __normalizeRoundAndPackFloat64(aSign, zExp - 11, zFrac0, zFrac1);
 724       return mix(retval_0, retval_1, zexp_normal);
 725    }
 726 }
 727
 728 /* Multiplies the 64-bit value formed by concatenating `a0' and `a1' to the
 729  * 64-bit value formed by concatenating `b0' and `b1' to obtain a 128-bit
 730  * product.  The product is broken into four 32-bit pieces which are stored at
 731  * the locations pointed to by `z0Ptr', `z1Ptr', `z2Ptr', and `z3Ptr'.
 732  */
 733 void
 734 __mul64To128(uint a0, uint a1, uint b0, uint b1,
 735              out uint z0Ptr,
 736              out uint z1Ptr,
 737              out uint z2Ptr,
 738              out uint z3Ptr)
 739 {
 740    uint z0 = 0u;
 741    uint z1 = 0u;
 742    uint z2 = 0u;
 743    uint z3 = 0u;
 744    uint more1 = 0u;
 745    uint more2 = 0u;
 746
 747    umulExtended(a1, b1, z2, z3);
 748    umulExtended(a1, b0, z1, more2);
 749    __add64(z1, more2, 0u, z2, z1, z2);
 750    umulExtended(a0, b0, z0, more1);
 751    __add64(z0, more1, 0u, z1, z0, z1);
 752    umulExtended(a0, b1, more1, more2);
 753    __add64(more1, more2, 0u, z2, more1, z2);
 754    __add64(z0, z1, 0u, more1, z0, z1);
 755    z3Ptr = z3;
 756    z2Ptr = z2;
 757    z1Ptr = z1;
 758    z0Ptr = z0;
 759 }
 760
 761 /* Normalizes the subnormal double-precision floating-point value represented
 762  * by the denormalized significand formed by the concatenation of `aFrac0' and
 763  * `aFrac1'.  The normalized exponent is stored at the location pointed to by
 764  * `zExpPtr'.  The most significant 21 bits of the normalized significand are
 765  * stored at the location pointed to by `zFrac0Ptr', and the least significant
 766  * 32 bits of the normalized significand are stored at the location pointed to
 767  * by `zFrac1Ptr'.
 768  */
 769 void
 770 __normalizeFloat64Subnormal(uint aFrac0, uint aFrac1,
 771                             out int zExpPtr,
 772                             out uint zFrac0Ptr,
 773                             out uint zFrac1Ptr)
 774 {
 775    int shiftCount;
 776    uint temp_zfrac0, temp_zfrac1;
 777    shiftCount = __countLeadingZeros32(mix(aFrac0, aFrac1, aFrac0 == 0u)) - 11;
 778    zExpPtr = mix(1 - shiftCount, -shiftCount - 31, aFrac0 == 0u);
 779
 780    temp_zfrac0 = mix(aFrac1<<shiftCount, aFrac1>>(-shiftCount), shiftCount < 0);
 781    temp_zfrac1 = mix(0u, aFrac1<<(shiftCount & 31), shiftCount < 0);
 782
 783    __shortShift64Left(aFrac0, aFrac1, shiftCount, zFrac0Ptr, zFrac1Ptr);
 784
 785    zFrac0Ptr = mix(zFrac0Ptr, temp_zfrac0, aFrac0 == 0);
 786    zFrac1Ptr = mix(zFrac1Ptr, temp_zfrac1, aFrac0 == 0);
 787 }
 788
 789 /* Returns the result of multiplying the double-precision floating-point values
 790  * `a' and `b'.  The operation is performed according to the IEEE Standard for
 791  * Floating-Point Arithmetic.
 792  */
 793 uint64_t
 794 __fmul64(uint64_t a, uint64_t b)
 795 {
 796    uint zFrac0 = 0u;
 797    uint zFrac1 = 0u;
 798    uint zFrac2 = 0u;
 799    uint zFrac3 = 0u;
 800    int zExp;
 801
 802    uint aFracLo = __extractFloat64FracLo(a);
 803    uint aFracHi = __extractFloat64FracHi(a);
 804    uint bFracLo = __extractFloat64FracLo(b);
 805    uint bFracHi = __extractFloat64FracHi(b);
 806    int aExp = __extractFloat64Exp(a);
 807    uint aSign = __extractFloat64Sign(a);
 808    int bExp = __extractFloat64Exp(b);
 809    uint bSign = __extractFloat64Sign(b);
 810    uint zSign = aSign ^ bSign;
 811    if (aExp == 0x7FF) {
 812       if (((aFracHi | aFracLo) != 0u) ||
 813          ((bExp == 0x7FF) && ((bFracHi | bFracLo) != 0u))) {
 814          return __propagateFloat64NaN(a, b);
 815       }
 816       if ((uint(bExp) | bFracHi | bFracLo) == 0u)
 817             return 0xFFFFFFFFFFFFFFFFUL;
 818       return __packFloat64(zSign, 0x7FF, 0u, 0u);
 819    }
 820    if (bExp == 0x7FF) {
 821       if ((bFracHi | bFracLo) != 0u)
 822          return __propagateFloat64NaN(a, b);
 823       if ((uint(aExp) | aFracHi | aFracLo) == 0u)
 824          return 0xFFFFFFFFFFFFFFFFUL;
 825       return __packFloat64(zSign, 0x7FF, 0u, 0u);
 826    }
 827    if (aExp == 0) {
 828       if ((aFracHi | aFracLo) == 0u)
 829          return __packFloat64(zSign, 0, 0u, 0u);
 830       __normalizeFloat64Subnormal(aFracHi, aFracLo, aExp, aFracHi, aFracLo);
 831    }
 832    if (bExp == 0) {
 833       if ((bFracHi | bFracLo) == 0u)
 834          return __packFloat64(zSign, 0, 0u, 0u);
 835       __normalizeFloat64Subnormal(bFracHi, bFracLo, bExp, bFracHi, bFracLo);
 836    }
 837    zExp = aExp + bExp - 0x400;
 838    aFracHi |= 0x00100000u;
 839    __shortShift64Left(bFracHi, bFracLo, 12, bFracHi, bFracLo);
 840    __mul64To128(
 841       aFracHi, aFracLo, bFracHi, bFracLo, zFrac0, zFrac1, zFrac2, zFrac3);
 842    __add64(zFrac0, zFrac1, aFracHi, aFracLo, zFrac0, zFrac1);
 843    zFrac2 |= uint(zFrac3 != 0u);
 844    if (0x00200000u <= zFrac0) {
 845       __shift64ExtraRightJamming(
 846          zFrac0, zFrac1, zFrac2, 1, zFrac0, zFrac1, zFrac2);
 847       ++zExp;
 848    }
 849    return __roundAndPackFloat64(zSign, zExp, zFrac0, zFrac1, zFrac2);
 850 }
 851
 852 uint64_t
 853 __ffma64(uint64_t a, uint64_t b, uint64_t c)
 854 {
 855    return __fadd64(__fmul64(a, b), c);
 856 }
 857
 858 /* Shifts the 64-bit value formed by concatenating `a0' and `a1' right by the
 859  * number of bits given in `count'.  Any bits shifted off are lost.  The value
 860  * of `count' can be arbitrarily large; in particular, if `count' is greater
 861  * than 64, the result will be 0.  The result is broken into two 32-bit pieces
 862  * which are stored at the locations pointed to by `z0Ptr' and `z1Ptr'.
 863  */
 864 void
 865 __shift64Right(uint a0, uint a1,
 866                int count,
 867                out uint z0Ptr,
 868                out uint z1Ptr)
 869 {
 870    uint z0;
 871    uint z1;
 872    int negCount = (-count) & 31;
 873
 874    z0 = 0u;
 875    z0 = mix(z0, (a0 >> count), count < 32);
 876    z0 = mix(z0, a0, count == 0);
 877
 878    z1 = mix(0u, (a0 >> (count & 31)), count < 64);
 879    z1 = mix(z1, (a0<<negCount) | (a1>>count), count < 32);
 880    z1 = mix(z1, a0, count == 0);
 881
 882    z1Ptr = z1;
 883    z0Ptr = z0;
 884 }
 885
 886 /* Returns the result of converting the double-precision floating-point value
 887  * `a' to the unsigned integer format.  The conversion is performed according
 888  * to the IEEE Standard for Floating-Point Arithmetic.
 889  */
 890 uint
 891 __fp64_to_uint(uint64_t a)
 892 {
 893    uint aFracLo = __extractFloat64FracLo(a);
 894    uint aFracHi = __extractFloat64FracHi(a);
 895    int aExp = __extractFloat64Exp(a);
 896    uint aSign = __extractFloat64Sign(a);
 897
 898    if ((aExp == 0x7FF) && ((aFracHi | aFracLo) != 0u))
 899       return 0xFFFFFFFFu;
 900
 901    aFracHi |= mix(0u, 0x00100000u, aExp != 0);
 902
 903    int shiftDist = 0x427 - aExp;
 904    if (0 < shiftDist)
 905       __shift64RightJamming(aFracHi, aFracLo, shiftDist, aFracHi, aFracLo);
 906
 907    if ((aFracHi & 0xFFFFF000u) != 0u)
 908       return mix(~0u, 0u, aSign != 0u);
 909
 910    uint z = 0u;
 911    uint zero = 0u;
 912    __shift64Right(aFracHi, aFracLo, 12, zero, z);
 913
 914    uint expt = mix(~0u, 0u, aSign != 0u);
 915
 916    return mix(z, expt, (aSign != 0u) && (z != 0u));
 917 }
 918
 919 uint64_t
 920 __uint_to_fp64(uint a)
 921 {
 922    if (a == 0u)
 923       return 0ul;
 924
 925    int shiftDist = __countLeadingZeros32(a) + 21;
 926
 927    uint aHigh = 0u;
 928    uint aLow = 0u;
 929    int negCount = (- shiftDist) & 31;
 930
 931    aHigh = mix(0u, a<< shiftDist - 32, shiftDist < 64);
 932    aLow = 0u;
 933    aHigh = mix(aHigh, 0u, shiftDist == 0);
 934    aLow = mix(aLow, a, shiftDist ==0);
 935    aHigh = mix(aHigh, a >> negCount, shiftDist < 32);
 936    aLow = mix(aLow, a << shiftDist, shiftDist < 32);
 937
 938    return __packFloat64(0u, 0x432 - shiftDist, aHigh, aLow);
 939 }
 940
 941 uint64_t
 942 __uint64_to_fp64(uint64_t a)
 943 {
 944    if (a == 0u)
 945       return 0ul;
 946
 947    uvec2 aFrac = unpackUint2x32(a);
 948    uint aFracLo = __extractFloat64FracLo(a);
 949    uint aFracHi = __extractFloat64FracHi(a);
 950
 951    if ((aFracHi & 0x80000000u) != 0u) {
 952       __shift64RightJamming(aFracHi, aFracLo, 1, aFracHi, aFracLo);
 953       return __roundAndPackFloat64(0, 0x433, aFracHi, aFracLo, 0u);
 954    } else {
 955       return __normalizeRoundAndPackFloat64(0, 0x432, aFrac.y, aFrac.x);
 956    }
 957 }
 958
 959 uint64_t
 960 __fp64_to_uint64(uint64_t a)
 961 {
 962    uint aFracLo = __extractFloat64FracLo(a);
 963    uint aFracHi = __extractFloat64FracHi(a);
 964    int aExp = __extractFloat64Exp(a);
 965    uint aSign = __extractFloat64Sign(a);
 966    uint zFrac2 = 0u;
 967    uint64_t default_nan = 0xFFFFFFFFFFFFFFFFUL;
 968
 969    aFracHi = mix(aFracHi, aFracHi | 0x00100000u, aExp != 0);
 970    int shiftCount = 0x433 - aExp;
 971
 972    if ( shiftCount <= 0 ) {
 973       if (shiftCount < -11 && aExp == 0x7FF) {
 974          if ((aFracHi | aFracLo) != 0u)
 975             return __propagateFloat64NaN(a, a);
 976          return mix(default_nan, a, aSign == 0u);
 977       }
 978       __shortShift64Left(aFracHi, aFracLo, -shiftCount, aFracHi, aFracLo);
 979    } else {
 980       __shift64ExtraRightJamming(aFracHi, aFracLo, zFrac2, shiftCount,
 981                                  aFracHi, aFracLo, zFrac2);
 982    }
 983    return __roundAndPackUInt64(aSign, aFracHi, aFracLo, zFrac2);
 984 }
 985
 986 int64_t
 987 __fp64_to_int64(uint64_t a)
 988 {
 989    uint zFrac2 = 0u;
 990    uint aFracLo = __extractFloat64FracLo(a);
 991    uint aFracHi = __extractFloat64FracHi(a);
 992    int aExp = __extractFloat64Exp(a);
 993    uint aSign = __extractFloat64Sign(a);
 994    int64_t default_NegNaN = -0x7FFFFFFFFFFFFFFEL;
 995    int64_t default_PosNaN = 0xFFFFFFFFFFFFFFFFL;
 996
 997    aFracHi = mix(aFracHi, aFracHi | 0x00100000u, aExp != 0);
 998    int shiftCount = 0x433 - aExp;
 999
1000    if (shiftCount <= 0) {
1001       if (shiftCount < -11 && aExp == 0x7FF) {
1002          if ((aFracHi | aFracLo) != 0u)
1003             return default_NegNaN;
1004          return mix(default_NegNaN, default_PosNaN, aSign == 0u);
1005       }
1006       __shortShift64Left(aFracHi, aFracLo, -shiftCount, aFracHi, aFracLo);
1007    } else {
1008       __shift64ExtraRightJamming(aFracHi, aFracLo, zFrac2, shiftCount,
1009                                  aFracHi, aFracLo, zFrac2);
1010    }
1011
1012    return __roundAndPackInt64(aSign, aFracHi, aFracLo, zFrac2);
1013 }
1014
1015 uint64_t
1016 __fp32_to_uint64(float f)
1017 {
1018    uint a = floatBitsToUint(f);
1019    uint aFrac = a & 0x007FFFFFu;
1020    int aExp = int((a>>23) & 0xFFu);
1021    uint aSign = a & 0x80000000u;
1022    uint zFrac0 = 0u;
1023    uint zFrac1 = 0u;
1024    uint zFrac2 = 0u;
1025    uint64_t default_nan = 0xFFFFFFFFFFFFFFFFUL;
1026    int shiftCount = 0xBE - aExp;
1027
1028    if (shiftCount <0) {
1029       if (aExp == 0xFF)
1030          return default_nan;
1031    }
1032
1033    aFrac = mix(aFrac, aFrac | 0x00800000u, aExp != 0);
1034    __shortShift64Left(aFrac, 0, 40, zFrac0, zFrac1);
1035
1036    if (shiftCount != 0) {
1037       __shift64ExtraRightJamming(zFrac0, zFrac1, zFrac2, shiftCount,
1038                                  zFrac0, zFrac1, zFrac2);
1039    }
1040
1041    return __roundAndPackUInt64(aSign, zFrac0, zFrac1, zFrac2);
1042 }
1043
1044 int64_t
1045 __fp32_to_int64(float f)
1046 {
1047    uint a = floatBitsToUint(f);
1048    uint aFrac = a & 0x007FFFFFu;
1049    int aExp = int((a>>23) & 0xFFu);
1050    uint aSign = a & 0x80000000u;
1051    uint zFrac0 = 0u;
1052    uint zFrac1 = 0u;
1053    uint zFrac2 = 0u;
1054    int64_t default_NegNaN = -0x7FFFFFFFFFFFFFFEL;
1055    int64_t default_PosNaN = 0xFFFFFFFFFFFFFFFFL;
1056    int shiftCount = 0xBE - aExp;
1057
1058    if (shiftCount <0) {
1059       if (aExp == 0xFF && aFrac != 0u)
1060          return default_NegNaN;
1061       return mix(default_NegNaN, default_PosNaN, aSign == 0u);
1062    }
1063
1064    aFrac = mix(aFrac, aFrac | 0x00800000u, aExp != 0);
1065    __shortShift64Left(aFrac, 0, 40, zFrac0, zFrac1);
1066
1067    if (shiftCount != 0) {
1068       __shift64ExtraRightJamming(zFrac0, zFrac1, zFrac2, shiftCount,
1069                                  zFrac0, zFrac1, zFrac2);
1070    }
1071
1072    return __roundAndPackInt64(aSign, zFrac0, zFrac1, zFrac2);
1073 }
1074
1075 uint64_t
1076 __int64_to_fp64(int64_t a)
1077 {
1078    if (a==0)
1079       return 0ul;
1080
1081    uint64_t absA = mix(uint64_t(a), uint64_t(-a), a < 0);
1082    uint aFracHi = __extractFloat64FracHi(absA);
1083    uvec2 aFrac = unpackUint2x32(absA);
1084    uint zSign = uint(unpackInt2x32(a).y) & 0x80000000u;
1085
1086    if ((aFracHi & 0x80000000u) != 0u) {
1087       return mix(0ul, __packFloat64(0x80000000u, 0x434, 0u, 0u), a < 0);
1088    }
1089
1090    return __normalizeRoundAndPackFloat64(zSign, 0x432, aFrac.y, aFrac.x);
1091 }
1092
1093 /* Returns the result of converting the double-precision floating-point value
1094  * `a' to the 32-bit two's complement integer format.  The conversion is
1095  * performed according to the IEEE Standard for Floating-Point Arithmetic---
1096  * which means in particular that the conversion is rounded according to the
1097  * current rounding mode.  If `a' is a NaN, the largest positive integer is
1098  * returned.  Otherwise, if the conversion overflows, the largest integer with
1099  * the same sign as `a' is returned.
1100  */
1101 int
1102 __fp64_to_int(uint64_t a)
1103 {
1104    uint aFracLo = __extractFloat64FracLo(a);
1105    uint aFracHi = __extractFloat64FracHi(a);
1106    int aExp = __extractFloat64Exp(a);
1107    uint aSign = __extractFloat64Sign(a);
1108
1109    uint absZ = 0u;
1110    uint aFracExtra = 0u;
1111    int shiftCount = aExp - 0x413;
1112
1113    if (0 <= shiftCount) {
1114       if (0x41E < aExp) {
1115          if ((aExp == 0x7FF) && bool(aFracHi | aFracLo))
1116             aSign = 0u;
1117          return mix(0x7FFFFFFF, 0x80000000, aSign != 0u);
1118       }
1119       __shortShift64Left(aFracHi | 0x00100000u, aFracLo, shiftCount, absZ, aFracExtra);
1120    } else {
1121       if (aExp < 0x3FF)
1122          return 0;
1123
1124       aFracHi |= 0x00100000u;
1125       aFracExtra = ( aFracHi << (shiftCount & 31)) | aFracLo;
1126       absZ = aFracHi >> (- shiftCount);
1127    }
1128
1129    int z = mix(int(absZ), -int(absZ), aSign != 0u);
1130    int nan = mix(0x7FFFFFFF, 0x80000000, aSign != 0u);
1131    return mix(z, nan, ((aSign != 0u) != (z < 0)) && bool(z));
1132 }
1133
1134 /* Returns the result of converting the 32-bit two's complement integer `a'
1135  * to the double-precision floating-point format.  The conversion is performed
1136  * according to the IEEE Standard for Floating-Point Arithmetic.
1137  */
1138 uint64_t
1139 __int_to_fp64(int a)
1140 {
1141    uint zFrac0 = 0u;
1142    uint zFrac1 = 0u;
1143    if (a==0)
1144       return __packFloat64(0u, 0, 0u, 0u);
1145    uint zSign = uint(a) & 0x80000000u;
1146    uint absA = mix(uint(a), uint(-a), a < 0);
1147    int shiftCount = __countLeadingZeros32(absA) - 11;
1148    if (0 <= shiftCount) {
1149       zFrac0 = absA << shiftCount;
1150       zFrac1 = 0u;
1151    } else {
1152       __shift64Right(absA, 0u, -shiftCount, zFrac0, zFrac1);
1153    }
1154    return __packFloat64(zSign, 0x412 - shiftCount, zFrac0, zFrac1);
1155 }
1156
1157 bool
1158 __fp64_to_bool(uint64_t a)
1159 {
1160    return !__feq64_nonnan(__fabs64(a), 0ul);
1161 }
1162
1163 uint64_t
1164 __bool_to_fp64(bool a)
1165 {
1166    return packUint2x32(uvec2(0x00000000u, uint(-int(a) & 0x3ff00000)));
1167 }
1168
1169 /* Packs the sign `zSign', exponent `zExp', and significand `zFrac' into a
1170  * single-precision floating-point value, returning the result.  After being
1171  * shifted into the proper positions, the three fields are simply added
1172  * together to form the result.  This means that any integer portion of `zSig'
1173  * will be added into the exponent.  Since a properly normalized significand
1174  * will have an integer portion equal to 1, the `zExp' input should be 1 less
1175  * than the desired result exponent whenever `zFrac' is a complete, normalized
1176  * significand.
1177  */
1178 float
1179 __packFloat32(uint zSign, int zExp, uint zFrac)
1180 {
1181    return uintBitsToFloat(zSign + (uint(zExp)<<23) + zFrac);
1182 }
1183
1184 /* Takes an abstract floating-point value having sign `zSign', exponent `zExp',
1185  * and significand `zFrac', and returns the proper single-precision floating-
1186  * point value corresponding to the abstract input.  Ordinarily, the abstract
1187  * value is simply rounded and packed into the single-precision format, with
1188  * the inexact exception raised if the abstract input cannot be represented
1189  * exactly.  However, if the abstract value is too large, the overflow and
1190  * inexact exceptions are raised and an infinity or maximal finite value is
1191  * returned.  If the abstract value is too small, the input value is rounded to
1192  * a subnormal number, and the underflow and inexact exceptions are raised if
1193  * the abstract input cannot be represented exactly as a subnormal single-
1194  * precision floating-point number.
1195  *     The input significand `zFrac' has its binary point between bits 30
1196  * and 29, which is 7 bits to the left of the usual location.  This shifted
1197  * significand must be normalized or smaller.  If `zFrac' is not normalized,
1198  * `zExp' must be 0; in that case, the result returned is a subnormal number,
1199  * and it must not require rounding.  In the usual case that `zFrac' is
1200  * normalized, `zExp' must be 1 less than the "true" floating-point exponent.
1201  * The handling of underflow and overflow follows the IEEE Standard for
1202  * Floating-Point Arithmetic.
1203  */
1204 float
1205 __roundAndPackFloat32(uint zSign, int zExp, uint zFrac)
1206 {
1207    bool roundNearestEven;
1208    int roundIncrement;
1209    int roundBits;
1210
1211    roundNearestEven = FLOAT_ROUNDING_MODE == FLOAT_ROUND_NEAREST_EVEN;
1212    roundIncrement = 0x40;
1213    if (!roundNearestEven) {
1214       if (FLOAT_ROUNDING_MODE == FLOAT_ROUND_TO_ZERO) {
1215          roundIncrement = 0;
1216       } else {
1217          roundIncrement = 0x7F;
1218          if (zSign != 0u) {
1219             if (FLOAT_ROUNDING_MODE == FLOAT_ROUND_UP)
1220                roundIncrement = 0;
1221          } else {
1222             if (FLOAT_ROUNDING_MODE == FLOAT_ROUND_DOWN)
1223                roundIncrement = 0;
1224          }
1225       }
1226    }
1227    roundBits = int(zFrac & 0x7Fu);
1228    if (0xFDu <= uint(zExp)) {
1229       if ((0xFD < zExp) || ((zExp == 0xFD) && (int(zFrac) + roundIncrement) < 0))
1230          return __packFloat32(zSign, 0xFF, 0u) - float(roundIncrement == 0);
1231       int count = -zExp;
1232       bool zexp_lt0 = zExp < 0;
1233       uint zFrac_lt0 = mix(uint(zFrac != 0u), (zFrac>>count) | uint((zFrac<<((-count) & 31)) != 0u), (-zExp) < 32);
1234       zFrac = mix(zFrac, zFrac_lt0, zexp_lt0);
1235       roundBits = mix(roundBits, int(zFrac) & 0x7f, zexp_lt0);
1236       zExp = mix(zExp, 0, zexp_lt0);
1237    }
1238    zFrac = (zFrac + uint(roundIncrement))>>7;
1239    zFrac &= ~uint(((roundBits ^ 0x40) == 0) && roundNearestEven);
1240
1241    return __packFloat32(zSign, mix(zExp, 0, zFrac == 0u), zFrac);
1242 }
1243
1244 /* Returns the result of converting the double-precision floating-point value
1245  * `a' to the single-precision floating-point format.  The conversion is
1246  * performed according to the IEEE Standard for Floating-Point Arithmetic.
1247  */
1248 float
1249 __fp64_to_fp32(uint64_t __a)
1250 {
1251    uvec2 a = unpackUint2x32(__a);
1252    uint zFrac = 0u;
1253    uint allZero = 0u;
1254
1255    uint aFracLo = __extractFloat64FracLo(__a);
1256    uint aFracHi = __extractFloat64FracHi(__a);
1257    int aExp = __extractFloat64Exp(__a);
1258    uint aSign = __extractFloat64Sign(__a);
1259    if (aExp == 0x7FF) {
1260       __shortShift64Left(a.y, a.x, 12, a.y, a.x);
1261       float rval = uintBitsToFloat(aSign | 0x7FC00000u | (a.y>>9));
1262       rval = mix(__packFloat32(aSign, 0xFF, 0u), rval, (aFracHi | aFracLo) != 0u);
1263       return rval;
1264    }
1265    __shift64RightJamming(aFracHi, aFracLo, 22, allZero, zFrac);
1266    zFrac = mix(zFrac, zFrac | 0x40000000u, aExp != 0);
1267    return __roundAndPackFloat32(aSign, aExp - 0x381, zFrac);
1268 }
1269
1270 float
1271 __uint64_to_fp32(uint64_t __a)
1272 {
1273    uvec2 aFrac = unpackUint2x32(__a);
1274    int shiftCount = mix(__countLeadingZeros32(aFrac.y) - 33,
1275                         __countLeadingZeros32(aFrac.x) - 1,
1276                         aFrac.y == 0u);
1277
1278    if (0 <= shiftCount)
1279       __shortShift64Left(aFrac.y, aFrac.x, shiftCount, aFrac.y, aFrac.x);
1280    else
1281       __shift64RightJamming(aFrac.y, aFrac.x, -shiftCount, aFrac.y, aFrac.x);
1282
1283    return __roundAndPackFloat32(0u, 0x9C - shiftCount, aFrac.x);
1284 }
1285
1286 float
1287 __int64_to_fp32(int64_t __a)
1288 {
1289    uint aSign = uint(unpackInt2x32(__a).y) & 0x80000000u;
1290    uint64_t absA = mix(uint64_t(__a), uint64_t(-__a), __a < 0);
1291    uvec2 aFrac = unpackUint2x32(absA);
1292    int shiftCount = mix(__countLeadingZeros32(aFrac.y) - 33,
1293                         __countLeadingZeros32(aFrac.x) - 1,
1294                         aFrac.y == 0u);
1295
1296    if (0 <= shiftCount)
1297       __shortShift64Left(aFrac.y, aFrac.x, shiftCount, aFrac.y, aFrac.x);
1298    else
1299       __shift64RightJamming(aFrac.y, aFrac.x, -shiftCount, aFrac.y, aFrac.x);
1300
1301    return __roundAndPackFloat32(aSign, 0x9C - shiftCount, aFrac.x);
1302 }
1303
1304 /* Returns the result of converting the single-precision floating-point value
1305  * `a' to the double-precision floating-point format.
1306  */
1307 uint64_t
1308 __fp32_to_fp64(float f)
1309 {
1310    uint a = floatBitsToUint(f);
1311    uint aFrac = a & 0x007FFFFFu;
1312    int aExp = int((a>>23) & 0xFFu);
1313    uint aSign = a & 0x80000000u;
1314    uint zFrac0 = 0u;
1315    uint zFrac1 = 0u;
1316
1317    if (aExp == 0xFF) {
1318       if (aFrac != 0u) {
1319          uint nanLo = 0u;
1320          uint nanHi = a<<9;
1321          __shift64Right(nanHi, nanLo, 12, nanHi, nanLo);
1322          nanHi |= aSign | 0x7FF80000u;
1323          return packUint2x32(uvec2(nanLo, nanHi));
1324       }
1325       return __packFloat64(aSign, 0x7FF, 0u, 0u);
1326     }
1327
1328    if (aExp == 0) {
1329       if (aFrac == 0u)
1330          return __packFloat64(aSign, 0, 0u, 0u);
1331       /* Normalize subnormal */
1332       int shiftCount = __countLeadingZeros32(aFrac) - 8;
1333       aFrac <<= shiftCount;
1334       aExp = 1 - shiftCount;
1335       --aExp;
1336    }
1337
1338    __shift64Right(aFrac, 0u, 3, zFrac0, zFrac1);
1339    return __packFloat64(aSign, aExp + 0x380, zFrac0, zFrac1);
1340 }
1341
1342 /* Adds the 96-bit value formed by concatenating `a0', `a1', and `a2' to the
1343  * 96-bit value formed by concatenating `b0', `b1', and `b2'.  Addition is
1344  * modulo 2^96, so any carry out is lost.  The result is broken into three
1345  * 32-bit pieces which are stored at the locations pointed to by `z0Ptr',
1346  * `z1Ptr', and `z2Ptr'.
1347  */
1348 void
1349 __add96(uint a0, uint a1, uint a2,
1350         uint b0, uint b1, uint b2,
1351         out uint z0Ptr,
1352         out uint z1Ptr,
1353         out uint z2Ptr)
1354 {
1355    uint z2 = a2 + b2;
1356    uint carry1 = uint(z2 < a2);
1357    uint z1 = a1 + b1;
1358    uint carry0 = uint(z1 < a1);
1359    uint z0 = a0 + b0;
1360    z1 += carry1;
1361    z0 += uint(z1 < carry1);
1362    z0 += carry0;
1363    z2Ptr = z2;
1364    z1Ptr = z1;
1365    z0Ptr = z0;
1366 }
1367
1368 /* Subtracts the 96-bit value formed by concatenating `b0', `b1', and `b2' from
1369  * the 96-bit value formed by concatenating `a0', `a1', and `a2'.  Subtraction
1370  * is modulo 2^96, so any borrow out (carry out) is lost.  The result is broken
1371  * into three 32-bit pieces which are stored at the locations pointed to by
1372  * `z0Ptr', `z1Ptr', and `z2Ptr'.
1373  */
1374 void
1375 __sub96(uint a0, uint a1, uint a2,
1376         uint b0, uint b1, uint b2,
1377         out uint z0Ptr,
1378         out uint z1Ptr,
1379         out uint z2Ptr)
1380 {
1381    uint z2 = a2 - b2;
1382    uint borrow1 = uint(a2 < b2);
1383    uint z1 = a1 - b1;
1384    uint borrow0 = uint(a1 < b1);
1385    uint z0 = a0 - b0;
1386    z0 -= uint(z1 < borrow1);
1387    z1 -= borrow1;
1388    z0 -= borrow0;
1389    z2Ptr = z2;
1390    z1Ptr = z1;
1391    z0Ptr = z0;
1392 }
1393
1394 /* Returns an approximation to the 32-bit integer quotient obtained by dividing
1395  * `b' into the 64-bit value formed by concatenating `a0' and `a1'.  The
1396  * divisor `b' must be at least 2^31.  If q is the exact quotient truncated
1397  * toward zero, the approximation returned lies between q and q + 2 inclusive.
1398  * If the exact quotient q is larger than 32 bits, the maximum positive 32-bit
1399  * unsigned integer is returned.
1400  */
1401 uint
1402 __estimateDiv64To32(uint a0, uint a1, uint b)
1403 {
1404    uint b0;
1405    uint b1;
1406    uint rem0 = 0u;
1407    uint rem1 = 0u;
1408    uint term0 = 0u;
1409    uint term1 = 0u;
1410    uint z;
1411
1412    if (b <= a0)
1413       return 0xFFFFFFFFu;
1414    b0 = b>>16;
1415    z = (b0<<16 <= a0) ? 0xFFFF0000u : (a0 / b0)<<16;
1416    umulExtended(b, z, term0, term1);
1417    __sub64(a0, a1, term0, term1, rem0, rem1);
1418    while (int(rem0) < 0) {
1419       z -= 0x10000u;
1420       b1 = b<<16;
1421       __add64(rem0, rem1, b0, b1, rem0, rem1);
1422    }
1423    rem0 = (rem0<<16) | (rem1>>16);
1424    z |= (b0<<16 <= rem0) ? 0xFFFFu : rem0 / b0;
1425    return z;
1426 }
1427
1428 uint
1429 __sqrtOddAdjustments(int index)
1430 {
1431    uint res = 0u;
1432    if (index == 0)
1433       res = 0x0004u;
1434    if (index == 1)
1435       res = 0x0022u;
1436    if (index == 2)
1437       res = 0x005Du;
1438    if (index == 3)
1439       res = 0x00B1u;
1440    if (index == 4)
1441       res = 0x011Du;
1442    if (index == 5)
1443       res = 0x019Fu;
1444    if (index == 6)
1445       res = 0x0236u;
1446    if (index == 7)
1447       res = 0x02E0u;
1448    if (index == 8)
1449       res = 0x039Cu;
1450    if (index == 9)
1451       res = 0x0468u;
1452    if (index == 10)
1453       res = 0x0545u;
1454    if (index == 11)
1455       res = 0x631u;
1456    if (index == 12)
1457       res = 0x072Bu;
1458    if (index == 13)
1459       res = 0x0832u;
1460    if (index == 14)
1461       res = 0x0946u;
1462    if (index == 15)
1463       res = 0x0A67u;
1464
1465    return res;
1466 }
1467
1468 uint
1469 __sqrtEvenAdjustments(int index)
1470 {
1471    uint res = 0u;
1472    if (index == 0)
1473       res = 0x0A2Du;
1474    if (index == 1)
1475       res = 0x08AFu;
1476    if (index == 2)
1477       res = 0x075Au;
1478    if (index == 3)
1479       res = 0x0629u;
1480    if (index == 4)
1481       res = 0x051Au;
1482    if (index == 5)
1483       res = 0x0429u;
1484    if (index == 6)
1485       res = 0x0356u;
1486    if (index == 7)
1487       res = 0x029Eu;
1488    if (index == 8)
1489       res = 0x0200u;
1490    if (index == 9)
1491       res = 0x0179u;
1492    if (index == 10)
1493       res = 0x0109u;
1494    if (index == 11)
1495       res = 0x00AFu;
1496    if (index == 12)
1497       res = 0x0068u;
1498    if (index == 13)
1499       res = 0x0034u;
1500    if (index == 14)
1501       res = 0x0012u;
1502    if (index == 15)
1503       res = 0x0002u;
1504
1505    return res;
1506 }
1507
1508 /* Returns an approximation to the square root of the 32-bit significand given
1509  * by `a'.  Considered as an integer, `a' must be at least 2^31.  If bit 0 of
1510  * `aExp' (the least significant bit) is 1, the integer returned approximates
1511  * 2^31*sqrt(`a'/2^31), where `a' is considered an integer.  If bit 0 of `aExp'
1512  * is 0, the integer returned approximates 2^31*sqrt(`a'/2^30).  In either
1513  * case, the approximation returned lies strictly within +/-2 of the exact
1514  * value.
1515  */
1516 uint
1517 __estimateSqrt32(int aExp, uint a)
1518 {
1519    uint z;
1520
1521    int index = int(a>>27 & 15u);
1522    if ((aExp & 1) != 0) {
1523       z = 0x4000u + (a>>17) - __sqrtOddAdjustments(index);
1524       z = ((a / z)<<14) + (z<<15);
1525       a >>= 1;
1526    } else {
1527       z = 0x8000u + (a>>17) - __sqrtEvenAdjustments(index);
1528       z = a / z + z;
1529       z = (0x20000u <= z) ? 0xFFFF8000u : (z<<15);
1530       if (z <= a)
1531          return uint(int(a)>>1);
1532    }
1533    return ((__estimateDiv64To32(a, 0u, z))>>1) + (z>>1);
1534 }
1535
1536 /* Returns the square root of the double-precision floating-point value `a'.
1537  * The operation is performed according to the IEEE Standard for Floating-Point
1538  * Arithmetic.
1539  */
1540 uint64_t
1541 __fsqrt64(uint64_t a)
1542 {
1543    uint zFrac0 = 0u;
1544    uint zFrac1 = 0u;
1545    uint zFrac2 = 0u;
1546    uint doubleZFrac0 = 0u;
1547    uint rem0 = 0u;
1548    uint rem1 = 0u;
1549    uint rem2 = 0u;
1550    uint rem3 = 0u;
1551    uint term0 = 0u;
1552    uint term1 = 0u;
1553    uint term2 = 0u;
1554    uint term3 = 0u;
1555    uint64_t default_nan = 0xFFFFFFFFFFFFFFFFUL;
1556
1557    uint aFracLo = __extractFloat64FracLo(a);
1558    uint aFracHi = __extractFloat64FracHi(a);
1559    int aExp = __extractFloat64Exp(a);
1560    uint aSign = __extractFloat64Sign(a);
1561    if (aExp == 0x7FF) {
1562       if ((aFracHi | aFracLo) != 0u)
1563          return __propagateFloat64NaN(a, a);
1564       if (aSign == 0u)
1565          return a;
1566       return default_nan;
1567    }
1568    if (aSign != 0u) {
1569       if ((uint(aExp) | aFracHi | aFracLo) == 0u)
1570          return a;
1571       return default_nan;
1572    }
1573    if (aExp == 0) {
1574       if ((aFracHi | aFracLo) == 0u)
1575          return __packFloat64(0u, 0, 0u, 0u);
1576       __normalizeFloat64Subnormal(aFracHi, aFracLo, aExp, aFracHi, aFracLo);
1577    }
1578    int zExp = ((aExp - 0x3FF)>>1) + 0x3FE;
1579    aFracHi |= 0x00100000u;
1580    __shortShift64Left(aFracHi, aFracLo, 11, term0, term1);
1581    zFrac0 = (__estimateSqrt32(aExp, term0)>>1) + 1u;
1582    if (zFrac0 == 0u)
1583       zFrac0 = 0x7FFFFFFFu;
1584    doubleZFrac0 = zFrac0 + zFrac0;
1585    __shortShift64Left(aFracHi, aFracLo, 9 - (aExp & 1), aFracHi, aFracLo);
1586    umulExtended(zFrac0, zFrac0, term0, term1);
1587    __sub64(aFracHi, aFracLo, term0, term1, rem0, rem1);
1588    while (int(rem0) < 0) {
1589       --zFrac0;
1590       doubleZFrac0 -= 2u;
1591       __add64(rem0, rem1, 0u, doubleZFrac0 | 1u, rem0, rem1);
1592    }
1593    zFrac1 = __estimateDiv64To32(rem1, 0u, doubleZFrac0);
1594    if ((zFrac1 & 0x1FFu) <= 5u) {
1595       if (zFrac1 == 0u)
1596          zFrac1 = 1u;
1597       umulExtended(doubleZFrac0, zFrac1, term1, term2);
1598       __sub64(rem1, 0u, term1, term2, rem1, rem2);
1599       umulExtended(zFrac1, zFrac1, term2, term3);
1600       __sub96(rem1, rem2, 0u, 0u, term2, term3, rem1, rem2, rem3);
1601       while (int(rem1) < 0) {
1602          --zFrac1;
1603          __shortShift64Left(0u, zFrac1, 1, term2, term3);
1604          term3 |= 1u;
1605          term2 |= doubleZFrac0;
1606          __add96(rem1, rem2, rem3, 0u, term2, term3, rem1, rem2, rem3);
1607       }
1608       zFrac1 |= uint((rem1 | rem2 | rem3) != 0u);
1609    }
1610    __shift64ExtraRightJamming(zFrac0, zFrac1, 0u, 10, zFrac0, zFrac1, zFrac2);
1611    return __roundAndPackFloat64(0u, zExp, zFrac0, zFrac1, zFrac2);
1612 }
1613
1614 uint64_t
1615 __ftrunc64(uint64_t __a)
1616 {
1617    uvec2 a = unpackUint2x32(__a);
1618    int aExp = __extractFloat64Exp(__a);
1619    uint zLo;
1620    uint zHi;
1621
1622    int unbiasedExp = aExp - 1023;
1623    int fracBits = 52 - unbiasedExp;
1624    uint maskLo = mix(~0u << fracBits, 0u, fracBits >= 32);
1625    uint maskHi = mix(~0u << (fracBits - 32), ~0u, fracBits < 33);
1626    zLo = maskLo & a.x;
1627    zHi = maskHi & a.y;
1628
1629    zLo = mix(zLo, 0u, unbiasedExp < 0);
1630    zHi = mix(zHi, 0u, unbiasedExp < 0);
1631    zLo = mix(zLo, a.x, unbiasedExp > 52);
1632    zHi = mix(zHi, a.y, unbiasedExp > 52);
1633    return packUint2x32(uvec2(zLo, zHi));
1634 }
1635
1636 uint64_t
1637 __ffloor64(uint64_t a)
1638 {
1639    bool is_positive = __fge64(a, 0ul);
1640    uint64_t tr = __ftrunc64(a);
1641
1642    if (is_positive || __feq64(tr, a)) {
1643       return tr;
1644    } else {
1645       return __fadd64(tr, 0xbff0000000000000ul /* -1.0 */);
1646    }
1647 }
1648
1649 uint64_t
1650 __fround64(uint64_t __a)
1651 {
1652    uvec2 a = unpackUint2x32(__a);
1653    int unbiasedExp = __extractFloat64Exp(__a) - 1023;
1654    uint aHi = a.y;
1655    uint aLo = a.x;
1656
1657    if (unbiasedExp < 20) {
1658       if (unbiasedExp < 0) {
1659          if ((aHi & 0x80000000u) != 0u && aLo == 0u) {
1660             return 0;
1661          }
1662          aHi &= 0x80000000u;
1663          if ((a.y & 0x000FFFFFu) == 0u && a.x == 0u) {
1664             aLo = 0u;
1665             return packUint2x32(uvec2(aLo, aHi));
1666          }
1667          aHi = mix(aHi, (aHi | 0x3FF00000u), unbiasedExp == -1);
1668          aLo = 0u;
1669       } else {
1670          uint maskExp = 0x000FFFFFu >> unbiasedExp;
1671          uint lastBit = maskExp + 1;
1672          aHi += 0x00080000u >> unbiasedExp;
1673          if ((aHi & maskExp) == 0u)
1674             aHi &= ~lastBit;
1675          aHi &= ~maskExp;
1676          aLo = 0u;
1677       }
1678    } else if (unbiasedExp > 51 || unbiasedExp == 1024) {
1679       return __a;
1680    } else {
1681       uint maskExp = 0xFFFFFFFFu >> (unbiasedExp - 20);
1682       if ((aLo & maskExp) == 0u)
1683          return __a;
1684       uint tmp = aLo + (1u << (51 - unbiasedExp));
1685       if(tmp < aLo)
1686          aHi += 1u;
1687       aLo = tmp;
1688       aLo &= ~maskExp;
1689    }
1690
1691    return packUint2x32(uvec2(aLo, aHi));
1692 }
1693
1694 uint64_t
1695 __fmin64(uint64_t a, uint64_t b)
1696 {
1697    if (__is_nan(a)) return b;
1698    if (__is_nan(b)) return a;
1699
1700    if (__flt64_nonnan(a, b)) return a;
1701    return b;
1702 }
1703
1704 uint64_t
1705 __fmax64(uint64_t a, uint64_t b)
1706 {
1707    if (__is_nan(a)) return b;
1708    if (__is_nan(b)) return a;
1709
1710    if (__flt64_nonnan(a, b)) return b;
1711    return a;
1712 }
1713
1714 uint64_t
1715 __ffract64(uint64_t a)
1716 {
1717    return __fadd64(a, __fneg64(__ffloor64(a)));
1718 }