title: How to refine floating-point division on FMA-capable GPUs?
body: <p>When writing computational code for GPUs using APIs where compute shaders are translated via SPIR-V (in particular, Vulkan), I am guaranteed that ULP error of floating-point division will be at most 3. Other basic arithmetic (addition, multiplication) is correctly rounded.</p>
  <p>How to efficiently implement correctly rounded division in these circumstances? Let's assume that FMA is available and correctly rounded.</p>
  <p>What happens with denormals will depend on the underlying hardware. Vulkan API allows to query whether the device can preserve denormals and whether it can flush them to zero (so a GPU that does not fully support denormals will have &quot;canPreserve: false, canFlush: true&quot;). So let's additionally assume that the GPU can produce and handle denormals, without flushing them to zero (otherwise trying to produce a correctly rounded result which is subnormal seems futile).</p>