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:
6 * License for Berkeley SoftFloat Release 3e
11 * The following applies to the whole of SoftFloat Release 3e as well as to
12 * each source file individually.
14 * Copyright 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018 The Regents of the
15 * University of California. All rights reserved.
17 * Redistribution and use in source and binary forms, with or without
18 * modification, are permitted provided that the following conditions are met:
20 * 1. Redistributions of source code must retain the above copyright notice,
21 * this list of conditions, and the following disclaimer.
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.
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.
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.
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
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.
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
62 /* Absolute value of a Float64 :
66 __fabs64(uint64_t __a)
68 uvec2 a = unpackUint2x32(__a);
70 return packUint2x32(a);
73 /* Returns 1 if the double-precision floating-point value `a' is a NaN;
74 * otherwise returns 0.
77 __is_nan(uint64_t __a)
79 uvec2 a = unpackUint2x32(__a);
80 return (0xFFE00000u <= (a.y<<1)) &&
81 ((a.x != 0u) || ((a.y & 0x000FFFFFu) != 0u));
84 /* Negate value of a Float64 :
88 __fneg64(uint64_t __a)
90 uvec2 a = unpackUint2x32(__a);
94 a.y = mix(t, a.y, __is_nan(__a));
95 return packUint2x32(a);
99 __fsign64(uint64_t __a)
101 uvec2 a = unpackUint2x32(__a);
104 retval.y = mix((a.y & 0x80000000u) | 0x3FF00000u, 0u, (a.y << 1 | a.x) == 0u);
105 return packUint2x32(retval);
108 /* Returns the fraction bits of the double-precision floating-point value `a'.*/
110 __extractFloat64FracLo(uint64_t a)
112 return unpackUint2x32(a).x;
116 __extractFloat64FracHi(uint64_t a)
118 return unpackUint2x32(a).y & 0x000FFFFFu;
121 /* Returns the exponent bits of the double-precision floating-point value `a'.*/
123 __extractFloat64Exp(uint64_t __a)
125 uvec2 a = unpackUint2x32(__a);
126 return int((a.y>>20) & 0x7FFu);
130 __feq64_nonnan(uint64_t __a, uint64_t __b)
132 uvec2 a = unpackUint2x32(__a);
133 uvec2 b = unpackUint2x32(__b);
134 return (a.x == b.x) &&
135 ((a.y == b.y) || ((a.x == 0u) && (((a.y | b.y)<<1) == 0u)));
138 /* Returns true if the double-precision floating-point value `a' is equal to the
139 * corresponding value `b', and false otherwise. The comparison is performed
140 * according to the IEEE Standard for Floating-Point Arithmetic.
143 __feq64(uint64_t a, uint64_t b)
145 if (__is_nan(a) || __is_nan(b))
148 return __feq64_nonnan(a, b);
151 /* Returns true if the double-precision floating-point value `a' is not equal
152 * to the corresponding value `b', and false otherwise. The comparison is
153 * performed according to the IEEE Standard for Floating-Point Arithmetic.
156 __fne64(uint64_t a, uint64_t b)
158 if (__is_nan(a) || __is_nan(b))
161 return !__feq64_nonnan(a, b);
164 /* Returns the sign bit of the double-precision floating-point value `a'.*/
166 __extractFloat64Sign(uint64_t a)
168 return unpackUint2x32(a).y >> 31;
171 /* Returns true if the 64-bit value formed by concatenating `a0' and `a1' is less
172 * than the 64-bit value formed by concatenating `b0' and `b1'. Otherwise,
176 lt64(uint a0, uint a1, uint b0, uint b1)
178 return (a0 < b0) || ((a0 == b0) && (a1 < b1));
182 __flt64_nonnan(uint64_t __a, uint64_t __b)
184 uvec2 a = unpackUint2x32(__a);
185 uvec2 b = unpackUint2x32(__b);
186 uint aSign = __extractFloat64Sign(__a);
187 uint bSign = __extractFloat64Sign(__b);
189 return (aSign != 0u) && ((((a.y | b.y)<<1) | a.x | b.x) != 0u);
191 return mix(lt64(a.y, a.x, b.y, b.x), lt64(b.y, b.x, a.y, a.x), aSign != 0u);
194 /* Returns true if the double-precision floating-point value `a' is less than
195 * the corresponding value `b', and false otherwise. The comparison is performed
196 * according to the IEEE Standard for Floating-Point Arithmetic.
199 __flt64(uint64_t a, uint64_t b)
201 if (__is_nan(a) || __is_nan(b))
204 return __flt64_nonnan(a, b);
207 /* Returns true if the double-precision floating-point value `a' is greater
208 * than or equal to * the corresponding value `b', and false otherwise. The
209 * comparison is performed * according to the IEEE Standard for Floating-Point
213 __fge64(uint64_t a, uint64_t b)
215 if (__is_nan(a) || __is_nan(b))
218 return !__flt64_nonnan(a, b);
221 /* Adds the 64-bit value formed by concatenating `a0' and `a1' to the 64-bit
222 * value formed by concatenating `b0' and `b1'. Addition is modulo 2^64, so
223 * any carry out is lost. The result is broken into two 32-bit pieces which
224 * are stored at the locations pointed to by `z0Ptr' and `z1Ptr'.
227 __add64(uint a0, uint a1, uint b0, uint b1,
233 z0Ptr = a0 + b0 + uint(z1 < a1);
237 /* Subtracts the 64-bit value formed by concatenating `b0' and `b1' from the
238 * 64-bit value formed by concatenating `a0' and `a1'. Subtraction is modulo
239 * 2^64, so any borrow out (carry out) is lost. The result is broken into two
240 * 32-bit pieces which are stored at the locations pointed to by `z0Ptr' and
244 __sub64(uint a0, uint a1, uint b0, uint b1,
249 z0Ptr = a0 - b0 - uint(a1 < b1);
252 /* Shifts the 64-bit value formed by concatenating `a0' and `a1' right by the
253 * number of bits given in `count'. If any nonzero bits are shifted off, they
254 * are "jammed" into the least significant bit of the result by setting the
255 * least significant bit to 1. The value of `count' can be arbitrarily large;
256 * in particular, if `count' is greater than 64, the result will be either 0
257 * or 1, depending on whether the concatenation of `a0' and `a1' is zero or
258 * nonzero. The result is broken into two 32-bit pieces which are stored at
259 * the locations pointed to by `z0Ptr' and `z1Ptr'.
262 __shift64RightJamming(uint a0,
270 int negCount = (-count) & 31;
272 z0 = mix(0u, a0, count == 0);
273 z0 = mix(z0, (a0 >> count), count < 32);
275 z1 = uint((a0 | a1) != 0u); /* count >= 64 */
276 uint z1_lt64 = (a0>>(count & 31)) | uint(((a0<<negCount) | a1) != 0u);
277 z1 = mix(z1, z1_lt64, count < 64);
278 z1 = mix(z1, (a0 | uint(a1 != 0u)), count == 32);
279 uint z1_lt32 = (a0<<negCount) | (a1>>count) | uint ((a1<<negCount) != 0u);
280 z1 = mix(z1, z1_lt32, count < 32);
281 z1 = mix(z1, a1, count == 0);
286 /* Shifts the 96-bit value formed by concatenating `a0', `a1', and `a2' right
287 * by 32 _plus_ the number of bits given in `count'. The shifted result is
288 * at most 64 nonzero bits; these are broken into two 32-bit pieces which are
289 * stored at the locations pointed to by `z0Ptr' and `z1Ptr'. The bits shifted
290 * off form a third 32-bit result as follows: The _last_ bit shifted off is
291 * the most-significant bit of the extra result, and the other 31 bits of the
292 * extra result are all zero if and only if _all_but_the_last_ bits shifted off
293 * were all zero. This extra result is stored in the location pointed to by
294 * `z2Ptr'. The value of `count' can be arbitrarily large.
295 * (This routine makes more sense if `a0', `a1', and `a2' are considered
296 * to form a fixed-point value with binary point between `a1' and `a2'. This
297 * fixed-point value is shifted right by the number of bits given in `count',
298 * and the integer part of the result is returned at the locations pointed to
299 * by `z0Ptr' and `z1Ptr'. The fractional part of the result may be slightly
300 * corrupted as described above, and is returned at the location pointed to by
304 __shift64ExtraRightJamming(uint a0, uint a1, uint a2,
313 int negCount = (-count) & 31;
315 z2 = mix(uint(a0 != 0u), a0, count == 64);
316 z2 = mix(z2, a0 << negCount, count < 64);
317 z2 = mix(z2, a1 << negCount, count < 32);
319 z1 = mix(0u, (a0 >> (count & 31)), count < 64);
320 z1 = mix(z1, (a0<<negCount) | (a1>>count), count < 32);
322 a2 = mix(a2 | a1, a2, count < 32);
323 z0 = mix(z0, a0 >> count, count < 32);
324 z2 |= uint(a2 != 0u);
326 z0 = mix(z0, 0u, (count == 32));
327 z1 = mix(z1, a0, (count == 32));
328 z2 = mix(z2, a1, (count == 32));
329 z0 = mix(z0, a0, (count == 0));
330 z1 = mix(z1, a1, (count == 0));
331 z2 = mix(z2, a2, (count == 0));
337 /* Shifts the 64-bit value formed by concatenating `a0' and `a1' left by the
338 * number of bits given in `count'. Any bits shifted off are lost. The value
339 * of `count' must be less than 32. The result is broken into two 32-bit
340 * pieces which are stored at the locations pointed to by `z0Ptr' and `z1Ptr'.
343 __shortShift64Left(uint a0, uint a1,
349 z0Ptr = mix((a0 << count | (a1 >> ((-count) & 31))), a0, count == 0);
352 /* Packs the sign `zSign', the exponent `zExp', and the significand formed by
353 * the concatenation of `zFrac0' and `zFrac1' into a double-precision floating-
354 * point value, returning the result. After being shifted into the proper
355 * positions, the three fields `zSign', `zExp', and `zFrac0' are simply added
356 * together to form the most significant 32 bits of the result. This means
357 * that any integer portion of `zFrac0' will be added into the exponent. Since
358 * a properly normalized significand will have an integer portion equal to 1,
359 * the `zExp' input should be 1 less than the desired result exponent whenever
360 * `zFrac0' and `zFrac1' concatenated form a complete, normalized significand.
363 __packFloat64(uint zSign, int zExp, uint zFrac0, uint zFrac1)
367 z.y = (zSign << 31) + (uint(zExp) << 20) + zFrac0;
369 return packUint2x32(z);
372 /* Takes an abstract floating-point value having sign `zSign', exponent `zExp',
373 * and extended significand formed by the concatenation of `zFrac0', `zFrac1',
374 * and `zFrac2', and returns the proper double-precision floating-point value
375 * corresponding to the abstract input. Ordinarily, the abstract value is
376 * simply rounded and packed into the double-precision format, with the inexact
377 * exception raised if the abstract input cannot be represented exactly.
378 * However, if the abstract value is too large, the overflow and inexact
379 * exceptions are raised and an infinity or maximal finite value is returned.
380 * If the abstract value is too small, the input value is rounded to a
381 * subnormal number, and the underflow and inexact exceptions are raised if the
382 * abstract input cannot be represented exactly as a subnormal double-precision
383 * floating-point number.
384 * The input significand must be normalized or smaller. If the input
385 * significand is not normalized, `zExp' must be 0; in that case, the result
386 * returned is a subnormal number, and it must not require rounding. In the
387 * usual case that the input significand is normalized, `zExp' must be 1 less
388 * than the "true" floating-point exponent. The handling of underflow and
389 * overflow follows the IEEE Standard for Floating-Point Arithmetic.
392 __roundAndPackFloat64(uint zSign,
398 bool roundNearestEven;
401 roundNearestEven = FLOAT_ROUNDING_MODE == FLOAT_ROUND_NEAREST_EVEN;
402 increment = int(zFrac2) < 0;
403 if (!roundNearestEven) {
404 if (FLOAT_ROUNDING_MODE == FLOAT_ROUND_TO_ZERO) {
408 increment = (FLOAT_ROUNDING_MODE == FLOAT_ROUND_DOWN) &&
411 increment = (FLOAT_ROUNDING_MODE == FLOAT_ROUND_UP) &&
417 if ((0x7FD < zExp) ||
419 (0x001FFFFFu == zFrac0 && 0xFFFFFFFFu == zFrac1) &&
421 if ((FLOAT_ROUNDING_MODE == FLOAT_ROUND_TO_ZERO) ||
422 ((zSign != 0u) && (FLOAT_ROUNDING_MODE == FLOAT_ROUND_UP)) ||
423 ((zSign == 0u) && (FLOAT_ROUNDING_MODE == FLOAT_ROUND_DOWN))) {
424 return __packFloat64(zSign, 0x7FE, 0x000FFFFFu, 0xFFFFFFFFu);
426 return __packFloat64(zSign, 0x7FF, 0u, 0u);
429 __shift64ExtraRightJamming(
430 zFrac0, zFrac1, zFrac2, -zExp, zFrac0, zFrac1, zFrac2);
432 if (roundNearestEven) {
433 increment = zFrac2 < 0u;
436 increment = (FLOAT_ROUNDING_MODE == FLOAT_ROUND_DOWN) &&
439 increment = (FLOAT_ROUNDING_MODE == FLOAT_ROUND_UP) &&
446 __add64(zFrac0, zFrac1, 0u, 1u, zFrac0, zFrac1);
447 zFrac1 &= ~((zFrac2 + uint(zFrac2 == 0u)) & uint(roundNearestEven));
449 zExp = mix(zExp, 0, (zFrac0 | zFrac1) == 0u);
451 return __packFloat64(zSign, zExp, zFrac0, zFrac1);
455 __roundAndPackUInt64(uint zSign, uint zFrac0, uint zFrac1, uint zFrac2)
457 bool roundNearestEven;
459 uint64_t default_nan = 0xFFFFFFFFFFFFFFFFUL;
461 roundNearestEven = FLOAT_ROUNDING_MODE == FLOAT_ROUND_NEAREST_EVEN;
463 if (zFrac2 >= 0x80000000u)
466 if (!roundNearestEven) {
468 if ((FLOAT_ROUNDING_MODE == FLOAT_ROUND_DOWN) && (zFrac2 != 0u)) {
472 increment = (FLOAT_ROUNDING_MODE == FLOAT_ROUND_UP) &&
478 __add64(zFrac0, zFrac1, 0u, 1u, zFrac0, zFrac1);
479 if ((zFrac0 | zFrac1) != 0u)
480 zFrac1 &= ~(1u) + uint(zFrac2 == 0u) & uint(roundNearestEven);
482 return mix(packUint2x32(uvec2(zFrac1, zFrac0)), default_nan,
483 (zSign !=0u && (zFrac0 | zFrac1) != 0u));
487 __roundAndPackInt64(uint zSign, uint zFrac0, uint zFrac1, uint zFrac2)
489 bool roundNearestEven;
491 int64_t default_NegNaN = -0x7FFFFFFFFFFFFFFEL;
492 int64_t default_PosNaN = 0xFFFFFFFFFFFFFFFFL;
494 roundNearestEven = FLOAT_ROUNDING_MODE == FLOAT_ROUND_NEAREST_EVEN;
496 if (zFrac2 >= 0x80000000u)
499 if (!roundNearestEven) {
501 increment = ((FLOAT_ROUNDING_MODE == FLOAT_ROUND_DOWN) &&
504 increment = (FLOAT_ROUNDING_MODE == FLOAT_ROUND_UP) &&
510 __add64(zFrac0, zFrac1, 0u, 1u, zFrac0, zFrac1);
511 if ((zFrac0 | zFrac1) != 0u)
512 zFrac1 &= ~(1u) + uint(zFrac2 == 0u) & uint(roundNearestEven);
515 int64_t absZ = mix(int64_t(packUint2x32(uvec2(zFrac1, zFrac0))),
516 -int64_t(packUint2x32(uvec2(zFrac1, zFrac0))),
518 int64_t nan = mix(default_PosNaN, default_NegNaN, bool(zSign));
519 return mix(absZ, nan, bool(zSign ^ uint(absZ < 0)) && bool(absZ));
522 /* Returns the number of leading 0 bits before the most-significant 1 bit of
523 * `a'. If `a' is zero, 32 is returned.
526 __countLeadingZeros32(uint a)
529 shiftCount = mix(31 - findMSB(a), 32, a == 0u);
533 /* Takes an abstract floating-point value having sign `zSign', exponent `zExp',
534 * and significand formed by the concatenation of `zSig0' and `zSig1', and
535 * returns the proper double-precision floating-point value corresponding
536 * to the abstract input. This routine is just like `__roundAndPackFloat64'
537 * except that the input significand has fewer bits and does not have to be
538 * normalized. In all cases, `zExp' must be 1 less than the "true" floating-
542 __normalizeRoundAndPackFloat64(uint zSign,
556 shiftCount = __countLeadingZeros32(zFrac0) - 11;
557 if (0 <= shiftCount) {
559 __shortShift64Left(zFrac0, zFrac1, shiftCount, zFrac0, zFrac1);
561 __shift64ExtraRightJamming(
562 zFrac0, zFrac1, 0u, -shiftCount, zFrac0, zFrac1, zFrac2);
565 return __roundAndPackFloat64(zSign, zExp, zFrac0, zFrac1, zFrac2);
568 /* Takes two double-precision floating-point values `a' and `b', one of which
569 * is a NaN, and returns the appropriate NaN result.
572 __propagateFloat64NaN(uint64_t __a, uint64_t __b)
574 bool aIsNaN = __is_nan(__a);
575 bool bIsNaN = __is_nan(__b);
576 uvec2 a = unpackUint2x32(__a);
577 uvec2 b = unpackUint2x32(__b);
581 return packUint2x32(mix(b, mix(a, b, bvec2(bIsNaN, bIsNaN)), bvec2(aIsNaN, aIsNaN)));
584 /* Returns the result of adding the double-precision floating-point values
585 * `a' and `b'. The operation is performed according to the IEEE Standard for
586 * Floating-Point Arithmetic.
589 __fadd64(uint64_t a, uint64_t b)
591 uint aSign = __extractFloat64Sign(a);
592 uint bSign = __extractFloat64Sign(b);
593 uint aFracLo = __extractFloat64FracLo(a);
594 uint aFracHi = __extractFloat64FracHi(a);
595 uint bFracLo = __extractFloat64FracLo(b);
596 uint bFracHi = __extractFloat64FracHi(b);
597 int aExp = __extractFloat64Exp(a);
598 int bExp = __extractFloat64Exp(b);
601 int expDiff = aExp - bExp;
602 if (aSign == bSign) {
605 bool orig_exp_diff_is_zero = (expDiff == 0);
607 if (orig_exp_diff_is_zero) {
609 bool propagate = (aFracHi | aFracLo | bFracHi | bFracLo) != 0u;
610 return mix(a, __propagateFloat64NaN(a, b), propagate);
612 __add64(aFracHi, aFracLo, bFracHi, bFracLo, zFrac0, zFrac1);
614 return __packFloat64(aSign, 0, zFrac0, zFrac1);
616 zFrac0 |= 0x00200000u;
618 __shift64ExtraRightJamming(
619 zFrac0, zFrac1, zFrac2, 1, zFrac0, zFrac1, zFrac2);
620 } else if (0 < expDiff) {
622 bool propagate = (aFracHi | aFracLo) != 0u;
623 return mix(a, __propagateFloat64NaN(a, b), propagate);
626 expDiff = mix(expDiff, expDiff - 1, bExp == 0);
627 bFracHi = mix(bFracHi | 0x00100000u, bFracHi, bExp == 0);
628 __shift64ExtraRightJamming(
629 bFracHi, bFracLo, 0u, expDiff, bFracHi, bFracLo, zFrac2);
631 } else if (expDiff < 0) {
633 bool propagate = (bFracHi | bFracLo) != 0u;
634 return mix(__packFloat64(aSign, 0x7ff, 0u, 0u), __propagateFloat64NaN(a, b), propagate);
636 expDiff = mix(expDiff, expDiff + 1, aExp == 0);
637 aFracHi = mix(aFracHi | 0x00100000u, aFracHi, aExp == 0);
638 __shift64ExtraRightJamming(
639 aFracHi, aFracLo, 0u, - expDiff, aFracHi, aFracLo, zFrac2);
642 if (!orig_exp_diff_is_zero) {
643 aFracHi |= 0x00100000u;
644 __add64(aFracHi, aFracLo, bFracHi, bFracLo, zFrac0, zFrac1);
646 if (!(zFrac0 < 0x00200000u)) {
647 __shift64ExtraRightJamming(zFrac0, zFrac1, zFrac2, 1, zFrac0, zFrac1, zFrac2);
651 return __roundAndPackFloat64(aSign, zExp, zFrac0, zFrac1, zFrac2);
656 __shortShift64Left(aFracHi, aFracLo, 10, aFracHi, aFracLo);
657 __shortShift64Left(bFracHi, bFracLo, 10, bFracHi, bFracLo);
660 bool propagate = (aFracHi | aFracLo) != 0u;
661 return mix(a, __propagateFloat64NaN(a, b), propagate);
663 expDiff = mix(expDiff, expDiff - 1, bExp == 0);
664 bFracHi = mix(bFracHi | 0x40000000u, bFracHi, bExp == 0);
665 __shift64RightJamming(bFracHi, bFracLo, expDiff, bFracHi, bFracLo);
666 aFracHi |= 0x40000000u;
667 __sub64(aFracHi, aFracLo, bFracHi, bFracLo, zFrac0, zFrac1);
670 return __normalizeRoundAndPackFloat64(aSign, zExp - 10, zFrac0, zFrac1);
674 bool propagate = (bFracHi | bFracLo) != 0u;
675 return mix(__packFloat64(aSign ^ 1u, 0x7ff, 0u, 0u), __propagateFloat64NaN(a, b), propagate);
677 expDiff = mix(expDiff, expDiff + 1, aExp == 0);
678 aFracHi = mix(aFracHi | 0x40000000u, aFracHi, aExp == 0);
679 __shift64RightJamming(aFracHi, aFracLo, - expDiff, aFracHi, aFracLo);
680 bFracHi |= 0x40000000u;
681 __sub64(bFracHi, bFracLo, aFracHi, aFracLo, zFrac0, zFrac1);
685 return __normalizeRoundAndPackFloat64(aSign, zExp - 10, zFrac0, zFrac1);
688 bool propagate = (aFracHi | aFracLo | bFracHi | bFracLo) != 0u;
689 return mix(0xFFFFFFFFFFFFFFFFUL, __propagateFloat64NaN(a, b), propagate);
691 bExp = mix(bExp, 1, aExp == 0);
692 aExp = mix(aExp, 1, aExp == 0);
693 bool zexp_normal = false;
695 if (bFracHi < aFracHi) {
696 __sub64(aFracHi, aFracLo, bFracHi, bFracLo, zFrac0, zFrac1);
699 else if (aFracHi < bFracHi) {
700 __sub64(bFracHi, bFracLo, aFracHi, aFracLo, zFrac0, zFrac1);
704 else if (bFracLo < aFracLo) {
705 __sub64(aFracHi, aFracLo, bFracHi, bFracLo, zFrac0, zFrac1);
708 else if (aFracLo < bFracLo) {
709 __sub64(bFracHi, bFracLo, aFracHi, aFracLo, zFrac0, zFrac1);
713 zExp = mix(bExp, aExp, blta);
714 aSign = mix(aSign ^ 1u, aSign, blta);
715 uint64_t retval_0 = __packFloat64(uint(FLOAT_ROUNDING_MODE == FLOAT_ROUND_DOWN), 0, 0u, 0u);
716 uint64_t retval_1 = __normalizeRoundAndPackFloat64(aSign, zExp - 11, zFrac0, zFrac1);
717 return mix(retval_0, retval_1, zexp_normal);
721 /* Multiplies `a' by `b' to obtain a 64-bit product. The product is broken
722 * into two 32-bit pieces which are stored at the locations pointed to by
723 * `z0Ptr' and `z1Ptr'.
726 __mul32To64(uint a, uint b, out uint z0Ptr, out uint z1Ptr)
728 uint aLow = a & 0x0000FFFFu;
730 uint bLow = b & 0x0000FFFFu;
732 uint z1 = aLow * bLow;
733 uint zMiddleA = aLow * bHigh;
734 uint zMiddleB = aHigh * bLow;
735 uint z0 = aHigh * bHigh;
736 zMiddleA += zMiddleB;
737 z0 += ((uint(zMiddleA < zMiddleB)) << 16) + (zMiddleA >> 16);
740 z0 += uint(z1 < zMiddleA);
745 /* Multiplies the 64-bit value formed by concatenating `a0' and `a1' to the
746 * 64-bit value formed by concatenating `b0' and `b1' to obtain a 128-bit
747 * product. The product is broken into four 32-bit pieces which are stored at
748 * the locations pointed to by `z0Ptr', `z1Ptr', `z2Ptr', and `z3Ptr'.
751 __mul64To128(uint a0, uint a1, uint b0, uint b1,
764 __mul32To64(a1, b1, z2, z3);
765 __mul32To64(a1, b0, z1, more2);
766 __add64(z1, more2, 0u, z2, z1, z2);
767 __mul32To64(a0, b0, z0, more1);
768 __add64(z0, more1, 0u, z1, z0, z1);
769 __mul32To64(a0, b1, more1, more2);
770 __add64(more1, more2, 0u, z2, more1, z2);
771 __add64(z0, z1, 0u, more1, z0, z1);
778 /* Normalizes the subnormal double-precision floating-point value represented
779 * by the denormalized significand formed by the concatenation of `aFrac0' and
780 * `aFrac1'. The normalized exponent is stored at the location pointed to by
781 * `zExpPtr'. The most significant 21 bits of the normalized significand are
782 * stored at the location pointed to by `zFrac0Ptr', and the least significant
783 * 32 bits of the normalized significand are stored at the location pointed to
787 __normalizeFloat64Subnormal(uint aFrac0, uint aFrac1,
793 uint temp_zfrac0, temp_zfrac1;
794 shiftCount = __countLeadingZeros32(mix(aFrac0, aFrac1, aFrac0 == 0u)) - 11;
795 zExpPtr = mix(1 - shiftCount, -shiftCount - 31, aFrac0 == 0u);
797 temp_zfrac0 = mix(aFrac1<<shiftCount, aFrac1>>(-shiftCount), shiftCount < 0);
798 temp_zfrac1 = mix(0u, aFrac1<<(shiftCount & 31), shiftCount < 0);
800 __shortShift64Left(aFrac0, aFrac1, shiftCount, zFrac0Ptr, zFrac1Ptr);
802 zFrac0Ptr = mix(zFrac0Ptr, temp_zfrac0, aFrac0 == 0);
803 zFrac1Ptr = mix(zFrac1Ptr, temp_zfrac1, aFrac0 == 0);
806 /* Returns the result of multiplying the double-precision floating-point values
807 * `a' and `b'. The operation is performed according to the IEEE Standard for
808 * Floating-Point Arithmetic.
811 __fmul64(uint64_t a, uint64_t b)
819 uint aFracLo = __extractFloat64FracLo(a);
820 uint aFracHi = __extractFloat64FracHi(a);
821 uint bFracLo = __extractFloat64FracLo(b);
822 uint bFracHi = __extractFloat64FracHi(b);
823 int aExp = __extractFloat64Exp(a);
824 uint aSign = __extractFloat64Sign(a);
825 int bExp = __extractFloat64Exp(b);
826 uint bSign = __extractFloat64Sign(b);
827 uint zSign = aSign ^ bSign;
829 if (((aFracHi | aFracLo) != 0u) ||
830 ((bExp == 0x7FF) && ((bFracHi | bFracLo) != 0u))) {
831 return __propagateFloat64NaN(a, b);
833 if ((uint(bExp) | bFracHi | bFracLo) == 0u)
834 return 0xFFFFFFFFFFFFFFFFUL;
835 return __packFloat64(zSign, 0x7FF, 0u, 0u);
838 if ((bFracHi | bFracLo) != 0u)
839 return __propagateFloat64NaN(a, b);
840 if ((uint(aExp) | aFracHi | aFracLo) == 0u)
841 return 0xFFFFFFFFFFFFFFFFUL;
842 return __packFloat64(zSign, 0x7FF, 0u, 0u);
845 if ((aFracHi | aFracLo) == 0u)
846 return __packFloat64(zSign, 0, 0u, 0u);
847 __normalizeFloat64Subnormal(aFracHi, aFracLo, aExp, aFracHi, aFracLo);
850 if ((bFracHi | bFracLo) == 0u)
851 return __packFloat64(zSign, 0, 0u, 0u);
852 __normalizeFloat64Subnormal(bFracHi, bFracLo, bExp, bFracHi, bFracLo);
854 zExp = aExp + bExp - 0x400;
855 aFracHi |= 0x00100000u;
856 __shortShift64Left(bFracHi, bFracLo, 12, bFracHi, bFracLo);
858 aFracHi, aFracLo, bFracHi, bFracLo, zFrac0, zFrac1, zFrac2, zFrac3);
859 __add64(zFrac0, zFrac1, aFracHi, aFracLo, zFrac0, zFrac1);
860 zFrac2 |= uint(zFrac3 != 0u);
861 if (0x00200000u <= zFrac0) {
862 __shift64ExtraRightJamming(
863 zFrac0, zFrac1, zFrac2, 1, zFrac0, zFrac1, zFrac2);
866 return __roundAndPackFloat64(zSign, zExp, zFrac0, zFrac1, zFrac2);
870 __ffma64(uint64_t a, uint64_t b, uint64_t c)
872 return __fadd64(__fmul64(a, b), c);
875 /* Shifts the 64-bit value formed by concatenating `a0' and `a1' right by the
876 * number of bits given in `count'. Any bits shifted off are lost. The value
877 * of `count' can be arbitrarily large; in particular, if `count' is greater
878 * than 64, the result will be 0. The result is broken into two 32-bit pieces
879 * which are stored at the locations pointed to by `z0Ptr' and `z1Ptr'.
882 __shift64Right(uint a0, uint a1,
889 int negCount = (-count) & 31;
892 z0 = mix(z0, (a0 >> count), count < 32);
893 z0 = mix(z0, a0, count == 0);
895 z1 = mix(0u, (a0 >> (count & 31)), count < 64);
896 z1 = mix(z1, (a0<<negCount) | (a1>>count), count < 32);
897 z1 = mix(z1, a0, count == 0);
903 /* Returns the result of converting the double-precision floating-point value
904 * `a' to the unsigned integer format. The conversion is performed according
905 * to the IEEE Standard for Floating-Point Arithmetic.
908 __fp64_to_uint(uint64_t a)
910 uint aFracLo = __extractFloat64FracLo(a);
911 uint aFracHi = __extractFloat64FracHi(a);
912 int aExp = __extractFloat64Exp(a);
913 uint aSign = __extractFloat64Sign(a);
915 if ((aExp == 0x7FF) && ((aFracHi | aFracLo) != 0u))
918 aFracHi |= mix(0u, 0x00100000u, aExp != 0);
920 int shiftDist = 0x427 - aExp;
922 __shift64RightJamming(aFracHi, aFracLo, shiftDist, aFracHi, aFracLo);
924 if ((aFracHi & 0xFFFFF000u) != 0u)
925 return mix(~0u, 0u, (aSign != 0u));
929 __shift64Right(aFracHi, aFracLo, 12, zero, z);
931 uint expt = mix(~0u, 0u, (aSign != 0u));
933 return mix(z, expt, (aSign != 0u) && (z != 0u));
937 __uint_to_fp64(uint a)
942 int shiftDist = __countLeadingZeros32(a) + 21;
946 int negCount = (- shiftDist) & 31;
948 aHigh = mix(0u, a<< shiftDist - 32, shiftDist < 64);
950 aHigh = mix(aHigh, 0u, shiftDist == 0);
951 aLow = mix(aLow, a, shiftDist ==0);
952 aHigh = mix(aHigh, a >> negCount, shiftDist < 32);
953 aLow = mix(aLow, a << shiftDist, shiftDist < 32);
955 return __packFloat64(0u, 0x432 - shiftDist, aHigh, aLow);
959 __uint64_to_fp64(uint64_t a)
964 uvec2 aFrac = unpackUint2x32(a);
965 uint aFracLo = __extractFloat64FracLo(a);
966 uint aFracHi = __extractFloat64FracHi(a);
968 if ((aFracHi & 0x80000000u) != 0u) {
969 __shift64RightJamming(aFracHi, aFracLo, 1, aFracHi, aFracLo);
970 return __roundAndPackFloat64(0, 0x433, aFracHi, aFracLo, 0u);
972 return __normalizeRoundAndPackFloat64(0, 0x432, aFrac.y, aFrac.x);
977 __fp64_to_uint64(uint64_t a)
979 uint aFracLo = __extractFloat64FracLo(a);
980 uint aFracHi = __extractFloat64FracHi(a);
981 int aExp = __extractFloat64Exp(a);
982 uint aSign = __extractFloat64Sign(a);
984 uint64_t default_nan = 0xFFFFFFFFFFFFFFFFUL;
986 aFracHi = mix(aFracHi, aFracHi | 0x00100000u, aExp != 0);
987 int shiftCount = 0x433 - aExp;
989 if ( shiftCount <= 0 ) {
990 if (shiftCount < -11 && aExp == 0x7FF) {
991 if ((aFracHi | aFracLo) != 0u)
992 return __propagateFloat64NaN(a, a);
993 return mix(default_nan, a, aSign == 0u);
995 __shortShift64Left(aFracHi, aFracLo, -shiftCount, aFracHi, aFracLo);
997 __shift64ExtraRightJamming(aFracHi, aFracLo, zFrac2, shiftCount,
998 aFracHi, aFracLo, zFrac2);
1000 return __roundAndPackUInt64(aSign, aFracHi, aFracLo, zFrac2);
1004 __fp64_to_int64(uint64_t a)
1007 uint aFracLo = __extractFloat64FracLo(a);
1008 uint aFracHi = __extractFloat64FracHi(a);
1009 int aExp = __extractFloat64Exp(a);
1010 uint aSign = __extractFloat64Sign(a);
1011 int64_t default_NegNaN = -0x7FFFFFFFFFFFFFFEL;
1012 int64_t default_PosNaN = 0xFFFFFFFFFFFFFFFFL;
1014 aFracHi = mix(aFracHi, aFracHi | 0x00100000u, aExp != 0);
1015 int shiftCount = 0x433 - aExp;
1017 if (shiftCount <= 0) {
1018 if (shiftCount < -11 && aExp == 0x7FF) {
1019 if ((aFracHi | aFracLo) != 0u)
1020 return default_NegNaN;
1021 return mix(default_NegNaN, default_PosNaN, aSign == 0u);
1023 __shortShift64Left(aFracHi, aFracLo, -shiftCount, aFracHi, aFracLo);
1025 __shift64ExtraRightJamming(aFracHi, aFracLo, zFrac2, shiftCount,
1026 aFracHi, aFracLo, zFrac2);
1029 return __roundAndPackInt64(aSign, aFracHi, aFracLo, zFrac2);
1033 __fp32_to_uint64(float f)
1035 uint a = floatBitsToUint(f);
1036 uint aFrac = a & 0x007FFFFFu;
1037 int aExp = int((a>>23) & 0xFFu);
1042 uint64_t default_nan = 0xFFFFFFFFFFFFFFFFUL;
1043 int shiftCount = 0xBE - aExp;
1045 if (shiftCount <0) {
1050 aFrac = mix(aFrac, aFrac | 0x00800000u, aExp != 0);
1051 __shortShift64Left(aFrac, 0, 40, zFrac0, zFrac1);
1053 if (shiftCount != 0) {
1054 __shift64ExtraRightJamming(zFrac0, zFrac1, zFrac2, shiftCount,
1055 zFrac0, zFrac1, zFrac2);
1058 return __roundAndPackUInt64(aSign, zFrac0, zFrac1, zFrac2);
1062 __int64_to_fp64(int64_t a)
1067 uint64_t absA = mix(uint64_t(a), uint64_t(-a), a < 0);
1068 uint aFracHi = __extractFloat64FracHi(absA);
1069 uvec2 aFrac = unpackUint2x32(absA);
1070 uint zSign = uint(a < 0);
1072 if ((aFracHi & 0x80000000u) != 0u) {
1073 return mix(0ul, __packFloat64(1, 0x434, 0u, 0u), a < 0);
1076 return __normalizeRoundAndPackFloat64(zSign, 0x432, aFrac.y, aFrac.x);
1079 /* Returns the result of converting the double-precision floating-point value
1080 * `a' to the 32-bit two's complement integer format. The conversion is
1081 * performed according to the IEEE Standard for Floating-Point Arithmetic---
1082 * which means in particular that the conversion is rounded according to the
1083 * current rounding mode. If `a' is a NaN, the largest positive integer is
1084 * returned. Otherwise, if the conversion overflows, the largest integer with
1085 * the same sign as `a' is returned.
1088 __fp64_to_int(uint64_t a)
1090 uint aFracLo = __extractFloat64FracLo(a);
1091 uint aFracHi = __extractFloat64FracHi(a);
1092 int aExp = __extractFloat64Exp(a);
1093 uint aSign = __extractFloat64Sign(a);
1096 uint aFracExtra = 0u;
1097 int shiftCount = aExp - 0x413;
1099 if (0 <= shiftCount) {
1101 if ((aExp == 0x7FF) && bool(aFracHi | aFracLo))
1103 return mix(0x7FFFFFFF, 0x80000000, bool(aSign));
1105 __shortShift64Left(aFracHi | 0x00100000u, aFracLo, shiftCount, absZ, aFracExtra);
1110 aFracHi |= 0x00100000u;
1111 aFracExtra = ( aFracHi << (shiftCount & 31)) | aFracLo;
1112 absZ = aFracHi >> (- shiftCount);
1115 int z = mix(int(absZ), -int(absZ), (aSign != 0u));
1116 int nan = mix(0x7FFFFFFF, 0x80000000, bool(aSign));
1117 return mix(z, nan, bool(aSign ^ uint(z < 0)) && bool(z));
1120 /* Returns the result of converting the 32-bit two's complement integer `a'
1121 * to the double-precision floating-point format. The conversion is performed
1122 * according to the IEEE Standard for Floating-Point Arithmetic.
1125 __int_to_fp64(int a)
1130 return __packFloat64(0u, 0, 0u, 0u);
1131 uint zSign = uint(a < 0);
1132 uint absA = mix(uint(a), uint(-a), a < 0);
1133 int shiftCount = __countLeadingZeros32(absA) - 11;
1134 if (0 <= shiftCount) {
1135 zFrac0 = absA << shiftCount;
1138 __shift64Right(absA, 0u, -shiftCount, zFrac0, zFrac1);
1140 return __packFloat64(zSign, 0x412 - shiftCount, zFrac0, zFrac1);
1144 __fp64_to_bool(uint64_t a)
1146 return !__feq64_nonnan(__fabs64(a), 0ul);
1150 __bool_to_fp64(bool a)
1152 return __int_to_fp64(int(a));
1155 /* Packs the sign `zSign', exponent `zExp', and significand `zFrac' into a
1156 * single-precision floating-point value, returning the result. After being
1157 * shifted into the proper positions, the three fields are simply added
1158 * together to form the result. This means that any integer portion of `zSig'
1159 * will be added into the exponent. Since a properly normalized significand
1160 * will have an integer portion equal to 1, the `zExp' input should be 1 less
1161 * than the desired result exponent whenever `zFrac' is a complete, normalized
1165 __packFloat32(uint zSign, int zExp, uint zFrac)
1167 return uintBitsToFloat((zSign<<31) + (uint(zExp)<<23) + zFrac);
1170 /* Takes an abstract floating-point value having sign `zSign', exponent `zExp',
1171 * and significand `zFrac', and returns the proper single-precision floating-
1172 * point value corresponding to the abstract input. Ordinarily, the abstract
1173 * value is simply rounded and packed into the single-precision format, with
1174 * the inexact exception raised if the abstract input cannot be represented
1175 * exactly. However, if the abstract value is too large, the overflow and
1176 * inexact exceptions are raised and an infinity or maximal finite value is
1177 * returned. If the abstract value is too small, the input value is rounded to
1178 * a subnormal number, and the underflow and inexact exceptions are raised if
1179 * the abstract input cannot be represented exactly as a subnormal single-
1180 * precision floating-point number.
1181 * The input significand `zFrac' has its binary point between bits 30
1182 * and 29, which is 7 bits to the left of the usual location. This shifted
1183 * significand must be normalized or smaller. If `zFrac' is not normalized,
1184 * `zExp' must be 0; in that case, the result returned is a subnormal number,
1185 * and it must not require rounding. In the usual case that `zFrac' is
1186 * normalized, `zExp' must be 1 less than the "true" floating-point exponent.
1187 * The handling of underflow and overflow follows the IEEE Standard for
1188 * Floating-Point Arithmetic.
1191 __roundAndPackFloat32(uint zSign, int zExp, uint zFrac)
1193 bool roundNearestEven;
1197 roundNearestEven = FLOAT_ROUNDING_MODE == FLOAT_ROUND_NEAREST_EVEN;
1198 roundIncrement = 0x40;
1199 if (!roundNearestEven) {
1200 if (FLOAT_ROUNDING_MODE == FLOAT_ROUND_TO_ZERO) {
1203 roundIncrement = 0x7F;
1205 if (FLOAT_ROUNDING_MODE == FLOAT_ROUND_UP)
1208 if (FLOAT_ROUNDING_MODE == FLOAT_ROUND_DOWN)
1213 roundBits = int(zFrac & 0x7Fu);
1214 if (0xFDu <= uint(zExp)) {
1215 if ((0xFD < zExp) || ((zExp == 0xFD) && (int(zFrac) + roundIncrement) < 0))
1216 return __packFloat32(zSign, 0xFF, 0u) - float(roundIncrement == 0);
1218 bool zexp_lt0 = zExp < 0;
1219 uint zFrac_lt0 = mix(uint(zFrac != 0u), (zFrac>>count) | uint((zFrac<<((-count) & 31)) != 0u), (-zExp) < 32);
1220 zFrac = mix(zFrac, zFrac_lt0, zexp_lt0);
1221 roundBits = mix(roundBits, int(zFrac) & 0x7f, zexp_lt0);
1222 zExp = mix(zExp, 0, zexp_lt0);
1224 zFrac = (zFrac + uint(roundIncrement))>>7;
1225 zFrac &= ~uint(((roundBits ^ 0x40) == 0) && roundNearestEven);
1227 return __packFloat32(zSign, mix(zExp, 0, zFrac == 0u), zFrac);
1230 /* Returns the result of converting the double-precision floating-point value
1231 * `a' to the single-precision floating-point format. The conversion is
1232 * performed according to the IEEE Standard for Floating-Point Arithmetic.
1235 __fp64_to_fp32(uint64_t __a)
1237 uvec2 a = unpackUint2x32(__a);
1241 uint aFracLo = __extractFloat64FracLo(__a);
1242 uint aFracHi = __extractFloat64FracHi(__a);
1243 int aExp = __extractFloat64Exp(__a);
1244 uint aSign = __extractFloat64Sign(__a);
1245 if (aExp == 0x7FF) {
1246 __shortShift64Left(a.y, a.x, 12, a.y, a.x);
1247 float rval = uintBitsToFloat((aSign<<31) | 0x7FC00000u | (a.y>>9));
1248 rval = mix(__packFloat32(aSign, 0xFF, 0u), rval, (aFracHi | aFracLo) != 0u);
1251 __shift64RightJamming(aFracHi, aFracLo, 22, allZero, zFrac);
1252 zFrac = mix(zFrac, zFrac | 0x40000000u, aExp != 0);
1253 return __roundAndPackFloat32(aSign, aExp - 0x381, zFrac);
1257 __uint64_to_fp32(uint64_t __a)
1260 uvec2 aFrac = unpackUint2x32(__a);
1261 int shiftCount = __countLeadingZeros32(mix(aFrac.y, aFrac.x, aFrac.y == 0u));
1262 shiftCount -= mix(40, 8, aFrac.y == 0u);
1264 if (0 <= shiftCount) {
1265 __shortShift64Left(aFrac.y, aFrac.x, shiftCount, aFrac.y, aFrac.x);
1266 bool is_zero = (aFrac.y | aFrac.x) == 0u;
1267 return mix(__packFloat32(0u, 0x95 - shiftCount, aFrac.x), 0, is_zero);
1271 __shift64RightJamming(aFrac.y, aFrac.x, -shiftCount, aFrac.y, aFrac.x);
1272 zFrac = mix(aFrac.x<<shiftCount, aFrac.x, shiftCount < 0);
1273 return __roundAndPackFloat32(0u, 0x9C - shiftCount, zFrac);
1277 __int64_to_fp32(int64_t __a)
1280 uint aSign = uint(__a < 0);
1281 uint64_t absA = mix(uint64_t(__a), uint64_t(-__a), __a < 0);
1282 uvec2 aFrac = unpackUint2x32(absA);
1283 int shiftCount = __countLeadingZeros32(mix(aFrac.y, aFrac.x, aFrac.y == 0u));
1284 shiftCount -= mix(40, 8, aFrac.y == 0u);
1286 if (0 <= shiftCount) {
1287 __shortShift64Left(aFrac.y, aFrac.x, shiftCount, aFrac.y, aFrac.x);
1288 bool is_zero = (aFrac.y | aFrac.x) == 0u;
1289 return mix(__packFloat32(aSign, 0x95 - shiftCount, aFrac.x), 0, absA == 0u);
1293 __shift64RightJamming(aFrac.y, aFrac.x, -shiftCount, aFrac.y, aFrac.x);
1294 zFrac = mix(aFrac.x<<shiftCount, aFrac.x, shiftCount < 0);
1295 return __roundAndPackFloat32(aSign, 0x9C - shiftCount, zFrac);
1298 /* Returns the result of converting the single-precision floating-point value
1299 * `a' to the double-precision floating-point format.
1302 __fp32_to_fp64(float f)
1304 uint a = floatBitsToUint(f);
1305 uint aFrac = a & 0x007FFFFFu;
1306 int aExp = int((a>>23) & 0xFFu);
1315 __shift64Right(nanHi, nanLo, 12, nanHi, nanLo);
1316 nanHi |= ((aSign<<31) | 0x7FF80000u);
1317 return packUint2x32(uvec2(nanLo, nanHi));
1319 return __packFloat64(aSign, 0x7FF, 0u, 0u);
1324 return __packFloat64(aSign, 0, 0u, 0u);
1325 /* Normalize subnormal */
1326 int shiftCount = __countLeadingZeros32(aFrac) - 8;
1327 aFrac <<= shiftCount;
1328 aExp = 1 - shiftCount;
1332 __shift64Right(aFrac, 0u, 3, zFrac0, zFrac1);
1333 return __packFloat64(aSign, aExp + 0x380, zFrac0, zFrac1);
1336 /* Adds the 96-bit value formed by concatenating `a0', `a1', and `a2' to the
1337 * 96-bit value formed by concatenating `b0', `b1', and `b2'. Addition is
1338 * modulo 2^96, so any carry out is lost. The result is broken into three
1339 * 32-bit pieces which are stored at the locations pointed to by `z0Ptr',
1340 * `z1Ptr', and `z2Ptr'.
1343 __add96(uint a0, uint a1, uint a2,
1344 uint b0, uint b1, uint b2,
1350 uint carry1 = uint(z2 < a2);
1352 uint carry0 = uint(z1 < a1);
1355 z0 += uint(z1 < carry1);
1362 /* Subtracts the 96-bit value formed by concatenating `b0', `b1', and `b2' from
1363 * the 96-bit value formed by concatenating `a0', `a1', and `a2'. Subtraction
1364 * is modulo 2^96, so any borrow out (carry out) is lost. The result is broken
1365 * into three 32-bit pieces which are stored at the locations pointed to by
1366 * `z0Ptr', `z1Ptr', and `z2Ptr'.
1369 __sub96(uint a0, uint a1, uint a2,
1370 uint b0, uint b1, uint b2,
1376 uint borrow1 = uint(a2 < b2);
1378 uint borrow0 = uint(a1 < b1);
1380 z0 -= uint(z1 < borrow1);
1388 /* Returns an approximation to the 32-bit integer quotient obtained by dividing
1389 * `b' into the 64-bit value formed by concatenating `a0' and `a1'. The
1390 * divisor `b' must be at least 2^31. If q is the exact quotient truncated
1391 * toward zero, the approximation returned lies between q and q + 2 inclusive.
1392 * If the exact quotient q is larger than 32 bits, the maximum positive 32-bit
1393 * unsigned integer is returned.
1396 __estimateDiv64To32(uint a0, uint a1, uint b)
1409 z = (b0<<16 <= a0) ? 0xFFFF0000u : (a0 / b0)<<16;
1410 __mul32To64(b, z, term0, term1);
1411 __sub64(a0, a1, term0, term1, rem0, rem1);
1412 while (int(rem0) < 0) {
1415 __add64(rem0, rem1, b0, b1, rem0, rem1);
1417 rem0 = (rem0<<16) | (rem1>>16);
1418 z |= (b0<<16 <= rem0) ? 0xFFFFu : rem0 / b0;
1423 __sqrtOddAdjustments(int index)
1463 __sqrtEvenAdjustments(int index)
1502 /* Returns an approximation to the square root of the 32-bit significand given
1503 * by `a'. Considered as an integer, `a' must be at least 2^31. If bit 0 of
1504 * `aExp' (the least significant bit) is 1, the integer returned approximates
1505 * 2^31*sqrt(`a'/2^31), where `a' is considered an integer. If bit 0 of `aExp'
1506 * is 0, the integer returned approximates 2^31*sqrt(`a'/2^30). In either
1507 * case, the approximation returned lies strictly within +/-2 of the exact
1511 __estimateSqrt32(int aExp, uint a)
1515 int index = int(a>>27 & 15u);
1516 if ((aExp & 1) != 0) {
1517 z = 0x4000u + (a>>17) - __sqrtOddAdjustments(index);
1518 z = ((a / z)<<14) + (z<<15);
1521 z = 0x8000u + (a>>17) - __sqrtEvenAdjustments(index);
1523 z = (0x20000u <= z) ? 0xFFFF8000u : (z<<15);
1525 return uint(int(a)>>1);
1527 return ((__estimateDiv64To32(a, 0u, z))>>1) + (z>>1);
1530 /* Returns the square root of the double-precision floating-point value `a'.
1531 * The operation is performed according to the IEEE Standard for Floating-Point
1535 __fsqrt64(uint64_t a)
1540 uint doubleZFrac0 = 0u;
1549 uint64_t default_nan = 0xFFFFFFFFFFFFFFFFUL;
1551 uint aFracLo = __extractFloat64FracLo(a);
1552 uint aFracHi = __extractFloat64FracHi(a);
1553 int aExp = __extractFloat64Exp(a);
1554 uint aSign = __extractFloat64Sign(a);
1555 if (aExp == 0x7FF) {
1556 if ((aFracHi | aFracLo) != 0u)
1557 return __propagateFloat64NaN(a, a);
1563 if ((uint(aExp) | aFracHi | aFracLo) == 0u)
1568 if ((aFracHi | aFracLo) == 0u)
1569 return __packFloat64(0u, 0, 0u, 0u);
1570 __normalizeFloat64Subnormal(aFracHi, aFracLo, aExp, aFracHi, aFracLo);
1572 int zExp = ((aExp - 0x3FF)>>1) + 0x3FE;
1573 aFracHi |= 0x00100000u;
1574 __shortShift64Left(aFracHi, aFracLo, 11, term0, term1);
1575 zFrac0 = (__estimateSqrt32(aExp, term0)>>1) + 1u;
1577 zFrac0 = 0x7FFFFFFFu;
1578 doubleZFrac0 = zFrac0 + zFrac0;
1579 __shortShift64Left(aFracHi, aFracLo, 9 - (aExp & 1), aFracHi, aFracLo);
1580 __mul32To64(zFrac0, zFrac0, term0, term1);
1581 __sub64(aFracHi, aFracLo, term0, term1, rem0, rem1);
1582 while (int(rem0) < 0) {
1585 __add64(rem0, rem1, 0u, doubleZFrac0 | 1u, rem0, rem1);
1587 zFrac1 = __estimateDiv64To32(rem1, 0u, doubleZFrac0);
1588 if ((zFrac1 & 0x1FFu) <= 5u) {
1591 __mul32To64(doubleZFrac0, zFrac1, term1, term2);
1592 __sub64(rem1, 0u, term1, term2, rem1, rem2);
1593 __mul32To64(zFrac1, zFrac1, term2, term3);
1594 __sub96(rem1, rem2, 0u, 0u, term2, term3, rem1, rem2, rem3);
1595 while (int(rem1) < 0) {
1597 __shortShift64Left(0u, zFrac1, 1, term2, term3);
1599 term2 |= doubleZFrac0;
1600 __add96(rem1, rem2, rem3, 0u, term2, term3, rem1, rem2, rem3);
1602 zFrac1 |= uint((rem1 | rem2 | rem3) != 0u);
1604 __shift64ExtraRightJamming(zFrac0, zFrac1, 0u, 10, zFrac0, zFrac1, zFrac2);
1605 return __roundAndPackFloat64(0u, zExp, zFrac0, zFrac1, zFrac2);
1609 __ftrunc64(uint64_t __a)
1611 uvec2 a = unpackUint2x32(__a);
1612 int aExp = __extractFloat64Exp(__a);
1616 int unbiasedExp = aExp - 1023;
1617 int fracBits = 52 - unbiasedExp;
1618 uint maskLo = mix(~0u << fracBits, 0u, fracBits >= 32);
1619 uint maskHi = mix(~0u << (fracBits - 32), ~0u, fracBits < 33);
1623 zLo = mix(zLo, 0u, unbiasedExp < 0);
1624 zHi = mix(zHi, 0u, unbiasedExp < 0);
1625 zLo = mix(zLo, a.x, unbiasedExp > 52);
1626 zHi = mix(zHi, a.y, unbiasedExp > 52);
1627 return packUint2x32(uvec2(zLo, zHi));
1631 __ffloor64(uint64_t a)
1633 bool is_positive = __fge64(a, 0ul);
1634 uint64_t tr = __ftrunc64(a);
1636 if (is_positive || __feq64(tr, a)) {
1639 return __fadd64(tr, 0xbff0000000000000ul /* -1.0 */);
1644 __fround64(uint64_t __a)
1646 uvec2 a = unpackUint2x32(__a);
1647 int unbiasedExp = __extractFloat64Exp(__a) - 1023;
1651 if (unbiasedExp < 20) {
1652 if (unbiasedExp < 0) {
1654 if (unbiasedExp == -1 && aLo != 0u)
1655 aHi |= (1023u << 20);
1658 uint maskExp = 0x000FFFFFu >> unbiasedExp;
1659 /* a is an integral value */
1660 if (((aHi & maskExp) == 0u) && (aLo == 0u))
1663 aHi += 0x00080000u >> unbiasedExp;
1667 } else if (unbiasedExp > 51 || unbiasedExp == 1024) {
1670 uint maskExp = 0xFFFFFFFFu >> (unbiasedExp - 20);
1671 if ((aLo & maskExp) == 0u)
1673 uint tmp = aLo + (1u << (51 - unbiasedExp));
1682 return packUint2x32(a);
1686 __fmin64(uint64_t a, uint64_t b)
1688 if (__is_nan(a)) return b;
1689 if (__is_nan(b)) return a;
1691 if (__flt64_nonnan(a, b)) return a;
1696 __fmax64(uint64_t a, uint64_t b)
1698 if (__is_nan(a)) return b;
1699 if (__is_nan(b)) return a;
1701 if (__flt64_nonnan(a, b)) return b;
1706 __ffract64(uint64_t a)
1708 return __fadd64(a, __fneg64(__ffloor64(a)));