[RFC V4] Enable libmvec support for RISC-V

[RFC V4] Enable libmvec support for RISC-V

[shiyulong]

From: yulong <shiyulong@iscas.ac.cn>

Diff: Chande the version from GLIBC_2.39 to GLIBC_2.40.
This patch tries to enable libmvec on RISC-V. I also have demonstrated
how this all fits together by adding implementations for vector cos.
This patch is a try and we hope to receive valuable comments.

Thanks,
yulong

---
 sysdeps/riscv/configure                       |   4 +
 sysdeps/riscv/configure.ac                    |   4 +
 sysdeps/riscv/rvd/Makefile                    |   5 +
 sysdeps/riscv/rvd/Versions                    |   5 +
 sysdeps/riscv/rvd/bits/math-vector.h          |  29 ++++
 sysdeps/riscv/rvd/cos.c                       |  94 ++++++++++++
 sysdeps/riscv/rvd/math_private.h              |  42 ++++++
 sysdeps/riscv/rvd/v_math.h                    | 139 ++++++++++++++++++
 sysdeps/riscv/rvd/vecmath_config.h            |  33 +++++
 sysdeps/unix/sysv/linux/riscv/libmvec.abilist |   1 +
 10 files changed, 356 insertions(+)
 mode change 100644 => 100755 sysdeps/riscv/configure
 create mode 100644 sysdeps/riscv/rvd/Makefile
 create mode 100644 sysdeps/riscv/rvd/Versions
 create mode 100644 sysdeps/riscv/rvd/bits/math-vector.h
 create mode 100644 sysdeps/riscv/rvd/cos.c
 create mode 100644 sysdeps/riscv/rvd/math_private.h
 create mode 100644 sysdeps/riscv/rvd/v_math.h
 create mode 100644 sysdeps/riscv/rvd/vecmath_config.h
 create mode 100644 sysdeps/unix/sysv/linux/riscv/libmvec.abilist

diff --git a/sysdeps/riscv/configure b/sysdeps/riscv/configure
old mode 100644
new mode 100755
index c8f01709f8..a6d0b4becb
--- a/sysdeps/riscv/configure
+++ b/sysdeps/riscv/configure
@@ -80,3 +80,7 @@ if test "$libc_cv_static_pie_on_riscv" = yes; then
   printf "%s\n" "#define SUPPORT_STATIC_PIE 1" >>confdefs.h
 
 fi
+
+if test x"$build_mathvec" = xnotset; then
+  build_mathvec=yes
+fi
diff --git a/sysdeps/riscv/configure.ac b/sysdeps/riscv/configure.ac
index ee3d1ed014..b1c1105baa 100644
--- a/sysdeps/riscv/configure.ac
+++ b/sysdeps/riscv/configure.ac
@@ -43,3 +43,7 @@ EOF
 if test "$libc_cv_static_pie_on_riscv" = yes; then
   AC_DEFINE(SUPPORT_STATIC_PIE)
 fi
