Add narrowing square root functions

This patch adds the narrowing square root functions from TS 18661-1 /
TS 18661-3 / C2X to glibc's libm: fsqrt, fsqrtl, dsqrtl, f32sqrtf64,
f32sqrtf32x, f32xsqrtf64 for all configurations; f32sqrtf64x,
f32sqrtf128, f64sqrtf64x, f64sqrtf128, f32xsqrtf64x, f32xsqrtf128,
f64xsqrtf128 for configurations with _Float64x and _Float128;
__f32sqrtieee128 and __f64sqrtieee128 aliases in the powerpc64le case
(for calls to fsqrtl and dsqrtl when long double is IEEE binary128).
Corresponding tgmath.h macro support is also added.

The changes are mostly similar to those for the other narrowing
functions previously added, so the description of those generally
applies to this patch as well.  However, the not-actually-narrowing
cases (where the two types involved in the function have the same
floating-point format) are aliased to sqrt, sqrtl or sqrtf128 rather
than needing a separately built not-actually-narrowing function such
as was needed for add / sub / mul / div.  Thus, there is no
__nldbl_dsqrtl name for ldbl-opt because no such name was needed
(whereas the other functions needed such a name since the only other
name for that entry point was e.g. f32xaddf64, not reserved by TS
18661-1); the headers are made to arrange for sqrt to be called in
that case instead.

The DIAG_* calls in sysdeps/ieee754/soft-fp/s_dsqrtl.c are because
they were observed to be needed in GCC 7 testing of
riscv32-linux-gnu-rv32imac-ilp32.  The other sysdeps/ieee754/soft-fp/
files added didn't need such DIAG_* in any configuration I tested with
build-many-glibcs.py, but if they do turn out to be needed in more
files with some other configuration / GCC version, they can always be
added there.

I reused the same test inputs in auto-libm-test-in as for
non-narrowing sqrt rather than adding extra or separate inputs for
narrowing sqrt.  The tests in libm-test-narrow-sqrt.inc also follow
those for non-narrowing sqrt.

Tested as followed: natively with the full glibc testsuite for x86_64
(GCC 11, 7, 6) and x86 (GCC 11); with build-many-glibcs.py with GCC
11, 7 and 6; cross testing of math/ tests for powerpc64le, powerpc32
hard float, mips64 (all three ABIs, both hard and soft float).  The
different GCC versions are to cover the different cases in tgmath.h
and tgmath.h tests properly (GCC 6 has _Float* only as typedefs in
glibc headers, GCC 7 has proper _Float* support, GCC 8 adds
__builtin_tgmath).

20 files changed, 523 insertions, 1 deletions

diff --git a/sysdeps/ieee754/dbl-64/s_f32xsqrtf64.c b/sysdeps/ieee754/dbl-64/s_f32xsqrtf64.c
new file mode 100644
index 0000000000..9af76a25b0
--- /dev/null
+++ b/sysdeps/ieee754/dbl-64/s_f32xsqrtf64.c
@@ -0,0 +1 @@
+/* Defined as an alias of sqrt.  */
diff --git a/sysdeps/ieee754/dbl-64/s_fsqrt.c b/sysdeps/ieee754/dbl-64/s_fsqrt.c
new file mode 100644
index 0000000000..0354ba9fbf
--- /dev/null
+++ b/sysdeps/ieee754/dbl-64/s_fsqrt.c
@@ -0,0 +1,34 @@
+/* Square root of double value, narrowing the result to float.
+   Copyright (C) 2021 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/>.  */
+
+#define f32sqrtf64 __hide_f32sqrtf64
+#define f32sqrtf32x __hide_f32sqrtf32x
+#define fsqrtl __hide_fsqrtl
+#include <math.h>
+#undef f32sqrtf64
+#undef f32sqrtf32x
+#undef fsqrtl
+
+#include <math-narrow.h>
+
+float
+__fsqrt (double x)
+{
+  NARROW_SQRT_ROUND_TO_ODD (x, float, union ieee754_double, , mantissa1);
+}
+libm_alias_float_double (sqrt)
diff --git a/sysdeps/ieee754/float128/float128_private.h b/sysdeps/ieee754/float128/float128_private.h
index b6b6d3d5fd..13d1d63baf 100644
--- a/sysdeps/ieee754/float128/float128_private.h
+++ b/sysdeps/ieee754/float128/float128_private.h
@@ -138,6 +138,9 @@
 #define libm_alias_float_ldouble(func) libm_alias_float32_float128 (func)
 #undef libm_alias_double_ldouble
 #define libm_alias_double_ldouble(func) libm_alias_float64_float128 (func)
