amdgcn: Improve FP division accuracy

GCN has a reciprocal-approximation instruction but no
hardware divide. This patch adjusts the open-coded reciprocal
approximation/Newton-Raphson refinement steps to use fused multiply-add
instructions as is necessary to obtain a properly-rounded result, and
adds further refinement steps to correctly round the full division result.

The patterns in question are still guarded by a flag_reciprocal_math
condition, and do not yet support denormals.

2021-01-13  Julian Brown  <julian@codesourcery.com>

gcc/
	* config/gcn/gcn-valu.md (recip<mode>2<exec>, recip<mode>2): Use unspec
	for reciprocal-approximation instructions.
	(div<mode>3): Use fused multiply-accumulate operations for reciprocal
	refinement and division result.
	* config/gcn/gcn.md (UNSPEC_RCP): New unspec constant.

gcc/testsuite/
	* gcc.target/gcn/fpdiv.c: New test.

3 files changed, 81 insertions, 20 deletions

diff --git a/gcc/config/gcn/gcn-valu.md b/gcc/config/gcn/gcn-valu.md
index 24c1721..beefcf7 100644
--- a/gcc/config/gcn/gcn-valu.md
+++ b/gcc/config/gcn/gcn-valu.md
@@ -2351,9 +2351,9 @@
 
 (define_insn "recip<mode>2<exec>"
   [(set (match_operand:V_FP 0 "register_operand"  "=  v")
-	(div:V_FP
-	  (vec_duplicate:V_FP (float:<SCALAR_MODE> (const_int 1)))
-	  (match_operand:V_FP 1 "gcn_alu_operand" "vSvB")))]
+	(unspec:V_FP
+	  [(match_operand:V_FP 1 "gcn_alu_operand" "vSvB")]
+	  UNSPEC_RCP))]
   ""
   "v_rcp%i0\t%0, %1"
   [(set_attr "type" "vop1")
@@ -2361,9 +2361,9 @@
 
 (define_insn "recip<mode>2"
   [(set (match_operand:FP 0 "register_operand"	 "=  v")
-	(div:FP
-	  (float:FP (const_int 1))
-	  (match_operand:FP 1 "gcn_alu_operand"	 "vSvB")))]
+	(unspec:FP
+	  [(match_operand:FP 1 "gcn_alu_operand" "vSvB")]
+	  UNSPEC_RCP))]
   ""
   "v_rcp%i0\t%0, %1"
   [(set_attr "type" "vop1")
@@ -2382,28 +2382,39 @@
    (match_operand:V_FP 2 "gcn_valu_src0_operand")]
   "flag_reciprocal_math"
   {
-    rtx two = gcn_vec_constant (<MODE>mode,
-		  const_double_from_real_value (dconst2, <SCALAR_MODE>mode));
+    rtx one = gcn_vec_constant (<MODE>mode,
+		  const_double_from_real_value (dconst1, <SCALAR_MODE>mode));
     rtx initrcp = gen_reg_rtx (<MODE>mode);
     rtx fma = gen_reg_rtx (<MODE>mode);
     rtx rcp;
+    rtx num = operands[1], denom = operands[2];
 
-    bool is_rcp = (GET_CODE (operands[1]) == CONST_VECTOR
+    bool is_rcp = (GET_CODE (num) == CONST_VECTOR
 		   && real_identical
 		        (CONST_DOUBLE_REAL_VALUE
-			  (CONST_VECTOR_ELT (operands[1], 0)), &dconstm1));
+			  (CONST_VECTOR_ELT (num, 0)), &dconstm1));
 
     if (is_rcp)
       rcp = operands[0];
     else
       rcp = gen_reg_rtx (<MODE>mode);
 
-    emit_insn (gen_recip<mode>2 (initrcp, operands[2]));
-    emit_insn (gen_fma<mode>4_negop2 (fma, initrcp, operands[2], two));
-    emit_insn (gen_mul<mode>3 (rcp, initrcp, fma));
+    emit_insn (gen_recip<mode>2 (initrcp, denom));
+    emit_insn (gen_fma<mode>4_negop2 (fma, initrcp, denom, one));
+    emit_insn (gen_fma<mode>4 (rcp, fma, initrcp, initrcp));
 
     if (!is_rcp)
