mirrors/glibc
1
0
Fork 0
mirror of https://sourceware.org/git/glibc.git synced 2024-11-23 09:43:32 +08:00
Code Issues Actions Packages Projects Releases Wiki Activity

aarch64: Add vector implementations of tan routines

Browse Source 
This includes some utility headers for evaluating polynomials using
various schemes.
This commit is contained in:
Joe Ramsay 2023-10-05 17:10:48 +01:00 committed by Szabolcs Nagy
parent 2aa0974d25
commit f554334c05
21 changed files with 1274 additions and 27 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
math/auto-libm-test-in
View File

@ -7655,7 +7655,7 @@ sqrt min
sqrt min_subnorm
tan 0
tan -0
tan -0 no-mathvec
tan pi/4
tan pi/2
tan -pi/2

50
math/auto-libm-test-out-tan
View File

@ -23,31 +23,31 @@ tan 0
= tan tonearest ibm128 0x0p+0 : 0x0p+0 : inexact-ok
= tan towardzero ibm128 0x0p+0 : 0x0p+0 : inexact-ok
= tan upward ibm128 0x0p+0 : 0x0p+0 : inexact-ok
tan -0
= tan downward binary32 -0x0p+0 : -0x0p+0 : inexact-ok
= tan tonearest binary32 -0x0p+0 : -0x0p+0 : inexact-ok
= tan towardzero binary32 -0x0p+0 : -0x0p+0 : inexact-ok
= tan upward binary32 -0x0p+0 : -0x0p+0 : inexact-ok
= tan downward binary64 -0x0p+0 : -0x0p+0 : inexact-ok
= tan tonearest binary64 -0x0p+0 : -0x0p+0 : inexact-ok
= tan towardzero binary64 -0x0p+0 : -0x0p+0 : inexact-ok
= tan upward binary64 -0x0p+0 : -0x0p+0 : inexact-ok
= tan downward intel96 -0x0p+0 : -0x0p+0 : inexact-ok
= tan tonearest intel96 -0x0p+0 : -0x0p+0 : inexact-ok
= tan towardzero intel96 -0x0p+0 : -0x0p+0 : inexact-ok
= tan upward intel96 -0x0p+0 : -0x0p+0 : inexact-ok
= tan downward m68k96 -0x0p+0 : -0x0p+0 : inexact-ok
= tan tonearest m68k96 -0x0p+0 : -0x0p+0 : inexact-ok
= tan towardzero m68k96 -0x0p+0 : -0x0p+0 : inexact-ok
= tan upward m68k96 -0x0p+0 : -0x0p+0 : inexact-ok
= tan downward binary128 -0x0p+0 : -0x0p+0 : inexact-ok
= tan tonearest binary128 -0x0p+0 : -0x0p+0 : inexact-ok
= tan towardzero binary128 -0x0p+0 : -0x0p+0 : inexact-ok
= tan upward binary128 -0x0p+0 : -0x0p+0 : inexact-ok
= tan downward ibm128 -0x0p+0 : -0x0p+0 : inexact-ok
= tan tonearest ibm128 -0x0p+0 : -0x0p+0 : inexact-ok
= tan towardzero ibm128 -0x0p+0 : -0x0p+0 : inexact-ok
= tan upward ibm128 -0x0p+0 : -0x0p+0 : inexact-ok
tan -0 no-mathvec
= tan downward binary32 -0x0p+0 : -0x0p+0 : no-mathvec inexact-ok
= tan tonearest binary32 -0x0p+0 : -0x0p+0 : no-mathvec inexact-ok
= tan towardzero binary32 -0x0p+0 : -0x0p+0 : no-mathvec inexact-ok
= tan upward binary32 -0x0p+0 : -0x0p+0 : no-mathvec inexact-ok
= tan downward binary64 -0x0p+0 : -0x0p+0 : no-mathvec inexact-ok
= tan tonearest binary64 -0x0p+0 : -0x0p+0 : no-mathvec inexact-ok
= tan towardzero binary64 -0x0p+0 : -0x0p+0 : no-mathvec inexact-ok
= tan upward binary64 -0x0p+0 : -0x0p+0 : no-mathvec inexact-ok
= tan downward intel96 -0x0p+0 : -0x0p+0 : no-mathvec inexact-ok
= tan tonearest intel96 -0x0p+0 : -0x0p+0 : no-mathvec inexact-ok
= tan towardzero intel96 -0x0p+0 : -0x0p+0 : no-mathvec inexact-ok
= tan upward intel96 -0x0p+0 : -0x0p+0 : no-mathvec inexact-ok
= tan downward m68k96 -0x0p+0 : -0x0p+0 : no-mathvec inexact-ok
= tan tonearest m68k96 -0x0p+0 : -0x0p+0 : no-mathvec inexact-ok
= tan towardzero m68k96 -0x0p+0 : -0x0p+0 : no-mathvec inexact-ok
= tan upward m68k96 -0x0p+0 : -0x0p+0 : no-mathvec inexact-ok
= tan downward binary128 -0x0p+0 : -0x0p+0 : no-mathvec inexact-ok
= tan tonearest binary128 -0x0p+0 : -0x0p+0 : no-mathvec inexact-ok
= tan towardzero binary128 -0x0p+0 : -0x0p+0 : no-mathvec inexact-ok
= tan upward binary128 -0x0p+0 : -0x0p+0 : no-mathvec inexact-ok
= tan downward ibm128 -0x0p+0 : -0x0p+0 : no-mathvec inexact-ok
= tan tonearest ibm128 -0x0p+0 : -0x0p+0 : no-mathvec inexact-ok
= tan towardzero ibm128 -0x0p+0 : -0x0p+0 : no-mathvec inexact-ok
= tan upward ibm128 -0x0p+0 : -0x0p+0 : no-mathvec inexact-ok
tan pi/4
= tan downward binary32 0xc.90fdbp-4 : 0x1p+0 : inexact-ok
= tan tonearest binary32 0xc.90fdbp-4 : 0x1p+0 : inexact-ok

3
sysdeps/aarch64/fpu/Makefile
View File

@ -1,7 +1,8 @@
libmvec-supported-funcs = cos \
exp \
log \
sin
sin \
tan
float-advsimd-funcs = $(libmvec-supported-funcs)
double-advsimd-funcs = $(libmvec-supported-funcs)

6
sysdeps/aarch64/fpu/Versions
View File

@ -17,4 +17,10 @@ libmvec {
_ZGVsMxv_sin;
_ZGVsMxv_sinf;
}
GLIBC_2.39 {
_ZGVnN4v_tanf;
_ZGVnN2v_tan;
_ZGVsMxv_tanf;
_ZGVsMxv_tan;
}
}

4
sysdeps/aarch64/fpu/bits/math-vector.h
View File

@ -53,11 +53,13 @@ __vpcs __f32x4_t _ZGVnN4v_cosf (__f32x4_t);
__vpcs __f32x4_t _ZGVnN4v_expf (__f32x4_t);
__vpcs __f32x4_t _ZGVnN4v_logf (__f32x4_t);
__vpcs __f32x4_t _ZGVnN4v_sinf (__f32x4_t);
__vpcs __f32x4_t _ZGVnN4v_tanf (__f32x4_t);
__vpcs __f64x2_t _ZGVnN2v_cos (__f64x2_t);
__vpcs __f64x2_t _ZGVnN2v_exp (__f64x2_t);
__vpcs __f64x2_t _ZGVnN2v_log (__f64x2_t);
__vpcs __f64x2_t _ZGVnN2v_sin (__f64x2_t);
__vpcs __f64x2_t _ZGVnN2v_tan (__f64x2_t);
# undef __ADVSIMD_VEC_MATH_SUPPORTED
#endif /* __ADVSIMD_VEC_MATH_SUPPORTED */
@ -68,11 +70,13 @@ __sv_f32_t _ZGVsMxv_cosf (__sv_f32_t, __sv_bool_t);
__sv_f32_t _ZGVsMxv_expf (__sv_f32_t, __sv_bool_t);
__sv_f32_t _ZGVsMxv_logf (__sv_f32_t, __sv_bool_t);
__sv_f32_t _ZGVsMxv_sinf (__sv_f32_t, __sv_bool_t);
__sv_f32_t _ZGVsMxv_tanf (__sv_f32_t, __sv_bool_t);
__sv_f64_t _ZGVsMxv_cos (__sv_f64_t, __sv_bool_t);
__sv_f64_t _ZGVsMxv_exp (__sv_f64_t, __sv_bool_t);
__sv_f64_t _ZGVsMxv_log (__sv_f64_t, __sv_bool_t);
__sv_f64_t _ZGVsMxv_sin (__sv_f64_t, __sv_bool_t);
__sv_f64_t _ZGVsMxv_tan (__sv_f64_t, __sv_bool_t);
# undef __SVE_VEC_MATH_SUPPORTED
#endif /* __SVE_VEC_MATH_SUPPORTED */