+
+if test x"$build_mathvec" = xnotset; then
+  build_mathvec=yes
+fi
diff --git a/sysdeps/riscv/rvd/Makefile b/sysdeps/riscv/rvd/Makefile
new file mode 100644
index 0000000000..1adb2ee582
--- /dev/null
+++ b/sysdeps/riscv/rvd/Makefile
@@ -0,0 +1,5 @@
+libmvec-supported-funcs = cos
+
+ifeq ($(subdir),mathvec)
+libmvec-support = $(addprefix d,$(libmvec-supported-funcs))
+endif
diff --git a/sysdeps/riscv/rvd/Versions b/sysdeps/riscv/rvd/Versions
new file mode 100644
index 0000000000..0fd283329c
--- /dev/null
+++ b/sysdeps/riscv/rvd/Versions
@@ -0,0 +1,5 @@
+libmvec {
+  GLIBC_2.40 {
+    _ZGVnN2v_cos;
+  }
+}
diff --git a/sysdeps/riscv/rvd/bits/math-vector.h b/sysdeps/riscv/rvd/bits/math-vector.h
new file mode 100644
index 0000000000..b34ffc9bc1
--- /dev/null
+++ b/sysdeps/riscv/rvd/bits/math-vector.h
@@ -0,0 +1,29 @@
+/* Platform-specific SIMD declarations of math functions.
+
+   Copyright (C) 2024 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/>.  */
+
+#ifndef _MATH_H
+#  error "Never include <bits/math-vector.h> directly;\
+ include <math.h> instead."
+#endif
+
+#if defined __riscv__
+# define __DECL_RVV_RISCV _Pragma
+# undef __DECL_RVV_cos
+# define __DECL_RVV_cos __DECL_RVV_RISCV
+#endif
diff --git a/sysdeps/riscv/rvd/cos.c b/sysdeps/riscv/rvd/cos.c
new file mode 100644
index 0000000000..1806acd629
--- /dev/null
+++ b/sysdeps/riscv/rvd/cos.c
@@ -0,0 +1,94 @@
+/* Double-precision vector cos function.
+
+   Copyright (C) 2024 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 "v_math.h"
+
+
+static const struct data
+{
+  vfloat64m2_t poly[7];
+  vfloat64m2_t range_val, shift, inv_pi, half_pi, pi_1, pi_2, pi_3;
+} data = {
+  /* Worst-case error is 3.3 ulp in [-pi/2, pi/2].  */
+  .poly = { V2 (-0x1.555555555547bp-3), V2 (0x1.1111111108a4dp-7),
+	    V2 (-0x1.a01a019936f27p-13), V2 (0x1.71de37a97d93ep-19),
+	    V2 (-0x1.ae633919987c6p-26), V2 (0x1.60e277ae07cecp-33),
+	    V2 (-0x1.9e9540300a1p-41) },
+  .inv_pi = V2 (0x1.45f306dc9c883p-2),
+  .half_pi = V2 (0x1.921fb54442d18p+0),
+  .pi_1 = V2 (0x1.921fb54442d18p+1),
+  .pi_2 = V2 (0x1.1a62633145c06p-53),
+  .pi_3 = V2 (0x1.c1cd129024e09p-106),
+  .shift = V2 (0x1.8p52),
+  .range_val = V2 (0x1p23)
+};
+
+#define C(i) d->poly[i]
+
+static vfloat64m2_t NOINLINE
+special_case (vfloat64m2_t x, vfloat64m2_t y, vuint64m2_t odd, vuint64m2_t cmp)
+{
+  y = vreinterpret_v_u64m2_f64m2 (vor (vreinterpret_v_f64m2_u64m2 (y), odd, 1));
+  return v_call_f64 (cos, x, y, cmp);
+}
+
+vfloat64m2_t V_NAME_D1 (cos) (vfloat64m2_t x)
+{
+  const struct data *d = ptr_barrier (&data);
+  vfloat64m2_t n, r, r2, r3, r4, t1, t2, t3, y;
+  vuint64m2_t odd, cmp;
+
+  r = vfabs_v_f64m2 (x, 2);
+  cmp = (vuint64m2_t) vmsgeu (vreinterpret_v_f64m2_u64m2 (r),
+		   vreinterpret_v_f64m2_u64m2 (d->range_val));
+  if (__glibc_unlikely (v_any_u64 (cmp)))
+    /* If fenv exceptions are to be triggered correctly, set any special lanes
+       to 1 (which is neutral w.r.t. fenv). These lanes will be fixed by
+       special-case handler later.  */
+    r = vmsltu (cmp, v_f64 (1.0), r);
+
+  /* n = rint((|x|+pi/2)/pi) - 0.5.  */
+  n = vfmadd (d->shift, d->inv_pi, vfadd (r, d->half_pi,2), 2);
+  odd = vshlq_n_u64 (vreinterpret_v_f64m2_u64m2 (n), 63);
+  n = vfsub (n, d->shift, 2);
+  n = vfsub (n, v_f64 (0.5), 2);
+
+  /* r = |x| - n*pi  (range reduction into -pi/2 .. pi/2).  */
+  r = vfmsub (r, d->pi_1, n, 2);
+  r = vfmsub (r, d->pi_2, n, 2);
+  r = vfmsub (r, d->pi_3, n, 2);
+
+  /* sin(r) poly approx.  */
+  r2 = vfmul (r, r, 2);
+  r3 = vfmul (r2, r, 2);
+  r4 = vfmul (r2, r2, 2);
+
+  t1 = vfmadd (C (4), C (5), r2, 2);
+  t2 = vfmadd (C (2), C (3), r2, 2);
+  t3 = vfmadd (C (0), C (1), r2, 2);
+
+  y = vfmadd (t1, C (6), r4, 2);
+  y = vfmadd (t2, y, r4, 2);
+  y = vfmadd (t3, y, r4, 2);
+  y = vfmadd (r, y, r3, 2);
+
+  if (__glibc_unlikely (v_any_u64 (cmp)))
+    return special_case (x, y, odd, cmp);
+  return vreinterpretq_f64_u64 (vor (vreinterpret_v_f64m2_u64m2 (y), odd, 2));
+}
diff --git a/sysdeps/riscv/rvd/math_private.h b/sysdeps/riscv/rvd/math_private.h
new file mode 100644
index 0000000000..655a4dcd55
--- /dev/null
+++ b/sysdeps/riscv/rvd/math_private.h
@@ -0,0 +1,42 @@
+/* Configure optimized libm functions.  RISC-V version.
+   Copyright (C) 2024 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/>.  */
+
+#ifndef RISCV_MATH_PRIVATE_H
+#define RISCV_MATH_PRIVATE_H 1
+
+#include <stdint.h>
+#include <math.h>
+
+/* Use inline round and lround instructions.  */
+#define TOINT_INTRINSICS 1
+
+static inline double_t
+roundtoint (double_t x)
+{
+  return round (x);
+}
+
+static inline int32_t
+converttoint (double_t x)
+{
+  return lround (x);
+}
+
+#include_next <math_private.h>
+
+#endif
diff --git a/sysdeps/riscv/rvd/v_math.h b/sysdeps/riscv/rvd/v_math.h
new file mode 100644
index 0000000000..d2e821aeb2
--- /dev/null
+++ b/sysdeps/riscv/rvd/v_math.h
@@ -0,0 +1,139 @@
+/* Utilities for Advanced SIMD libmvec routines.
+   Copyright (C) 2024 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/>.  */
+
+#ifndef _V_MATH_H
+#define _V_MATH_H
+
+#include <riscv_vector.h>
+#include "vecmath_config.h"
+
+#define V_NAME_D1(fun) _ZGVnN2v_##fun
+
+/* Shorthand helpers for declaring constants.  */
+#define V2(X) { X, X }
+#define V4(X) { X, X, X, X }
+#define V8(X) { X, X, X, X, X, X, X, X }
+
+static inline vfloat32m4_t
+v_f32 (float x)
+{
+  return (vfloat32m4_t) V4 (x);
+}
+static inline vuint32m4_t
+v_u32 (uint32_t x)
+{
+  return (vuint32m4_t) V4 (x);
+}
+static inline vint32m4_t
+v_s32 (int32_t x)
+{
+  return (vint32m4_t) V4 (x);
+}
+
+/* true if any elements of a vector compare result is non-zero.  */
+static inline int
+v_any_u32 (vuint32m4_t x)
+{
+  /* assume elements in x are either 0 or -1u.  */
+  return vpaddd_u64 (vreinterpret_v_u64m2_u32m2 (x)) != 0;
+}
+static inline int
+v_any_u32h (vuint32m2_t x)
+{
+  return vget_lane_u64 (vreinterpret_v_u32m2_u64m2 (x), 0) != 0;
+}
+static inline vfloat32m4_t
+v_lookup_f32 (const float *tab, vuint32m4_t idx)
+{
+  return (vfloat32m4_t){ tab[idx[0]], tab[idx[1]], tab[idx[2]], tab[idx[3]] };
+}
+static inline vuint32m4_t
+v_lookup_u32 (const uint32_t *tab, vuint32m4_t idx)
+{
+  return (vuint32m4_t){ tab[idx[0]], tab[idx[1]], tab[idx[2]], tab[idx[3]] };
+}
+static inline vfloat32m4_t
+v_call_f32 (float (*f) (float), vfloat32m4_t x, vfloat32m4_t y, vuint32m4_t p)
+{
+  return (vfloat32m4_t){ p[0] ? f (x[0]) : y[0], p[1] ? f (x[1]) : y[1],
+			p[2] ? f (x[2]) : y[2], p[3] ? f (x[3]) : y[3] };
+}
+static inline vfloat32m4_t
+v_call2_f32 (float (*f) (float, float), vfloat32m4_t x1, vfloat32m4_t x2,
+	     vfloat32m4_t y, vuint32m4_t p)
+{
+  return (vfloat32m4_t){ p[0] ? f (x1[0], x2[0]) : y[0],
+			p[1] ? f (x1[1], x2[1]) : y[1],
+			p[2] ? f (x1[2], x2[2]) : y[2],
+			p[3] ? f (x1[3], x2[3]) : y[3] };
+}
+
+static inline vfloat64m2_t
+v_f64 (double x)
+{
+  return (vfloat64m2_t) V2 (x);
+}
+static inline vuint64m2_t
+v_u64 (uint64_t x)
+{
+  return (vuint64m2_t) V2 (x);
+}
+static inline vint64m2_t
+v_s64 (int64_t x)
+{
+  return (vint64m2_t) V2 (x);
+}
+
+/* true if any elements of a vector compare result is non-zero.  */
+static inline int
+v_any_u64 (vuint64m1_t x)
+{
+  /* assume elements in x are either 0 or -1u.  */
+  return vpaddd_u64 (x) != 0;
+}
+/* true if all elements of a vector compare result is 1.  */
+static inline int
+v_all_u64 (vuint64m1_t x)
+{
+  /* assume elements in x are either 0 or -1u.  */
+  return vpaddd_s64 (vreinterpretq_s64_u64 (x)) == -2;
+}
+static inline vfloat64m1_t
+v_lookup_f64 (const double *tab, vuint64m1_t idx)
+{
+  return (vfloat64m1_t){ tab[idx[0]], tab[idx[1]] };
+}
+static inline vuint64m1_t
+v_lookup_u64 (const uint64_t *tab, vuint64m1_t idx)
+{
+  return (vuint64m1_t){ tab[idx[0]], tab[idx[1]] };
+}
+static inline vfloat64m1_t
+v_call_f64 (double (*f) (double), vfloat64m1_t x, vfloat64m1_t y, vuint64m1_t p)
+{
+  return (vfloat64m1_t){ p[0] ? f (x[0]) : y[0], p[1] ? f (x[1]) : y[1] };
+}
+static inline vfloat64m1_t
+v_call2_f64 (double (*f) (double, double), vfloat64m1_t x1, vfloat64m1_t x2,
+	     vfloat64m1_t y, vuint64m1_t p)
+{
+  return (vfloat64m1_t){ p[0] ? f (x1[0], x2[0]) : y[0],
+			p[1] ? f (x1[1], x2[1]) : y[1] };
+}
+
+#endif
diff --git a/sysdeps/riscv/rvd/vecmath_config.h b/sysdeps/riscv/rvd/vecmath_config.h
new file mode 100644
index 0000000000..290ea1e33c
--- /dev/null
+++ b/sysdeps/riscv/rvd/vecmath_config.h
@@ -0,0 +1,33 @@
+/* Configuration for libmvec routines.
+   Copyright (C) 2023 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/>.  */
+
+#ifndef _VECMATH_CONFIG_H
+#define _VECMATH_CONFIG_H
+
+#include <math_private.h>
+
+/* Return ptr but hide its value from the compiler so accesses through it
+   cannot be optimized based on the contents.  */
+#define ptr_barrier(ptr)                                                      \
+  ({                                                                          \
+    __typeof (ptr) __ptr = (ptr);                                             \
+    __asm("" : "+r"(__ptr));                                                  \
+    __ptr;                                                                    \
+  })
+
+#endif
diff --git a/sysdeps/unix/sysv/linux/riscv/libmvec.abilist b/sysdeps/unix/sysv/linux/riscv/libmvec.abilist
new file mode 100644
index 0000000000..fe8141b189
--- /dev/null
+++ b/sysdeps/unix/sysv/linux/riscv/libmvec.abilist
@@ -0,0 +1 @@
+GLIBC_2.40 _ZGVnN2v_cos F
-- 
2.34.1