+#undef libm_alias_ldouble_narrow
+#define libm_alias_ldouble_narrow(from, to)	\
+  libm_alias_float128_narrow (from, to)
 
 #include <math-use-builtins.h>
 #undef USE_NEARBYINTL_BUILTIN
@@ -284,6 +287,8 @@
 #define __ddivl __f64divf128
 #define __fmull __f32mulf128
 #define __dmull __f64mulf128
+#define __fsqrtl __f32sqrtf128
+#define __dsqrtl __f64sqrtf128
 #define __fsubl __f32subf128
 #define __dsubl __f64subf128
 
diff --git a/sysdeps/ieee754/float128/s_f32sqrtf128.c b/sysdeps/ieee754/float128/s_f32sqrtf128.c
new file mode 100644
index 0000000000..74e91abd4d
--- /dev/null
+++ b/sysdeps/ieee754/float128/s_f32sqrtf128.c
@@ -0,0 +1,6 @@
+#define f32sqrtf64x __hide_f32sqrtf64x
+#define f32sqrtf128 __hide_f32sqrtf128
+#include <float128_private.h>
+#undef f32sqrtf64x
+#undef f32sqrtf128
+#include "../ldbl-128/s_fsqrtl.c"
diff --git a/sysdeps/ieee754/float128/s_f64sqrtf128.c b/sysdeps/ieee754/float128/s_f64sqrtf128.c
new file mode 100644
index 0000000000..8373879d42
--- /dev/null
+++ b/sysdeps/ieee754/float128/s_f64sqrtf128.c
@@ -0,0 +1,10 @@
+#define f32xsqrtf64x __hide_f32xsqrtf64x
+#define f32xsqrtf128 __hide_f32xsqrtf128
+#define f64sqrtf64x __hide_f64sqrtf64x
+#define f64sqrtf128 __hide_f64sqrtf128
+#include <float128_private.h>
+#undef f32xsqrtf64x
+#undef f32xsqrtf128
+#undef f64sqrtf64x
+#undef f64sqrtf128
+#include "../ldbl-128/s_dsqrtl.c"
diff --git a/sysdeps/ieee754/float128/s_f64xsqrtf128.c b/sysdeps/ieee754/float128/s_f64xsqrtf128.c
new file mode 100644
index 0000000000..383f92a94e
--- /dev/null
+++ b/sysdeps/ieee754/float128/s_f64xsqrtf128.c
@@ -0,0 +1,2 @@
+#include <float128_private.h>
+#include "../ldbl-128/s_f64xsqrtf128.c"
diff --git a/sysdeps/ieee754/ldbl-128/s_dsqrtl.c b/sysdeps/ieee754/ldbl-128/s_dsqrtl.c
new file mode 100644
index 0000000000..b0d0162467
--- /dev/null
+++ b/sysdeps/ieee754/ldbl-128/s_dsqrtl.c
@@ -0,0 +1,37 @@
+/* Square root of long double (ldbl-128) value, narrowing the result to double.
+   Copyright (C) 2021 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/>.  */
+
+#define f32xsqrtf64x __hide_f32xsqrtf64x
+#define f32xsqrtf128 __hide_f32xsqrtf128
+#define f64sqrtf64x __hide_f64sqrtf64x
+#define f64sqrtf128 __hide_f64sqrtf128
+#include <math.h>
+#undef f32xsqrtf64x
+#undef f32xsqrtf128
+#undef f64sqrtf64x
+#undef f64sqrtf128
+
+#include <math-narrow.h>
+
+double
+__dsqrtl (_Float128 x)
+{
+  NARROW_SQRT_ROUND_TO_ODD (x, double, union ieee854_long_double, l,
+			    mantissa3);
+}
+libm_alias_double_ldouble (sqrt)
diff --git a/sysdeps/ieee754/ldbl-128/s_f64xsqrtf128.c b/sysdeps/ieee754/ldbl-128/s_f64xsqrtf128.c
new file mode 100644
index 0000000000..2f1898d48c
--- /dev/null
+++ b/sysdeps/ieee754/ldbl-128/s_f64xsqrtf128.c
@@ -0,0 +1,38 @@
+/* Square root of _Float128 value, converting the result to _Float64x.
+   Copyright (C) 2021 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 <math.h>
+#include <math-narrow.h>
+
+/* math_ldbl.h defines _Float128 to long double for this directory,
+   but when they are different, this function must be defined with
+   _Float128 arguments to avoid defining an alias with an incompatible
+   type.  */
+#undef _Float128
+
+#if __HAVE_FLOAT64X_LONG_DOUBLE && __HAVE_DISTINCT_FLOAT128
+_Float64x
+__f64xsqrtf128 (_Float128 x)
+{
+  NARROW_SQRT_ROUND_TO_ODD (x, _Float64x, union ieee854_long_double, l,
+			    mantissa3);
+}
+libm_alias_float64x_float128 (sqrt)
+#else
+/* Defined as an alias of sqrtl.  */
+#endif
diff --git a/sysdeps/ieee754/ldbl-128/s_fsqrtl.c b/sysdeps/ieee754/ldbl-128/s_fsqrtl.c
new file mode 100644
index 0000000000..95afe7b7ad
--- /dev/null
+++ b/sysdeps/ieee754/ldbl-128/s_fsqrtl.c
@@ -0,0 +1,33 @@
+/* Square root of long double (ldbl-128) value, narrowing the result to float.
+   Copyright (C) 2021 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/>.  */
+
+#define f32sqrtf64x __hide_f32sqrtf64x
+#define f32sqrtf128 __hide_f32sqrtf128
+#include <math.h>
+#undef f32sqrtf64x
+#undef f32sqrtf128
+
+#include <math-narrow.h>
+
+float
+__fsqrtl (_Float128 x)
+{
+  NARROW_SQRT_ROUND_TO_ODD (x, float, union ieee854_long_double, l,
+			    mantissa3);
+}
+libm_alias_float_ldouble (sqrt)
diff --git a/sysdeps/ieee754/ldbl-128ibm-compat/Versions b/sysdeps/ieee754/ldbl-128ibm-compat/Versions
index 07e6a5b08a..885ffb977f 100644
--- a/sysdeps/ieee754/ldbl-128ibm-compat/Versions
+++ b/sysdeps/ieee754/ldbl-128ibm-compat/Versions
@@ -123,6 +123,10 @@ libm {
     __y1ieee128;
     __ynieee128;
   }