36
sysdeps/aarch64/fpu/poly_advsimd_f32.h Normal file
View File

@ -0,0 +1,36 @@
/* Helpers for evaluating polynomials on single-precision AdvSIMD input, using
various schemes.
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 AARCH64_FPU_POLY_ADVSIMD_F32_H
#define AARCH64_FPU_POLY_ADVSIMD_F32_H
#include <arm_neon.h>
/* Wrap AdvSIMD f32 helpers: evaluation of some scheme/order has form:
v_[scheme]_[order]_f32. */
#define VTYPE float32x4_t
#define FMA(x, y, z) vfmaq_f32 (z, x, y)
#define VWRAP(f) v_##f##_f32
#include "poly_generic.h"
#undef VWRAP
#undef FMA
#undef VTYPE
#endif

36
sysdeps/aarch64/fpu/poly_advsimd_f64.h Normal file
View File

@ -0,0 +1,36 @@
/* Helpers for evaluating polynomials on double-precision AdvSIMD input, using
various schemes.
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 AARCH64_FPU_POLY_ADVSIMD_F64_H
#define AARCH64_FPU_POLY_ADVSIMD_F64_H
#include <arm_neon.h>
/* Wrap AdvSIMD f64 helpers: evaluation of some scheme/order has form:
v_[scheme]_[order]_f64. */
#define VTYPE float64x2_t
#define FMA(x, y, z) vfmaq_f64 (z, x, y)
#define VWRAP(f) v_##f##_f64
#include "poly_generic.h"
#undef VWRAP
#undef FMA
#undef VTYPE
#endif

285
sysdeps/aarch64/fpu/poly_generic.h Normal file
View File

