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