ftz shaders 0.0.1
HLSL 2021 interface
Loading...
Searching...
No Matches
math.h
Go to the documentation of this file.
1#pragma once
2#include "ftz/config.h"
3#include "ftz/math/float.h"
4#include "ftz/math/exp_coefficients.h"
5
6namespace ftz { namespace detail { namespace math {
7
8 // The host and shader use the same endpoints and selected-degree graph.
9 // The active interval has n in [-126,128], and excludes the minimum-normal
10 // rounding strip. Only n=128 requires two normal scaling factors.
11 template <unsigned int Degree = 6>
12 uint2 exp_value(unsigned int bits) {
13 if ((bits & 0x7f800000u) == 0x7f800000u) return uint2(0u, 0u);
14 precise float x = asfloat(bits);
15 bool active = !(x < asfloat(0xc2aeac4fu));
16 bool overflow = x > asfloat(0x42b17217u);
17 precise float product = x * asfloat(0x3fb8aa3bu);
18 precise float n = round(product);
19 precise float first = fp32_fma<false>(n, asfloat(0xbf317200u), x);
20 precise float r = fp32_fma<false>(n, asfloat(0xb5bfbe8eu), first);
21 // Degree is a template constant: unused Horner stages disappear at compile time.
22 precise float y = fp32_fma<false>(r, asfloat(exp_coefficients<Degree>::leading), asfloat(exp_coefficients<Degree>::next));
23 if (Degree >= 7) y = fp32_fma<false>(r, y, asfloat(exp_coefficients<Degree>::c5));
24 if (Degree >= 6) y = fp32_fma<false>(r, y, asfloat(exp_coefficients<Degree>::c4));
25 if (Degree >= 5) y = fp32_fma<false>(r, y, asfloat(exp_coefficients<Degree>::c3));
26 if (Degree >= 4) y = fp32_fma<false>(r, y, asfloat(exp_coefficients<Degree>::c2));
27 if (Degree >= 3) y = fp32_fma<false>(r, y, asfloat(exp_coefficients<Degree>::c1));
28 if (Degree >= 2) y = fp32_fma<false>(r, y, asfloat(exp_coefficients<Degree>::c0));
29 // HLSL float-to-uint conversion has no ARM FCVTZU saturation contract.
30 // Keep only this conversion bounded; range masks do not feed the reducer.
31 bool high = n > 127.0f;
32 precise float biased = active && !overflow ? n + (high ? 126.0f : 127.0f) : 0.0f;
33 float factor = asfloat((unsigned int)biased << 23);
34 precise float first_scaled = y * factor;
35 precise float high_scaled = (high ? first_scaled : 0.0f) * 2.0f;
36 precise float scaled = high ? high_scaled : first_scaled;
37 precise float value = overflow ? asfloat(0x7f800000u) : active ? scaled : 0.0f;
38 return uint2(asuint(value), 1u);
39 }
40
41}}}
42
43namespace ftz { namespace detail { namespace math {
44 struct expm1_result {
45 float value;
46 unsigned int valid;
47 };
48 // Internal finite argument <=1; clamped below. Callers apply signed input FTZ.
49 float expm1_core(float input) {
50 precise float x = max(input, -18.0f);
51 precise float product = x * asfloat(0x3fb8aa3bu);
52 precise float n = round(product);
53 precise float r0 = fp32_fma<false>(n, asfloat(0xbf317200u), x);
54 precise float r1 = fp32_fma<false>(n, asfloat(0xb5bfbe8eu), r0);
55 precise float r = n == 0.0f ? x : r1;
56 bool tiny = (asuint(r) & 0x7fffffffu) <= 0x33000000u;
57#if FTZ_FP32_HARDWARE_FTZ
58 precise float t = r;
59#else
60 precise float t = tiny ? 0.0f : r;
61#endif
62 precise float z = t * t;
63 precise float h0 = asfloat(0x3493f27eu);
64 precise float h1 = fp32_fma<false>(t, h0, asfloat(0x3638ef1du));
65 precise float h2 = fp32_fma<false>(t, h1, asfloat(0x37d00d01u));
66 precise float h3 = fp32_fma<false>(t, h2, asfloat(0x39500d01u));
67 precise float h4 = fp32_fma<false>(t, h3, asfloat(0x3ab60b61u));
68 precise float h5 = fp32_fma<false>(t, h4, asfloat(0x3c088889u));
69 precise float h6 = fp32_fma<false>(t, h5, asfloat(0x3d2aaaabu));
70 precise float h7 = fp32_fma<false>(t, h6, asfloat(0x3e2aaaabu));
71 precise float h8 = fp32_fma<false>(t, h7, asfloat(0x3f000000u));
72 precise float fused = fp32_fma<false>(z, h8, r);
73 precise float p = tiny ? r : fused;
74 int index = (int)max(n, -24.0f);
75 precise float scale = asfloat((unsigned int)(index + 127) << 23u);
76 precise float offset = scale - 1.0f;
77 precise float scaled = fp32_fma<false>(scale, p, offset);
78 precise float result = n == 0.0f ? p : scaled;
79 result = n == -25.0f ? (p > 0.0f ? asfloat(0xbf7fffffu) : -1.0f) : result;
80 result = n < -25.0f ? -1.0f : result;
81 // Entry-canonicalized graph outputs are already normal or signed zero.
82 return fp32_ftz(result);
83 }
84 expm1_result expm1_checked(float input) {
85 unsigned int word = asuint(input), magnitude = word & 0x7fffffffu;
86 bool valid = magnitude < 0x7f800000u &&
87 ((word & 0x80000000u) != 0u || magnitude <= 0x3f800000u);
88 word = valid ? word : 0u;
89 word = (word & 0x7f800000u) == 0u ? (word & 0x80000000u) : word;
90 precise float evaluated = expm1_core(asfloat(word));
91 expm1_result result;
92 result.value = valid ? evaluated : 0.0f;
93 result.valid = valid ? 1u : 0u;
94 return result;
95 }
96 expm1_result damping_gain_checked(float input) {
97 unsigned int word = asuint(input), magnitude = word & 0x7fffffffu;
98 bool valid = magnitude < 0x7f800000u &&
99 ((word & 0x80000000u) == 0u || magnitude == 0u);
100 word = valid ? word : 0u;
101 word = (word & 0x7f800000u) == 0u ? (word & 0x80000000u) : word;
102 precise float evaluated = expm1_core(asfloat(word ^ 0x80000000u));
103 expm1_result result;
104 result.value = asfloat(valid ? (asuint(evaluated) ^ 0x80000000u) : 0u);
105 result.valid = valid ? 1u : 0u;
106 return result;
107 }
108
109}}}
110
111
112
113// Altered source: Pommier coefficients and reducer, signed input FTZ,
114// mask-before-square and bitwise quadrant reconstruction. Original notices below.
115// Require RNE, fused mad, and precise non-fused stages.
116namespace ftz { namespace detail { namespace math {
117 struct sincos_result { float sine, cosine; };
118
119 sincos_result trig_ftz_graph(float v, bool bounded) {
120 uint bits = asuint(v);
121 uint canonical = (bits & 0x7f800000u) == 0 ? bits & 0x80000000u : bits;
122 uint magnitude = canonical & 0x7fffffffu;
123 uint index = 0u;
124 precise float x;
125 if (bounded) {
126 precise float product = asfloat(magnitude) * asfloat(0x3fa2f983u);
127 index = ((uint)product + 1u) & 0xfffffffeu;
128 precise float multiple = (float)index;
129 precise float r0 = fp32_fma<false>(multiple, asfloat(0xbf490000u), asfloat(magnitude));
130 precise float r1 = fp32_fma<false>(multiple, asfloat(0xb97da000u), r0);
131 x = fp32_fma<false>(multiple, asfloat(0xb3222169u), r1);
132 } else {
133 x = asfloat(canonical);
134 }
135 bool active = (asuint(x) & 0x7fffffffu) > 0x39800000u;
136#if FTZ_FP32_HARDWARE_FTZ
137 precise float t = x;
138#else
139 precise float t = asfloat(active ? asuint(x) : 0u);
140#endif
141 precise float z = t * t;
142 precise float s0 = fp32_fma<false>(asfloat(0xb94ca1f9u), z, asfloat(0x3c08839eu));
143 precise float c0 = fp32_fma<false>(asfloat(0x37ccf5ceu), z, asfloat(0xbab6061au));
144 precise float s1 = fp32_fma<false>(s0, z, asfloat(0xbe2aaaa3u));
145 precise float c1 = fp32_fma<false>(c0, z, asfloat(0x3d2aaaa5u));
146 precise float s2 = s1 * z;
147 precise float c2 = c1 * z;
148 precise float c3 = c2 * z;
149 precise float half_z = z * 0.5f;
150 precise float c4 = c3 - half_z;
151 precise float s3 = fp32_fma<false>(s2, t, x);
152 precise float c5 = c4 + 1.0f;
153 uint sine = active ? asuint(s3) : asuint(x);
154 uint cosine = asuint(c5);
155 sincos_result result;
156 if (bounded) {
157 bool swap_pair = (index & 2u) != 0;
158 result.sine = asfloat((swap_pair ? cosine : sine) ^
159 (canonical & 0x80000000u) ^ ((index & 4u) << 29));
160 result.cosine = asfloat((swap_pair ? sine : cosine) ^
161 ((~(index - 2u) & 4u) << 29));
162 } else {
163 result.sine = asfloat(sine);
164 result.cosine = asfloat(cosine);
165 }
166 return result;
167 }
168
169 // Finite |x| <= 1; subnormal inputs become signed zero.
170 sincos_result sincos_reduced_ftz(float x) { return trig_ftz_graph(x, false); }
171 // Finite |x| < 8192. One shared three-FMA reducer.
172 sincos_result sincos_ftz(float x) { return trig_ftz_graph(x, true); }
173 float sin_reduced_ftz(float x) { return sincos_reduced_ftz(x).sine; }
174 float cos_reduced_ftz(float x) { return sincos_reduced_ftz(x).cosine; }
175 float sin_ftz(float x) { return sincos_ftz(x).sine; }
176 float cos_ftz(float x) { return sincos_ftz(x).cosine; }
177}}}
178
179/*
180 AVX implementation of sin, cos, sincos, exp and log
181
182 Based on "sse_mathfun.h", by Julien Pommier
183 http://gruntthepeon.free.fr/ssemath/
184
185 Copyright (C) 2012 Giovanni Garberoglio
186 Interdisciplinary Laboratory for Computational Science (LISC)
187 Fondazione Bruno Kessler and University of Trento
188 via Sommarive, 18
189 I-38123 Trento (Italy)
190
191 This software is provided 'as-is', without any express or implied
192 warranty. In no event will the authors be held liable for any damages
193 arising from the use of this software.
194
195 Permission is granted to anyone to use this software for any purpose,
196 including commercial applications, and to alter it and redistribute it
197 freely, subject to the following restrictions:
198
199 1. The origin of this software must not be misrepresented; you must not
200 claim that you wrote the original software. If you use this software
201 in a product, an acknowledgment in the product documentation would be
202 appreciated but is not required.
203 2. Altered source versions must be plainly marked as such, and must not be
204 misrepresented as being the original software.
205 3. This notice may not be removed or altered from any source distribution.
206*/
207
208/* RTS repository license (retained verbatim):
209Software License Agreement (BSD 2-Clause License)
210========================================
211
212Copyright 2017 Edward Kmett
213
214Redistribution and use in source and binary forms, with or without
215modification, are permitted provided that the following conditions are met:
216
217 * Redistributions of source code must retain the above copyright
218 notice, this list of conditions and the following disclaimer.
219
220 * Redistributions in binary form must reproduce the above copyright
221 notice, this list of conditions and the following disclaimer in the
222 documentation and/or other materials provided with the distribution.
223
224THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
225ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
226WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
227DISCLAIMED. IN NO EVENT SHALL YAHOO! INC. BE LIABLE FOR ANY
228DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
229(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
230LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
231ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
232(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
233SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
234*/
235