1 /* Single-precision vector (SVE) erf function
3 Copyright (C) 2024 Free Software Foundation, Inc.
4 This file is part of the GNU C Library.
6 The GNU C Library is free software; you can redistribute it and/or
7 modify it under the terms of the GNU Lesser General Public
8 License as published by the Free Software Foundation; either
9 version 2.1 of the License, or (at your option) any later version.
11 The GNU C Library is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 Lesser General Public License for more details.
16 You should have received a copy of the GNU Lesser General Public
17 License along with the GNU C Library; if not, see
18 <https://www.gnu.org/licenses/>. */
22 static const struct data
24 float min
, max
, scale
, shift
, third
;
26 .min
= 0x1.cp
-7f
, /* 1/64 - 1/512. */
27 .max
= 3.9375, /* 4 - 8/128. */
28 .scale
= 0x1.20dd76p
+0f
, /* 2/sqrt(pi). */
30 .third
= 0x1.555556p
-2f
, /* 1/3. */
33 #define SignMask (0x80000000)
35 /* Single-precision implementation of vector erf(x).
36 Approximation based on series expansion near x rounded to
37 nearest multiple of 1/128.
38 Let d = x - r, and scale = 2 / sqrt(pi) * exp(-r^2). For x near r,
40 erf(x) ~ erf(r) + scale * d * [1 - r * d - 1/3 * d^2]
42 Values of erf(r) and scale are read from lookup tables.
43 For |x| < 0x1.cp-7, the algorithm sets r = 0, erf(r) = 0, and scale = 2 /
44 sqrt(pi), so it simply boils down to a Taylor series expansion near 0. For
45 |x| > 3.9375, erf(|x|) rounds to 1.0f.
47 Maximum error on each interval:
48 - [0, 0x1.cp-7]: 1.93 ULP
49 _ZGVsMxv_erff(0x1.c373e6p-9) got 0x1.fd686cp-9 want 0x1.fd6868p-9
50 - [0x1.cp-7, 4.0]: 1.26 ULP
51 _ZGVsMxv_erff(0x1.1d002ep+0) got 0x1.c4eb9ap-1 want 0x1.c4eb98p-1. */
52 svfloat32_t
SV_NAME_F1 (erf
) (svfloat32_t x
, const svbool_t pg
)
54 const struct data
*dat
= ptr_barrier (&data
);
56 /* |x| > 1/64 - 1/512. */
57 svbool_t a_gt_min
= svacgt (pg
, x
, dat
->min
);
59 /* |x| >= 4.0 - 8/128. */
60 svbool_t a_ge_max
= svacge (pg
, x
, dat
->max
);
61 svfloat32_t a
= svabs_x (pg
, x
);
63 svfloat32_t shift
= sv_f32 (dat
->shift
);
64 svfloat32_t z
= svadd_x (pg
, a
, shift
);
66 = svsub_x (pg
, svreinterpret_u32 (z
), svreinterpret_u32 (shift
));
68 /* Saturate lookup index. */
69 i
= svsel (a_ge_max
, sv_u32 (512), i
);
71 /* r and erf(r) set to 0 for |x| below min. */
72 svfloat32_t r
= svsub_z (a_gt_min
, z
, shift
);
73 svfloat32_t erfr
= svld1_gather_index (a_gt_min
, __sv_erff_data
.erf
, i
);
75 /* scale set to 2/sqrt(pi) for |x| below min. */
76 svfloat32_t scale
= svld1_gather_index (a_gt_min
, __sv_erff_data
.scale
, i
);
77 scale
= svsel (a_gt_min
, scale
, sv_f32 (dat
->scale
));
79 /* erf(x) ~ erf(r) + scale * d * (1 - r * d + 1/3 * d^2). */
80 svfloat32_t d
= svsub_x (pg
, a
, r
);
81 svfloat32_t d2
= svmul_x (pg
, d
, d
);
82 svfloat32_t y
= svmla_x (pg
, r
, d
, dat
->third
);
83 y
= svmla_x (pg
, erfr
, scale
, svmls_x (pg
, d
, d2
, y
));
85 /* Solves the |x| = inf case. */
86 y
= svsel (a_ge_max
, sv_f32 (1.0f
), y
);
89 svuint32_t ix
= svreinterpret_u32 (x
);
90 svuint32_t iy
= svreinterpret_u32 (y
);
91 svuint32_t sign
= svand_x (pg
, ix
, SignMask
);
92 return svreinterpret_f32 (svorr_x (pg
, sign
, iy
));