@ -0,0 +1,285 @@
/* Generic helpers for evaluating polynomials with various schemes.
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 VTYPE
# error Cannot use poly_generic without defining VTYPE
#endif
#ifndef VWRAP
# error Cannot use poly_generic without defining VWRAP
#endif
#ifndef FMA
# error Cannot use poly_generic without defining FMA
#endif
static inline VTYPE VWRAP (pairwise_poly_3) (VTYPE x, VTYPE x2,
const VTYPE *poly)
{
/* At order 3, Estrin and Pairwise Horner are identical. */
VTYPE p01 = FMA (poly[1], x, poly[0]);
VTYPE p23 = FMA (poly[3], x, poly[2]);
return FMA (p23, x2, p01);
}
static inline VTYPE VWRAP (estrin_4) (VTYPE x, VTYPE x2, VTYPE x4,
const VTYPE *poly)
{
VTYPE p03 = VWRAP (pairwise_poly_3) (x, x2, poly);
return FMA (poly[4], x4, p03);
}
static inline VTYPE VWRAP (estrin_5) (VTYPE x, VTYPE x2, VTYPE x4,
const VTYPE *poly)
{
VTYPE p03 = VWRAP (pairwise_poly_3) (x, x2, poly);
VTYPE p45 = FMA (poly[5], x, poly[4]);
return FMA (p45, x4, p03);
}
static inline VTYPE VWRAP (estrin_6) (VTYPE x, VTYPE x2, VTYPE x4,
const VTYPE *poly)
{
VTYPE p03 = VWRAP (pairwise_poly_3) (x, x2, poly);
VTYPE p45 = FMA (poly[5], x, poly[4]);
VTYPE p46 = FMA (poly[6], x2, p45);
return FMA (p46, x4, p03);
}
static inline VTYPE VWRAP (estrin_7) (VTYPE x, VTYPE x2, VTYPE x4,
const VTYPE *poly)
{
VTYPE p03 = VWRAP (pairwise_poly_3) (x, x2, poly);
VTYPE p47 = VWRAP (pairwise_poly_3) (x, x2, poly + 4);
return FMA (p47, x4, p03);
}
static inline VTYPE VWRAP (estrin_8) (VTYPE x, VTYPE x2, VTYPE x4, VTYPE x8,
const VTYPE *poly)
{
return FMA (poly[8], x8, VWRAP (estrin_7) (x, x2, x4, poly));
}
static inline VTYPE VWRAP (estrin_9) (VTYPE x, VTYPE x2, VTYPE x4, VTYPE x8,
const VTYPE *poly)
{
VTYPE p89 = FMA (poly[9], x, poly[8]);
return FMA (p89, x8, VWRAP (estrin_7) (x, x2, x4, poly));
}
static inline VTYPE VWRAP (estrin_10) (VTYPE x, VTYPE x2, VTYPE x4, VTYPE x8,
const VTYPE *poly)
{
VTYPE p89 = FMA (poly[9], x, poly[8]);
VTYPE p8_10 = FMA (poly[10], x2, p89);
return FMA (p8_10, x8, VWRAP (estrin_7) (x, x2, x4, poly));
}
static inline VTYPE VWRAP (estrin_11) (VTYPE x, VTYPE x2, VTYPE x4, VTYPE x8,
const VTYPE *poly)
{
VTYPE p8_11 = VWRAP (pairwise_poly_3) (x, x2, poly + 8);
return FMA (p8_11, x8, VWRAP (estrin_7) (x, x2, x4, poly));
}
static inline VTYPE VWRAP (estrin_12) (VTYPE x, VTYPE x2, VTYPE x4, VTYPE x8,
const VTYPE *poly)
{
return FMA (VWRAP (estrin_4) (x, x2, x4, poly + 8), x8,
VWRAP (estrin_7) (x, x2, x4, poly));
}
static inline VTYPE VWRAP (estrin_13) (VTYPE x, VTYPE x2, VTYPE x4, VTYPE x8,
const VTYPE *poly)
{
return FMA (VWRAP (estrin_5) (x, x2, x4, poly + 8), x8,
VWRAP (estrin_7) (x, x2, x4, poly));
}
static inline VTYPE VWRAP (estrin_14) (VTYPE x, VTYPE x2, VTYPE x4, VTYPE x8,
const VTYPE *poly)
{
return FMA (VWRAP (estrin_6) (x, x2, x4, poly + 8), x8,
VWRAP (estrin_7) (x, x2, x4, poly));
}
static inline VTYPE VWRAP (estrin_15) (VTYPE x, VTYPE x2, VTYPE x4, VTYPE x8,
const VTYPE *poly)
{
return FMA (VWRAP (estrin_7) (x, x2, x4, poly + 8), x8,
VWRAP (estrin_7) (x, x2, x4, poly));
}
static inline VTYPE VWRAP (estrin_16) (VTYPE x, VTYPE x2, VTYPE x4, VTYPE x8,
VTYPE x16, const VTYPE *poly)
{
return FMA (poly[16], x16, VWRAP (estrin_15) (x, x2, x4, x8, poly));
}
static inline VTYPE VWRAP (estrin_17) (VTYPE x, VTYPE x2, VTYPE x4, VTYPE x8,
VTYPE x16, const VTYPE *poly)
{
VTYPE p16_17 = FMA (poly[17], x, poly[16]);
return FMA (p16_17, x16, VWRAP (estrin_15) (x, x2, x4, x8, poly));
}
static inline VTYPE VWRAP (estrin_18) (VTYPE x, VTYPE x2, VTYPE x4, VTYPE x8,
VTYPE x16, const VTYPE *poly)
{
VTYPE p16_17 = FMA (poly[17], x, poly[16]);
VTYPE p16_18 = FMA (poly[18], x2, p16_17);
return FMA (p16_18, x16, VWRAP (estrin_15) (x, x2, x4, x8, poly));
}
static inline VTYPE VWRAP (estrin_19) (VTYPE x, VTYPE x2, VTYPE x4, VTYPE x8,
VTYPE x16, const VTYPE *poly)
{
VTYPE p16_19 = VWRAP (pairwise_poly_3) (x, x2, poly + 16);
return FMA (p16_19, x16, VWRAP (estrin_15) (x, x2, x4, x8, poly));
}
static inline VTYPE VWRAP (horner_3) (VTYPE x, const VTYPE *poly)
{
VTYPE p = FMA (poly[3], x, poly[2]);
p = FMA (x, p, poly[1]);
p = FMA (x, p, poly[0]);
return p;
}
static inline VTYPE VWRAP (horner_4) (VTYPE x, const VTYPE *poly)
{
VTYPE p = FMA (poly[4], x, poly[3]);
p = FMA (x, p, poly[2]);
p = FMA (x, p, poly[1]);
p = FMA (x, p, poly[0]);
return p;
}
static inline VTYPE VWRAP (horner_5) (VTYPE x, const VTYPE *poly)
{
return FMA (x, VWRAP (horner_4) (x, poly + 1), poly[0]);
}
static inline VTYPE VWRAP (horner_6) (VTYPE x, const VTYPE *poly)
{
return FMA (x, VWRAP (horner_5) (x, poly + 1), poly[0]);
}
static inline VTYPE VWRAP (horner_7) (VTYPE x, const VTYPE *poly)
{
return FMA (x, VWRAP (horner_6) (x, poly + 1), poly[0]);
}
static inline VTYPE VWRAP (horner_8) (VTYPE x, const VTYPE *poly)
{
return FMA (x, VWRAP (horner_7) (x, poly + 1), poly[0]);
}
static inline VTYPE VWRAP (horner_9) (VTYPE x, const VTYPE *poly)
{
return FMA (x, VWRAP (horner_8) (x, poly + 1), poly[0]);
}
static inline VTYPE VWRAP (horner_10) (VTYPE x, const VTYPE *poly)
{
return FMA (x, VWRAP (horner_9) (x, poly + 1), poly[0]);
}
static inline VTYPE VWRAP (horner_11) (VTYPE x, const VTYPE *poly)
{
return FMA (x, VWRAP (horner_10) (x, poly + 1), poly[0]);
}
static inline VTYPE VWRAP (horner_12) (VTYPE x, const VTYPE *poly)
{
return FMA (x, VWRAP (horner_11) (x, poly + 1), poly[0]);
}
static inline VTYPE VWRAP (pw_horner_4) (VTYPE x, VTYPE x2, const VTYPE *poly)
{
VTYPE p01 = FMA (poly[1], x, poly[0]);
VTYPE p23 = FMA (poly[3], x, poly[2]);
VTYPE p;
p = FMA (x2, poly[4], p23);
p = FMA (x2, p, p01);
return p;
}
static inline VTYPE VWRAP (pw_horner_5) (VTYPE x, VTYPE x2, const VTYPE *poly)
{
VTYPE p01 = FMA (poly[1], x, poly[0]);
VTYPE p23 = FMA (poly[3], x, poly[2]);
VTYPE p45 = FMA (poly[5], x, poly[4]);
VTYPE p;
p = FMA (x2, p45, p23);
p = FMA (x2, p, p01);
return p;
}
static inline VTYPE VWRAP (pw_horner_6) (VTYPE x, VTYPE x2, const VTYPE *poly)
{
VTYPE p26 = VWRAP (pw_horner_4) (x, x2, poly + 2);
VTYPE p01 = FMA (poly[1], x, poly[0]);
return FMA (x2, p26, p01);
}
static inline VTYPE VWRAP (pw_horner_7) (VTYPE x, VTYPE x2, const VTYPE *poly)
{
VTYPE p27 = VWRAP (pw_horner_5) (x, x2, poly + 2);
VTYPE p01 = FMA (poly[1], x, poly[0]);
return FMA (x2, p27, p01);
}
static inline VTYPE VWRAP (pw_horner_8) (VTYPE x, VTYPE x2, const VTYPE *poly)
{
VTYPE p28 = VWRAP (pw_horner_6) (x, x2, poly + 2);
VTYPE p01 = FMA (poly[1], x, poly[0]);
return FMA (x2, p28, p01);
}
static inline VTYPE VWRAP (pw_horner_9) (VTYPE x, VTYPE x2, const VTYPE *poly)
{
VTYPE p29 = VWRAP (pw_horner_7) (x, x2, poly + 2);
VTYPE p01 = FMA (poly[1], x, poly[0]);
return FMA (x2, p29, p01);
}
static inline VTYPE VWRAP (pw_horner_10) (VTYPE x, VTYPE x2, const VTYPE *poly)
{
VTYPE p2_10 = VWRAP (pw_horner_8) (x, x2, poly + 2);
VTYPE p01 = FMA (poly[1], x, poly[0]);
return FMA (x2, p2_10, p01);
}
static inline VTYPE VWRAP (pw_horner_11) (VTYPE x, VTYPE x2, const VTYPE *poly)
{
VTYPE p2_11 = VWRAP (pw_horner_9) (x, x2, poly + 2);
VTYPE p01 = FMA (poly[1], x, poly[0]);
return FMA (x2, p2_11, p01);
}
static inline VTYPE VWRAP (pw_horner_12) (VTYPE x, VTYPE x2, const VTYPE *poly)
{
VTYPE p2_12 = VWRAP (pw_horner_10) (x, x2, poly + 2);
VTYPE p01 = FMA (poly[1], x, poly[0]);
return FMA (x2, p2_12, p01);
}
static inline VTYPE VWRAP (pw_horner_13) (VTYPE x, VTYPE x2, const VTYPE *poly)
{
VTYPE p2_13 = VWRAP (pw_horner_11) (x, x2, poly + 2);
VTYPE p01 = FMA (poly[1], x, poly[0]);
return FMA (x2, p2_13, p01);
}
static inline VTYPE VWRAP (pw_horner_14) (VTYPE x, VTYPE x2, const VTYPE *poly)
{
VTYPE p2_14 = VWRAP (pw_horner_12) (x, x2, poly + 2);
VTYPE p01 = FMA (poly[1], x, poly[0]);
return FMA (x2, p2_14, p01);
}
static inline VTYPE VWRAP (pw_horner_15) (VTYPE x, VTYPE x2, const VTYPE *poly)
{
VTYPE p2_15 = VWRAP (pw_horner_13) (x, x2, poly + 2);
VTYPE p01 = FMA (poly[1], x, poly[0]);
return FMA (x2, p2_15, p01);
}
static inline VTYPE VWRAP (pw_horner_16) (VTYPE x, VTYPE x2, const VTYPE *poly)
{
VTYPE p2_16 = VWRAP (pw_horner_14) (x, x2, poly + 2);
VTYPE p01 = FMA (poly[1], x, poly[0]);
return FMA (x2, p2_16, p01);
}
static inline VTYPE VWRAP (pw_horner_17) (VTYPE x, VTYPE x2, const VTYPE *poly)
{
VTYPE p2_17 = VWRAP (pw_horner_15) (x, x2, poly + 2);
VTYPE p01 = FMA (poly[1], x, poly[0]);
return FMA (x2, p2_17, p01);
}
static inline VTYPE VWRAP (pw_horner_18) (VTYPE x, VTYPE x2, const VTYPE *poly)
{
VTYPE p2_18 = VWRAP (pw_horner_16) (x, x2, poly + 2);
VTYPE p01 = FMA (poly[1], x, poly[0]);
return FMA (x2, p2_18, p01);
}

38
sysdeps/aarch64/fpu/poly_sve_f32.h Normal file
View File

@ -0,0 +1,38 @@
/* Helpers for evaluating polynomials on single-precision SVE input, using
various schemes.
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 AARCH64_FPU_POLY_SVE_F32_H
#define AARCH64_FPU_POLY_SVE_F32_H
#include <arm_sve.h>
/* Wrap SVE f32 helpers: evaluation of some scheme/order has form:
sv_[scheme]_[order]_f32_x. */
#define VTYPE svfloat32_t
#define STYPE float
#define VWRAP(f) sv_##f##_f32_x
#define DUP svdup_n_f32
#include "poly_sve_generic.h"
#undef DUP
#undef VWRAP
#undef STYPE
#undef VTYPE
#endif

