[PATCH 3/4] math: Make rsqrtf correctly rounded

[Adhemerval Zanella]

The generic implementation computes 1.0f / sqrtf (x), which rounds twice
and shows up to 1 ulp for rounding to nearest (for about 26% of the
inputs) and up to 2 ulp for the directed roundings.

Using double precision is a straightforward way to get the correctly
rounded resuls, and it was exhaustively tested for all possible inputs
on different machines (x86_64, aarch64, powerpc, loongarch64).

The performacne results shows a performance hit, since to use the
double version requires extra single to double rounding.  Benchtests on
x86_64 (Ryzen 9 5900X), aarch64 (Neoverse-N1), powerpc (POWER10, built
for power8), and loongarch64 (Loongson-3C5000L-LL) shows:

* latency
- input: [0.5,2)
                   master       patched   improvement
x86_64            23.8349       35.2275       -32.34%
x86_64-v2         24.0356       35.2532       -31.82%
x86_64-v3         22.7667       34.0389       -33.12%
aarch64            7.7400       13.4675       -42.53%
powerpc            7.1457        9.7876       -26.99%
loongarch64        2.1518        3.4530       -37.68%

- input: normal
                   master       patched   improvement
x86_64            23.6777       35.1879       -32.71%
x86_64-v2         24.0492       35.7858       -32.80%
x86_64-v3         22.7498       33.9989       -33.09%
aarch64            7.7401       13.4679       -42.53%
powerpc            7.1465        9.7657       -26.82%
loongarch64        2.1519        3.4527       -37.68%

- input: subnormal
                   master       patched   improvement
x86_64            35.6423       35.1984         1.26%
x86_64-v2         35.5350       35.4863         0.14%
x86_64-v3         32.1197       33.9955        -5.52%
aarch64            8.0830       13.4669       -39.98%
powerpc            7.1431        9.7794       -26.96%
loongarch64        2.2020        3.4524       -36.22%

* reciprocal-throughput
- input: [0.5,2)
                   master       patched   improvement
x86_64             8.4902       10.8043       -21.42%
x86_64-v2          9.3151       10.5439       -11.65%
x86_64-v3          8.4450       10.8503       -22.17%
aarch64            3.8338        5.8337       -34.28%
powerpc            1.7494        2.4880       -29.69%
loongarch64        1.7027        1.8515        -8.04%

- input: normal
                   master       patched   improvement
x86_64             8.4228       10.7892       -21.93%
x86_64-v2          9.2843       10.6977       -13.21%
x86_64-v3          8.4090       10.8325       -22.37%
aarch64            3.8345        5.8413       -34.36%
powerpc            1.7454        2.4812       -29.66%
loongarch64        1.7026        1.8514        -8.04%

- input: subnormal
                   master       patched   improvement
x86_64            15.0495       10.7866        39.52%
x86_64-v2         14.2217       10.6250        33.85%
x86_64-v3         14.2447       10.8250        31.59%
aarch64            3.9800        5.8460       -31.92%
powerpc            1.7488        2.4808       -29.51%
loongarch64        1.7026        1.8514        -8.04%

Checked on aarch64-linux-gnu, x86_64-linux-gnu, powerpc64le-linux-gnu,
and loongarch64-linux-gnuf64.
---
 sysdeps/ieee754/flt-32/libm-test-ulps | 12 +++++++++++
 sysdeps/ieee754/flt-32/s_rsqrtf.c     | 30 +++++++++++++++++++++++++++
 2 files changed, 42 insertions(+)
 create mode 100644 sysdeps/ieee754/flt-32/s_rsqrtf.c

diff --git a/sysdeps/ieee754/flt-32/libm-test-ulps b/sysdeps/ieee754/flt-32/libm-test-ulps
index 15a7248d8b4..52a9238c828 100644
--- a/sysdeps/ieee754/flt-32/libm-test-ulps
+++ b/sysdeps/ieee754/flt-32/libm-test-ulps
@@ -251,6 +251,18 @@ float: 0
 Function: "log2p1_upward":
 float: 0
 
+Function: "rsqrt":
+float: 0
+
+Function: "rsqrt_downward":
+float: 0
+
+Function: "rsqrt_towardzero":
+float: 0
+
+Function: "rsqrt_upward":
+float: 0
+
 Function: "sinh":
 float: 0
 
diff --git a/sysdeps/ieee754/flt-32/s_rsqrtf.c b/sysdeps/ieee754/flt-32/s_rsqrtf.c
new file mode 100644
index 00000000000..112b024794d
--- /dev/null
+++ b/sysdeps/ieee754/flt-32/s_rsqrtf.c
@@ -0,0 +1,30 @@
+/* Correctly rounded reciprocal square root of binary32 value.
+   Copyright (C) 2026 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+
+   The GNU C Library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Lesser General Public
+   License as published by the Free Software Foundation; either
+   version 2.1 of the License, or (at your option) any later version.
+
+   The GNU C Library is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+   Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public
+   License along with the GNU C Library; if not, see
+   <https://www.gnu.org/licenses/>.  */
+
+#include <errno.h>
+#include <math.h>
+#include <libm-alias-float.h>
+
+float
+__rsqrtf (float x)
+{
+  if (islessequal (x, 0.0f))
+    __set_errno (x < 0 ? EDOM : ERANGE);
+  return 1.0 / sqrt ((double) x);
+}
+libm_alias_float (__rsqrt, rsqrt)
-- 
2.53.0