+  GLIBC_2.35 {
+    __f32sqrtieee128;
+    __f64sqrtieee128;
+  }
 }
 libc {
   LDBL_IBM128_VERSION {
diff --git a/sysdeps/ieee754/ldbl-128ibm/s_dsqrtl.c b/sysdeps/ieee754/ldbl-128ibm/s_dsqrtl.c
new file mode 100644
index 0000000000..786ea7001b
--- /dev/null
+++ b/sysdeps/ieee754/ldbl-128ibm/s_dsqrtl.c
@@ -0,0 +1,28 @@
+/* Square root of long double (ldbl-128ibm) value, narrowing the result to
+   double.
+   Copyright (C) 2021 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 <math.h>
+#include <math-narrow.h>
+
+double
+__dsqrtl (long double x)
+{
+  NARROW_SQRT_TRIVIAL (x, double, l);
+}
+libm_alias_double_ldouble (sqrt)
diff --git a/sysdeps/ieee754/ldbl-128ibm/s_fsqrtl.c b/sysdeps/ieee754/ldbl-128ibm/s_fsqrtl.c
new file mode 100644
index 0000000000..abfc0d69fe
--- /dev/null
+++ b/sysdeps/ieee754/ldbl-128ibm/s_fsqrtl.c
@@ -0,0 +1,28 @@
+/* Square root of long double (ldbl-128ibm) value, narrowing the result to
+   float.
+   Copyright (C) 2021 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 <math.h>
+#include <math-narrow.h>
+
+float
+__fsqrtl (long double x)
+{
+  NARROW_SQRT_TRIVIAL (x, float, l);
+}
+libm_alias_float_ldouble (sqrt)
diff --git a/sysdeps/ieee754/ldbl-96/s_dsqrtl.c b/sysdeps/ieee754/ldbl-96/s_dsqrtl.c
new file mode 100644
index 0000000000..8bfc9624c2
--- /dev/null
+++ b/sysdeps/ieee754/ldbl-96/s_dsqrtl.c
@@ -0,0 +1,33 @@
+/* Square root of long double (ldbl-96) value, narrowing the result to double.
+   Copyright (C) 2021 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/>.  */
+
+#define f32xsqrtf64x __hide_f32xsqrtf64x
+#define f64sqrtf64x __hide_f64sqrtf64x
+#include <math.h>
+#undef f32xsqrtf64x
+#undef f64sqrtf64x
+
+#include <math-narrow.h>
+
+double
+__dsqrtl (long double x)
+{
+  NARROW_SQRT_ROUND_TO_ODD (x, double, union ieee854_long_double, l,
+			    mantissa1);
+}
+libm_alias_double_ldouble (sqrt)
diff --git a/sysdeps/ieee754/ldbl-96/s_fsqrtl.c b/sysdeps/ieee754/ldbl-96/s_fsqrtl.c
new file mode 100644
index 0000000000..026add5d03
--- /dev/null
+++ b/sysdeps/ieee754/ldbl-96/s_fsqrtl.c
@@ -0,0 +1,31 @@
+/* Square root of long double (ldbl-96) value, narrowing the result to float.
+   Copyright (C) 2021 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/>.  */
+
+#define f32sqrtf64x __hide_f32sqrtf64x
+#include <math.h>
+#undef f32sqrtf64x
+
+#include <math-narrow.h>
+
+float
+__fsqrtl (long double x)
+{
+  NARROW_SQRT_ROUND_TO_ODD (x, float, union ieee854_long_double, l,
+			    mantissa1);
+}
+libm_alias_float_ldouble (sqrt)
diff --git a/sysdeps/ieee754/ldbl-opt/Makefile b/sysdeps/ieee754/ldbl-opt/Makefile
index 4e41a27e80..dff53cc0d8 100644
--- a/sysdeps/ieee754/ldbl-opt/Makefile
+++ b/sysdeps/ieee754/ldbl-opt/Makefile
@@ -45,7 +45,7 @@ libnldbl-calls = asprintf dprintf fprintf fscanf fwprintf fwscanf iovfscanf \
 		 nextup nextdown totalorder totalordermag getpayload \
 		 canonicalize setpayload setpayloadsig llogb fmaxmag fminmag \
 		 roundeven fromfp ufromfp fromfpx ufromfpx fadd dadd \
-		 fdiv ddiv fmul dmul fsub dsub
+		 fdiv ddiv fmul dmul fsqrt dsqrt fsub dsub
 libnldbl-routines = $(libnldbl-calls:%=nldbl-%)
 libnldbl-inhibit-o = $(object-suffixes)
 libnldbl-static-only-routines = $(libnldbl-routines)
@@ -90,6 +90,7 @@ CFLAGS-nldbl-ctanh.c = -fno-builtin-ctanhl
 CFLAGS-nldbl-dadd.c = -fno-builtin-daddl
 CFLAGS-nldbl-ddiv.c = -fno-builtin-ddivl
 CFLAGS-nldbl-dmul.c = -fno-builtin-dmull
+CFLAGS-nldbl-dsqrt.c = -fno-builtin-dsqrtl
 CFLAGS-nldbl-dsub.c = -fno-builtin-dsubl
 CFLAGS-nldbl-erf.c = -fno-builtin-erfl
 CFLAGS-nldbl-erfc.c = -fno-builtin-erfcl
@@ -113,6 +114,7 @@ CFLAGS-nldbl-fmul.c = -fno-builtin-fmull
 CFLAGS-nldbl-frexp.c = -fno-builtin-frexpl
 CFLAGS-nldbl-fromfp.c = -fno-builtin-fromfpl
 CFLAGS-nldbl-fromfpx.c = -fno-builtin-fromfpxl
+CFLAGS-nldbl-fsqrt.c = -fno-builtin-fsqrtl
 CFLAGS-nldbl-fsub.c = -fno-builtin-fsubl
 CFLAGS-nldbl-gamma.c = -fno-builtin-gammal
 CFLAGS-nldbl-getpayload.c = -fno-builtin-getpayloadl
diff --git a/sysdeps/ieee754/ldbl-opt/nldbl-dsqrt.c b/sysdeps/ieee754/ldbl-opt/nldbl-dsqrt.c
new file mode 100644
index 0000000000..1d131f083d
--- /dev/null
+++ b/sysdeps/ieee754/ldbl-opt/nldbl-dsqrt.c
@@ -0,0 +1,28 @@
+/* Compatibility routine for IEEE double as long double for dsqrt.
+   Copyright (C) 2021 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/>.  */
+
+#define dsqrtl __hide_dsqrtl
+#include "nldbl-compat.h"
+#undef dsqrtl
+
+double
+attribute_hidden
+dsqrtl (double x)
+{
+  return sqrt (x);
+}
diff --git a/sysdeps/ieee754/ldbl-opt/nldbl-fsqrt.c b/sysdeps/ieee754/ldbl-opt/nldbl-fsqrt.c
new file mode 100644
index 0000000000..bcdf1aa556
--- /dev/null
+++ b/sysdeps/ieee754/ldbl-opt/nldbl-fsqrt.c
@@ -0,0 +1,28 @@
+/* Compatibility routine for IEEE double as long double for fsqrt.
+   Copyright (C) 2021 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/>.  */
+
+#define fsqrtl __hide_fsqrtl
+#include "nldbl-compat.h"
+#undef fsqrtl
+
+float
+attribute_hidden
+fsqrtl (double x)
+{
+  return fsqrt (x);
+}
diff --git a/sysdeps/ieee754/soft-fp/s_dsqrtl.c b/sysdeps/ieee754/soft-fp/s_dsqrtl.c
new file mode 100644
index 0000000000..85ea8538a5
--- /dev/null
+++ b/sysdeps/ieee754/soft-fp/s_dsqrtl.c
@@ -0,0 +1,67 @@
+/* Square root of long double (ldbl-128) value, narrowing the result to
+   double, using soft-fp.
+   Copyright (C) 2021 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/>.  */
+
+#define f32xsqrtf64x __hide_f32xsqrtf64x
+#define f32xsqrtf128 __hide_f32xsqrtf128
+#define f64sqrtf64x __hide_f64sqrtf64x
+#define f64sqrtf128 __hide_f64sqrtf128
+#include <math.h>
+#undef f32xsqrtf64x
+#undef f32xsqrtf128
+#undef f64sqrtf64x
+#undef f64sqrtf128
+
+#include <math-narrow.h>
+#include <libc-diag.h>
+
+/* Some components of the result's significand and exponent are not
+   set in cases where they are not used in packing, but the compiler
+   does not see that they are set in all cases where they are used,
+   resulting in warnings that they may be used uninitialized.  */
+DIAG_PUSH_NEEDS_COMMENT;
+DIAG_IGNORE_NEEDS_COMMENT (7, "-Wmaybe-uninitialized");
+#include <soft-fp.h>
+#include <double.h>
+#include <quad.h>
+
+double
+__dsqrtl (_Float128 x)
+{
+  FP_DECL_EX;
+  FP_DECL_Q (X);
+  FP_DECL_Q (R);
+  FP_DECL_D (RN);
+  double ret;
+
+  FP_INIT_ROUNDMODE;
+  FP_UNPACK_Q (X, x);
+  FP_SQRT_Q (R, X);
+#if _FP_W_TYPE_SIZE < 64
+  FP_TRUNC_COOKED (D, Q, 2, 4, RN, R);
+#else
+  FP_TRUNC_COOKED (D, Q, 1, 2, RN, R);
+#endif
+  FP_PACK_D (ret, RN);
+  FP_HANDLE_EXCEPTIONS;
+  CHECK_NARROW_SQRT (ret, x);
+  return ret;
+}
+DIAG_POP_NEEDS_COMMENT;
+
+libm_alias_double_ldouble (sqrt)
diff --git a/sysdeps/ieee754/soft-fp/s_fsqrt.c b/sysdeps/ieee754/soft-fp/s_fsqrt.c
new file mode 100644
index 0000000000..b6f7c4e962
--- /dev/null
+++ b/sysdeps/ieee754/soft-fp/s_fsqrt.c
@@ -0,0 +1,54 @@
+/* Square root of double value, narrowing the result to float, using soft-fp.
+   Copyright (C) 2021 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/>.  */
+
+#define f32sqrtf64 __hide_f32sqrtf64
+#define f32sqrtf32x __hide_f32sqrtf32x
+#define fsqrtl __hide_fsqrtl
+#include <math.h>
+#undef f32sqrtf64
+#undef f32sqrtf32x
+#undef fsqrtl
+
+#include <math-narrow.h>
+#include <soft-fp.h>
+#include <single.h>
+#include <double.h>
+
+float
+__fsqrt (double x)
+{
+  FP_DECL_EX;
+  FP_DECL_D (X);
+  FP_DECL_D (R);
+  FP_DECL_S (RN);
+  float ret;
+
+  FP_INIT_ROUNDMODE;
+  FP_UNPACK_D (X, x);
+  FP_SQRT_D (R, X);
+#if _FP_W_TYPE_SIZE < _FP_FRACBITS_D
+  FP_TRUNC_COOKED (S, D, 1, 2, RN, R);
+#else
+  FP_TRUNC_COOKED (S, D, 1, 1, RN, R);
+#endif
+  FP_PACK_S (ret, RN);
+  FP_HANDLE_EXCEPTIONS;
+  CHECK_NARROW_SQRT (ret, x);
+  return ret;
+}
+libm_alias_float_double (sqrt)
diff --git a/sysdeps/ieee754/soft-fp/s_fsqrtl.c b/sysdeps/ieee754/soft-fp/s_fsqrtl.c
new file mode 100644
index 0000000000..618e734744
--- /dev/null
+++ b/sysdeps/ieee754/soft-fp/s_fsqrtl.c
@@ -0,0 +1,53 @@
+/* Square root of long double (ldbl-128) value, narrowing the result to
+   float, using soft-fp.
+   Copyright (C) 2021 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/>.  */
+
+#define f32sqrtf64x __hide_f32sqrtf64x
+#define f32sqrtf128 __hide_f32sqrtf128
+#include <math.h>
+#undef f32sqrtf64x
+#undef f32sqrtf128
+
+#include <math-narrow.h>
+#include <soft-fp.h>
+#include <single.h>
+#include <quad.h>
+
+float
+__fsqrtl (_Float128 x)
+{
+  FP_DECL_EX;
+  FP_DECL_Q (X);
+  FP_DECL_Q (R);
+  FP_DECL_S (RN);
+  float ret;
+
+  FP_INIT_ROUNDMODE;
+  FP_UNPACK_Q (X, x);
+  FP_SQRT_Q (R, X);
+#if _FP_W_TYPE_SIZE < 64
+  FP_TRUNC_COOKED (S, Q, 1, 4, RN, R);
+#else
+  FP_TRUNC_COOKED (S, Q, 1, 2, RN, R);
+#endif
+  FP_PACK_S (ret, RN);
+  FP_HANDLE_EXCEPTIONS;
+  CHECK_NARROW_SQRT (ret, x);
+  return ret;
+}
+libm_alias_float_ldouble (sqrt)