38
sysdeps/aarch64/fpu/poly_sve_f64.h Normal file
View File

@ -0,0 +1,38 @@
/* Helpers for evaluating polynomials on double-precision SVE input, using
various schemes.
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 AARCH64_FPU_POLY_SVE_F64_H
#define AARCH64_FPU_POLY_SVE_F64_H
#include <arm_sve.h>
/* Wrap SVE f64 helpers: evaluation of some scheme/order has form:
sv_[scheme]_[order]_f64_x. */
#define VTYPE svfloat64_t
#define STYPE double
#define VWRAP(f) sv_##f##_f64_x
#define DUP svdup_n_f64
#include "poly_sve_generic.h"
#undef DUP
#undef VWRAP
#undef STYPE
#undef VTYPE
#endif

313
sysdeps/aarch64/fpu/poly_sve_generic.h Normal file
View File

@ -0,0 +1,313 @@
/* Helpers for evaluating polynomials with various schemes - specific to SVE
but precision-agnostic.
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 VTYPE
# error Cannot use poly_generic without defining VTYPE
#endif
#ifndef STYPE
# error Cannot use poly_generic without defining STYPE
#endif
#ifndef VWRAP
# error Cannot use poly_generic without defining VWRAP
#endif
#ifndef DUP
# error Cannot use poly_generic without defining DUP
#endif
static inline VTYPE VWRAP (pairwise_poly_3) (svbool_t pg, VTYPE x, VTYPE x2,
const STYPE *poly)
{
/* At order 3, Estrin and Pairwise Horner are identical. */
VTYPE p01 = svmla_x (pg, DUP (poly[0]), x, poly[1]);
VTYPE p23 = svmla_x (pg, DUP (poly[2]), x, poly[3]);
return svmla_x (pg, p01, p23, x2);
}
static inline VTYPE VWRAP (estrin_4) (svbool_t pg, VTYPE x, VTYPE x2, VTYPE x4,
const STYPE *poly)
{
VTYPE p03 = VWRAP (pairwise_poly_3) (pg, x, x2, poly);
return svmla_x (pg, p03, x4, poly[4]);
}
static inline VTYPE VWRAP (estrin_5) (svbool_t pg, VTYPE x, VTYPE x2, VTYPE x4,
const STYPE *poly)
{
VTYPE p03 = VWRAP (pairwise_poly_3) (pg, x, x2, poly);
VTYPE p45 = svmla_x (pg, DUP (poly[4]), x, poly[5]);
return svmla_x (pg, p03, p45, x4);
}
static inline VTYPE VWRAP (estrin_6) (svbool_t pg, VTYPE x, VTYPE x2, VTYPE x4,
const STYPE *poly)
{
VTYPE p03 = VWRAP (pairwise_poly_3) (pg, x, x2, poly);
VTYPE p45 = svmla_x (pg, DUP (poly[4]), x, poly[5]);
VTYPE p46 = svmla_x (pg, p45, x, poly[6]);
return svmla_x (pg, p03, p46, x4);
}
static inline VTYPE VWRAP (estrin_7) (svbool_t pg, VTYPE x, VTYPE x2, VTYPE x4,
const STYPE *poly)
{
VTYPE p03 = VWRAP (pairwise_poly_3) (pg, x, x2, poly);
VTYPE p47 = VWRAP (pairwise_poly_3) (pg, x, x2, poly + 4);
return svmla_x (pg, p03, p47, x4);
}
static inline VTYPE VWRAP (estrin_8) (svbool_t pg, VTYPE x, VTYPE x2, VTYPE x4,
VTYPE x8, const STYPE *poly)
{
return svmla_x (pg, VWRAP (estrin_7) (pg, x, x2, x4, poly), x8, poly[8]);
}
static inline VTYPE VWRAP (estrin_9) (svbool_t pg, VTYPE x, VTYPE x2, VTYPE x4,
VTYPE x8, const STYPE *poly)
{
VTYPE p89 = svmla_x (pg, DUP (poly[8]), x, poly[9]);
return svmla_x (pg, VWRAP (estrin_7) (pg, x, x2, x4, poly), p89, x8);
}
static inline VTYPE VWRAP (estrin_10) (svbool_t pg, VTYPE x, VTYPE x2,
VTYPE x4, VTYPE x8, const STYPE *poly)
{
VTYPE p89 = svmla_x (pg, DUP (poly[8]), x, poly[9]);
VTYPE p8_10 = svmla_x (pg, p89, x2, poly[10]);
return svmla_x (pg, VWRAP (estrin_7) (pg, x, x2, x4, poly), p8_10, x8);
}
static inline VTYPE VWRAP (estrin_11) (svbool_t pg, VTYPE x, VTYPE x2,
VTYPE x4, VTYPE x8, const STYPE *poly)
{
VTYPE p8_11 = VWRAP (pairwise_poly_3) (pg, x, x2, poly + 8);
return svmla_x (pg, VWRAP (estrin_7) (pg, x, x2, x4, poly), p8_11, x8);
}
static inline VTYPE VWRAP (estrin_12) (svbool_t pg, VTYPE x, VTYPE x2,
VTYPE x4, VTYPE x8, const STYPE *poly)
{
return svmla_x (pg, VWRAP (estrin_7) (pg, x, x2, x4, poly),
VWRAP (estrin_4) (pg, x, x2, x4, poly + 8), x8);
}
static inline VTYPE VWRAP (estrin_13) (svbool_t pg, VTYPE x, VTYPE x2,
VTYPE x4, VTYPE x8, const STYPE *poly)
{
return svmla_x (pg, VWRAP (estrin_7) (pg, x, x2, x4, poly),
VWRAP (estrin_5) (pg, x, x2, x4, poly + 8), x8);
}
static inline VTYPE VWRAP (estrin_14) (svbool_t pg, VTYPE x, VTYPE x2,
VTYPE x4, VTYPE x8, const STYPE *poly)
{
return svmla_x (pg, VWRAP (estrin_7) (pg, x, x2, x4, poly),
VWRAP (estrin_6) (pg, x, x2, x4, poly + 8), x8);
}
static inline VTYPE VWRAP (estrin_15) (svbool_t pg, VTYPE x, VTYPE x2,
VTYPE x4, VTYPE x8, const STYPE *poly)
{
return svmla_x (pg, VWRAP (estrin_7) (pg, x, x2, x4, poly),
VWRAP (estrin_7) (pg, x, x2, x4, poly + 8), x8);
}
static inline VTYPE VWRAP (estrin_16) (svbool_t pg, VTYPE x, VTYPE x2,
VTYPE x4, VTYPE x8, VTYPE x16,
const STYPE *poly)
{
return svmla_x (pg, VWRAP (estrin_15) (pg, x, x2, x4, x8, poly), x16,
poly[16]);
}
static inline VTYPE VWRAP (estrin_17) (svbool_t pg, VTYPE x, VTYPE x2,
VTYPE x4, VTYPE x8, VTYPE x16,
const STYPE *poly)
{
VTYPE p16_17 = svmla_x (pg, DUP (poly[16]), x, poly[17]);
return svmla_x (pg, VWRAP (estrin_15) (pg, x, x2, x4, x8, poly), p16_17,
x16);
}
static inline VTYPE VWRAP (estrin_18) (svbool_t pg, VTYPE x, VTYPE x2,
VTYPE x4, VTYPE x8, VTYPE x16,
const STYPE *poly)
{
VTYPE p16_17 = svmla_x (pg, DUP (poly[16]), x, poly[17]);
VTYPE p16_18 = svmla_x (pg, p16_17, x2, poly[18]);
return svmla_x (pg, VWRAP (estrin_15) (pg, x, x2, x4, x8, poly), p16_18,
x16);
}
static inline VTYPE VWRAP (estrin_19) (svbool_t pg, VTYPE x, VTYPE x2,
VTYPE x4, VTYPE x8, VTYPE x16,
const STYPE *poly)
{
return svmla_x (pg, VWRAP (estrin_15) (pg, x, x2, x4, x8, poly),
VWRAP (pairwise_poly_3) (pg, x, x2, poly + 16), x16);
}
static inline VTYPE VWRAP (horner_3) (svbool_t pg, VTYPE x, const STYPE *poly)
{
VTYPE p = svmla_x (pg, DUP (poly[2]), x, poly[3]);
p = svmad_x (pg, x, p, poly[1]);
p = svmad_x (pg, x, p, poly[0]);
return p;
}
static inline VTYPE VWRAP (horner_4) (svbool_t pg, VTYPE x, const STYPE *poly)
{
VTYPE p = svmla_x (pg, DUP (poly[3]), x, poly[4]);
p = svmad_x (pg, x, p, poly[2]);
p = svmad_x (pg, x, p, poly[1]);
p = svmad_x (pg, x, p, poly[0]);
return p;
}
static inline VTYPE VWRAP (horner_5) (svbool_t pg, VTYPE x, const STYPE *poly)
{
return svmad_x (pg, x, VWRAP (horner_4) (pg, x, poly + 1), poly[0]);
}
static inline VTYPE VWRAP (horner_6) (svbool_t pg, VTYPE x, const STYPE *poly)
{
return svmad_x (pg, x, VWRAP (horner_5) (pg, x, poly + 1), poly[0]);
}
static inline VTYPE VWRAP (horner_7) (svbool_t pg, VTYPE x, const STYPE *poly)
{
return svmad_x (pg, x, VWRAP (horner_6) (pg, x, poly + 1), poly[0]);
}
static inline VTYPE VWRAP (horner_8) (svbool_t pg, VTYPE x, const STYPE *poly)
{
return svmad_x (pg, x, VWRAP (horner_7) (pg, x, poly + 1), poly[0]);
}
static inline VTYPE VWRAP (horner_9) (svbool_t pg, VTYPE x, const STYPE *poly)
{
return svmad_x (pg, x, VWRAP (horner_8) (pg, x, poly + 1), poly[0]);
}
static inline VTYPE
sv_horner_10_f32_x (svbool_t pg, VTYPE x, const STYPE *poly)
{
return svmad_x (pg, x, VWRAP (horner_9) (pg, x, poly + 1), poly[0]);
}
static inline VTYPE
sv_horner_11_f32_x (svbool_t pg, VTYPE x, const STYPE *poly)
{
return svmad_x (pg, x, sv_horner_10_f32_x (pg, x, poly + 1), poly[0]);
}
static inline VTYPE
sv_horner_12_f32_x (svbool_t pg, VTYPE x, const STYPE *poly)
{
return svmad_x (pg, x, sv_horner_11_f32_x (pg, x, poly + 1), poly[0]);
}
static inline VTYPE VWRAP (pw_horner_4) (svbool_t pg, VTYPE x, VTYPE x2,
const STYPE *poly)
{
VTYPE p01 = svmla_x (pg, DUP (poly[0]), x, poly[1]);
VTYPE p23 = svmla_x (pg, DUP (poly[2]), x, poly[3]);
VTYPE p;
p = svmla_x (pg, p23, x2, poly[4]);
p = svmla_x (pg, p01, x2, p);
return p;
}
static inline VTYPE VWRAP (pw_horner_5) (svbool_t pg, VTYPE x, VTYPE x2,
const STYPE *poly)
{
VTYPE p01 = svmla_x (pg, DUP (poly[0]), x, poly[1]);
VTYPE p23 = svmla_x (pg, DUP (poly[2]), x, poly[3]);
VTYPE p45 = svmla_x (pg, DUP (poly[4]), x, poly[5]);
VTYPE p;
p = svmla_x (pg, p23, x2, p45);
p = svmla_x (pg, p01, x2, p);
return p;
}
static inline VTYPE VWRAP (pw_horner_6) (svbool_t pg, VTYPE x, VTYPE x2,
const STYPE *poly)
{
VTYPE p26 = VWRAP (pw_horner_4) (pg, x, x2, poly + 2);
VTYPE p01 = svmla_x (pg, DUP (poly[0]), x, poly[1]);
return svmla_x (pg, p01, x2, p26);
}
static inline VTYPE VWRAP (pw_horner_7) (svbool_t pg, VTYPE x, VTYPE x2,
const STYPE *poly)
{
VTYPE p27 = VWRAP (pw_horner_5) (pg, x, x2, poly + 2);
VTYPE p01 = svmla_x (pg, DUP (poly[0]), x, poly[1]);
return svmla_x (pg, p01, x2, p27);
}
static inline VTYPE VWRAP (pw_horner_8) (svbool_t pg, VTYPE x, VTYPE x2,
const STYPE *poly)
{
VTYPE p28 = VWRAP (pw_horner_6) (pg, x, x2, poly + 2);
VTYPE p01 = svmla_x (pg, DUP (poly[0]), x, poly[1]);
return svmla_x (pg, p01, x2, p28);
}
static inline VTYPE VWRAP (pw_horner_9) (svbool_t pg, VTYPE x, VTYPE x2,
const STYPE *poly)
{
VTYPE p29 = VWRAP (pw_horner_7) (pg, x, x2, poly + 2);
VTYPE p01 = svmla_x (pg, DUP (poly[0]), x, poly[1]);
return svmla_x (pg, p01, x2, p29);
}
static inline VTYPE VWRAP (pw_horner_10) (svbool_t pg, VTYPE x, VTYPE x2,
const STYPE *poly)
{
VTYPE p2_10 = VWRAP (pw_horner_8) (pg, x, x2, poly + 2);
VTYPE p01 = svmla_x (pg, DUP (poly[0]), x, poly[1]);
return svmla_x (pg, p01, x2, p2_10);
}
static inline VTYPE VWRAP (pw_horner_11) (svbool_t pg, VTYPE x, VTYPE x2,
const STYPE *poly)
{
VTYPE p2_11 = VWRAP (pw_horner_9) (pg, x, x2, poly + 2);
VTYPE p01 = svmla_x (pg, DUP (poly[0]), x, poly[1]);
return svmla_x (pg, p01, x2, p2_11);
}
static inline VTYPE VWRAP (pw_horner_12) (svbool_t pg, VTYPE x, VTYPE x2,
const STYPE *poly)
{
VTYPE p2_12 = VWRAP (pw_horner_10) (pg, x, x2, poly + 2);
VTYPE p01 = svmla_x (pg, DUP (poly[0]), x, poly[1]);
return svmla_x (pg, p01, x2, p2_12);
}
static inline VTYPE VWRAP (pw_horner_13) (svbool_t pg, VTYPE x, VTYPE x2,
const STYPE *poly)
{
VTYPE p2_13 = VWRAP (pw_horner_11) (pg, x, x2, poly + 2);
VTYPE p01 = svmla_x (pg, DUP (poly[0]), x, poly[1]);
return svmla_x (pg, p01, x2, p2_13);
}
static inline VTYPE VWRAP (pw_horner_14) (svbool_t pg, VTYPE x, VTYPE x2,
const STYPE *poly)
{
VTYPE p2_14 = VWRAP (pw_horner_12) (pg, x, x2, poly + 2);
VTYPE p01 = svmla_x (pg, DUP (poly[0]), x, poly[1]);
return svmla_x (pg, p01, x2, p2_14);
}
static inline VTYPE VWRAP (pw_horner_15) (svbool_t pg, VTYPE x, VTYPE x2,
const STYPE *poly)
{
VTYPE p2_15 = VWRAP (pw_horner_13) (pg, x, x2, poly + 2);
VTYPE p01 = svmla_x (pg, DUP (poly[0]), x, poly[1]);
return svmla_x (pg, p01, x2, p2_15);
}
static inline VTYPE VWRAP (pw_horner_16) (svbool_t pg, VTYPE x, VTYPE x2,
const STYPE *poly)
{
VTYPE p2_16 = VWRAP (pw_horner_14) (pg, x, x2, poly + 2);
VTYPE p01 = svmla_x (pg, DUP (poly[0]), x, poly[1]);
return svmla_x (pg, p01, x2, p2_16);
}
static inline VTYPE VWRAP (pw_horner_17) (svbool_t pg, VTYPE x, VTYPE x2,
const STYPE *poly)
{
VTYPE p2_17 = VWRAP (pw_horner_15) (pg, x, x2, poly + 2);
VTYPE p01 = svmla_x (pg, DUP (poly[0]), x, poly[1]);
return svmla_x (pg, p01, x2, p2_17);
}
static inline VTYPE VWRAP (pw_horner_18) (svbool_t pg, VTYPE x, VTYPE x2,
const STYPE *poly)
{
VTYPE p2_18 = VWRAP (pw_horner_16) (pg, x, x2, poly + 2);
VTYPE p01 = svmla_x (pg, DUP (poly[0]), x, poly[1]);
return svmla_x (pg, p01, x2, p2_18);
}

