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 /* Relax propagation of NaN. Binary operations with a NaN source will still
63 * produce a NaN result, but it won't follow strict IEEE rules.
65 #define RELAXED_NAN_PROPAGATION
67 /* Absolute value of a Float64 :
71 __fabs64(uint64_t __a)
73 uvec2 a = unpackUint2x32(__a);
75 return packUint2x32(a);
78 /* Returns 1 if the double-precision floating-point value `a' is a NaN;
79 * otherwise returns 0.
82 __is_nan(uint64_t __a)
84 uvec2 a = unpackUint2x32(__a);
85 return (0xFFE00000u <= (a.y<<1)) &&
86 ((a.x != 0u) || ((a.y & 0x000FFFFFu) != 0u));
89 /* Negate value of a Float64 :
93 __fneg64(uint64_t __a)
95 uvec2 a = unpackUint2x32(__a);
97 return packUint2x32(a);
101 __fsign64(uint64_t __a)
103 uvec2 a = unpackUint2x32(__a);
106 retval.y = mix((a.y & 0x80000000u) | 0x3FF00000u, 0u, (a.y << 1 | a.x) == 0u);
107 return packUint2x32(retval);
110 /* Returns the fraction bits of the double-precision floating-point value `a'.*/
112 __extractFloat64FracLo(uint64_t a)
114 return unpackUint2x32(a).x;
118 __extractFloat64FracHi(uint64_t a)
120 return unpackUint2x32(a).y & 0x000FFFFFu;
123 /* Returns the exponent bits of the double-precision floating-point value `a'.*/
125 __extractFloat64Exp(uint64_t __a)
127 uvec2 a = unpackUint2x32(__a);
128 return int((a.y>>20) & 0x7FFu);
132 __feq64_nonnan(uint64_t __a, uint64_t __b)
134 uvec2 a = unpackUint2x32(__a);
135 uvec2 b = unpackUint2x32(__b);
136 return (a.x == b.x) &&
137 ((a.y == b.y) || ((a.x == 0u) && (((a.y | b.y)<<1) == 0u)));
140 /* Returns true if the double-precision floating-point value `a' is equal to the
141 * corresponding value `b', and false otherwise. The comparison is performed
142 * according to the IEEE Standard for Floating-Point Arithmetic.
145 __feq64(uint64_t a, uint64_t b)
147 if (__is_nan(a) || __is_nan(b))
150 return __feq64_nonnan(a, b);
153 /* Returns true if the double-precision floating-point value `a' is not equal
154 * to the corresponding value `b', and false otherwise. The comparison is
155 * performed according to the IEEE Standard for Floating-Point Arithmetic.
158 __fne64(uint64_t a, uint64_t b)
160 if (__is_nan(a) || __is_nan(b))
163 return !__feq64_nonnan(a, b);
166 /* Returns the sign bit of the double-precision floating-point value `a'.*/
168 __extractFloat64Sign(uint64_t a)
170 return unpackUint2x32(a).y & 0x80000000u;
173 /* Returns true if the signed 64-bit value formed by concatenating `a0' and
174 * `a1' is less than the signed 64-bit value formed by concatenating `b0' and
175 * `b1'. Otherwise, returns false.
178 ilt64(uint a0, uint a1, uint b0, uint b1)
180 return (int(a0) < int(b0)) || ((a0 == b0) && (a1 < b1));
184 __flt64_nonnan(uint64_t __a, uint64_t __b)
186 uvec2 a = unpackUint2x32(__a);
187 uvec2 b = unpackUint2x32(__b);
189 /* IEEE 754 floating point numbers are specifically designed so that, with
190 * two exceptions, values can be compared by bit-casting to signed integers
191 * with the same number of bits.
193 * From https://en.wikipedia.org/wiki/IEEE_754-1985#Comparing_floating-point_numbers:
195 * When comparing as 2's-complement integers: If the sign bits differ,
196 * the negative number precedes the positive number, so 2's complement
197 * gives the correct result (except that negative zero and positive zero
198 * should be considered equal). If both values are positive, the 2's
199 * complement comparison again gives the correct result. Otherwise (two
200 * negative numbers), the correct FP ordering is the opposite of the 2's
201 * complement ordering.
203 * The logic implied by the above quotation is:
205 * !both_are_zero(a, b) && (both_negative(a, b) ? a > b : a < b)
207 * This is equivalent to
209 * fne(a, b) && (both_negative(a, b) ? a >= b : a < b)
211 * fne(a, b) && (both_negative(a, b) ? !(a < b) : a < b)
213 * fne(a, b) && ((both_negative(a, b) && !(a < b)) ||
214 * (!both_negative(a, b) && (a < b)))
216 * (A!|B)&(A|!B) is (A xor B) which is implemented here using !=.
218 * fne(a, b) && (both_negative(a, b) != (a < b))
220 bool lt = ilt64(a.y, a.x, b.y, b.x);
221 bool both_negative = (a.y & b.y & 0x80000000u) != 0;
223 return !__feq64_nonnan(__a, __b) && (lt != both_negative);
226 /* Returns true if the double-precision floating-point value `a' is less than
227 * the corresponding value `b', and false otherwise. The comparison is performed
228 * according to the IEEE Standard for Floating-Point Arithmetic.
231 __flt64(uint64_t a, uint64_t b)
233 /* This weird layout matters. Doing the "obvious" thing results in extra
234 * flow control being inserted to implement the short-circuit evaluation
235 * rules. Flow control is bad!
237 bool x = !__is_nan(a);
238 bool y = !__is_nan(b);
239 bool z = __flt64_nonnan(a, b);
241 return (x && y && z);
244 /* Returns true if the double-precision floating-point value `a' is greater
245 * than or equal to * the corresponding value `b', and false otherwise. The
246 * comparison is performed * according to the IEEE Standard for Floating-Point
250 __fge64(uint64_t a, uint64_t b)
252 /* This weird layout matters. Doing the "obvious" thing results in extra
253 * flow control being inserted to implement the short-circuit evaluation
254 * rules. Flow control is bad!
256 bool x = !__is_nan(a);
257 bool y = !__is_nan(b);
258 bool z = !__flt64_nonnan(a, b);
260 return (x && y && z);
264 __fsat64(uint64_t __a)
266 uvec2 a = unpackUint2x32(__a);
268 /* fsat(NaN) should be zero. */
269 if (__is_nan(__a) || int(a.y) < 0)
272 /* IEEE 754 floating point numbers are specifically designed so that, with
273 * two exceptions, values can be compared by bit-casting to signed integers
274 * with the same number of bits.
276 * From https://en.wikipedia.org/wiki/IEEE_754-1985#Comparing_floating-point_numbers:
278 * When comparing as 2's-complement integers: If the sign bits differ,
279 * the negative number precedes the positive number, so 2's complement
280 * gives the correct result (except that negative zero and positive zero
281 * should be considered equal). If both values are positive, the 2's
282 * complement comparison again gives the correct result. Otherwise (two
283 * negative numbers), the correct FP ordering is the opposite of the 2's
284 * complement ordering.
286 * We know that both values are not negative, and we know that at least one
287 * value is not zero. Therefore, we can just use the 2's complement
288 * comparison ordering.
290 if (ilt64(0x3FF00000, 0x00000000, a.y, a.x))
291 return 0x3FF0000000000000ul;
296 /* Adds the 64-bit value formed by concatenating `a0' and `a1' to the 64-bit
297 * value formed by concatenating `b0' and `b1'. Addition is modulo 2^64, so
298 * any carry out is lost. The result is broken into two 32-bit pieces which
299 * are stored at the locations pointed to by `z0Ptr' and `z1Ptr'.
302 __add64(uint a0, uint a1, uint b0, uint b1,
308 z0Ptr = a0 + b0 + uint(z1 < a1);
312 /* Subtracts the 64-bit value formed by concatenating `b0' and `b1' from the
313 * 64-bit value formed by concatenating `a0' and `a1'. Subtraction is modulo
314 * 2^64, so any borrow out (carry out) is lost. The result is broken into two
315 * 32-bit pieces which are stored at the locations pointed to by `z0Ptr' and
319 __sub64(uint a0, uint a1, uint b0, uint b1,
324 z0Ptr = a0 - b0 - uint(a1 < b1);
327 /* Shifts the 64-bit value formed by concatenating `a0' and `a1' right by the
328 * number of bits given in `count'. If any nonzero bits are shifted off, they
329 * are "jammed" into the least significant bit of the result by setting the
330 * least significant bit to 1. The value of `count' can be arbitrarily large;
331 * in particular, if `count' is greater than 64, the result will be either 0
332 * or 1, depending on whether the concatenation of `a0' and `a1' is zero or
333 * nonzero. The result is broken into two 32-bit pieces which are stored at
334 * the locations pointed to by `z0Ptr' and `z1Ptr'.
337 __shift64RightJamming(uint a0,
345 int negCount = (-count) & 31;
347 z0 = mix(0u, a0, count == 0);
348 z0 = mix(z0, (a0 >> count), count < 32);
350 z1 = uint((a0 | a1) != 0u); /* count >= 64 */
351 uint z1_lt64 = (a0>>(count & 31)) | uint(((a0<<negCount) | a1) != 0u);
352 z1 = mix(z1, z1_lt64, count < 64);
353 z1 = mix(z1, (a0 | uint(a1 != 0u)), count == 32);
354 uint z1_lt32 = (a0<<negCount) | (a1>>count) | uint ((a1<<negCount) != 0u);
355 z1 = mix(z1, z1_lt32, count < 32);
356 z1 = mix(z1, a1, count == 0);
361 /* Shifts the 96-bit value formed by concatenating `a0', `a1', and `a2' right
362 * by 32 _plus_ the number of bits given in `count'. The shifted result is
363 * at most 64 nonzero bits; these are broken into two 32-bit pieces which are
364 * stored at the locations pointed to by `z0Ptr' and `z1Ptr'. The bits shifted
365 * off form a third 32-bit result as follows: The _last_ bit shifted off is
366 * the most-significant bit of the extra result, and the other 31 bits of the
367 * extra result are all zero if and only if _all_but_the_last_ bits shifted off
368 * were all zero. This extra result is stored in the location pointed to by
369 * `z2Ptr'. The value of `count' can be arbitrarily large.
370 * (This routine makes more sense if `a0', `a1', and `a2' are considered
371 * to form a fixed-point value with binary point between `a1' and `a2'. This
372 * fixed-point value is shifted right by the number of bits given in `count',
373 * and the integer part of the result is returned at the locations pointed to
374 * by `z0Ptr' and `z1Ptr'. The fractional part of the result may be slightly
375 * corrupted as described above, and is returned at the location pointed to by
379 __shift64ExtraRightJamming(uint a0, uint a1, uint a2,
388 int negCount = (-count) & 31;
390 z2 = mix(uint(a0 != 0u), a0, count == 64);
391 z2 = mix(z2, a0 << negCount, count < 64);
392 z2 = mix(z2, a1 << negCount, count < 32);
394 z1 = mix(0u, (a0 >> (count & 31)), count < 64);
395 z1 = mix(z1, (a0<<negCount) | (a1>>count), count < 32);
397 a2 = mix(a2 | a1, a2, count < 32);
398 z0 = mix(z0, a0 >> count, count < 32);
399 z2 |= uint(a2 != 0u);
401 z0 = mix(z0, 0u, (count == 32));
402 z1 = mix(z1, a0, (count == 32));
403 z2 = mix(z2, a1, (count == 32));
404 z0 = mix(z0, a0, (count == 0));
405 z1 = mix(z1, a1, (count == 0));
406 z2 = mix(z2, a2, (count == 0));
412 /* Shifts the 64-bit value formed by concatenating `a0' and `a1' left by the
413 * number of bits given in `count'. Any bits shifted off are lost. The value
414 * of `count' must be less than 32. The result is broken into two 32-bit
415 * pieces which are stored at the locations pointed to by `z0Ptr' and `z1Ptr'.
418 __shortShift64Left(uint a0, uint a1,
424 z0Ptr = mix((a0 << count | (a1 >> ((-count) & 31))), a0, count == 0);
427 /* Packs the sign `zSign', the exponent `zExp', and the significand formed by
428 * the concatenation of `zFrac0' and `zFrac1' into a double-precision floating-
429 * point value, returning the result. After being shifted into the proper
430 * positions, the three fields `zSign', `zExp', and `zFrac0' are simply added
431 * together to form the most significant 32 bits of the result. This means
432 * that any integer portion of `zFrac0' will be added into the exponent. Since
433 * a properly normalized significand will have an integer portion equal to 1,
434 * the `zExp' input should be 1 less than the desired result exponent whenever
435 * `zFrac0' and `zFrac1' concatenated form a complete, normalized significand.
438 __packFloat64(uint zSign, int zExp, uint zFrac0, uint zFrac1)
442 z.y = zSign + (uint(zExp) << 20) + zFrac0;
444 return packUint2x32(z);
447 /* Takes an abstract floating-point value having sign `zSign', exponent `zExp',
448 * and extended significand formed by the concatenation of `zFrac0', `zFrac1',
449 * and `zFrac2', and returns the proper double-precision floating-point value
450 * corresponding to the abstract input. Ordinarily, the abstract value is
451 * simply rounded and packed into the double-precision format, with the inexact
452 * exception raised if the abstract input cannot be represented exactly.
453 * However, if the abstract value is too large, the overflow and inexact
454 * exceptions are raised and an infinity or maximal finite value is returned.
455 * If the abstract value is too small, the input value is rounded to a
456 * subnormal number, and the underflow and inexact exceptions are raised if the
457 * abstract input cannot be represented exactly as a subnormal double-precision
458 * floating-point number.
459 * The input significand must be normalized or smaller. If the input
460 * significand is not normalized, `zExp' must be 0; in that case, the result
461 * returned is a subnormal number, and it must not require rounding. In the
462 * usual case that the input significand is normalized, `zExp' must be 1 less
463 * than the "true" floating-point exponent. The handling of underflow and
464 * overflow follows the IEEE Standard for Floating-Point Arithmetic.
467 __roundAndPackFloat64(uint zSign,
473 bool roundNearestEven;
476 roundNearestEven = FLOAT_ROUNDING_MODE == FLOAT_ROUND_NEAREST_EVEN;
477 increment = int(zFrac2) < 0;
478 if (!roundNearestEven) {
479 if (FLOAT_ROUNDING_MODE == FLOAT_ROUND_TO_ZERO) {
483 increment = (FLOAT_ROUNDING_MODE == FLOAT_ROUND_DOWN) &&
486 increment = (FLOAT_ROUNDING_MODE == FLOAT_ROUND_UP) &&
492 if ((0x7FD < zExp) ||
494 (0x001FFFFFu == zFrac0 && 0xFFFFFFFFu == zFrac1) &&
496 if ((FLOAT_ROUNDING_MODE == FLOAT_ROUND_TO_ZERO) ||
497 ((zSign != 0u) && (FLOAT_ROUNDING_MODE == FLOAT_ROUND_UP)) ||
498 ((zSign == 0u) && (FLOAT_ROUNDING_MODE == FLOAT_ROUND_DOWN))) {
499 return __packFloat64(zSign, 0x7FE, 0x000FFFFFu, 0xFFFFFFFFu);
501 return __packFloat64(zSign, 0x7FF, 0u, 0u);
504 __shift64ExtraRightJamming(
505 zFrac0, zFrac1, zFrac2, -zExp, zFrac0, zFrac1, zFrac2);
507 if (roundNearestEven) {
508 increment = zFrac2 < 0u;
511 increment = (FLOAT_ROUNDING_MODE == FLOAT_ROUND_DOWN) &&
514 increment = (FLOAT_ROUNDING_MODE == FLOAT_ROUND_UP) &&
521 __add64(zFrac0, zFrac1, 0u, 1u, zFrac0, zFrac1);
522 zFrac1 &= ~((zFrac2 + uint(zFrac2 == 0u)) & uint(roundNearestEven));
524 zExp = mix(zExp, 0, (zFrac0 | zFrac1) == 0u);
526 return __packFloat64(zSign, zExp, zFrac0, zFrac1);
530 __roundAndPackUInt64(uint zSign, uint zFrac0, uint zFrac1, uint zFrac2)
532 bool roundNearestEven;
534 uint64_t default_nan = 0xFFFFFFFFFFFFFFFFUL;
536 roundNearestEven = FLOAT_ROUNDING_MODE == FLOAT_ROUND_NEAREST_EVEN;
538 if (zFrac2 >= 0x80000000u)
541 if (!roundNearestEven) {
543 if ((FLOAT_ROUNDING_MODE == FLOAT_ROUND_DOWN) && (zFrac2 != 0u)) {
547 increment = (FLOAT_ROUNDING_MODE == FLOAT_ROUND_UP) &&
553 __add64(zFrac0, zFrac1, 0u, 1u, zFrac0, zFrac1);
554 if ((zFrac0 | zFrac1) != 0u)
555 zFrac1 &= ~(1u) + uint(zFrac2 == 0u) & uint(roundNearestEven);
557 return mix(packUint2x32(uvec2(zFrac1, zFrac0)), default_nan,
558 (zSign != 0u && (zFrac0 | zFrac1) != 0u));
562 __roundAndPackInt64(uint zSign, uint zFrac0, uint zFrac1, uint zFrac2)
564 bool roundNearestEven;
566 int64_t default_NegNaN = -0x7FFFFFFFFFFFFFFEL;
567 int64_t default_PosNaN = 0xFFFFFFFFFFFFFFFFL;
569 roundNearestEven = FLOAT_ROUNDING_MODE == FLOAT_ROUND_NEAREST_EVEN;
571 if (zFrac2 >= 0x80000000u)
574 if (!roundNearestEven) {
576 increment = ((FLOAT_ROUNDING_MODE == FLOAT_ROUND_DOWN) &&
579 increment = (FLOAT_ROUNDING_MODE == FLOAT_ROUND_UP) &&
585 __add64(zFrac0, zFrac1, 0u, 1u, zFrac0, zFrac1);
586 if ((zFrac0 | zFrac1) != 0u)
587 zFrac1 &= ~(1u) + uint(zFrac2 == 0u) & uint(roundNearestEven);
590 int64_t absZ = mix(int64_t(packUint2x32(uvec2(zFrac1, zFrac0))),
591 -int64_t(packUint2x32(uvec2(zFrac1, zFrac0))),
593 int64_t nan = mix(default_PosNaN, default_NegNaN, zSign != 0u);
594 return mix(absZ, nan, ((zSign != 0u) != (absZ < 0)) && bool(absZ));
597 /* Returns the number of leading 0 bits before the most-significant 1 bit of
598 * `a'. If `a' is zero, 32 is returned.
601 __countLeadingZeros32(uint a)
603 return 31 - findMSB(a);
606 /* Takes an abstract floating-point value having sign `zSign', exponent `zExp',
607 * and significand formed by the concatenation of `zSig0' and `zSig1', and
608 * returns the proper double-precision floating-point value corresponding
609 * to the abstract input. This routine is just like `__roundAndPackFloat64'
610 * except that the input significand has fewer bits and does not have to be
611 * normalized. In all cases, `zExp' must be 1 less than the "true" floating-
615 __normalizeRoundAndPackFloat64(uint zSign,
629 shiftCount = __countLeadingZeros32(zFrac0) - 11;
630 if (0 <= shiftCount) {
632 __shortShift64Left(zFrac0, zFrac1, shiftCount, zFrac0, zFrac1);
634 __shift64ExtraRightJamming(
635 zFrac0, zFrac1, 0u, -shiftCount, zFrac0, zFrac1, zFrac2);
638 return __roundAndPackFloat64(zSign, zExp, zFrac0, zFrac1, zFrac2);
641 /* Takes two double-precision floating-point values `a' and `b', one of which
642 * is a NaN, and returns the appropriate NaN result.
645 __propagateFloat64NaN(uint64_t __a, uint64_t __b)
647 #if defined RELAXED_NAN_PROPAGATION
648 uvec2 a = unpackUint2x32(__a);
649 uvec2 b = unpackUint2x32(__b);
651 return packUint2x32(uvec2(a.x | b.x, a.y | b.y));
653 bool aIsNaN = __is_nan(__a);
654 bool bIsNaN = __is_nan(__b);
655 uvec2 a = unpackUint2x32(__a);
656 uvec2 b = unpackUint2x32(__b);
660 return packUint2x32(mix(b, mix(a, b, bvec2(bIsNaN, bIsNaN)), bvec2(aIsNaN, aIsNaN)));
664 /* Returns the result of adding the double-precision floating-point values
665 * `a' and `b'. The operation is performed according to the IEEE Standard for
666 * Floating-Point Arithmetic.
669 __fadd64(uint64_t a, uint64_t b)
671 uint aSign = __extractFloat64Sign(a);
672 uint bSign = __extractFloat64Sign(b);
673 uint aFracLo = __extractFloat64FracLo(a);
674 uint aFracHi = __extractFloat64FracHi(a);
675 uint bFracLo = __extractFloat64FracLo(b);
676 uint bFracHi = __extractFloat64FracHi(b);
677 int aExp = __extractFloat64Exp(a);
678 int bExp = __extractFloat64Exp(b);
679 int expDiff = aExp - bExp;
680 if (aSign == bSign) {
685 bool orig_exp_diff_is_zero = (expDiff == 0);
687 if (orig_exp_diff_is_zero) {
689 bool propagate = ((aFracHi | bFracHi) | (aFracLo| bFracLo)) != 0u;
690 return mix(a, __propagateFloat64NaN(a, b), propagate);
692 __add64(aFracHi, aFracLo, bFracHi, bFracLo, zFrac0, zFrac1);
694 return __packFloat64(aSign, 0, zFrac0, zFrac1);
696 zFrac0 |= 0x00200000u;
698 __shift64ExtraRightJamming(
699 zFrac0, zFrac1, zFrac2, 1, zFrac0, zFrac1, zFrac2);
700 } else if (0 < expDiff) {
702 bool propagate = (aFracHi | aFracLo) != 0u;
703 return mix(a, __propagateFloat64NaN(a, b), propagate);
706 expDiff = mix(expDiff, expDiff - 1, bExp == 0);
707 bFracHi = mix(bFracHi | 0x00100000u, bFracHi, bExp == 0);
708 __shift64ExtraRightJamming(
709 bFracHi, bFracLo, 0u, expDiff, bFracHi, bFracLo, zFrac2);
711 } else if (expDiff < 0) {
713 bool propagate = (bFracHi | bFracLo) != 0u;
714 return mix(__packFloat64(aSign, 0x7ff, 0u, 0u), __propagateFloat64NaN(a, b), propagate);
716 expDiff = mix(expDiff, expDiff + 1, aExp == 0);
717 aFracHi = mix(aFracHi | 0x00100000u, aFracHi, aExp == 0);
718 __shift64ExtraRightJamming(
719 aFracHi, aFracLo, 0u, - expDiff, aFracHi, aFracLo, zFrac2);
722 if (!orig_exp_diff_is_zero) {
723 aFracHi |= 0x00100000u;
724 __add64(aFracHi, aFracLo, bFracHi, bFracLo, zFrac0, zFrac1);
726 if (!(zFrac0 < 0x00200000u)) {
727 __shift64ExtraRightJamming(zFrac0, zFrac1, zFrac2, 1, zFrac0, zFrac1, zFrac2);
731 return __roundAndPackFloat64(aSign, zExp, zFrac0, zFrac1, zFrac2);
736 __shortShift64Left(aFracHi, aFracLo, 10, aFracHi, aFracLo);
737 __shortShift64Left(bFracHi, bFracLo, 10, bFracHi, bFracLo);
743 bool propagate = (aFracHi | aFracLo) != 0u;
744 return mix(a, __propagateFloat64NaN(a, b), propagate);
746 expDiff = mix(expDiff, expDiff - 1, bExp == 0);
747 bFracHi = mix(bFracHi | 0x40000000u, bFracHi, bExp == 0);
748 __shift64RightJamming(bFracHi, bFracLo, expDiff, bFracHi, bFracLo);
749 aFracHi |= 0x40000000u;
750 __sub64(aFracHi, aFracLo, bFracHi, bFracLo, zFrac0, zFrac1);
753 return __normalizeRoundAndPackFloat64(aSign, zExp - 10, zFrac0, zFrac1);
760 bool propagate = (bFracHi | bFracLo) != 0u;
761 return mix(__packFloat64(aSign ^ 0x80000000u, 0x7ff, 0u, 0u), __propagateFloat64NaN(a, b), propagate);
763 expDiff = mix(expDiff, expDiff + 1, aExp == 0);
764 aFracHi = mix(aFracHi | 0x40000000u, aFracHi, aExp == 0);
765 __shift64RightJamming(aFracHi, aFracLo, - expDiff, aFracHi, aFracLo);
766 bFracHi |= 0x40000000u;
767 __sub64(bFracHi, bFracLo, aFracHi, aFracLo, zFrac0, zFrac1);
769 aSign ^= 0x80000000u;
771 return __normalizeRoundAndPackFloat64(aSign, zExp - 10, zFrac0, zFrac1);
774 bool propagate = ((aFracHi | bFracHi) | (aFracLo | bFracLo)) != 0u;
775 return mix(0xFFFFFFFFFFFFFFFFUL, __propagateFloat64NaN(a, b), propagate);
777 bExp = mix(bExp, 1, aExp == 0);
778 aExp = mix(aExp, 1, aExp == 0);
782 uint sign_of_difference = 0;
783 if (bFracHi < aFracHi) {
784 __sub64(aFracHi, aFracLo, bFracHi, bFracLo, zFrac0, zFrac1);
786 else if (aFracHi < bFracHi) {
787 __sub64(bFracHi, bFracLo, aFracHi, aFracLo, zFrac0, zFrac1);
788 sign_of_difference = 0x80000000;
790 else if (bFracLo < aFracLo) {
791 __sub64(aFracHi, aFracLo, bFracHi, bFracLo, zFrac0, zFrac1);
793 else if (aFracLo < bFracLo) {
794 __sub64(bFracHi, bFracLo, aFracHi, aFracLo, zFrac0, zFrac1);
795 sign_of_difference = 0x80000000;
797 zExp = mix(bExp, aExp, sign_of_difference == 0u);
798 aSign ^= sign_of_difference;
799 uint64_t retval_0 = __packFloat64(uint(FLOAT_ROUNDING_MODE == FLOAT_ROUND_DOWN) << 31, 0, 0u, 0u);
800 uint64_t retval_1 = __normalizeRoundAndPackFloat64(aSign, zExp - 11, zFrac0, zFrac1);
801 return mix(retval_0, retval_1, zFrac0 != 0u || zFrac1 != 0u);
805 /* Multiplies the 64-bit value formed by concatenating `a0' and `a1' to the
806 * 64-bit value formed by concatenating `b0' and `b1' to obtain a 128-bit
807 * product. The product is broken into four 32-bit pieces which are stored at
808 * the locations pointed to by `z0Ptr', `z1Ptr', `z2Ptr', and `z3Ptr'.
811 __mul64To128(uint a0, uint a1, uint b0, uint b1,
824 umulExtended(a1, b1, z2, z3);
825 umulExtended(a1, b0, z1, more2);
826 __add64(z1, more2, 0u, z2, z1, z2);
827 umulExtended(a0, b0, z0, more1);
828 __add64(z0, more1, 0u, z1, z0, z1);
829 umulExtended(a0, b1, more1, more2);
830 __add64(more1, more2, 0u, z2, more1, z2);
831 __add64(z0, z1, 0u, more1, z0, z1);
838 /* Normalizes the subnormal double-precision floating-point value represented
839 * by the denormalized significand formed by the concatenation of `aFrac0' and
840 * `aFrac1'. The normalized exponent is stored at the location pointed to by
841 * `zExpPtr'. The most significant 21 bits of the normalized significand are
842 * stored at the location pointed to by `zFrac0Ptr', and the least significant
843 * 32 bits of the normalized significand are stored at the location pointed to
847 __normalizeFloat64Subnormal(uint aFrac0, uint aFrac1,
853 uint temp_zfrac0, temp_zfrac1;
854 shiftCount = __countLeadingZeros32(mix(aFrac0, aFrac1, aFrac0 == 0u)) - 11;
855 zExpPtr = mix(1 - shiftCount, -shiftCount - 31, aFrac0 == 0u);
857 temp_zfrac0 = mix(aFrac1<<shiftCount, aFrac1>>(-shiftCount), shiftCount < 0);
858 temp_zfrac1 = mix(0u, aFrac1<<(shiftCount & 31), shiftCount < 0);
860 __shortShift64Left(aFrac0, aFrac1, shiftCount, zFrac0Ptr, zFrac1Ptr);
862 zFrac0Ptr = mix(zFrac0Ptr, temp_zfrac0, aFrac0 == 0);
863 zFrac1Ptr = mix(zFrac1Ptr, temp_zfrac1, aFrac0 == 0);
866 /* Returns the result of multiplying the double-precision floating-point values
867 * `a' and `b'. The operation is performed according to the IEEE Standard for
868 * Floating-Point Arithmetic.
871 __fmul64(uint64_t a, uint64_t b)
879 uint aFracLo = __extractFloat64FracLo(a);
880 uint aFracHi = __extractFloat64FracHi(a);
881 uint bFracLo = __extractFloat64FracLo(b);
882 uint bFracHi = __extractFloat64FracHi(b);
883 int aExp = __extractFloat64Exp(a);
884 uint aSign = __extractFloat64Sign(a);
885 int bExp = __extractFloat64Exp(b);
886 uint bSign = __extractFloat64Sign(b);
887 uint zSign = aSign ^ bSign;
889 if (((aFracHi | aFracLo) != 0u) ||
890 ((bExp == 0x7FF) && ((bFracHi | bFracLo) != 0u))) {
891 return __propagateFloat64NaN(a, b);
893 if ((uint(bExp) | bFracHi | bFracLo) == 0u)
894 return 0xFFFFFFFFFFFFFFFFUL;
895 return __packFloat64(zSign, 0x7FF, 0u, 0u);
898 /* a cannot be NaN, but is b NaN? */
899 if ((bFracHi | bFracLo) != 0u)
900 #if defined RELAXED_NAN_PROPAGATION
903 return __propagateFloat64NaN(a, b);
905 if ((uint(aExp) | aFracHi | aFracLo) == 0u)
906 return 0xFFFFFFFFFFFFFFFFUL;
907 return __packFloat64(zSign, 0x7FF, 0u, 0u);
910 if ((aFracHi | aFracLo) == 0u)
911 return __packFloat64(zSign, 0, 0u, 0u);
912 __normalizeFloat64Subnormal(aFracHi, aFracLo, aExp, aFracHi, aFracLo);
915 if ((bFracHi | bFracLo) == 0u)
916 return __packFloat64(zSign, 0, 0u, 0u);
917 __normalizeFloat64Subnormal(bFracHi, bFracLo, bExp, bFracHi, bFracLo);
919 zExp = aExp + bExp - 0x400;
920 aFracHi |= 0x00100000u;
921 __shortShift64Left(bFracHi, bFracLo, 12, bFracHi, bFracLo);
923 aFracHi, aFracLo, bFracHi, bFracLo, zFrac0, zFrac1, zFrac2, zFrac3);
924 __add64(zFrac0, zFrac1, aFracHi, aFracLo, zFrac0, zFrac1);
925 zFrac2 |= uint(zFrac3 != 0u);
926 if (0x00200000u <= zFrac0) {
927 __shift64ExtraRightJamming(
928 zFrac0, zFrac1, zFrac2, 1, zFrac0, zFrac1, zFrac2);
931 return __roundAndPackFloat64(zSign, zExp, zFrac0, zFrac1, zFrac2);
935 __ffma64(uint64_t a, uint64_t b, uint64_t c)
937 return __fadd64(__fmul64(a, b), c);
940 /* Shifts the 64-bit value formed by concatenating `a0' and `a1' right by the
941 * number of bits given in `count'. Any bits shifted off are lost. The value
942 * of `count' can be arbitrarily large; in particular, if `count' is greater
943 * than 64, the result will be 0. The result is broken into two 32-bit pieces
944 * which are stored at the locations pointed to by `z0Ptr' and `z1Ptr'.
947 __shift64Right(uint a0, uint a1,
954 int negCount = (-count) & 31;
957 z0 = mix(z0, (a0 >> count), count < 32);
958 z0 = mix(z0, a0, count == 0);
960 z1 = mix(0u, (a0 >> (count & 31)), count < 64);
961 z1 = mix(z1, (a0<<negCount) | (a1>>count), count < 32);
962 z1 = mix(z1, a0, count == 0);
968 /* Returns the result of converting the double-precision floating-point value
969 * `a' to the unsigned integer format. The conversion is performed according
970 * to the IEEE Standard for Floating-Point Arithmetic.
973 __fp64_to_uint(uint64_t a)
975 uint aFracLo = __extractFloat64FracLo(a);
976 uint aFracHi = __extractFloat64FracHi(a);
977 int aExp = __extractFloat64Exp(a);
978 uint aSign = __extractFloat64Sign(a);
980 if ((aExp == 0x7FF) && ((aFracHi | aFracLo) != 0u))
983 aFracHi |= mix(0u, 0x00100000u, aExp != 0);
985 int shiftDist = 0x427 - aExp;
987 __shift64RightJamming(aFracHi, aFracLo, shiftDist, aFracHi, aFracLo);
989 if ((aFracHi & 0xFFFFF000u) != 0u)
990 return mix(~0u, 0u, aSign != 0u);
994 __shift64Right(aFracHi, aFracLo, 12, zero, z);
996 uint expt = mix(~0u, 0u, aSign != 0u);
998 return mix(z, expt, (aSign != 0u) && (z != 0u));
1002 __uint_to_fp64(uint a)
1007 int shiftDist = __countLeadingZeros32(a) + 21;
1011 int negCount = (- shiftDist) & 31;
1013 aHigh = mix(0u, a<< shiftDist - 32, shiftDist < 64);
1015 aHigh = mix(aHigh, 0u, shiftDist == 0);
1016 aLow = mix(aLow, a, shiftDist ==0);
1017 aHigh = mix(aHigh, a >> negCount, shiftDist < 32);
1018 aLow = mix(aLow, a << shiftDist, shiftDist < 32);
1020 return __packFloat64(0u, 0x432 - shiftDist, aHigh, aLow);
1024 __uint64_to_fp64(uint64_t a)
1029 uvec2 aFrac = unpackUint2x32(a);
1030 uint aFracLo = __extractFloat64FracLo(a);
1031 uint aFracHi = __extractFloat64FracHi(a);
1033 if ((aFracHi & 0x80000000u) != 0u) {
1034 __shift64RightJamming(aFracHi, aFracLo, 1, aFracHi, aFracLo);
1035 return __roundAndPackFloat64(0, 0x433, aFracHi, aFracLo, 0u);
1037 return __normalizeRoundAndPackFloat64(0, 0x432, aFrac.y, aFrac.x);
1042 __fp64_to_uint64(uint64_t a)
1044 uint aFracLo = __extractFloat64FracLo(a);
1045 uint aFracHi = __extractFloat64FracHi(a);
1046 int aExp = __extractFloat64Exp(a);
1047 uint aSign = __extractFloat64Sign(a);
1049 uint64_t default_nan = 0xFFFFFFFFFFFFFFFFUL;
1051 aFracHi = mix(aFracHi, aFracHi | 0x00100000u, aExp != 0);
1052 int shiftCount = 0x433 - aExp;
1054 if ( shiftCount <= 0 ) {
1055 if (shiftCount < -11 && aExp == 0x7FF) {
1056 if ((aFracHi | aFracLo) != 0u)
1057 return __propagateFloat64NaN(a, a);
1058 return mix(default_nan, a, aSign == 0u);
1060 __shortShift64Left(aFracHi, aFracLo, -shiftCount, aFracHi, aFracLo);
1062 __shift64ExtraRightJamming(aFracHi, aFracLo, zFrac2, shiftCount,
1063 aFracHi, aFracLo, zFrac2);
1065 return __roundAndPackUInt64(aSign, aFracHi, aFracLo, zFrac2);
1069 __fp64_to_int64(uint64_t a)
1072 uint aFracLo = __extractFloat64FracLo(a);
1073 uint aFracHi = __extractFloat64FracHi(a);
1074 int aExp = __extractFloat64Exp(a);
1075 uint aSign = __extractFloat64Sign(a);
1076 int64_t default_NegNaN = -0x7FFFFFFFFFFFFFFEL;
1077 int64_t default_PosNaN = 0xFFFFFFFFFFFFFFFFL;
1079 aFracHi = mix(aFracHi, aFracHi | 0x00100000u, aExp != 0);
1080 int shiftCount = 0x433 - aExp;
1082 if (shiftCount <= 0) {
1083 if (shiftCount < -11 && aExp == 0x7FF) {
1084 if ((aFracHi | aFracLo) != 0u)
1085 return default_NegNaN;
1086 return mix(default_NegNaN, default_PosNaN, aSign == 0u);
1088 __shortShift64Left(aFracHi, aFracLo, -shiftCount, aFracHi, aFracLo);
1090 __shift64ExtraRightJamming(aFracHi, aFracLo, zFrac2, shiftCount,
1091 aFracHi, aFracLo, zFrac2);
1094 return __roundAndPackInt64(aSign, aFracHi, aFracLo, zFrac2);
1098 __fp32_to_uint64(float f)
1100 uint a = floatBitsToUint(f);
1101 uint aFrac = a & 0x007FFFFFu;
1102 int aExp = int((a>>23) & 0xFFu);
1103 uint aSign = a & 0x80000000u;
1107 uint64_t default_nan = 0xFFFFFFFFFFFFFFFFUL;
1108 int shiftCount = 0xBE - aExp;
1110 if (shiftCount <0) {
1115 aFrac = mix(aFrac, aFrac | 0x00800000u, aExp != 0);
1116 __shortShift64Left(aFrac, 0, 40, zFrac0, zFrac1);
1118 if (shiftCount != 0) {
1119 __shift64ExtraRightJamming(zFrac0, zFrac1, zFrac2, shiftCount,
1120 zFrac0, zFrac1, zFrac2);
1123 return __roundAndPackUInt64(aSign, zFrac0, zFrac1, zFrac2);
1127 __fp32_to_int64(float f)
1129 uint a = floatBitsToUint(f);
1130 uint aFrac = a & 0x007FFFFFu;
1131 int aExp = int((a>>23) & 0xFFu);
1132 uint aSign = a & 0x80000000u;
1136 int64_t default_NegNaN = -0x7FFFFFFFFFFFFFFEL;
1137 int64_t default_PosNaN = 0xFFFFFFFFFFFFFFFFL;
1138 int shiftCount = 0xBE - aExp;
1140 if (shiftCount <0) {
1141 if (aExp == 0xFF && aFrac != 0u)
1142 return default_NegNaN;
1143 return mix(default_NegNaN, default_PosNaN, aSign == 0u);
1146 aFrac = mix(aFrac, aFrac | 0x00800000u, aExp != 0);
1147 __shortShift64Left(aFrac, 0, 40, zFrac0, zFrac1);
1149 if (shiftCount != 0) {
1150 __shift64ExtraRightJamming(zFrac0, zFrac1, zFrac2, shiftCount,
1151 zFrac0, zFrac1, zFrac2);
1154 return __roundAndPackInt64(aSign, zFrac0, zFrac1, zFrac2);
1158 __int64_to_fp64(int64_t a)
1163 uint64_t absA = mix(uint64_t(a), uint64_t(-a), a < 0);
1164 uint aFracHi = __extractFloat64FracHi(absA);
1165 uvec2 aFrac = unpackUint2x32(absA);
1166 uint zSign = uint(unpackInt2x32(a).y) & 0x80000000u;
1168 if ((aFracHi & 0x80000000u) != 0u) {
1169 return mix(0ul, __packFloat64(0x80000000u, 0x434, 0u, 0u), a < 0);
1172 return __normalizeRoundAndPackFloat64(zSign, 0x432, aFrac.y, aFrac.x);
1175 /* Returns the result of converting the double-precision floating-point value
1176 * `a' to the 32-bit two's complement integer format. The conversion is
1177 * performed according to the IEEE Standard for Floating-Point Arithmetic---
1178 * which means in particular that the conversion is rounded according to the
1179 * current rounding mode. If `a' is a NaN, the largest positive integer is
1180 * returned. Otherwise, if the conversion overflows, the largest integer with
1181 * the same sign as `a' is returned.
1184 __fp64_to_int(uint64_t a)
1186 uint aFracLo = __extractFloat64FracLo(a);
1187 uint aFracHi = __extractFloat64FracHi(a);
1188 int aExp = __extractFloat64Exp(a);
1189 uint aSign = __extractFloat64Sign(a);
1192 uint aFracExtra = 0u;
1193 int shiftCount = aExp - 0x413;
1195 if (0 <= shiftCount) {
1197 if ((aExp == 0x7FF) && bool(aFracHi | aFracLo))
1199 return mix(0x7FFFFFFF, 0x80000000, aSign != 0u);
1201 __shortShift64Left(aFracHi | 0x00100000u, aFracLo, shiftCount, absZ, aFracExtra);
1206 aFracHi |= 0x00100000u;
1207 aFracExtra = ( aFracHi << (shiftCount & 31)) | aFracLo;
1208 absZ = aFracHi >> (- shiftCount);
1211 int z = mix(int(absZ), -int(absZ), aSign != 0u);
1212 int nan = mix(0x7FFFFFFF, 0x80000000, aSign != 0u);
1213 return mix(z, nan, ((aSign != 0u) != (z < 0)) && bool(z));
1216 /* Returns the result of converting the 32-bit two's complement integer `a'
1217 * to the double-precision floating-point format. The conversion is performed
1218 * according to the IEEE Standard for Floating-Point Arithmetic.
1221 __int_to_fp64(int a)
1226 return __packFloat64(0u, 0, 0u, 0u);
1227 uint zSign = uint(a) & 0x80000000u;
1228 uint absA = mix(uint(a), uint(-a), a < 0);
1229 int shiftCount = __countLeadingZeros32(absA) - 11;
1230 if (0 <= shiftCount) {
1231 zFrac0 = absA << shiftCount;
1234 __shift64Right(absA, 0u, -shiftCount, zFrac0, zFrac1);
1236 return __packFloat64(zSign, 0x412 - shiftCount, zFrac0, zFrac1);
1240 __fp64_to_bool(uint64_t a)
1242 return !__feq64_nonnan(__fabs64(a), 0ul);
1246 __bool_to_fp64(bool a)
1248 return packUint2x32(uvec2(0x00000000u, uint(-int(a) & 0x3ff00000)));
1251 /* Packs the sign `zSign', exponent `zExp', and significand `zFrac' into a
1252 * single-precision floating-point value, returning the result. After being
1253 * shifted into the proper positions, the three fields are simply added
1254 * together to form the result. This means that any integer portion of `zSig'
1255 * will be added into the exponent. Since a properly normalized significand
1256 * will have an integer portion equal to 1, the `zExp' input should be 1 less
1257 * than the desired result exponent whenever `zFrac' is a complete, normalized
1261 __packFloat32(uint zSign, int zExp, uint zFrac)
1263 return uintBitsToFloat(zSign + (uint(zExp)<<23) + zFrac);
1266 /* Takes an abstract floating-point value having sign `zSign', exponent `zExp',
1267 * and significand `zFrac', and returns the proper single-precision floating-
1268 * point value corresponding to the abstract input. Ordinarily, the abstract
1269 * value is simply rounded and packed into the single-precision format, with
1270 * the inexact exception raised if the abstract input cannot be represented
1271 * exactly. However, if the abstract value is too large, the overflow and
1272 * inexact exceptions are raised and an infinity or maximal finite value is
1273 * returned. If the abstract value is too small, the input value is rounded to
1274 * a subnormal number, and the underflow and inexact exceptions are raised if
1275 * the abstract input cannot be represented exactly as a subnormal single-
1276 * precision floating-point number.
1277 * The input significand `zFrac' has its binary point between bits 30
1278 * and 29, which is 7 bits to the left of the usual location. This shifted
1279 * significand must be normalized or smaller. If `zFrac' is not normalized,
1280 * `zExp' must be 0; in that case, the result returned is a subnormal number,
1281 * and it must not require rounding. In the usual case that `zFrac' is
1282 * normalized, `zExp' must be 1 less than the "true" floating-point exponent.
1283 * The handling of underflow and overflow follows the IEEE Standard for
1284 * Floating-Point Arithmetic.
1287 __roundAndPackFloat32(uint zSign, int zExp, uint zFrac)
1289 bool roundNearestEven;
1293 roundNearestEven = FLOAT_ROUNDING_MODE == FLOAT_ROUND_NEAREST_EVEN;
1294 roundIncrement = 0x40;
1295 if (!roundNearestEven) {
1296 if (FLOAT_ROUNDING_MODE == FLOAT_ROUND_TO_ZERO) {
1299 roundIncrement = 0x7F;
1301 if (FLOAT_ROUNDING_MODE == FLOAT_ROUND_UP)
1304 if (FLOAT_ROUNDING_MODE == FLOAT_ROUND_DOWN)
1309 roundBits = int(zFrac & 0x7Fu);
1310 if (0xFDu <= uint(zExp)) {
1311 if ((0xFD < zExp) || ((zExp == 0xFD) && (int(zFrac) + roundIncrement) < 0))
1312 return __packFloat32(zSign, 0xFF, 0u) - float(roundIncrement == 0);
1314 bool zexp_lt0 = zExp < 0;
1315 uint zFrac_lt0 = mix(uint(zFrac != 0u), (zFrac>>count) | uint((zFrac<<((-count) & 31)) != 0u), (-zExp) < 32);
1316 zFrac = mix(zFrac, zFrac_lt0, zexp_lt0);
1317 roundBits = mix(roundBits, int(zFrac) & 0x7f, zexp_lt0);
1318 zExp = mix(zExp, 0, zexp_lt0);
1320 zFrac = (zFrac + uint(roundIncrement))>>7;
1321 zFrac &= ~uint(((roundBits ^ 0x40) == 0) && roundNearestEven);
1323 return __packFloat32(zSign, mix(zExp, 0, zFrac == 0u), zFrac);
1326 /* Returns the result of converting the double-precision floating-point value
1327 * `a' to the single-precision floating-point format. The conversion is
1328 * performed according to the IEEE Standard for Floating-Point Arithmetic.
1331 __fp64_to_fp32(uint64_t __a)
1333 uvec2 a = unpackUint2x32(__a);
1337 uint aFracLo = __extractFloat64FracLo(__a);
1338 uint aFracHi = __extractFloat64FracHi(__a);
1339 int aExp = __extractFloat64Exp(__a);
1340 uint aSign = __extractFloat64Sign(__a);
1341 if (aExp == 0x7FF) {
1342 __shortShift64Left(a.y, a.x, 12, a.y, a.x);
1343 float rval = uintBitsToFloat(aSign | 0x7FC00000u | (a.y>>9));
1344 rval = mix(__packFloat32(aSign, 0xFF, 0u), rval, (aFracHi | aFracLo) != 0u);
1347 __shift64RightJamming(aFracHi, aFracLo, 22, allZero, zFrac);
1348 zFrac = mix(zFrac, zFrac | 0x40000000u, aExp != 0);
1349 return __roundAndPackFloat32(aSign, aExp - 0x381, zFrac);
1353 __uint64_to_fp32(uint64_t __a)
1355 uvec2 aFrac = unpackUint2x32(__a);
1356 int shiftCount = mix(__countLeadingZeros32(aFrac.y) - 33,
1357 __countLeadingZeros32(aFrac.x) - 1,
1360 if (0 <= shiftCount)
1361 __shortShift64Left(aFrac.y, aFrac.x, shiftCount, aFrac.y, aFrac.x);
1363 __shift64RightJamming(aFrac.y, aFrac.x, -shiftCount, aFrac.y, aFrac.x);
1365 return __roundAndPackFloat32(0u, 0x9C - shiftCount, aFrac.x);
1369 __int64_to_fp32(int64_t __a)
1371 uint aSign = uint(unpackInt2x32(__a).y) & 0x80000000u;
1372 uint64_t absA = mix(uint64_t(__a), uint64_t(-__a), __a < 0);
1373 uvec2 aFrac = unpackUint2x32(absA);
1374 int shiftCount = mix(__countLeadingZeros32(aFrac.y) - 33,
1375 __countLeadingZeros32(aFrac.x) - 1,
1378 if (0 <= shiftCount)
1379 __shortShift64Left(aFrac.y, aFrac.x, shiftCount, aFrac.y, aFrac.x);
1381 __shift64RightJamming(aFrac.y, aFrac.x, -shiftCount, aFrac.y, aFrac.x);
1383 return __roundAndPackFloat32(aSign, 0x9C - shiftCount, aFrac.x);
1386 /* Returns the result of converting the single-precision floating-point value
1387 * `a' to the double-precision floating-point format.
1390 __fp32_to_fp64(float f)
1392 uint a = floatBitsToUint(f);
1393 uint aFrac = a & 0x007FFFFFu;
1394 int aExp = int((a>>23) & 0xFFu);
1395 uint aSign = a & 0x80000000u;
1403 __shift64Right(nanHi, nanLo, 12, nanHi, nanLo);
1404 nanHi |= aSign | 0x7FF80000u;
1405 return packUint2x32(uvec2(nanLo, nanHi));
1407 return __packFloat64(aSign, 0x7FF, 0u, 0u);
1412 return __packFloat64(aSign, 0, 0u, 0u);
1413 /* Normalize subnormal */
1414 int shiftCount = __countLeadingZeros32(aFrac) - 8;
1415 aFrac <<= shiftCount;
1416 aExp = 1 - shiftCount;
1420 __shift64Right(aFrac, 0u, 3, zFrac0, zFrac1);
1421 return __packFloat64(aSign, aExp + 0x380, zFrac0, zFrac1);
1424 /* Adds the 96-bit value formed by concatenating `a0', `a1', and `a2' to the
1425 * 96-bit value formed by concatenating `b0', `b1', and `b2'. Addition is
1426 * modulo 2^96, so any carry out is lost. The result is broken into three
1427 * 32-bit pieces which are stored at the locations pointed to by `z0Ptr',
1428 * `z1Ptr', and `z2Ptr'.
1431 __add96(uint a0, uint a1, uint a2,
1432 uint b0, uint b1, uint b2,
1438 uint carry1 = uint(z2 < a2);
1440 uint carry0 = uint(z1 < a1);
1443 z0 += uint(z1 < carry1);
1450 /* Subtracts the 96-bit value formed by concatenating `b0', `b1', and `b2' from
1451 * the 96-bit value formed by concatenating `a0', `a1', and `a2'. Subtraction
1452 * is modulo 2^96, so any borrow out (carry out) is lost. The result is broken
1453 * into three 32-bit pieces which are stored at the locations pointed to by
1454 * `z0Ptr', `z1Ptr', and `z2Ptr'.
1457 __sub96(uint a0, uint a1, uint a2,
1458 uint b0, uint b1, uint b2,
1464 uint borrow1 = uint(a2 < b2);
1466 uint borrow0 = uint(a1 < b1);
1468 z0 -= uint(z1 < borrow1);
1476 /* Returns an approximation to the 32-bit integer quotient obtained by dividing
1477 * `b' into the 64-bit value formed by concatenating `a0' and `a1'. The
1478 * divisor `b' must be at least 2^31. If q is the exact quotient truncated
1479 * toward zero, the approximation returned lies between q and q + 2 inclusive.
1480 * If the exact quotient q is larger than 32 bits, the maximum positive 32-bit
1481 * unsigned integer is returned.
1484 __estimateDiv64To32(uint a0, uint a1, uint b)
1497 z = (b0<<16 <= a0) ? 0xFFFF0000u : (a0 / b0)<<16;
1498 umulExtended(b, z, term0, term1);
1499 __sub64(a0, a1, term0, term1, rem0, rem1);
1500 while (int(rem0) < 0) {
1503 __add64(rem0, rem1, b0, b1, rem0, rem1);
1505 rem0 = (rem0<<16) | (rem1>>16);
1506 z |= (b0<<16 <= rem0) ? 0xFFFFu : rem0 / b0;
1511 __sqrtOddAdjustments(int index)
1551 __sqrtEvenAdjustments(int index)
1590 /* Returns an approximation to the square root of the 32-bit significand given
1591 * by `a'. Considered as an integer, `a' must be at least 2^31. If bit 0 of
1592 * `aExp' (the least significant bit) is 1, the integer returned approximates
1593 * 2^31*sqrt(`a'/2^31), where `a' is considered an integer. If bit 0 of `aExp'
1594 * is 0, the integer returned approximates 2^31*sqrt(`a'/2^30). In either
1595 * case, the approximation returned lies strictly within +/-2 of the exact
1599 __estimateSqrt32(int aExp, uint a)
1603 int index = int(a>>27 & 15u);
1604 if ((aExp & 1) != 0) {
1605 z = 0x4000u + (a>>17) - __sqrtOddAdjustments(index);
1606 z = ((a / z)<<14) + (z<<15);
1609 z = 0x8000u + (a>>17) - __sqrtEvenAdjustments(index);
1611 z = (0x20000u <= z) ? 0xFFFF8000u : (z<<15);
1613 return uint(int(a)>>1);
1615 return ((__estimateDiv64To32(a, 0u, z))>>1) + (z>>1);
1618 /* Returns the square root of the double-precision floating-point value `a'.
1619 * The operation is performed according to the IEEE Standard for Floating-Point
1623 __fsqrt64(uint64_t a)
1628 uint doubleZFrac0 = 0u;
1637 uint64_t default_nan = 0xFFFFFFFFFFFFFFFFUL;
1639 uint aFracLo = __extractFloat64FracLo(a);
1640 uint aFracHi = __extractFloat64FracHi(a);
1641 int aExp = __extractFloat64Exp(a);
1642 uint aSign = __extractFloat64Sign(a);
1643 if (aExp == 0x7FF) {
1644 if ((aFracHi | aFracLo) != 0u)
1645 return __propagateFloat64NaN(a, a);
1651 if ((uint(aExp) | aFracHi | aFracLo) == 0u)
1656 if ((aFracHi | aFracLo) == 0u)
1657 return __packFloat64(0u, 0, 0u, 0u);
1658 __normalizeFloat64Subnormal(aFracHi, aFracLo, aExp, aFracHi, aFracLo);
1660 int zExp = ((aExp - 0x3FF)>>1) + 0x3FE;
1661 aFracHi |= 0x00100000u;
1662 __shortShift64Left(aFracHi, aFracLo, 11, term0, term1);
1663 zFrac0 = (__estimateSqrt32(aExp, term0)>>1) + 1u;
1665 zFrac0 = 0x7FFFFFFFu;
1666 doubleZFrac0 = zFrac0 + zFrac0;
1667 __shortShift64Left(aFracHi, aFracLo, 9 - (aExp & 1), aFracHi, aFracLo);
1668 umulExtended(zFrac0, zFrac0, term0, term1);
1669 __sub64(aFracHi, aFracLo, term0, term1, rem0, rem1);
1670 while (int(rem0) < 0) {
1673 __add64(rem0, rem1, 0u, doubleZFrac0 | 1u, rem0, rem1);
1675 zFrac1 = __estimateDiv64To32(rem1, 0u, doubleZFrac0);
1676 if ((zFrac1 & 0x1FFu) <= 5u) {
1679 umulExtended(doubleZFrac0, zFrac1, term1, term2);
1680 __sub64(rem1, 0u, term1, term2, rem1, rem2);
1681 umulExtended(zFrac1, zFrac1, term2, term3);
1682 __sub96(rem1, rem2, 0u, 0u, term2, term3, rem1, rem2, rem3);
1683 while (int(rem1) < 0) {
1685 __shortShift64Left(0u, zFrac1, 1, term2, term3);
1687 term2 |= doubleZFrac0;
1688 __add96(rem1, rem2, rem3, 0u, term2, term3, rem1, rem2, rem3);
1690 zFrac1 |= uint((rem1 | rem2 | rem3) != 0u);
1692 __shift64ExtraRightJamming(zFrac0, zFrac1, 0u, 10, zFrac0, zFrac1, zFrac2);
1693 return __roundAndPackFloat64(0u, zExp, zFrac0, zFrac1, zFrac2);
1697 __ftrunc64(uint64_t __a)
1699 uvec2 a = unpackUint2x32(__a);
1700 int aExp = __extractFloat64Exp(__a);
1704 int unbiasedExp = aExp - 1023;
1705 int fracBits = 52 - unbiasedExp;
1706 uint maskLo = mix(~0u << fracBits, 0u, fracBits >= 32);
1707 uint maskHi = mix(~0u << (fracBits - 32), ~0u, fracBits < 33);
1711 zLo = mix(zLo, 0u, unbiasedExp < 0);
1712 zHi = mix(zHi, 0u, unbiasedExp < 0);
1713 zLo = mix(zLo, a.x, unbiasedExp > 52);
1714 zHi = mix(zHi, a.y, unbiasedExp > 52);
1715 return packUint2x32(uvec2(zLo, zHi));
1719 __ffloor64(uint64_t a)
1721 /* The big assumtion is that when 'a' is NaN, __ftrunc(a) returns a. Based
1722 * on that assumption, NaN values that don't have the sign bit will safely
1723 * return NaN (identity). This is guarded by RELAXED_NAN_PROPAGATION
1724 * because otherwise the NaN should have the "signal" bit set. The
1725 * __fadd64 will ensure that occurs.
1728 #if defined RELAXED_NAN_PROPAGATION
1729 int(unpackUint2x32(a).y) >= 0
1734 uint64_t tr = __ftrunc64(a);
1736 if (is_positive || __feq64(tr, a)) {
1739 return __fadd64(tr, 0xbff0000000000000ul /* -1.0 */);
1744 __fround64(uint64_t __a)
1746 uvec2 a = unpackUint2x32(__a);
1747 int unbiasedExp = __extractFloat64Exp(__a) - 1023;
1751 if (unbiasedExp < 20) {
1752 if (unbiasedExp < 0) {
1753 if ((aHi & 0x80000000u) != 0u && aLo == 0u) {
1757 if ((a.y & 0x000FFFFFu) == 0u && a.x == 0u) {
1759 return packUint2x32(uvec2(aLo, aHi));
1761 aHi = mix(aHi, (aHi | 0x3FF00000u), unbiasedExp == -1);
1764 uint maskExp = 0x000FFFFFu >> unbiasedExp;
1765 uint lastBit = maskExp + 1;
1766 aHi += 0x00080000u >> unbiasedExp;
1767 if ((aHi & maskExp) == 0u)
1772 } else if (unbiasedExp > 51 || unbiasedExp == 1024) {
1775 uint maskExp = 0xFFFFFFFFu >> (unbiasedExp - 20);
1776 if ((aLo & maskExp) == 0u)
1778 uint tmp = aLo + (1u << (51 - unbiasedExp));
1785 return packUint2x32(uvec2(aLo, aHi));
1789 __fmin64(uint64_t a, uint64_t b)
1791 /* This weird layout matters. Doing the "obvious" thing results in extra
1792 * flow control being inserted to implement the short-circuit evaluation
1793 * rules. Flow control is bad!
1795 bool b_nan = __is_nan(b);
1796 bool a_lt_b = __flt64_nonnan(a, b);
1797 bool a_nan = __is_nan(a);
1799 return (b_nan || a_lt_b) && !a_nan ? a : b;
1803 __fmax64(uint64_t a, uint64_t b)
1805 /* This weird layout matters. Doing the "obvious" thing results in extra
1806 * flow control being inserted to implement the short-circuit evaluation
1807 * rules. Flow control is bad!
1809 bool b_nan = __is_nan(b);
1810 bool a_lt_b = __flt64_nonnan(a, b);
1811 bool a_nan = __is_nan(a);
1813 return (b_nan || a_lt_b) && !a_nan ? b : a;
1817 __ffract64(uint64_t a)
1819 return __fadd64(a, __fneg64(__ffloor64(a)));