@@ -134,7 +134,8 @@ type-double-routines := branred doasin dosincos mpa mpatan2 \
# float support
type-float-suffix := f
type-float-routines := math_errf e_exp2f_data e_logf_data \
- e_log2f_data e_powf_log2_data s_sincosf_data
+ e_log2f_data e_powf_log2_data s_sincosf_data \
+ e_reduce_pi
# _Float128 support
type-float128-suffix := f128
new file mode 100644
@@ -0,0 +1,186 @@
+/* Argument reduction modulo 2pi.
+ 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_config.h"
+
+/* Argument reduction mod 2pi: T[i] ~ 2^i mod (2*pi). */
+static const double T[128] = {
+ 0x1p+0,
+ 0x2p+0,
+ -0x2.487ed5110b462p+0,
+ 0x1.b7812aeef4b9fp+0,
+ -0x2.d97c7f3321d24p+0,
+ 0x9.585d6aac7a1a8p-4,
+ 0x1.2b0bad558f435p+0,
+ 0x2.56175aab1e86ap+0,
+ -0x1.9c501fbace38dp+0,
+ 0x3.0fde959b6ed46p+0,
+ -0x2.8c1a9da2d9d3cp-4,
+ -0x5.18353b45b3a78p-4,
+ -0xa.306a768b674fp-4,
+ -0x1.460d4ed16ce9ep+0,
+ -0x2.8c1a9da2d9d3cp+0,
+ 0x1.304999cb579e8p+0,
+ 0x2.60933396af3dp+0,
+ -0x1.87586de3accc3p+0,
+ -0x3.0eb0dbc759986p+0,
+ 0x2.b1d1d8258155p-4,
+ 0x5.63a3b04b02aap-4,
+ 0xa.c74760960554p-4,
+ 0x1.58e8ec12c0aa8p+0,
+ 0x2.b1d1d8258155p+0,
+ -0xe.4db24c6089bf8p-4,
+ -0x1.c9b6498c1137fp+0,
+ 0x2.b51241f8e8d64p+0,
+ -0xd.e5a511f3999bp-4,
+ -0x1.bcb4a23e73336p+0,
+ 0x2.cf15909424df4p+0,
+ -0xa.a53b3e8c1877p-4,
+ -0x1.54a767d1830eep+0,
+ -0x2.a94ecfa3061dcp+0,
+ 0xf.5e135caff0a78p-4,
+ 0x1.ebc26b95fe14fp+0,
+ -0x2.70f9fde50f1c4p+0,
+ 0x1.668ad946ed0dap+0,
+ 0x2.cd15b28dda1b4p+0,
+ -0xa.e536ff5570fbp-4,
+ -0x1.5ca6dfeaae1f6p+0,
+ -0x2.b94dbfd55c3ecp+0,
+ 0xd.5e3556652c89p-4,
+ 0x1.abc6aacca5912p+0,
+ -0x2.f0f17f77c023ep+0,
+ 0x6.69bd6218afe64p-4,
+ 0xc.d37ac4315fcc8p-4,
+ 0x1.9a6f58862bf99p+0,
+ -0x3.13a02404b352ep+0,
+ 0x2.13e8d07a4a046p-4,
+ 0x4.27d1a0f49408cp-4,
+ 0x8.4fa341e928118p-4,
+ 0x1.09f4683d25023p+0,
+ 0x2.13e8d07a4a046p+0,
+ -0x2.20ad341c773d4p+0,
+ 0x2.07246cd81ccb8p+0,
+ -0x2.3a35fb60d1af2p+0,
+ 0x1.d412de4f67e7ep+0,
+ -0x2.a05918723b764p+0,
+ 0x1.07cca42c94599p+0,
+ 0x2.0f99485928b32p+0,
+ -0x2.294c445eb9dfep+0,
+ 0x1.f5e64c5397865p+0,
+ -0x2.5cb23c69dc396p+0,
+ 0x1.8f1a5c3d52d34p+0,
+ 0x3.1e34b87aa5a68p+0,
+ -0xc.15641bbff90bp-8,
+ -0x1.82ac8377ff216p-4,
+ -0x3.055906effe42cp-4,
+ -0x6.0ab20ddffc858p-4,
+ -0xc.15641bbff90bp-4,
+ -0x1.82ac8377ff216p+0,
+ -0x3.055906effe42cp+0,
+ 0x3.dccc7310ec09ap-4,
+ 0x7.b998e621d8134p-4,
+ 0xf.7331cc43b0268p-4,
+ 0x1.ee6639887604dp+0,
+ -0x2.6bb262001f3c6p+0,
+ 0x1.711a1110cccd4p+0,
+ 0x2.e2342221999a8p+0,
+ -0x8.41690cdd81128p-4,
+ -0x1.082d219bb0225p+0,
+ -0x2.105a43376044ap+0,
+ 0x2.27ca4ea24abcep+0,
+ -0x1.f8ea37cc75cc6p+0,
+ 0x2.56aa65781fad6p+0,
+ -0x1.9b2a0a20cbeb6p+0,
+ 0x3.122ac0cf736f6p+0,
+ -0x2.4295372246764p-4,
+ -0x4.852a6e448cec8p-4,
+ -0x9.0a54dc8919d9p-4,
+ -0x1.214a9b91233b2p+0,
+ -0x2.4295372246764p+0,
+ 0x1.c35466cc7e59ap+0,
+ -0x2.c1d607780e92ep+0,
+ 0xc.4d2c620ee206p-4,
+ 0x1.89a58c41dc40cp+0,
+ 0x3.134b1883b8818p+0,
+ -0x2.1e8a4099a4324p-4,
+ -0x4.3d14813348648p-4,
+ -0x8.7a29026690c9p-4,
+ -0x1.0f45204cd2192p+0,
+ -0x2.1e8a4099a4324p+0,
+ 0x2.0b6a53ddc2e18p+0,
+ -0x2.31aa2d558583p+0,
+ 0x1.e52a7a6600402p+0,
+ -0x2.7e29e0450ac5cp+0,
+ 0x1.4c2b1486f5ba8p+0,
+ 0x2.9856290deb75p+0,
+ -0x1.17d282f5345bfp+0,
+ -0x2.2fa505ea68b7ep+0,
+ 0x1.e934c93c39d64p+0,
+ -0x2.761542989799ap+0,
+ 0x1.5c544fdfdc12fp+0,
+ 0x2.b8a89fbfb825ep+0,
+ -0xd.72d95919afa58p-4,
+ -0x1.ae5b2b2335f4bp+0,
+ 0x2.ebc87eca9f5cap+0,
+ -0x7.0edd77bcc8ccp-4,
+ -0xe.1dbaef799198p-4,
+ -0x1.c3b75def3233p+0,
+ 0x2.c1101932a6e02p+0,
+ -0xc.65ea2abbd85fp-4,
+ -0x1.8cbd45577b0bep+0,
+ -0x3.197a8aaef617cp+0,
+ 0x1.589bfb31f1687p-4,
+ 0x2.b137f663e2d0ep-4,
+ 0x5.626fecc7c5a1cp-4,
+ 0xa.c4dfd98f8b438p-4
+};
+
+/* Return h such that x - pi/4 mod (2*pi) ~ h,
+ with maximal error about 2^-45. */
+attribute_hidden double
+__math_reduce_mod_twopi (double x)
+{
+ double sign = 1, h;
+ if (x < 0)
+ {
+ x = -x;
+ sign = -1;
+ }
+ h = 0; /* Invariant is x+h mod (2*pi). */
+ while (x >= 4)
+ {
+ int i = ilogb (x);
+ x -= ldexp (1.0f, i);
+ h += T[i];
+ }
+ /* Add the remainder x. */
+ h += x;
+ /* Subtract pi/4. */
+ double piover4 = 0xc.90fdaa22168cp-4;
+ h -= piover4;
+ /* Reduce mod 2*pi. */
+ double twopi = 0x6.487ed5110b46p+0;
+ while (fabs (h) > twopi * 0.5)
+ {
+ if (h > 0)
+ h -= twopi;
+ else
+ h += twopi;
+ }
+ return sign * h;
+}
new file mode 100644
@@ -0,0 +1,19 @@
+/* Argument reduction modulo 2pi.
+ 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/>. */
+
+attribute_hidden double __math_reduce_mod_twopi (double);