123
sysdeps/aarch64/fpu/tan_advsimd.c Normal file
View File

@ -0,0 +1,123 @@
/* Double-precision vector (Advanced SIMD) tan function
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/>. */
#include "v_math.h"
#include "poly_advsimd_f64.h"
static const struct data
{
float64x2_t poly[9];
float64x2_t half_pi_hi, half_pi_lo, two_over_pi, shift;
#if !WANT_SIMD_EXCEPT
float64x2_t range_val;
#endif
} data = {
/* Coefficients generated using FPMinimax. */
.poly = { V2 (0x1.5555555555556p-2), V2 (0x1.1111111110a63p-3),
V2 (0x1.ba1ba1bb46414p-5), V2 (0x1.664f47e5b5445p-6),
V2 (0x1.226e5e5ecdfa3p-7), V2 (0x1.d6c7ddbf87047p-9),
V2 (0x1.7ea75d05b583ep-10), V2 (0x1.289f22964a03cp-11),
V2 (0x1.4e4fd14147622p-12) },
.half_pi_hi = V2 (0x1.921fb54442d18p0),
.half_pi_lo = V2 (0x1.1a62633145c07p-54),
.two_over_pi = V2 (0x1.45f306dc9c883p-1),
.shift = V2 (0x1.8p52),
#if !WANT_SIMD_EXCEPT
.range_val = V2 (0x1p23),
#endif
};
#define RangeVal 0x4160000000000000 /* asuint64(0x1p23). */
#define TinyBound 0x3e50000000000000 /* asuint64(2^-26). */
#define Thresh 0x310000000000000 /* RangeVal - TinyBound. */
/* Special cases (fall back to scalar calls). */
static float64x2_t VPCS_ATTR NOINLINE
special_case (float64x2_t x)
{
return v_call_f64 (tan, x, x, v_u64 (-1));
}
/* Vector approximation for double-precision tan.
Maximum measured error is 3.48 ULP:
__v_tan(0x1.4457047ef78d8p+20) got -0x1.f6ccd8ecf7dedp+37
want -0x1.f6ccd8ecf7deap+37. */
float64x2_t VPCS_ATTR V_NAME_D1 (tan) (float64x2_t x)
{
const struct data *dat = ptr_barrier (&data);
/* Our argument reduction cannot calculate q with sufficient accuracy for very
large inputs. Fall back to scalar routine for all lanes if any are too
large, or Inf/NaN. If fenv exceptions are expected, also fall back for tiny
input to avoid underflow. */
#if WANT_SIMD_EXCEPT
uint64x2_t iax = vreinterpretq_u64_f64 (vabsq_f64 (x));
/* iax - tiny_bound > range_val - tiny_bound. */
uint64x2_t special
= vcgtq_u64 (vsubq_u64 (iax, v_u64 (TinyBound)), v_u64 (Thresh));
if (__glibc_unlikely (v_any_u64 (special)))
return special_case (x);
#endif
/* q = nearest integer to 2 * x / pi. */
float64x2_t q
= vsubq_f64 (vfmaq_f64 (dat->shift, x, dat->two_over_pi), dat->shift);
int64x2_t qi = vcvtq_s64_f64 (q);
/* Use q to reduce x to r in [-pi/4, pi/4], by:
r = x - q * pi/2, in extended precision. */
float64x2_t r = x;
r = vfmsq_f64 (r, q, dat->half_pi_hi);
r = vfmsq_f64 (r, q, dat->half_pi_lo);
/* Further reduce r to [-pi/8, pi/8], to be reconstructed using double angle
formula. */
r = vmulq_n_f64 (r, 0.5);
/* Approximate tan(r) using order 8 polynomial.
tan(x) is odd, so polynomial has the form:
tan(x) ~= x + C0 * x^3 + C1 * x^5 + C3 * x^7 + ...
Hence we first approximate P(r) = C1 + C2 * r^2 + C3 * r^4 + ...
Then compute the approximation by:
tan(r) ~= r + r^3 * (C0 + r^2 * P(r)). */
float64x2_t r2 = vmulq_f64 (r, r), r4 = vmulq_f64 (r2, r2),
r8 = vmulq_f64 (r4, r4);
/* Offset coefficients to evaluate from C1 onwards. */
float64x2_t p = v_estrin_7_f64 (r2, r4, r8, dat->poly + 1);
p = vfmaq_f64 (dat->poly[0], p, r2);
p = vfmaq_f64 (r, r2, vmulq_f64 (p, r));
/* Recombination uses double-angle formula:
tan(2x) = 2 * tan(x) / (1 - (tan(x))^2)
and reciprocity around pi/2:
tan(x) = 1 / (tan(pi/2 - x))
to assemble result using change-of-sign and conditional selection of
numerator/denominator, dependent on odd/even-ness of q (hence quadrant). */
float64x2_t n = vfmaq_f64 (v_f64 (-1), p, p);
float64x2_t d = vaddq_f64 (p, p);
uint64x2_t no_recip = vtstq_u64 (vreinterpretq_u64_s64 (qi), v_u64 (1));
#if !WANT_SIMD_EXCEPT
uint64x2_t special = vceqzq_u64 (vcaleq_f64 (x, dat->range_val));
if (__glibc_unlikely (v_any_u64 (special)))
return special_case (x);
#endif
return vdivq_f64 (vbslq_f64 (no_recip, n, vnegq_f64 (d)),
vbslq_f64 (no_recip, d, n));
}