-      emit_insn (gen_mul<mode>3 (operands[0], operands[1], rcp));
+      {
+	rtx div_est = gen_reg_rtx (<MODE>mode);
+	rtx fma2 = gen_reg_rtx (<MODE>mode);
+	rtx fma3 = gen_reg_rtx (<MODE>mode);
+	rtx fma4 = gen_reg_rtx (<MODE>mode);
+	emit_insn (gen_mul<mode>3 (div_est, num, rcp));
+	emit_insn (gen_fma<mode>4_negop2 (fma2, div_est, denom, num));
+	emit_insn (gen_fma<mode>4 (fma3, fma2, rcp, div_est));
+	emit_insn (gen_fma<mode>4_negop2 (fma4, fma3, denom, num));
+	emit_insn (gen_fma<mode>4 (operands[0], fma4, rcp, fma3));
+      }
 
     DONE;
   })
@@ -2414,10 +2425,11 @@
    (match_operand:FP 2 "gcn_valu_src0_operand")]
   "flag_reciprocal_math"
   {
-    rtx two = const_double_from_real_value (dconst2, <MODE>mode);
+    rtx one = const_double_from_real_value (dconst1, <MODE>mode);
     rtx initrcp = gen_reg_rtx (<MODE>mode);
     rtx fma = gen_reg_rtx (<MODE>mode);
     rtx rcp;
+    rtx num = operands[1], denom = operands[2];
 
     bool is_rcp = (GET_CODE (operands[1]) == CONST_DOUBLE
 		   && real_identical (CONST_DOUBLE_REAL_VALUE (operands[1]),
@@ -2428,12 +2440,22 @@
     else
       rcp = gen_reg_rtx (<MODE>mode);
 
-    emit_insn (gen_recip<mode>2 (initrcp, operands[2]));
-    emit_insn (gen_fma<mode>4_negop2 (fma, initrcp, operands[2], two));
-    emit_insn (gen_mul<mode>3 (rcp, initrcp, fma));
+    emit_insn (gen_recip<mode>2 (initrcp, denom));
+    emit_insn (gen_fma<mode>4_negop2 (fma, initrcp, denom, one));
+    emit_insn (gen_fma<mode>4 (rcp, fma, initrcp, initrcp));
 
     if (!is_rcp)
-      emit_insn (gen_mul<mode>3 (operands[0], operands[1], rcp));
+      {
+	rtx div_est = gen_reg_rtx (<MODE>mode);
+	rtx fma2 = gen_reg_rtx (<MODE>mode);
+	rtx fma3 = gen_reg_rtx (<MODE>mode);
+	rtx fma4 = gen_reg_rtx (<MODE>mode);
+	emit_insn (gen_mul<mode>3 (div_est, num, rcp));
+	emit_insn (gen_fma<mode>4_negop2 (fma2, div_est, denom, num));
+	emit_insn (gen_fma<mode>4 (fma3, fma2, rcp, div_est));
+	emit_insn (gen_fma<mode>4_negop2 (fma4, fma3, denom, num));
+	emit_insn (gen_fma<mode>4 (operands[0], fma4, rcp, fma3));
+      }
 
     DONE;
   })
diff --git a/gcc/config/gcn/gcn.md b/gcc/config/gcn/gcn.md
index 757e575..b5f895a 100644
--- a/gcc/config/gcn/gcn.md
+++ b/gcc/config/gcn/gcn.md
@@ -80,7 +80,8 @@
   UNSPEC_MOV_DPP_SHR
   UNSPEC_MOV_FROM_LANE63
   UNSPEC_GATHER
-  UNSPEC_SCATTER])
+  UNSPEC_SCATTER
+  UNSPEC_RCP])
 
 ;; }}}
 ;; {{{ Attributes
diff --git a/gcc/testsuite/gcc.target/gcn/fpdiv.c b/gcc/testsuite/gcc.target/gcn/fpdiv.c
new file mode 100644
index 0000000..7125b6f
--- /dev/null
+++ b/gcc/testsuite/gcc.target/gcn/fpdiv.c
@@ -0,0 +1,38 @@
+/* { dg-do run } */
+/* { dg-options "-ffast-math" } */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <stdint.h>
+
+double num = 100;
+double denom = 50;
+
+union bits {
+  double fp;
+  uint64_t ull;
+};
+
+int main (int argc, char* argv[])
+{
+  union bits recip;
+  union bits res;
+
+  recip.fp = 1.0 / denom;
+
+  if (recip.ull != 0x3f947ae147ae147b)
+    {
+      fprintf (stderr, "incorrectly-rounded reciprocal: %llx", recip.ull);
+      exit (1);
+    }
+
+  res.fp = num / denom;
+
+  if (res.ull != 0x4000000000000000)
+    {
+      fprintf (stderr, "incorrectly-rounded quotient: %llx", res.ull);
+      exit (1);
+    }
+
+  return 0;
+}