2 * Copyright © 2015 Intel Corporation
4 * Permission is hereby granted, free of charge, to any person obtaining a
5 * copy of this software and associated documentation files (the "Software"),
6 * to deal in the Software without restriction, including without limitation
7 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
8 * and/or sell copies of the Software, and to permit persons to whom the
9 * Software is furnished to do so, subject to the following conditions:
11 * The above copyright notice and this permission notice (including the next
12 * paragraph) shall be included in all copies or substantial portions of the
15 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
16 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
17 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
18 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
19 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
20 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
26 #include "nir_builder.h"
30 * Lowers some unsupported double operations, using only:
32 * - pack/unpackDouble2x32
33 * - conversion to/from single-precision
34 * - double add, mul, and fma
35 * - conditional select
36 * - 32-bit integer and floating point arithmetic
39 /* Creates a double with the exponent bits set to a given integer value */
41 set_exponent(nir_builder
*b
, nir_ssa_def
*src
, nir_ssa_def
*exp
)
43 /* Split into bits 0-31 and 32-63 */
44 nir_ssa_def
*lo
= nir_unpack_double_2x32_split_x(b
, src
);
45 nir_ssa_def
*hi
= nir_unpack_double_2x32_split_y(b
, src
);
47 /* The exponent is bits 52-62, or 20-30 of the high word, so set the exponent
50 nir_ssa_def
*new_hi
= nir_bfi(b
, nir_imm_int(b
, 0x7ff00000), exp
, hi
);
52 return nir_pack_double_2x32_split(b
, lo
, new_hi
);
56 get_exponent(nir_builder
*b
, nir_ssa_def
*src
)
59 nir_ssa_def
*hi
= nir_unpack_double_2x32_split_y(b
, src
);
61 /* extract bits 20-30 of the high word */
62 return nir_ubitfield_extract(b
, hi
, nir_imm_int(b
, 20), nir_imm_int(b
, 11));
65 /* Return infinity with the sign of the given source which is +/-0 */
68 get_signed_inf(nir_builder
*b
, nir_ssa_def
*zero
)
70 nir_ssa_def
*zero_hi
= nir_unpack_double_2x32_split_y(b
, zero
);
72 /* The bit pattern for infinity is 0x7ff0000000000000, where the sign bit
73 * is the highest bit. Only the sign bit can be non-zero in the passed in
74 * source. So we essentially need to OR the infinity and the zero, except
75 * the low 32 bits are always 0 so we can construct the correct high 32
76 * bits and then pack it together with zero low 32 bits.
78 nir_ssa_def
*inf_hi
= nir_ior(b
, nir_imm_int(b
, 0x7ff00000), zero_hi
);
79 return nir_pack_double_2x32_split(b
, nir_imm_int(b
, 0), inf_hi
);
83 * Generates the correctly-signed infinity if the source was zero, and flushes
84 * the result to 0 if the source was infinity or the calculated exponent was
85 * too small to be representable.
89 fix_inv_result(nir_builder
*b
, nir_ssa_def
*res
, nir_ssa_def
*src
,
92 /* If the exponent is too small or the original input was infinity/NaN,
93 * force the result to 0 (flush denorms) to avoid the work of handling
94 * denorms properly. Note that this doesn't preserve positive/negative
95 * zeros, but GLSL doesn't require it.
97 res
= nir_bcsel(b
, nir_ior(b
, nir_ige(b
, nir_imm_int(b
, 0), exp
),
98 nir_feq(b
, nir_fabs(b
, src
),
99 nir_imm_double(b
, INFINITY
))),
100 nir_imm_double(b
, 0.0f
), res
);
102 /* If the original input was 0, generate the correctly-signed infinity */
103 res
= nir_bcsel(b
, nir_fne(b
, src
, nir_imm_double(b
, 0.0f
)),
104 res
, get_signed_inf(b
, src
));
111 lower_rcp(nir_builder
*b
, nir_ssa_def
*src
)
113 /* normalize the input to avoid range issues */
114 nir_ssa_def
*src_norm
= set_exponent(b
, src
, nir_imm_int(b
, 1023));
116 /* cast to float, do an rcp, and then cast back to get an approximate
119 nir_ssa_def
*ra
= nir_f2d(b
, nir_frcp(b
, nir_d2f(b
, src_norm
)));
121 /* Fixup the exponent of the result - note that we check if this is too
124 nir_ssa_def
*new_exp
= nir_isub(b
, get_exponent(b
, ra
),
125 nir_isub(b
, get_exponent(b
, src
),
126 nir_imm_int(b
, 1023)));
128 ra
= set_exponent(b
, ra
, new_exp
);
130 /* Do a few Newton-Raphson steps to improve precision.
132 * Each step doubles the precision, and we started off with around 24 bits,
133 * so we only need to do 2 steps to get to full precision. The step is:
135 * x_new = x * (2 - x*src)
137 * But we can re-arrange this to improve precision by using another fused
140 * x_new = x + x * (1 - x*src)
142 * See https://en.wikipedia.org/wiki/Division_algorithm for more details.
145 ra
= nir_ffma(b
, ra
, nir_ffma(b
, ra
, src
, nir_imm_double(b
, -1)), ra
);
146 ra
= nir_ffma(b
, ra
, nir_ffma(b
, ra
, src
, nir_imm_double(b
, -1)), ra
);
148 return fix_inv_result(b
, ra
, src
, new_exp
);
152 lower_sqrt_rsq(nir_builder
*b
, nir_ssa_def
*src
, bool sqrt
)
154 /* We want to compute:
158 * When the exponent is even, this is equivalent to:
160 * 1/sqrt(m) * 2^(-e/2)
162 * and then the exponent is odd, this is equal to:
164 * 1/sqrt(m * 2) * 2^(-(e - 1)/2)
166 * where the m * 2 is absorbed into the exponent. So we want the exponent
167 * inside the square root to be 1 if e is odd and 0 if e is even, and we
168 * want to subtract off e/2 from the final exponent, rounded to negative
169 * infinity. We can do the former by first computing the unbiased exponent,
170 * and then AND'ing it with 1 to get 0 or 1, and we can do the latter by
171 * shifting right by 1.
174 nir_ssa_def
*unbiased_exp
= nir_isub(b
, get_exponent(b
, src
),
175 nir_imm_int(b
, 1023));
176 nir_ssa_def
*even
= nir_iand(b
, unbiased_exp
, nir_imm_int(b
, 1));
177 nir_ssa_def
*half
= nir_ishr(b
, unbiased_exp
, nir_imm_int(b
, 1));
179 nir_ssa_def
*src_norm
= set_exponent(b
, src
,
180 nir_iadd(b
, nir_imm_int(b
, 1023),
183 nir_ssa_def
*ra
= nir_f2d(b
, nir_frsq(b
, nir_d2f(b
, src_norm
)));
184 nir_ssa_def
*new_exp
= nir_isub(b
, get_exponent(b
, ra
), half
);
185 ra
= set_exponent(b
, ra
, new_exp
);
188 * The following implements an iterative algorithm that's very similar
189 * between sqrt and rsqrt. We start with an iteration of Goldschmit's
190 * algorithm, which looks like:
193 * y_0 = initial (single-precision) rsqrt estimate
197 * r_0 = .5 - h_0 * g_0
198 * g_1 = g_0 * r_0 + g_0
199 * h_1 = h_0 * r_0 + h_0
201 * Now g_1 ~= sqrt(a), and h_1 ~= 1/(2 * sqrt(a)). We could continue
202 * applying another round of Goldschmit, but since we would never refer
203 * back to a (the original source), we would add too much rounding error.
204 * So instead, we do one last round of Newton-Raphson, which has better
205 * rounding characteristics, to get the final rounding correct. This is
206 * split into two cases:
210 * Normally, doing a round of Newton-Raphson for sqrt involves taking a
211 * reciprocal of the original estimate, which is slow since it isn't
212 * supported in HW. But we can take advantage of the fact that we already
213 * computed a good estimate of 1/(2 * g_1) by rearranging it like so:
215 * g_2 = .5 * (g_1 + a / g_1)
216 * = g_1 + .5 * (a / g_1 - g_1)
217 * = g_1 + (.5 / g_1) * (a - g_1^2)
218 * = g_1 + h_1 * (a - g_1^2)
220 * The second term represents the error, and by splitting it out we can get
221 * better precision by computing it as part of a fused multiply-add. Since
222 * both Newton-Raphson and Goldschmit approximately double the precision of
223 * the result, these two steps should be enough.
227 * First off, note that the first round of the Goldschmit algorithm is
228 * really just a Newton-Raphson step in disguise:
230 * h_1 = h_0 * (.5 - h_0 * g_0) + h_0
231 * = h_0 * (1.5 - h_0 * g_0)
232 * = h_0 * (1.5 - .5 * a * y_0^2)
233 * = (.5 * y_0) * (1.5 - .5 * a * y_0^2)
235 * which is the standard formula multiplied by .5. Unlike in the sqrt case,
236 * we don't need the inverse to do a Newton-Raphson step; we just need h_1,
237 * so we can skip the calculation of g_1. Instead, we simply do another
238 * Newton-Raphson step:
241 * r_1 = .5 - h_1 * y_1 * a
242 * y_2 = y_1 * r_1 + y_1
244 * Where the difference from Goldschmit is that we calculate y_1 * a
245 * instead of using g_1. Doing it this way should be as fast as computing
246 * y_1 up front instead of h_1, and it lets us share the code for the
247 * initial Goldschmit step with the sqrt case.
249 * Putting it together, the computations are:
253 * r_0 = .5 - h_0 * g_0
254 * h_1 = h_0 * r_0 + h_0
256 * g_1 = g_0 * r_0 + g_0
257 * r_1 = a - g_1 * g_1
258 * g_2 = h_1 * r_1 + g_1
261 * r_1 = .5 - y_1 * (h_1 * a)
262 * y_2 = y_1 * r_1 + y_1
264 * For more on the ideas behind this, see "Software Division and Square
265 * Root Using Goldschmit's Algorithms" by Markstein and the Wikipedia page
267 * (https://en.wikipedia.org/wiki/Methods_of_computing_square_roots).
270 nir_ssa_def
*one_half
= nir_imm_double(b
, 0.5);
271 nir_ssa_def
*h_0
= nir_fmul(b
, one_half
, ra
);
272 nir_ssa_def
*g_0
= nir_fmul(b
, src
, ra
);
273 nir_ssa_def
*r_0
= nir_ffma(b
, nir_fneg(b
, h_0
), g_0
, one_half
);
274 nir_ssa_def
*h_1
= nir_ffma(b
, h_0
, r_0
, h_0
);
277 nir_ssa_def
*g_1
= nir_ffma(b
, g_0
, r_0
, g_0
);
278 nir_ssa_def
*r_1
= nir_ffma(b
, nir_fneg(b
, g_1
), g_1
, src
);
279 res
= nir_ffma(b
, h_1
, r_1
, g_1
);
281 nir_ssa_def
*y_1
= nir_fmul(b
, nir_imm_double(b
, 2.0), h_1
);
282 nir_ssa_def
*r_1
= nir_ffma(b
, nir_fneg(b
, y_1
), nir_fmul(b
, h_1
, src
),
284 res
= nir_ffma(b
, y_1
, r_1
, y_1
);
288 /* Here, the special cases we need to handle are
292 res
= nir_bcsel(b
, nir_ior(b
, nir_feq(b
, src
, nir_imm_double(b
, 0.0)),
293 nir_feq(b
, src
, nir_imm_double(b
, INFINITY
))),
296 res
= fix_inv_result(b
, res
, src
, new_exp
);
303 lower_trunc(nir_builder
*b
, nir_ssa_def
*src
)
305 nir_ssa_def
*unbiased_exp
= nir_isub(b
, get_exponent(b
, src
),
306 nir_imm_int(b
, 1023));
308 nir_ssa_def
*frac_bits
= nir_isub(b
, nir_imm_int(b
, 52), unbiased_exp
);
311 * Decide the operation to apply depending on the unbiased exponent:
313 * if (unbiased_exp < 0)
315 * else if (unbiased_exp > 52)
318 * return src & (~0 << frac_bits)
320 * Notice that the else branch is a 64-bit integer operation that we need
321 * to implement in terms of 32-bit integer arithmetics (at least until we
322 * support 64-bit integer arithmetics).
325 /* Compute "~0 << frac_bits" in terms of hi/lo 32-bit integer math */
326 nir_ssa_def
*mask_lo
=
328 nir_ige(b
, frac_bits
, nir_imm_int(b
, 32)),
330 nir_ishl(b
, nir_imm_int(b
, ~0), frac_bits
));
332 nir_ssa_def
*mask_hi
=
334 nir_ilt(b
, frac_bits
, nir_imm_int(b
, 33)),
338 nir_isub(b
, frac_bits
, nir_imm_int(b
, 32))));
340 nir_ssa_def
*src_lo
= nir_unpack_double_2x32_split_x(b
, src
);
341 nir_ssa_def
*src_hi
= nir_unpack_double_2x32_split_y(b
, src
);
345 nir_ilt(b
, unbiased_exp
, nir_imm_int(b
, 0)),
346 nir_imm_double(b
, 0.0),
347 nir_bcsel(b
, nir_ige(b
, unbiased_exp
, nir_imm_int(b
, 53)),
349 nir_pack_double_2x32_split(b
,
350 nir_iand(b
, mask_lo
, src_lo
),
351 nir_iand(b
, mask_hi
, src_hi
))));
355 lower_floor(nir_builder
*b
, nir_ssa_def
*src
)
358 * For x >= 0, floor(x) = trunc(x)
360 * - if x is integer, floor(x) = x
361 * - otherwise, floor(x) = trunc(x) - 1
363 nir_ssa_def
*tr
= nir_ftrunc(b
, src
);
364 nir_ssa_def
*positive
= nir_fge(b
, src
, nir_imm_double(b
, 0.0));
366 nir_ior(b
, positive
, nir_feq(b
, src
, tr
)),
368 nir_fsub(b
, tr
, nir_imm_double(b
, 1.0)));
372 lower_ceil(nir_builder
*b
, nir_ssa_def
*src
)
374 /* if x < 0, ceil(x) = trunc(x)
375 * else if (x - trunc(x) == 0), ceil(x) = x
376 * else, ceil(x) = trunc(x) + 1
378 nir_ssa_def
*tr
= nir_ftrunc(b
, src
);
379 nir_ssa_def
*negative
= nir_flt(b
, src
, nir_imm_double(b
, 0.0));
381 nir_ior(b
, negative
, nir_feq(b
, src
, tr
)),
383 nir_fadd(b
, tr
, nir_imm_double(b
, 1.0)));
387 lower_fract(nir_builder
*b
, nir_ssa_def
*src
)
389 return nir_fsub(b
, src
, nir_ffloor(b
, src
));
393 lower_round_even(nir_builder
*b
, nir_ssa_def
*src
)
395 /* If fract(src) == 0.5, then we will have to decide the rounding direction.
396 * We will do this by computing the mod(abs(src), 2) and testing if it
399 * We compute mod(abs(src), 2) as:
400 * abs(src) - 2.0 * floor(abs(src) / 2.0)
402 nir_ssa_def
*two
= nir_imm_double(b
, 2.0);
403 nir_ssa_def
*abs_src
= nir_fabs(b
, src
);
412 nir_imm_double(b
, 0.5)))));
415 * If fract(src) != 0.5, then we round as floor(src + 0.5)
417 * If fract(src) == 0.5, then we have to check the modulo:
419 * if it is < 1 we need a trunc operation so we get:
420 * 0.5 -> 0, -0.5 -> -0
421 * 2.5 -> 2, -2.5 -> -2
423 * otherwise we need to check if src >= 0, in which case we need to round
424 * upwards, or not, in which case we need to round downwards so we get:
425 * 1.5 -> 2, -1.5 -> -2
426 * 3.5 -> 4, -3.5 -> -4
428 nir_ssa_def
*fract
= nir_ffract(b
, src
);
430 nir_fne(b
, fract
, nir_imm_double(b
, 0.5)),
431 nir_ffloor(b
, nir_fadd(b
, src
, nir_imm_double(b
, 0.5))),
433 nir_flt(b
, mod
, nir_imm_double(b
, 1.0)),
436 nir_fge(b
, src
, nir_imm_double(b
, 0.0)),
437 nir_fadd(b
, src
, nir_imm_double(b
, 0.5)),
438 nir_fsub(b
, src
, nir_imm_double(b
, 0.5)))));
442 lower_mod(nir_builder
*b
, nir_ssa_def
*src0
, nir_ssa_def
*src1
)
444 /* mod(x,y) = x - y * floor(x/y)
446 * If the division is lowered, it could add some rounding errors that make
447 * floor() to return the quotient minus one when x = N * y. If this is the
448 * case, we return zero because mod(x, y) output value is [0, y).
450 nir_ssa_def
*floor
= nir_ffloor(b
, nir_fdiv(b
, src0
, src1
));
451 nir_ssa_def
*mod
= nir_fsub(b
, src0
, nir_fmul(b
, src1
, floor
));
454 nir_fne(b
, mod
, src1
),
456 nir_imm_double(b
, 0.0));
460 lower_doubles_instr(nir_alu_instr
*instr
, nir_lower_doubles_options options
)
462 assert(instr
->dest
.dest
.is_ssa
);
463 if (instr
->dest
.dest
.ssa
.bit_size
!= 64)
468 if (!(options
& nir_lower_drcp
))
473 if (!(options
& nir_lower_dsqrt
))
478 if (!(options
& nir_lower_drsq
))
483 if (!(options
& nir_lower_dtrunc
))
488 if (!(options
& nir_lower_dfloor
))
493 if (!(options
& nir_lower_dceil
))
498 if (!(options
& nir_lower_dfract
))
502 case nir_op_fround_even
:
503 if (!(options
& nir_lower_dround_even
))
508 if (!(options
& nir_lower_dmod
))
517 nir_builder_init(&bld
, nir_cf_node_get_function(&instr
->instr
.block
->cf_node
));
518 bld
.cursor
= nir_before_instr(&instr
->instr
);
520 nir_ssa_def
*src
= nir_fmov_alu(&bld
, instr
->src
[0],
521 instr
->dest
.dest
.ssa
.num_components
);
527 result
= lower_rcp(&bld
, src
);
530 result
= lower_sqrt_rsq(&bld
, src
, true);
533 result
= lower_sqrt_rsq(&bld
, src
, false);
536 result
= lower_trunc(&bld
, src
);
539 result
= lower_floor(&bld
, src
);
542 result
= lower_ceil(&bld
, src
);
545 result
= lower_fract(&bld
, src
);
547 case nir_op_fround_even
:
548 result
= lower_round_even(&bld
, src
);
552 nir_ssa_def
*src1
= nir_fmov_alu(&bld
, instr
->src
[1],
553 instr
->dest
.dest
.ssa
.num_components
);
554 result
= lower_mod(&bld
, src
, src1
);
558 unreachable("unhandled opcode");
561 nir_ssa_def_rewrite_uses(&instr
->dest
.dest
.ssa
, nir_src_for_ssa(result
));
562 nir_instr_remove(&instr
->instr
);
566 nir_lower_doubles(nir_shader
*shader
, nir_lower_doubles_options options
)
568 nir_foreach_function(function
, shader
) {
572 nir_foreach_block(block
, function
->impl
) {
573 nir_foreach_instr_safe(instr
, block
) {
574 if (instr
->type
== nir_instr_type_alu
)
575 lower_doubles_instr(nir_instr_as_alu(instr
), options
);