104
sysdeps/aarch64/fpu/tan_sve.c Normal file
View File

@ -0,0 +1,104 @@
/* Double-precision vector (SVE) tan function
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/>. */
#include "sv_math.h"
#include "poly_sve_f64.h"
static const struct data
{
double poly[9];
double half_pi_hi, half_pi_lo, inv_half_pi, range_val, shift;
} data = {
/* Polynomial generated with FPMinimax. */
.poly = { 0x1.5555555555556p-2, 0x1.1111111110a63p-3, 0x1.ba1ba1bb46414p-5,
0x1.664f47e5b5445p-6, 0x1.226e5e5ecdfa3p-7, 0x1.d6c7ddbf87047p-9,
0x1.7ea75d05b583ep-10, 0x1.289f22964a03cp-11,
0x1.4e4fd14147622p-12, },
.half_pi_hi = 0x1.921fb54442d18p0,
.half_pi_lo = 0x1.1a62633145c07p-54,
.inv_half_pi = 0x1.45f306dc9c883p-1,
.range_val = 0x1p23,
.shift = 0x1.8p52,
};
static svfloat64_t NOINLINE
special_case (svfloat64_t x, svfloat64_t y, svbool_t special)
{
return sv_call_f64 (tan, x, y, special);
}
/* Vector approximation for double-precision tan.
Maximum measured error is 3.48 ULP:
_ZGVsMxv_tan(0x1.4457047ef78d8p+20) got -0x1.f6ccd8ecf7dedp+37
want -0x1.f6ccd8ecf7deap+37. */
svfloat64_t SV_NAME_D1 (tan) (svfloat64_t x, svbool_t pg)
{
const struct data *dat = ptr_barrier (&data);
/* Invert condition to catch NaNs and Infs as well as large values. */
svbool_t special = svnot_z (pg, svaclt (pg, x, dat->range_val));
/* q = nearest integer to 2 * x / pi. */
svfloat64_t shift = sv_f64 (dat->shift);
svfloat64_t q = svmla_x (pg, shift, x, dat->inv_half_pi);
q = svsub_x (pg, q, shift);
svint64_t qi = svcvt_s64_x (pg, q);
/* Use q to reduce x to r in [-pi/4, pi/4], by:
r = x - q * pi/2, in extended precision. */
svfloat64_t r = x;
svfloat64_t half_pi = svld1rq (svptrue_b64 (), &dat->half_pi_hi);
r = svmls_lane (r, q, half_pi, 0);
r = svmls_lane (r, q, half_pi, 1);
/* Further reduce r to [-pi/8, pi/8], to be reconstructed using double angle
formula. */
r = svmul_x (pg, r, 0.5);
/* Approximate tan(r) using order 8 polynomial.
tan(x) is odd, so polynomial has the form:
tan(x) ~= x + C0 * x^3 + C1 * x^5 + C3 * x^7 + ...
Hence we first approximate P(r) = C1 + C2 * r^2 + C3 * r^4 + ...
Then compute the approximation by:
tan(r) ~= r + r^3 * (C0 + r^2 * P(r)). */
svfloat64_t r2 = svmul_x (pg, r, r);
svfloat64_t r4 = svmul_x (pg, r2, r2);
svfloat64_t r8 = svmul_x (pg, r4, r4);
/* Use offset version coeff array by 1 to evaluate from C1 onwards. */
svfloat64_t p = sv_estrin_7_f64_x (pg, r2, r4, r8, dat->poly + 1);
p = svmad_x (pg, p, r2, dat->poly[0]);
p = svmla_x (pg, r, r2, svmul_x (pg, p, r));
/* Recombination uses double-angle formula:
tan(2x) = 2 * tan(x) / (1 - (tan(x))^2)
and reciprocity around pi/2:
tan(x) = 1 / (tan(pi/2 - x))
to assemble result using change-of-sign and conditional selection of
numerator/denominator dependent on odd/even-ness of q (hence quadrant). */
svbool_t use_recip
= svcmpeq (pg, svand_x (pg, svreinterpret_u64 (qi), 1), 0);
svfloat64_t n = svmad_x (pg, p, p, -1);
svfloat64_t d = svmul_x (pg, p, 2);
svfloat64_t swap = n;
n = svneg_m (n, use_recip, d);
d = svsel (use_recip, swap, d);
if (__glibc_unlikely (svptest_any (pg, special)))
return special_case (x, svdiv_x (svnot_z (pg, special), n, d), special);
return svdiv_x (pg, n, d);
}

