ftz 0.0.1
Fast, reproducible floating-point arithmetic
Loading...
Searching...
No Matches
webgpu_words.h
Go to the documentation of this file.
1#pragma once
2#include "ftz/config.h"
3
4
5namespace ftz { namespace detail { namespace math {
6 struct unsigned_word_product {
7 unsigned int low, high;
8 };
9 // Independently derived base-2^16 long multiplication. Each limb product and
10 // carry addition fits uint32; no native64 type, shift or intermediate is used.
11 inline unsigned_word_product unsigned_multiply_words(unsigned int a, unsigned int b) {
12 unsigned int a0 = a & 65535u, a1 = a >> 16;
13 unsigned int b0 = b & 65535u, b1 = b >> 16;
14 unsigned int low = a0 * b0;
15 unsigned int middle0 = a1 * b0 + (low >> 16);
16 unsigned int middle1 = a0 * b1 + (middle0 & 65535u);
17 unsigned_word_product result;
18 result.low = (middle1 << 16) | (low & 65535u);
19 result.high = a1 * b1 + (middle0 >> 16) + (middle1 >> 16);
20 return result;
21 }
22
23
24}}}
25
26namespace ftz { namespace detail { namespace math {
27 // Exact binary32 word operations for rare rounding-boundary repair and
28 // large-angle reduction. RNE with signed post-round FTZ retains the minimum
29 // normal result. Invalid/nonfinite input or overflow returns {0,0}.
30 struct fp32_result { unsigned int bits, valid; };
31 inline fp32_result fp32_result_of(unsigned int bits, unsigned int valid) {
32 fp32_result r; r.bits = bits; r.valid = valid; return r;
33 }
34 inline unsigned int fp32_flush_word(unsigned int bits) {
35 return (bits & 0x7fffffffu) < 0x00800000u ? bits & 0x80000000u : bits;
36 }
37 inline bool fp32_finite(unsigned int bits) {
38 return (bits & 0x7f800000u) != 0x7f800000u;
39 }
40 inline unsigned_word_product word_pair(unsigned int low, unsigned int high) {
41 unsigned_word_product r; r.low = low; r.high = high; return r;
42 }
43 inline bool word_pair_nonzero(unsigned_word_product a) { return (a.low | a.high) != 0u; }
44 inline unsigned_word_product word_pair_add(unsigned_word_product a, unsigned_word_product b) {
45 unsigned int low = a.low + b.low;
46 return word_pair(low, a.high + b.high + (low < a.low ? 1u : 0u));
47 }
48 inline unsigned_word_product word_pair_sub(unsigned_word_product a, unsigned_word_product b) {
49 return word_pair(a.low - b.low, a.high - b.high - (a.low < b.low ? 1u : 0u));
50 }
51 inline bool word_pair_less(unsigned_word_product a, unsigned_word_product b) {
52 return a.high < b.high || (a.high == b.high && a.low < b.low);
53 }
54 inline unsigned_word_product word_pair_left(unsigned_word_product a, unsigned int shift) {
55 if (shift == 0u) return a;
56 if (shift < 32u) return word_pair(a.low << shift, (a.high << shift) | (a.low >> (32u - shift)));
57 if (shift < 64u) return word_pair(0u, a.low << (shift - 32u));
58 return word_pair(0u, 0u);
59 }
60 // Shift with a sticky low bit: discarded information affects rounding, never
61 // the retained high value. Every variable shift is guarded against 32/64.
62 inline unsigned_word_product word_pair_right_jam(unsigned_word_product a, unsigned int shift) {
63 if (shift == 0u) return a;
64 if (shift < 32u) {
65 unsigned int low = (a.low >> shift) | (a.high << (32u - shift));
66 return word_pair(low | ((a.low << (32u - shift)) != 0u ? 1u : 0u), a.high >> shift);
67 }
68 if (shift == 32u) return word_pair(a.high | (a.low != 0u ? 1u : 0u), 0u);
69 if (shift < 64u) {
70 unsigned int tail = a.low | (a.high << (64u - shift));
71 return word_pair((a.high >> (shift - 32u)) | (tail != 0u ? 1u : 0u), 0u);
72 }
73 return word_pair(word_pair_nonzero(a) ? 1u : 0u, 0u);
74 }
75 inline unsigned int word_leading_zeros(unsigned int word) {
76 if (word == 0u) return 32u;
77 unsigned int count = 0u;
78 if ((word & 0xffff0000u) == 0u) { word <<= 16; count += 16u; }
79 if ((word & 0xff000000u) == 0u) { word <<= 8; count += 8u; }
80 if ((word & 0xf0000000u) == 0u) { word <<= 4; count += 4u; }
81 if ((word & 0xc0000000u) == 0u) { word <<= 2; count += 2u; }
82 if ((word & 0x80000000u) == 0u) ++count;
83 return count;
84 }
85 inline unsigned int word_pair_top(unsigned_word_product a) {
86 return a.high != 0u ? 63u - word_leading_zeros(a.high) : 31u - word_leading_zeros(a.low);
87 }
88 // Magnitude represents value / 2^(exponent-61), with sufficient low rounding
89 // information. Alignment can lose bits only when the operand exponents differ
90 // by >14, which excludes deep cancellation; the final rounding shift then
91 // discards the jam bit and keeps its sticky meaning.
92 inline fp32_result fp32_pack(unsigned_word_product magnitude, int exponent, unsigned int sign) {
93 if (!word_pair_nonzero(magnitude)) return fp32_result_of(sign, 1u);
94 int top = (int)word_pair_top(magnitude);
95 int biased = exponent - 61 + top + 127;
96 if (biased >= 255) return fp32_result_of(0u, 0u);
97 int shift = top - 23;
98 int denormal_shift = -88 - exponent;
99 if (shift < denormal_shift) shift = denormal_shift;
100 unsigned int quotient;
101 if (shift <= 0) {
102 quotient = word_pair_left(magnitude, (unsigned int)(-shift)).low;
103 } else if (shift == 1) {
104 quotient = (magnitude.low >> 1) | (magnitude.high << 31);
105 quotient += (magnitude.low & quotient & 1u);
106 } else {
107 unsigned int window = word_pair_right_jam(magnitude, (unsigned int)(shift - 2)).low;
108 quotient = window >> 2;
109 quotient += ((window & 2u) != 0u && ((window & 1u) != 0u || (quotient & 1u) != 0u)) ? 1u : 0u;
110 }
111 if (biased <= 0)
112 return fp32_result_of(sign | (quotient >= 0x00800000u ? 0x00800000u : 0u), 1u);
113 if (quotient == 0x01000000u) { quotient >>= 1; ++biased; }
114 if (biased >= 255) return fp32_result_of(0u, 0u);
115 return fp32_result_of(sign | ((unsigned int)biased << 23) | (quotient & 0x007fffffu), 1u);
116 }
117
118 inline fp32_result fp32_fma_words(unsigned int a, unsigned int b, unsigned int c) {
119 if (!fp32_finite(a) || !fp32_finite(b) || !fp32_finite(c))
120 return fp32_result_of(0u, 0u);
121 a = fp32_flush_word(a); b = fp32_flush_word(b); c = fp32_flush_word(c);
122 unsigned int sign = (a ^ b) & 0x80000000u;
123 unsigned int csign = c & 0x80000000u;
124 unsigned int ea = (a >> 23) & 255u, eb = (b >> 23) & 255u, ec = (c >> 23) & 255u;
125 if (ea == 0u || eb == 0u)
126 return fp32_result_of(ec == 0u ? sign & csign : c, 1u);
127 unsigned int ma = (a & 0x7fffffu) | 0x800000u, mb = (b & 0x7fffffu) | 0x800000u;
128 unsigned_word_product product = unsigned_multiply_words(ma, mb);
129 unsigned int top = word_pair_top(product);
130 int exponent = (int)ea + (int)eb - 300 + (int)top;
131 product = word_pair_left(product, 61u - top);
132 if (ec == 0u) return fp32_pack(product, exponent, sign);
133 unsigned_word_product addend = word_pair_left(word_pair((c & 0x7fffffu) | 0x800000u, 0u), 38u);
134 int ce = (int)ec - 127;
135 if (exponent < ce) {
136 product = word_pair_right_jam(product, (unsigned int)(ce - exponent));
137 exponent = ce;
138 } else {
139 addend = word_pair_right_jam(addend, (unsigned int)(exponent - ce));
140 }
141 unsigned_word_product magnitude;
142 if (sign == csign) {
143 magnitude = word_pair_add(product, addend);
144 } else if (word_pair_less(product, addend)) {
145 magnitude = word_pair_sub(addend, product); sign = csign;
146 } else {
147 magnitude = word_pair_sub(product, addend);
148 if (!word_pair_nonzero(magnitude)) sign = 0u;
149 }
150 return fp32_pack(magnitude, exponent, sign);
151 }
152
153}}}
154