Re: [RFC V4] Enable libmvec support for RISC-V

[Jeff Law]

On 4/15/24 1:21 AM, shiyulong@iscas.ac.cn wrote:
> From: yulong <shiyulong@iscas.ac.cn>
> 
> Diff: Chande the version from GLIBC_2.39 to GLIBC_2.40.
> This patch tries to enable libmvec on RISC-V. I also have demonstrated
> how this all fits together by adding implementations for vector cos.
> This patch is a try and we hope to receive valuable comments.
Just an FYI -- Palmer's team over at Rivos have implementations for a 
number of routines that would fit into libmvec.  You might reach out to 
Ping Tak Peter Tang  <ptpt@rivosinc.com> for information in his 
implementation.

> https://github.com/rivosinc/veclibm/


THeir implementations may provide good guidance on performant 
implementations of various routines that libmvec typically provides.

jeff

Re: [RFC V4] Enable libmvec support for RISC-V

[yulong]

在 2024/4/25 13:07, Jeff Law 写道:
>
>
> On 4/15/24 1:21 AM, shiyulong@iscas.ac.cn wrote:
>> From: yulong <shiyulong@iscas.ac.cn>
>>
>> Diff: Chande the version from GLIBC_2.39 to GLIBC_2.40.
>> This patch tries to enable libmvec on RISC-V. I also have demonstrated
>> how this all fits together by adding implementations for vector cos.
>> This patch is a try and we hope to receive valuable comments.
> Just an FYI -- Palmer's team over at Rivos have implementations for a 
> number of routines that would fit into libmvec.  You might reach out 
> to Ping Tak Peter Tang  <ptpt@rivosinc.com> for information in his 
> implementation.
>
>> https://github.com/rivosinc/veclibm/
>
>
> THeir implementations may provide good guidance on performant 
> implementations of various routines that libmvec typically provides.
>
> jeff
Thanks Jeff for your advice, I'm working on a new implementation after 
reading the above code.