129
sysdeps/aarch64/fpu/tanf_advsimd.c Normal file
View File

@ -0,0 +1,129 @@
/* Single-precision vector (Advanced SIMD) tan function
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/>. */
#include "v_math.h"
#include "poly_advsimd_f32.h"
static const struct data
{
float32x4_t poly[6];
float32x4_t neg_half_pi_1, neg_half_pi_2, neg_half_pi_3, two_over_pi, shift;
#if !WANT_SIMD_EXCEPT
float32x4_t range_val;
#endif
} data = {
/* Coefficients generated using FPMinimax. */
.poly = { V4 (0x1.55555p-2f), V4 (0x1.11166p-3f), V4 (0x1.b88a78p-5f),
V4 (0x1.7b5756p-6f), V4 (0x1.4ef4cep-8f), V4 (0x1.0e1e74p-7f) },
.neg_half_pi_1 = V4 (-0x1.921fb6p+0f),
.neg_half_pi_2 = V4 (0x1.777a5cp-25f),
.neg_half_pi_3 = V4 (0x1.ee59dap-50f),
.two_over_pi = V4 (0x1.45f306p-1f),
.shift = V4 (0x1.8p+23f),
#if !WANT_SIMD_EXCEPT
.range_val = V4 (0x1p15f),
#endif
};
#define RangeVal v_u32 (0x47000000) /* asuint32(0x1p15f). */
#define TinyBound v_u32 (0x30000000) /* asuint32 (0x1p-31f). */
#define Thresh v_u32 (0x16000000) /* asuint32(RangeVal) - TinyBound. */
/* Special cases (fall back to scalar calls). */
static float32x4_t VPCS_ATTR NOINLINE
special_case (float32x4_t x, float32x4_t y, uint32x4_t cmp)
{
return v_call_f32 (tanf, x, y, cmp);
}
/* Use a full Estrin scheme to evaluate polynomial. */
static inline float32x4_t
eval_poly (float32x4_t z, const struct data *d)
{
float32x4_t z2 = vmulq_f32 (z, z);
#if WANT_SIMD_EXCEPT
/* Tiny z (<= 0x1p-31) will underflow when calculating z^4. If fp exceptions
are to be triggered correctly, sidestep this by fixing such lanes to 0. */
uint32x4_t will_uflow
= vcleq_u32 (vreinterpretq_u32_f32 (vabsq_f32 (z)), TinyBound);
if (__glibc_unlikely (v_any_u32 (will_uflow)))
z2 = vbslq_f32 (will_uflow, v_f32 (0), z2);
#endif
float32x4_t z4 = vmulq_f32 (z2, z2);
return v_estrin_5_f32 (z, z2, z4, d->poly);
}
/* Fast implementation of AdvSIMD tanf.
Maximum error is 3.45 ULP:
__v_tanf(-0x1.e5f0cap+13) got 0x1.ff9856p-1
want 0x1.ff9850p-1. */
float32x4_t VPCS_ATTR V_NAME_F1 (tan) (float32x4_t x)
{
const struct data *d = ptr_barrier (&data);
float32x4_t special_arg = x;
/* iax >= RangeVal means x, if not inf or NaN, is too large to perform fast
regression. */
#if WANT_SIMD_EXCEPT
uint32x4_t iax = vreinterpretq_u32_f32 (vabsq_f32 (x));
/* If fp exceptions are to be triggered correctly, also special-case tiny
input, as this will load to overflow later. Fix any special lanes to 1 to
prevent any exceptions being triggered. */
uint32x4_t special = vcgeq_u32 (vsubq_u32 (iax, TinyBound), Thresh);
if (__glibc_unlikely (v_any_u32 (special)))
x = vbslq_f32 (special, v_f32 (1.0f), x);
#else
/* Otherwise, special-case large and special values. */
uint32x4_t special = vcageq_f32 (x, d->range_val);
#endif
/* n = rint(x/(pi/2)). */
float32x4_t q = vfmaq_f32 (d->shift, d->two_over_pi, x);
float32x4_t n = vsubq_f32 (q, d->shift);
/* Determine if x lives in an interval, where |tan(x)| grows to infinity. */
uint32x4_t pred_alt = vtstq_u32 (vreinterpretq_u32_f32 (q), v_u32 (1));
/* r = x - n * (pi/2) (range reduction into -pi./4 .. pi/4). */
float32x4_t r;
r = vfmaq_f32 (x, d->neg_half_pi_1, n);
r = vfmaq_f32 (r, d->neg_half_pi_2, n);
r = vfmaq_f32 (r, d->neg_half_pi_3, n);
/* If x lives in an interval, where |tan(x)|
- is finite, then use a polynomial approximation of the form
tan(r) ~ r + r^3 * P(r^2) = r + r * r^2 * P(r^2).
- grows to infinity then use symmetries of tangent and the identity
tan(r) = cotan(pi/2 - r) to express tan(x) as 1/tan(-r). Finally, use
the same polynomial approximation of tan as above. */
/* Invert sign of r if odd quadrant. */
float32x4_t z = vmulq_f32 (r, vbslq_f32 (pred_alt, v_f32 (-1), v_f32 (1)));
/* Evaluate polynomial approximation of tangent on [-pi/4, pi/4]. */
float32x4_t z2 = vmulq_f32 (r, r);
float32x4_t p = eval_poly (z2, d);
float32x4_t y = vfmaq_f32 (z, vmulq_f32 (z, z2), p);
/* Compute reciprocal and apply if required. */
float32x4_t inv_y = vdivq_f32 (v_f32 (1.0f), y);
if (__glibc_unlikely (v_any_u32 (special)))
return special_case (special_arg, vbslq_f32 (pred_alt, inv_y, y), special);
return vbslq_f32 (pred_alt, inv_y, y);
}

118
sysdeps/aarch64/fpu/tanf_sve.c Normal file
View File

