src/compiler/nir/nir_lower_double_ops.c - third_party/mesa - Git at Google

 /*
  * Copyright © 2015 Intel Corporation
  *
  * Permission is hereby granted, free of charge, to any person obtaining a
  * copy of this software and associated documentation files (the "Software"),
  * to deal in the Software without restriction, including without limitation
  * the rights to use, copy, modify, merge, publish, distribute, sublicense,
  * and/or sell copies of the Software, and to permit persons to whom the
  * Software is furnished to do so, subject to the following conditions:
  *
  * The above copyright notice and this permission notice (including the next
  * paragraph) shall be included in all copies or substantial portions of the
  * Software.
  *
  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
  * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.  IN NO EVENT SHALL
  * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
  * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
  * IN THE SOFTWARE.
  *
  */

 #include "nir.h"
 #include "nir_builder.h"
 #include "c99_math.h"

 /*
  * Lowers some unsupported double operations, using only:
  *
  * - pack/unpackDouble2x32
  * - conversion to/from single-precision
  * - double add, mul, and fma
  * - conditional select
  * - 32-bit integer and floating point arithmetic
  */

 /* Creates a double with the exponent bits set to a given integer value */
 static nir_ssa_def *
 set_exponent(nir_builder *b, nir_ssa_def *src, nir_ssa_def *exp)
 {
    /* Split into bits 0-31 and 32-63 */
    nir_ssa_def *lo = nir_unpack_64_2x32_split_x(b, src);
    nir_ssa_def *hi = nir_unpack_64_2x32_split_y(b, src);

    /* The exponent is bits 52-62, or 20-30 of the high word, so set the exponent
     * to 1023
     */
    nir_ssa_def *new_hi = nir_bfi(b, nir_imm_int(b, 0x7ff00000), exp, hi);
    /* recombine */
    return nir_pack_64_2x32_split(b, lo, new_hi);
 }

 static nir_ssa_def *
 get_exponent(nir_builder *b, nir_ssa_def *src)
 {
    /* get bits 32-63 */
    nir_ssa_def *hi = nir_unpack_64_2x32_split_y(b, src);

    /* extract bits 20-30 of the high word */
    return nir_ubitfield_extract(b, hi, nir_imm_int(b, 20), nir_imm_int(b, 11));
 }

 /* Return infinity with the sign of the given source which is +/-0 */

 static nir_ssa_def *
 get_signed_inf(nir_builder *b, nir_ssa_def *zero)
 {
    nir_ssa_def *zero_hi = nir_unpack_64_2x32_split_y(b, zero);

    /* The bit pattern for infinity is 0x7ff0000000000000, where the sign bit
     * is the highest bit. Only the sign bit can be non-zero in the passed in
     * source. So we essentially need to OR the infinity and the zero, except
     * the low 32 bits are always 0 so we can construct the correct high 32
     * bits and then pack it together with zero low 32 bits.
     */
    nir_ssa_def *inf_hi = nir_ior(b, nir_imm_int(b, 0x7ff00000), zero_hi);
    return nir_pack_64_2x32_split(b, nir_imm_int(b, 0), inf_hi);
 }

 /*
  * Generates the correctly-signed infinity if the source was zero, and flushes
  * the result to 0 if the source was infinity or the calculated exponent was
  * too small to be representable.
  */

 static nir_ssa_def *
 fix_inv_result(nir_builder *b, nir_ssa_def *res, nir_ssa_def *src,
                nir_ssa_def *exp)
 {
    /* If the exponent is too small or the original input was infinity/NaN,
     * force the result to 0 (flush denorms) to avoid the work of handling
     * denorms properly. Note that this doesn't preserve positive/negative
     * zeros, but GLSL doesn't require it.
     */
    res = nir_bcsel(b, nir_ior(b, nir_ige(b, nir_imm_int(b, 0), exp),
                               nir_feq(b, nir_fabs(b, src),
                                       nir_imm_double(b, INFINITY))),
                    nir_imm_double(b, 0.0f), res);

    /* If the original input was 0, generate the correctly-signed infinity */
    res = nir_bcsel(b, nir_fne(b, src, nir_imm_double(b, 0.0f)),
                    res, get_signed_inf(b, src));

    return res;

 }

 static nir_ssa_def *
 lower_rcp(nir_builder *b, nir_ssa_def *src)
 {
    /* normalize the input to avoid range issues */
    nir_ssa_def *src_norm = set_exponent(b, src, nir_imm_int(b, 1023));

    /* cast to float, do an rcp, and then cast back to get an approximate
     * result
     */
    nir_ssa_def *ra = nir_f2f64(b, nir_frcp(b, nir_f2f32(b, src_norm)));

    /* Fixup the exponent of the result - note that we check if this is too
     * small below.
     */
    nir_ssa_def *new_exp = nir_isub(b, get_exponent(b, ra),
                                    nir_isub(b, get_exponent(b, src),
                                             nir_imm_int(b, 1023)));

    ra = set_exponent(b, ra, new_exp);

    /* Do a few Newton-Raphson steps to improve precision.
     *
     * Each step doubles the precision, and we started off with around 24 bits,
     * so we only need to do 2 steps to get to full precision. The step is:
     *
     * x_new = x * (2 - x*src)
     *
     * But we can re-arrange this to improve precision by using another fused
     * multiply-add:
     *
     * x_new = x + x * (1 - x*src)
     *
     * See https://en.wikipedia.org/wiki/Division_algorithm for more details.
     */

    ra = nir_ffma(b, ra, nir_ffma(b, ra, src, nir_imm_double(b, -1)), ra);
    ra = nir_ffma(b, ra, nir_ffma(b, ra, src, nir_imm_double(b, -1)), ra);

    return fix_inv_result(b, ra, src, new_exp);
 }

 static nir_ssa_def *
 lower_sqrt_rsq(nir_builder *b, nir_ssa_def *src, bool sqrt)
 {
    /* We want to compute:
     *
     * 1/sqrt(m * 2^e)
     *
     * When the exponent is even, this is equivalent to:
     *
     * 1/sqrt(m) * 2^(-e/2)
     *
     * and then the exponent is odd, this is equal to:
     *
     * 1/sqrt(m * 2) * 2^(-(e - 1)/2)
     *
     * where the m * 2 is absorbed into the exponent. So we want the exponent
     * inside the square root to be 1 if e is odd and 0 if e is even, and we
     * want to subtract off e/2 from the final exponent, rounded to negative
     * infinity. We can do the former by first computing the unbiased exponent,
     * and then AND'ing it with 1 to get 0 or 1, and we can do the latter by
     * shifting right by 1.
     */

    nir_ssa_def *unbiased_exp = nir_isub(b, get_exponent(b, src),
                                         nir_imm_int(b, 1023));
    nir_ssa_def *even = nir_iand(b, unbiased_exp, nir_imm_int(b, 1));
    nir_ssa_def *half = nir_ishr(b, unbiased_exp, nir_imm_int(b, 1));

    nir_ssa_def *src_norm = set_exponent(b, src,
                                         nir_iadd(b, nir_imm_int(b, 1023),
                                                  even));

    nir_ssa_def *ra = nir_f2f64(b, nir_frsq(b, nir_f2f32(b, src_norm)));
    nir_ssa_def *new_exp = nir_isub(b, get_exponent(b, ra), half);
    ra = set_exponent(b, ra, new_exp);

    /*
     * The following implements an iterative algorithm that's very similar
     * between sqrt and rsqrt. We start with an iteration of Goldschmit's
     * algorithm, which looks like:
     *
     * a = the source
     * y_0 = initial (single-precision) rsqrt estimate
     *
     * h_0 = .5 * y_0
     * g_0 = a * y_0
     * r_0 = .5 - h_0 * g_0
     * g_1 = g_0 * r_0 + g_0
     * h_1 = h_0 * r_0 + h_0
     *
     * Now g_1 ~= sqrt(a), and h_1 ~= 1/(2 * sqrt(a)). We could continue
     * applying another round of Goldschmit, but since we would never refer
     * back to a (the original source), we would add too much rounding error.
     * So instead, we do one last round of Newton-Raphson, which has better
     * rounding characteristics, to get the final rounding correct. This is
     * split into two cases:
     *
     * 1. sqrt
     *
     * Normally, doing a round of Newton-Raphson for sqrt involves taking a
     * reciprocal of the original estimate, which is slow since it isn't
     * supported in HW. But we can take advantage of the fact that we already
     * computed a good estimate of 1/(2 * g_1) by rearranging it like so:
     *
     * g_2 = .5 * (g_1 + a / g_1)
     *     = g_1 + .5 * (a / g_1 - g_1)
     *     = g_1 + (.5 / g_1) * (a - g_1^2)
     *     = g_1 + h_1 * (a - g_1^2)
     *
     * The second term represents the error, and by splitting it out we can get
     * better precision by computing it as part of a fused multiply-add. Since
     * both Newton-Raphson and Goldschmit approximately double the precision of
     * the result, these two steps should be enough.
     *
     * 2. rsqrt
     *
     * First off, note that the first round of the Goldschmit algorithm is
     * really just a Newton-Raphson step in disguise:
     *
     * h_1 = h_0 * (.5 - h_0 * g_0) + h_0
     *     = h_0 * (1.5 - h_0 * g_0)
     *     = h_0 * (1.5 - .5 * a * y_0^2)
     *     = (.5 * y_0) * (1.5 - .5 * a * y_0^2)
     *
     * which is the standard formula multiplied by .5. Unlike in the sqrt case,
     * we don't need the inverse to do a Newton-Raphson step; we just need h_1,
     * so we can skip the calculation of g_1. Instead, we simply do another
     * Newton-Raphson step:
     *
     * y_1 = 2 * h_1
     * r_1 = .5 - h_1 * y_1 * a
     * y_2 = y_1 * r_1 + y_1
     *
     * Where the difference from Goldschmit is that we calculate y_1 * a
     * instead of using g_1. Doing it this way should be as fast as computing
     * y_1 up front instead of h_1, and it lets us share the code for the
     * initial Goldschmit step with the sqrt case.
     *
     * Putting it together, the computations are:
     *
     * h_0 = .5 * y_0
     * g_0 = a * y_0
     * r_0 = .5 - h_0 * g_0
     * h_1 = h_0 * r_0 + h_0
     * if sqrt:
     *    g_1 = g_0 * r_0 + g_0
     *    r_1 = a - g_1 * g_1
     *    g_2 = h_1 * r_1 + g_1
     * else:
     *    y_1 = 2 * h_1
     *    r_1 = .5 - y_1 * (h_1 * a)
     *    y_2 = y_1 * r_1 + y_1
     *
     * For more on the ideas behind this, see "Software Division and Square
     * Root Using Goldschmit's Algorithms" by Markstein and the Wikipedia page
     * on square roots
     * (https://en.wikipedia.org/wiki/Methods_of_computing_square_roots).
     */

    nir_ssa_def *one_half = nir_imm_double(b, 0.5);
    nir_ssa_def *h_0 = nir_fmul(b, one_half, ra);
    nir_ssa_def *g_0 = nir_fmul(b, src, ra);
    nir_ssa_def *r_0 = nir_ffma(b, nir_fneg(b, h_0), g_0, one_half);
    nir_ssa_def *h_1 = nir_ffma(b, h_0, r_0, h_0);
    nir_ssa_def *res;
    if (sqrt) {
       nir_ssa_def *g_1 = nir_ffma(b, g_0, r_0, g_0);
       nir_ssa_def *r_1 = nir_ffma(b, nir_fneg(b, g_1), g_1, src);
       res = nir_ffma(b, h_1, r_1, g_1);
    } else {
       nir_ssa_def *y_1 = nir_fmul(b, nir_imm_double(b, 2.0), h_1);
       nir_ssa_def *r_1 = nir_ffma(b, nir_fneg(b, y_1), nir_fmul(b, h_1, src),
                                   one_half);
       res = nir_ffma(b, y_1, r_1, y_1);
    }

    if (sqrt) {
       /* Here, the special cases we need to handle are
        * 0 -> 0 and
        * +inf -> +inf
        */
       res = nir_bcsel(b, nir_ior(b, nir_feq(b, src, nir_imm_double(b, 0.0)),
                                  nir_feq(b, src, nir_imm_double(b, INFINITY))),
                                  src, res);
    } else {
       res = fix_inv_result(b, res, src, new_exp);
    }

    return res;
 }

 static nir_ssa_def *
 lower_trunc(nir_builder *b, nir_ssa_def *src)
 {
    nir_ssa_def *unbiased_exp = nir_isub(b, get_exponent(b, src),
                                         nir_imm_int(b, 1023));

    nir_ssa_def *frac_bits = nir_isub(b, nir_imm_int(b, 52), unbiased_exp);

    /*
     * Decide the operation to apply depending on the unbiased exponent:
     *
     * if (unbiased_exp < 0)
     *    return 0
     * else if (unbiased_exp > 52)
     *    return src
     * else
     *    return src & (~0 << frac_bits)
     *
     * Notice that the else branch is a 64-bit integer operation that we need
     * to implement in terms of 32-bit integer arithmetics (at least until we
     * support 64-bit integer arithmetics).
     */

    /* Compute "~0 << frac_bits" in terms of hi/lo 32-bit integer math */
    nir_ssa_def *mask_lo =
       nir_bcsel(b,
                 nir_ige(b, frac_bits, nir_imm_int(b, 32)),
                 nir_imm_int(b, 0),
                 nir_ishl(b, nir_imm_int(b, ~0), frac_bits));

    nir_ssa_def *mask_hi =
       nir_bcsel(b,
                 nir_ilt(b, frac_bits, nir_imm_int(b, 33)),
                 nir_imm_int(b, ~0),
                 nir_ishl(b,
                          nir_imm_int(b, ~0),
                          nir_isub(b, frac_bits, nir_imm_int(b, 32))));

    nir_ssa_def *src_lo = nir_unpack_64_2x32_split_x(b, src);
    nir_ssa_def *src_hi = nir_unpack_64_2x32_split_y(b, src);

    return
       nir_bcsel(b,
                 nir_ilt(b, unbiased_exp, nir_imm_int(b, 0)),
                 nir_imm_double(b, 0.0),
                 nir_bcsel(b, nir_ige(b, unbiased_exp, nir_imm_int(b, 53)),
                           src,
                           nir_pack_64_2x32_split(b,
                                                  nir_iand(b, mask_lo, src_lo),
                                                  nir_iand(b, mask_hi, src_hi))));
 }

 static nir_ssa_def *
 lower_floor(nir_builder *b, nir_ssa_def *src)
 {
    /*
     * For x >= 0, floor(x) = trunc(x)
     * For x < 0,
     *    - if x is integer, floor(x) = x
     *    - otherwise, floor(x) = trunc(x) - 1
     */
    nir_ssa_def *tr = nir_ftrunc(b, src);
    nir_ssa_def *positive = nir_fge(b, src, nir_imm_double(b, 0.0));
    return nir_bcsel(b,
                     nir_ior(b, positive, nir_feq(b, src, tr)),
                     tr,
                     nir_fsub(b, tr, nir_imm_double(b, 1.0)));
 }

 static nir_ssa_def *
 lower_ceil(nir_builder *b, nir_ssa_def *src)
 {
    /* if x < 0,                    ceil(x) = trunc(x)
     * else if (x - trunc(x) == 0), ceil(x) = x
     * else,                        ceil(x) = trunc(x) + 1
     */
    nir_ssa_def *tr = nir_ftrunc(b, src);
    nir_ssa_def *negative = nir_flt(b, src, nir_imm_double(b, 0.0));
    return nir_bcsel(b,
                     nir_ior(b, negative, nir_feq(b, src, tr)),
                     tr,
                     nir_fadd(b, tr, nir_imm_double(b, 1.0)));
 }

 static nir_ssa_def *
 lower_fract(nir_builder *b, nir_ssa_def *src)
 {
    return nir_fsub(b, src, nir_ffloor(b, src));
 }

 static nir_ssa_def *
 lower_round_even(nir_builder *b, nir_ssa_def *src)
 {
    /* If fract(src) == 0.5, then we will have to decide the rounding direction.
     * We will do this by computing the mod(abs(src), 2) and testing if it
     * is < 1 or not.
     *
     * We compute mod(abs(src), 2) as:
     * abs(src) - 2.0 * floor(abs(src) / 2.0)
     */
    nir_ssa_def *two = nir_imm_double(b, 2.0);
    nir_ssa_def *abs_src = nir_fabs(b, src);
    nir_ssa_def *mod =
       nir_fsub(b,
                abs_src,
                nir_fmul(b,
                         two,
                         nir_ffloor(b,
                                    nir_fmul(b,
                                             abs_src,
                                             nir_imm_double(b, 0.5)))));

    /*
     * If fract(src) != 0.5, then we round as floor(src + 0.5)
     *
     * If fract(src) == 0.5, then we have to check the modulo:
     *
     *   if it is < 1 we need a trunc operation so we get:
     *      0.5 -> 0,   -0.5 -> -0
     *      2.5 -> 2,   -2.5 -> -2
     *
     *   otherwise we need to check if src >= 0, in which case we need to round
     *   upwards, or not, in which case we need to round downwards so we get:
     *      1.5 -> 2,   -1.5 -> -2
     *      3.5 -> 4,   -3.5 -> -4
     */
    nir_ssa_def *fract = nir_ffract(b, src);
    return nir_bcsel(b,
                     nir_fne(b, fract, nir_imm_double(b, 0.5)),
                     nir_ffloor(b, nir_fadd(b, src, nir_imm_double(b, 0.5))),
                     nir_bcsel(b,
                               nir_flt(b, mod, nir_imm_double(b, 1.0)),
                               nir_ftrunc(b, src),
                               nir_bcsel(b,
                                         nir_fge(b, src, nir_imm_double(b, 0.0)),
                                         nir_fadd(b, src, nir_imm_double(b, 0.5)),
                                         nir_fsub(b, src, nir_imm_double(b, 0.5)))));
 }

 static nir_ssa_def *
 lower_mod(nir_builder *b, nir_ssa_def *src0, nir_ssa_def *src1)
 {
    /* mod(x,y) = x - y * floor(x/y)
     *
     * If the division is lowered, it could add some rounding errors that make
     * floor() to return the quotient minus one when x = N * y. If this is the
     * case, we return zero because mod(x, y) output value is [0, y).
     */
    nir_ssa_def *floor = nir_ffloor(b, nir_fdiv(b, src0, src1));
    nir_ssa_def *mod = nir_fsub(b, src0, nir_fmul(b, src1, floor));

    return nir_bcsel(b,
                     nir_fne(b, mod, src1),
                     mod,
                     nir_imm_double(b, 0.0));
 }

 static bool
 lower_doubles_instr(nir_alu_instr *instr, nir_lower_doubles_options options)
 {
    assert(instr->dest.dest.is_ssa);
    if (instr->dest.dest.ssa.bit_size != 64)
       return false;

    switch (instr->op) {
    case nir_op_frcp:
       if (!(options & nir_lower_drcp))
          return false;
       break;

    case nir_op_fsqrt:
       if (!(options & nir_lower_dsqrt))
          return false;
       break;

    case nir_op_frsq:
       if (!(options & nir_lower_drsq))
          return false;
       break;

    case nir_op_ftrunc:
       if (!(options & nir_lower_dtrunc))
          return false;
       break;

    case nir_op_ffloor:
       if (!(options & nir_lower_dfloor))
          return false;
       break;

    case nir_op_fceil:
       if (!(options & nir_lower_dceil))
          return false;
       break;

    case nir_op_ffract:
       if (!(options & nir_lower_dfract))
          return false;
       break;

    case nir_op_fround_even:
       if (!(options & nir_lower_dround_even))
          return false;
       break;

    case nir_op_fmod:
       if (!(options & nir_lower_dmod))
          return false;
       break;

    default:
       return false;
    }

    nir_builder bld;
    nir_builder_init(&bld, nir_cf_node_get_function(&instr->instr.block->cf_node));
    bld.cursor = nir_before_instr(&instr->instr);

    nir_ssa_def *src = nir_fmov_alu(&bld, instr->src[0],
                                    instr->dest.dest.ssa.num_components);

    nir_ssa_def *result;

    switch (instr->op) {
    case nir_op_frcp:
       result = lower_rcp(&bld, src);
       break;
    case nir_op_fsqrt:
       result = lower_sqrt_rsq(&bld, src, true);
       break;
    case nir_op_frsq:
       result = lower_sqrt_rsq(&bld, src, false);
       break;
    case nir_op_ftrunc:
       result = lower_trunc(&bld, src);
       break;
    case nir_op_ffloor:
       result = lower_floor(&bld, src);
       break;
    case nir_op_fceil:
       result = lower_ceil(&bld, src);
       break;
    case nir_op_ffract:
       result = lower_fract(&bld, src);
       break;
    case nir_op_fround_even:
       result = lower_round_even(&bld, src);
       break;

    case nir_op_fmod: {
       nir_ssa_def *src1 = nir_fmov_alu(&bld, instr->src[1],
                                       instr->dest.dest.ssa.num_components);
       result = lower_mod(&bld, src, src1);
    }
       break;
    default:
       unreachable("unhandled opcode");
    }

    nir_ssa_def_rewrite_uses(&instr->dest.dest.ssa, nir_src_for_ssa(result));
    nir_instr_remove(&instr->instr);
    return true;
 }

 static bool
 nir_lower_doubles_impl(nir_function_impl *impl,
                        nir_lower_doubles_options options)
 {
    bool progress = false;

    nir_foreach_block(block, impl) {
       nir_foreach_instr_safe(instr, block) {
          if (instr->type == nir_instr_type_alu)
             progress |= lower_doubles_instr(nir_instr_as_alu(instr),
                                             options);
       }
    }

    if (progress)
       nir_metadata_preserve(impl, nir_metadata_block_index |
                                   nir_metadata_dominance);

    return progress;
 }

 bool
 nir_lower_doubles(nir_shader *shader, nir_lower_doubles_options options)
 {
    bool progress = false;

    nir_foreach_function(function, shader) {
       if (function->impl) {
          progress |= nir_lower_doubles_impl(function->impl, options);
       }
    }

    return progress;
 }
	/*
	* Copyright © 2015 Intel Corporation
	*
	* Permission is hereby granted, free of charge, to any person obtaining a
	* copy of this software and associated documentation files (the "Software"),
	* to deal in the Software without restriction, including without limitation
	* the rights to use, copy, modify, merge, publish, distribute, sublicense,
	* and/or sell copies of the Software, and to permit persons to whom the
	* Software is furnished to do so, subject to the following conditions:
	*
	* The above copyright notice and this permission notice (including the next
	* paragraph) shall be included in all copies or substantial portions of the
	* Software.
	*
	* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
	* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
	* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
	* THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
	* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
	* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
	* IN THE SOFTWARE.
	*
	*/

	#include "nir.h"
	#include "nir_builder.h"
	#include "c99_math.h"

	/*
	* Lowers some unsupported double operations, using only:
	*
	* - pack/unpackDouble2x32
	* - conversion to/from single-precision
	* - double add, mul, and fma
	* - conditional select
	* - 32-bit integer and floating point arithmetic
	*/

	/* Creates a double with the exponent bits set to a given integer value */
	static nir_ssa_def *
	set_exponent(nir_builder b, nir_ssa_def src, nir_ssa_def *exp)
	{
	/* Split into bits 0-31 and 32-63 */
	nir_ssa_def *lo = nir_unpack_64_2x32_split_x(b, src);
	nir_ssa_def *hi = nir_unpack_64_2x32_split_y(b, src);

	/* The exponent is bits 52-62, or 20-30 of the high word, so set the exponent
	* to 1023
	*/
	nir_ssa_def *new_hi = nir_bfi(b, nir_imm_int(b, 0x7ff00000), exp, hi);
	/* recombine */
	return nir_pack_64_2x32_split(b, lo, new_hi);
	}

	static nir_ssa_def *
	get_exponent(nir_builder b, nir_ssa_def src)
	{
	/* get bits 32-63 */
	nir_ssa_def *hi = nir_unpack_64_2x32_split_y(b, src);

	/* extract bits 20-30 of the high word */
	return nir_ubitfield_extract(b, hi, nir_imm_int(b, 20), nir_imm_int(b, 11));
	}

	/* Return infinity with the sign of the given source which is +/-0 */

	static nir_ssa_def *
	get_signed_inf(nir_builder b, nir_ssa_def zero)
	{
	nir_ssa_def *zero_hi = nir_unpack_64_2x32_split_y(b, zero);

	/* The bit pattern for infinity is 0x7ff0000000000000, where the sign bit
	* is the highest bit. Only the sign bit can be non-zero in the passed in
	* source. So we essentially need to OR the infinity and the zero, except
	* the low 32 bits are always 0 so we can construct the correct high 32
	* bits and then pack it together with zero low 32 bits.
	*/
	nir_ssa_def *inf_hi = nir_ior(b, nir_imm_int(b, 0x7ff00000), zero_hi);
	return nir_pack_64_2x32_split(b, nir_imm_int(b, 0), inf_hi);
	}

	/*
	* Generates the correctly-signed infinity if the source was zero, and flushes
	* the result to 0 if the source was infinity or the calculated exponent was
	* too small to be representable.
	*/

	static nir_ssa_def *
	fix_inv_result(nir_builder b, nir_ssa_def res, nir_ssa_def *src,
	nir_ssa_def *exp)
	{
	/* If the exponent is too small or the original input was infinity/NaN,
	* force the result to 0 (flush denorms) to avoid the work of handling
	* denorms properly. Note that this doesn't preserve positive/negative
	* zeros, but GLSL doesn't require it.
	*/
	res = nir_bcsel(b, nir_ior(b, nir_ige(b, nir_imm_int(b, 0), exp),
	nir_feq(b, nir_fabs(b, src),
	nir_imm_double(b, INFINITY))),
	nir_imm_double(b, 0.0f), res);

	/* If the original input was 0, generate the correctly-signed infinity */
	res = nir_bcsel(b, nir_fne(b, src, nir_imm_double(b, 0.0f)),
	res, get_signed_inf(b, src));

	return res;

	}

	static nir_ssa_def *
	lower_rcp(nir_builder b, nir_ssa_def src)
	{
	/* normalize the input to avoid range issues */
	nir_ssa_def *src_norm = set_exponent(b, src, nir_imm_int(b, 1023));

	/* cast to float, do an rcp, and then cast back to get an approximate
	* result
	*/
	nir_ssa_def *ra = nir_f2f64(b, nir_frcp(b, nir_f2f32(b, src_norm)));

	/* Fixup the exponent of the result - note that we check if this is too
	* small below.
	*/
	nir_ssa_def *new_exp = nir_isub(b, get_exponent(b, ra),
	nir_isub(b, get_exponent(b, src),
	nir_imm_int(b, 1023)));

	ra = set_exponent(b, ra, new_exp);

	/* Do a few Newton-Raphson steps to improve precision.
	*
	* Each step doubles the precision, and we started off with around 24 bits,
	* so we only need to do 2 steps to get to full precision. The step is:
	*
	* x_new = x * (2 - x*src)
	*
	* But we can re-arrange this to improve precision by using another fused
	* multiply-add:
	*
	* x_new = x + x * (1 - x*src)
	*
	* See https://en.wikipedia.org/wiki/Division_algorithm for more details.
	*/

	ra = nir_ffma(b, ra, nir_ffma(b, ra, src, nir_imm_double(b, -1)), ra);
	ra = nir_ffma(b, ra, nir_ffma(b, ra, src, nir_imm_double(b, -1)), ra);

	return fix_inv_result(b, ra, src, new_exp);
	}

	static nir_ssa_def *
	lower_sqrt_rsq(nir_builder b, nir_ssa_def src, bool sqrt)
	{
	/* We want to compute:
	*
	* 1/sqrt(m * 2^e)
	*
	* When the exponent is even, this is equivalent to:
	*
	* 1/sqrt(m) * 2^(-e/2)
	*
	* and then the exponent is odd, this is equal to:
	*
	* 1/sqrt(m * 2) * 2^(-(e - 1)/2)
	*
	* where the m * 2 is absorbed into the exponent. So we want the exponent
	* inside the square root to be 1 if e is odd and 0 if e is even, and we
	* want to subtract off e/2 from the final exponent, rounded to negative
	* infinity. We can do the former by first computing the unbiased exponent,
	* and then AND'ing it with 1 to get 0 or 1, and we can do the latter by
	* shifting right by 1.
	*/

	nir_ssa_def *unbiased_exp = nir_isub(b, get_exponent(b, src),
	nir_imm_int(b, 1023));
	nir_ssa_def *even = nir_iand(b, unbiased_exp, nir_imm_int(b, 1));
	nir_ssa_def *half = nir_ishr(b, unbiased_exp, nir_imm_int(b, 1));

	nir_ssa_def *src_norm = set_exponent(b, src,
	nir_iadd(b, nir_imm_int(b, 1023),
	even));

	nir_ssa_def *ra = nir_f2f64(b, nir_frsq(b, nir_f2f32(b, src_norm)));
	nir_ssa_def *new_exp = nir_isub(b, get_exponent(b, ra), half);
	ra = set_exponent(b, ra, new_exp);

	/*
	* The following implements an iterative algorithm that's very similar
	* between sqrt and rsqrt. We start with an iteration of Goldschmit's
	* algorithm, which looks like:
	*
	* a = the source
	* y_0 = initial (single-precision) rsqrt estimate
	*
	* h_0 = .5 * y_0
	* g_0 = a * y_0
	* r_0 = .5 - h_0 * g_0
	* g_1 = g_0 * r_0 + g_0
	* h_1 = h_0 * r_0 + h_0
	*
	* Now g_1 ~= sqrt(a), and h_1 ~= 1/(2 * sqrt(a)). We could continue
	* applying another round of Goldschmit, but since we would never refer
	* back to a (the original source), we would add too much rounding error.
	* So instead, we do one last round of Newton-Raphson, which has better
	* rounding characteristics, to get the final rounding correct. This is
	* split into two cases:
	*
	* 1. sqrt
	*
	* Normally, doing a round of Newton-Raphson for sqrt involves taking a
	* reciprocal of the original estimate, which is slow since it isn't
	* supported in HW. But we can take advantage of the fact that we already
	* computed a good estimate of 1/(2 * g_1) by rearranging it like so:
	*
	* g_2 = .5 * (g_1 + a / g_1)
	* = g_1 + .5 * (a / g_1 - g_1)
	* = g_1 + (.5 / g_1) * (a - g_1^2)
	* = g_1 + h_1 * (a - g_1^2)
	*
	* The second term represents the error, and by splitting it out we can get
	* better precision by computing it as part of a fused multiply-add. Since
	* both Newton-Raphson and Goldschmit approximately double the precision of
	* the result, these two steps should be enough.
	*
	* 2. rsqrt
	*
	* First off, note that the first round of the Goldschmit algorithm is
	* really just a Newton-Raphson step in disguise:
	*
	* h_1 = h_0 * (.5 - h_0 * g_0) + h_0
	* = h_0 * (1.5 - h_0 * g_0)
	* = h_0 * (1.5 - .5 * a * y_0^2)
	* = (.5 * y_0) * (1.5 - .5 * a * y_0^2)
	*
	* which is the standard formula multiplied by .5. Unlike in the sqrt case,
	* we don't need the inverse to do a Newton-Raphson step; we just need h_1,
	* so we can skip the calculation of g_1. Instead, we simply do another
	* Newton-Raphson step:
	*
	* y_1 = 2 * h_1
	* r_1 = .5 - h_1 * y_1 * a
	* y_2 = y_1 * r_1 + y_1
	*
	* Where the difference from Goldschmit is that we calculate y_1 * a
	* instead of using g_1. Doing it this way should be as fast as computing
	* y_1 up front instead of h_1, and it lets us share the code for the
	* initial Goldschmit step with the sqrt case.
	*
	* Putting it together, the computations are:
	*
	* h_0 = .5 * y_0
	* g_0 = a * y_0
	* r_0 = .5 - h_0 * g_0
	* h_1 = h_0 * r_0 + h_0
	* if sqrt:
	* g_1 = g_0 * r_0 + g_0
	* r_1 = a - g_1 * g_1
	* g_2 = h_1 * r_1 + g_1
	* else:
	* y_1 = 2 * h_1
	* r_1 = .5 - y_1 * (h_1 * a)
	* y_2 = y_1 * r_1 + y_1
	*
	* For more on the ideas behind this, see "Software Division and Square
	* Root Using Goldschmit's Algorithms" by Markstein and the Wikipedia page
	* on square roots
	* (https://en.wikipedia.org/wiki/Methods_of_computing_square_roots).
	*/

	nir_ssa_def *one_half = nir_imm_double(b, 0.5);
	nir_ssa_def *h_0 = nir_fmul(b, one_half, ra);
	nir_ssa_def *g_0 = nir_fmul(b, src, ra);
	nir_ssa_def *r_0 = nir_ffma(b, nir_fneg(b, h_0), g_0, one_half);
	nir_ssa_def *h_1 = nir_ffma(b, h_0, r_0, h_0);
	nir_ssa_def *res;
	if (sqrt) {
	nir_ssa_def *g_1 = nir_ffma(b, g_0, r_0, g_0);
	nir_ssa_def *r_1 = nir_ffma(b, nir_fneg(b, g_1), g_1, src);
	res = nir_ffma(b, h_1, r_1, g_1);
	} else {
	nir_ssa_def *y_1 = nir_fmul(b, nir_imm_double(b, 2.0), h_1);
	nir_ssa_def *r_1 = nir_ffma(b, nir_fneg(b, y_1), nir_fmul(b, h_1, src),
	one_half);
	res = nir_ffma(b, y_1, r_1, y_1);
	}

	if (sqrt) {
	/* Here, the special cases we need to handle are
	* 0 -> 0 and
	* +inf -> +inf
	*/
	res = nir_bcsel(b, nir_ior(b, nir_feq(b, src, nir_imm_double(b, 0.0)),
	nir_feq(b, src, nir_imm_double(b, INFINITY))),
	src, res);
	} else {
	res = fix_inv_result(b, res, src, new_exp);
	}

	return res;
	}

	static nir_ssa_def *
	lower_trunc(nir_builder b, nir_ssa_def src)
	{
	nir_ssa_def *unbiased_exp = nir_isub(b, get_exponent(b, src),
	nir_imm_int(b, 1023));

	nir_ssa_def *frac_bits = nir_isub(b, nir_imm_int(b, 52), unbiased_exp);

	/*
	* Decide the operation to apply depending on the unbiased exponent:
	*
	* if (unbiased_exp < 0)
	* return 0
	* else if (unbiased_exp > 52)
	* return src
	* else
	* return src & (~0 << frac_bits)
	*
	* Notice that the else branch is a 64-bit integer operation that we need
	* to implement in terms of 32-bit integer arithmetics (at least until we
	* support 64-bit integer arithmetics).
	*/

	/* Compute "~0 << frac_bits" in terms of hi/lo 32-bit integer math */
	nir_ssa_def *mask_lo =
	nir_bcsel(b,
	nir_ige(b, frac_bits, nir_imm_int(b, 32)),
	nir_imm_int(b, 0),
	nir_ishl(b, nir_imm_int(b, ~0), frac_bits));

	nir_ssa_def *mask_hi =
	nir_bcsel(b,
	nir_ilt(b, frac_bits, nir_imm_int(b, 33)),
	nir_imm_int(b, ~0),
	nir_ishl(b,
	nir_imm_int(b, ~0),
	nir_isub(b, frac_bits, nir_imm_int(b, 32))));

	nir_ssa_def *src_lo = nir_unpack_64_2x32_split_x(b, src);
	nir_ssa_def *src_hi = nir_unpack_64_2x32_split_y(b, src);

	return
	nir_bcsel(b,
	nir_ilt(b, unbiased_exp, nir_imm_int(b, 0)),
	nir_imm_double(b, 0.0),
	nir_bcsel(b, nir_ige(b, unbiased_exp, nir_imm_int(b, 53)),
	src,
	nir_pack_64_2x32_split(b,
	nir_iand(b, mask_lo, src_lo),
	nir_iand(b, mask_hi, src_hi))));
	}

	static nir_ssa_def *
	lower_floor(nir_builder b, nir_ssa_def src)
	{
	/*
	* For x >= 0, floor(x) = trunc(x)
	* For x < 0,
	* - if x is integer, floor(x) = x
	* - otherwise, floor(x) = trunc(x) - 1
	*/
	nir_ssa_def *tr = nir_ftrunc(b, src);
	nir_ssa_def *positive = nir_fge(b, src, nir_imm_double(b, 0.0));
	return nir_bcsel(b,
	nir_ior(b, positive, nir_feq(b, src, tr)),
	tr,
	nir_fsub(b, tr, nir_imm_double(b, 1.0)));
	}

	static nir_ssa_def *
	lower_ceil(nir_builder b, nir_ssa_def src)
	{
	/* if x < 0, ceil(x) = trunc(x)
	* else if (x - trunc(x) == 0), ceil(x) = x
	* else, ceil(x) = trunc(x) + 1
	*/
	nir_ssa_def *tr = nir_ftrunc(b, src);
	nir_ssa_def *negative = nir_flt(b, src, nir_imm_double(b, 0.0));
	return nir_bcsel(b,
	nir_ior(b, negative, nir_feq(b, src, tr)),
	tr,
	nir_fadd(b, tr, nir_imm_double(b, 1.0)));
	}

	static nir_ssa_def *
	lower_fract(nir_builder b, nir_ssa_def src)
	{
	return nir_fsub(b, src, nir_ffloor(b, src));
	}

	static nir_ssa_def *
	lower_round_even(nir_builder b, nir_ssa_def src)
	{
	/* If fract(src) == 0.5, then we will have to decide the rounding direction.
	* We will do this by computing the mod(abs(src), 2) and testing if it
	* is < 1 or not.
	*
	* We compute mod(abs(src), 2) as:
	* abs(src) - 2.0 * floor(abs(src) / 2.0)
	*/
	nir_ssa_def *two = nir_imm_double(b, 2.0);
	nir_ssa_def *abs_src = nir_fabs(b, src);
	nir_ssa_def *mod =
	nir_fsub(b,
	abs_src,
	nir_fmul(b,
	two,
	nir_ffloor(b,
	nir_fmul(b,
	abs_src,
	nir_imm_double(b, 0.5)))));

	/*
	* If fract(src) != 0.5, then we round as floor(src + 0.5)
	*
	* If fract(src) == 0.5, then we have to check the modulo:
	*
	* if it is < 1 we need a trunc operation so we get:
	* 0.5 -> 0, -0.5 -> -0
	* 2.5 -> 2, -2.5 -> -2
	*
	* otherwise we need to check if src >= 0, in which case we need to round
	* upwards, or not, in which case we need to round downwards so we get:
	* 1.5 -> 2, -1.5 -> -2
	* 3.5 -> 4, -3.5 -> -4
	*/
	nir_ssa_def *fract = nir_ffract(b, src);
	return nir_bcsel(b,
	nir_fne(b, fract, nir_imm_double(b, 0.5)),
	nir_ffloor(b, nir_fadd(b, src, nir_imm_double(b, 0.5))),
	nir_bcsel(b,
	nir_flt(b, mod, nir_imm_double(b, 1.0)),
	nir_ftrunc(b, src),
	nir_bcsel(b,
	nir_fge(b, src, nir_imm_double(b, 0.0)),
	nir_fadd(b, src, nir_imm_double(b, 0.5)),
	nir_fsub(b, src, nir_imm_double(b, 0.5)))));
	}

	static nir_ssa_def *
	lower_mod(nir_builder b, nir_ssa_def src0, nir_ssa_def *src1)
	{
	/* mod(x,y) = x - y * floor(x/y)
	*
	* If the division is lowered, it could add some rounding errors that make
	* floor() to return the quotient minus one when x = N * y. If this is the
	* case, we return zero because mod(x, y) output value is [0, y).
	*/
	nir_ssa_def *floor = nir_ffloor(b, nir_fdiv(b, src0, src1));
	nir_ssa_def *mod = nir_fsub(b, src0, nir_fmul(b, src1, floor));

	return nir_bcsel(b,
	nir_fne(b, mod, src1),
	mod,
	nir_imm_double(b, 0.0));
	}

	static bool
	lower_doubles_instr(nir_alu_instr *instr, nir_lower_doubles_options options)
	{
	assert(instr->dest.dest.is_ssa);
	if (instr->dest.dest.ssa.bit_size != 64)
	return false;

	switch (instr->op) {
	case nir_op_frcp:
	if (!(options & nir_lower_drcp))
	return false;
	break;

	case nir_op_fsqrt:
	if (!(options & nir_lower_dsqrt))
	return false;
	break;

	case nir_op_frsq:
	if (!(options & nir_lower_drsq))
	return false;
	break;

	case nir_op_ftrunc:
	if (!(options & nir_lower_dtrunc))
	return false;
	break;

	case nir_op_ffloor:
	if (!(options & nir_lower_dfloor))
	return false;
	break;

	case nir_op_fceil:
	if (!(options & nir_lower_dceil))
	return false;
	break;

	case nir_op_ffract:
	if (!(options & nir_lower_dfract))
	return false;
	break;

	case nir_op_fround_even:
	if (!(options & nir_lower_dround_even))
	return false;
	break;

	case nir_op_fmod:
	if (!(options & nir_lower_dmod))
	return false;
	break;

	default:
	return false;
	}

	nir_builder bld;
	nir_builder_init(&bld, nir_cf_node_get_function(&instr->instr.block->cf_node));
	bld.cursor = nir_before_instr(&instr->instr);

	nir_ssa_def *src = nir_fmov_alu(&bld, instr->src[0],
	instr->dest.dest.ssa.num_components);

	nir_ssa_def *result;

	switch (instr->op) {
	case nir_op_frcp:
	result = lower_rcp(&bld, src);
	break;
	case nir_op_fsqrt:
	result = lower_sqrt_rsq(&bld, src, true);
	break;
	case nir_op_frsq:
	result = lower_sqrt_rsq(&bld, src, false);
	break;
	case nir_op_ftrunc:
	result = lower_trunc(&bld, src);
	break;
	case nir_op_ffloor:
	result = lower_floor(&bld, src);
	break;
	case nir_op_fceil:
	result = lower_ceil(&bld, src);
	break;
	case nir_op_ffract:
	result = lower_fract(&bld, src);
	break;
	case nir_op_fround_even:
	result = lower_round_even(&bld, src);
	break;

	case nir_op_fmod: {
	nir_ssa_def *src1 = nir_fmov_alu(&bld, instr->src[1],
	instr->dest.dest.ssa.num_components);
	result = lower_mod(&bld, src, src1);
	}
	break;
	default:
	unreachable("unhandled opcode");
	}

	nir_ssa_def_rewrite_uses(&instr->dest.dest.ssa, nir_src_for_ssa(result));
	nir_instr_remove(&instr->instr);
	return true;
	}

	static bool
	nir_lower_doubles_impl(nir_function_impl *impl,
	nir_lower_doubles_options options)
	{
	bool progress = false;

	nir_foreach_block(block, impl) {
	nir_foreach_instr_safe(instr, block) {
	if (instr->type == nir_instr_type_alu)
	progress \|= lower_doubles_instr(nir_instr_as_alu(instr),
	options);
	}
	}

	if (progress)
	nir_metadata_preserve(impl, nir_metadata_block_index \|
	nir_metadata_dominance);

	return progress;
	}

	bool
	nir_lower_doubles(nir_shader *shader, nir_lower_doubles_options options)
	{
	bool progress = false;

	nir_foreach_function(function, shader) {
	if (function->impl) {
	progress \|= nir_lower_doubles_impl(function->impl, options);
	}
	}

	return progress;
	}