Re: [RFC V4] Enable libmvec support for RISC-V

[Palmer Dabbelt]

On Wed, 24 Apr 2024 22:07:31 PDT (-0700), jeffreyalaw@gmail.com wrote:
>
>
> On 4/15/24 1:21 AM, shiyulong@iscas.ac.cn wrote:
>> From: yulong <shiyulong@iscas.ac.cn>
>>
>> Diff: Chande the version from GLIBC_2.39 to GLIBC_2.40.
>> This patch tries to enable libmvec on RISC-V. I also have demonstrated
>> how this all fits together by adding implementations for vector cos.
>> This patch is a try and we hope to receive valuable comments.
> Just an FYI -- Palmer's team over at Rivos have implementations for a
> number of routines that would fit into libmvec.  You might reach out to
> Ping Tak Peter Tang  <ptpt@rivosinc.com> for information in his
> implementation.
>
>> https://github.com/rivosinc/veclibm/
>
>
> THeir implementations may provide good guidance on performant
> implementations of various routines that libmvec typically provides.

Ya, that's the idea of veclibm.  The actual functions are written in a 
way that's more suitable for some other libraries, but the core 
computational implemenations should be the same.  A few of us had 
briefly talked internally about getting these into glibc, IIUC all the 
code was written at Rivos and thus could be copyright assigned to the 
FSF and used in glibc.  We don't have time to do that right now, but if 
you're interested in helping that'd be awesome.  We'll need to be 
careful with the copyright/licensing, though.

That said, I've never really quite managed to figure out how all the 
libmvec stuff is supposed to fit together.  I'm more worried about the 
ABI side of things than the implementation, so I think starting with 
just one function to get the ABI template figure out is a reasonable way 
to go and we can get the rest of the implementations ported over next.  
The first thing that jumps out on the ABI side of things is cos() taking 
EMUL=2 types, I'm not sure if there's a reason for that but it seems 
we'd want EMUL=1 to fit more data in the argument registers?

Also, I think some of this can be split out: the roundtoint/converttoint 
isn't really a libmvec thing (see 
https://inbox.sourceware.org/libc-alpha/20220803174258.4235-1-palmer@rivosinc.com/, 
which fails some test), and ptr_barrier() can probably be pulled out to 
something generic as it's the same as arm64's version.

I'm also only seeing draft versions of the vector intrinsics.  I know we 
merged them into GCC and usually that means things are stable, but we 
merged these pre-freeze (based on some assertions things wouldn't 
change) and things have drifted around a bit it the spec.  I think we're 
probably safe just depending on the types, if there's no frozen version 
we should at least write down exactly which version we're following 
though.

Also: are there GCC patches for these?  It'd be great to be able to test 
things through the whole codegen stack so we can make sure it works.

>
> jeff