@ -0,0 +1,118 @@
/* Single-precision vector (SVE) tan function
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/>. */
#include "sv_math.h"
static const struct data
{
float pio2_1, pio2_2, pio2_3, invpio2;
float c1, c3, c5;
float c0, c2, c4, range_val, shift;
} data = {
/* Coefficients generated using:
poly = fpminimax((tan(sqrt(x))-sqrt(x))/x^(3/2),
deg,
[|single ...|],
[a*a;b*b]);
optimize relative error
final prec : 23 bits
deg : 5
a : 0x1p-126 ^ 2
b : ((pi) / 0x1p2) ^ 2
dirty rel error: 0x1.f7c2e4p-25
dirty abs error: 0x1.f7c2ecp-25. */
.c0 = 0x1.55555p-2, .c1 = 0x1.11166p-3,
.c2 = 0x1.b88a78p-5, .c3 = 0x1.7b5756p-6,
.c4 = 0x1.4ef4cep-8, .c5 = 0x1.0e1e74p-7,
.pio2_1 = 0x1.921fb6p+0f, .pio2_2 = -0x1.777a5cp-25f,
.pio2_3 = -0x1.ee59dap-50f, .invpio2 = 0x1.45f306p-1f,
.range_val = 0x1p15f, .shift = 0x1.8p+23f
};
static svfloat32_t NOINLINE
special_case (svfloat32_t x, svfloat32_t y, svbool_t cmp)
{
return sv_call_f32 (tanf, x, y, cmp);
}
/* Fast implementation of SVE tanf.
Maximum error is 3.45 ULP:
SV_NAME_F1 (tan)(-0x1.e5f0cap+13) got 0x1.ff9856p-1
want 0x1.ff9850p-1. */
svfloat32_t SV_NAME_F1 (tan) (svfloat32_t x, const svbool_t pg)
{
const struct data *d = ptr_barrier (&data);
/* Determine whether input is too large to perform fast regression. */
svbool_t cmp = svacge (pg, x, d->range_val);
svfloat32_t odd_coeffs = svld1rq (svptrue_b32 (), &d->c1);
svfloat32_t pi_vals = svld1rq (svptrue_b32 (), &d->pio2_1);
/* n = rint(x/(pi/2)). */
svfloat32_t q = svmla_lane (sv_f32 (d->shift), x, pi_vals, 3);
svfloat32_t n = svsub_x (pg, q, d->shift);
/* n is already a signed integer, simply convert it. */
svint32_t in = svcvt_s32_x (pg, n);
/* Determine if x lives in an interval, where |tan(x)| grows to infinity. */
svint32_t alt = svand_x (pg, in, 1);
svbool_t pred_alt = svcmpne (pg, alt, 0);
/* r = x - n * (pi/2) (range reduction into 0 .. pi/4). */
svfloat32_t r;
r = svmls_lane (x, n, pi_vals, 0);
r = svmls_lane (r, n, pi_vals, 1);
r = svmls_lane (r, n, pi_vals, 2);
/* If x lives in an interval, where |tan(x)|
- is finite, then use a polynomial approximation of the form
tan(r) ~ r + r^3 * P(r^2) = r + r * r^2 * P(r^2).
- grows to infinity then use symmetries of tangent and the identity
tan(r) = cotan(pi/2 - r) to express tan(x) as 1/tan(-r). Finally, use
the same polynomial approximation of tan as above. */
/* Perform additional reduction if required. */
svfloat32_t z = svneg_m (r, pred_alt, r);
/* Evaluate polynomial approximation of tangent on [-pi/4, pi/4],
using Estrin on z^2. */
svfloat32_t z2 = svmul_x (pg, z, z);
svfloat32_t p01 = svmla_lane (sv_f32 (d->c0), z2, odd_coeffs, 0);
svfloat32_t p23 = svmla_lane (sv_f32 (d->c2), z2, odd_coeffs, 1);
svfloat32_t p45 = svmla_lane (sv_f32 (d->c4), z2, odd_coeffs, 2);
svfloat32_t z4 = svmul_x (pg, z2, z2);
svfloat32_t p = svmla_x (pg, p01, z4, p23);
svfloat32_t z8 = svmul_x (pg, z4, z4);
p = svmla_x (pg, p, z8, p45);
svfloat32_t y = svmla_x (pg, z, p, svmul_x (pg, z, z2));
/* Transform result back, if necessary. */
svfloat32_t inv_y = svdivr_x (pg, y, 1.0f);
/* No need to pass pg to specialcase here since cmp is a strict subset,
guaranteed by the cmpge above. */
if (__glibc_unlikely (svptest_any (pg, cmp)))
return special_case (x, svsel (pred_alt, inv_y, y), cmp);
return svsel (pred_alt, inv_y, y);
}

1
sysdeps/aarch64/fpu/test-double-advsimd-wrappers.c
View File

@ -27,3 +27,4 @@ VPCS_VECTOR_WRAPPER (cos_advsimd, _ZGVnN2v_cos)
VPCS_VECTOR_WRAPPER (exp_advsimd, _ZGVnN2v_exp)
VPCS_VECTOR_WRAPPER (log_advsimd, _ZGVnN2v_log)
VPCS_VECTOR_WRAPPER (sin_advsimd, _ZGVnN2v_sin)
VPCS_VECTOR_WRAPPER (tan_advsimd, _ZGVnN2v_tan)

1
sysdeps/aarch64/fpu/test-double-sve-wrappers.c
View File

@ -36,3 +36,4 @@ SVE_VECTOR_WRAPPER (cos_sve, _ZGVsMxv_cos)
SVE_VECTOR_WRAPPER (exp_sve, _ZGVsMxv_exp)
SVE_VECTOR_WRAPPER (log_sve, _ZGVsMxv_log)
SVE_VECTOR_WRAPPER (sin_sve, _ZGVsMxv_sin)
SVE_VECTOR_WRAPPER (tan_sve, _ZGVsMxv_tan)

1
sysdeps/aarch64/fpu/test-float-advsimd-wrappers.c
View File

@ -27,3 +27,4 @@ VPCS_VECTOR_WRAPPER (cosf_advsimd, _ZGVnN4v_cosf)
VPCS_VECTOR_WRAPPER (expf_advsimd, _ZGVnN4v_expf)
VPCS_VECTOR_WRAPPER (logf_advsimd, _ZGVnN4v_logf)
VPCS_VECTOR_WRAPPER (sinf_advsimd, _ZGVnN4v_sinf)
VPCS_VECTOR_WRAPPER (tanf_advsimd, _ZGVnN4v_tanf)

1
sysdeps/aarch64/fpu/test-float-sve-wrappers.c
View File

@ -36,3 +36,4 @@ SVE_VECTOR_WRAPPER (cosf_sve, _ZGVsMxv_cosf)
SVE_VECTOR_WRAPPER (expf_sve, _ZGVsMxv_expf)
SVE_VECTOR_WRAPPER (logf_sve, _ZGVsMxv_logf)
SVE_VECTOR_WRAPPER (sinf_sve, _ZGVsMxv_sinf)
SVE_VECTOR_WRAPPER (tanf_sve, _ZGVsMxv_tanf)

8
sysdeps/aarch64/libm-test-ulps
View File

@ -1340,11 +1340,19 @@ Function: "tan":
float: 1
ldouble: 1
Function: "tan_advsimd":
double: 2
float: 2
Function: "tan_downward":
double: 1
float: 2
ldouble: 1
Function: "tan_sve":
double: 2
float: 2
Function: "tan_towardzero":
double: 1
float: 1

4
sysdeps/unix/sysv/linux/aarch64/libmvec.abilist
View File

@ -14,3 +14,7 @@ GLIBC_2.38 _ZGVsMxv_log F
GLIBC_2.38 _ZGVsMxv_logf F
GLIBC_2.38 _ZGVsMxv_sin F
GLIBC_2.38 _ZGVsMxv_sinf F
GLIBC_2.39 _ZGVnN2v_tan F
GLIBC_2.39 _ZGVnN4v_tanf F
GLIBC_2.39 _ZGVsMxv_tan F
GLIBC_2.39 _ZGVsMxv_tanf F