| Message ID | 1513143801-3290-2-git-send-email-raji@linux.vnet.ibm.com |
|---|---|
| State | New |
| Headers | show |
| Series | New generic sincosf | expand |
On Wed, 13 Dec 2017, Rajalakshmi Srinivasaraghavan wrote: > This implementation is based on generic s_sinf.c and s_cosf.c. > Tested on s390x, powerpc64le and powerpc32. > > 2017-12-12 Rajalakshmi Srinivasaraghavan <raji@linux.vnet.ibm.com> > > * sysdeps/ieee754/flt-32/s_cosf.c: Move reduced() and > constants to s_sincos.h file. > * sysdeps/ieee754/flt-32/s_sinf.c: Likewise. > * sysdeps/ieee754/flt-32/s_sincosf.c: New implementation. > * sysdeps/ieee754/flt-32/s_sincos.h: New file. OK.
On Thu, Dec 14, 2017 at 4:29 PM, Joseph Myers <joseph@codesourcery.com> wrote: > On Wed, 13 Dec 2017, Rajalakshmi Srinivasaraghavan wrote: > >> This implementation is based on generic s_sinf.c and s_cosf.c. >> Tested on s390x, powerpc64le and powerpc32. >> >> 2017-12-12 Rajalakshmi Srinivasaraghavan <raji@linux.vnet.ibm.com> >> >> * sysdeps/ieee754/flt-32/s_cosf.c: Move reduced() and >> constants to s_sincos.h file. >> * sysdeps/ieee754/flt-32/s_sinf.c: Likewise. >> * sysdeps/ieee754/flt-32/s_sincosf.c: New implementation. >> * sysdeps/ieee754/flt-32/s_sincos.h: New file. > > OK. > Please put those static constants into a different file and mark them hidden so that they can be shared when they are compiled twice for FMA. Thanks.
On 12/15/2017 07:49 AM, H.J. Lu wrote: > On Thu, Dec 14, 2017 at 4:29 PM, Joseph Myers <joseph@codesourcery.com> wrote: >> On Wed, 13 Dec 2017, Rajalakshmi Srinivasaraghavan wrote: >> >>> This implementation is based on generic s_sinf.c and s_cosf.c. >>> Tested on s390x, powerpc64le and powerpc32. >>> >>> 2017-12-12 Rajalakshmi Srinivasaraghavan <raji@linux.vnet.ibm.com> >>> >>> * sysdeps/ieee754/flt-32/s_cosf.c: Move reduced() and >>> constants to s_sincos.h file. >>> * sysdeps/ieee754/flt-32/s_sinf.c: Likewise. >>> * sysdeps/ieee754/flt-32/s_sincosf.c: New implementation. >>> * sysdeps/ieee754/flt-32/s_sincos.h: New file. >> >> OK. >> > > Please put those static constants into a different file and mark them > hidden so that they can be shared when they are compiled twice > for FMA. > I assume you are asking changes similar to sysdeps/ieee754/flt-32/e_exp2f_data.c. Could you please explain the advantage of doing this or suggest any previous mail thread which I can refer to? > Thanks. >
On Fri, 15 Dec 2017, H.J. Lu wrote: > On Thu, Dec 14, 2017 at 4:29 PM, Joseph Myers <joseph@codesourcery.com> wrote: > > On Wed, 13 Dec 2017, Rajalakshmi Srinivasaraghavan wrote: > > > >> This implementation is based on generic s_sinf.c and s_cosf.c. > >> Tested on s390x, powerpc64le and powerpc32. > >> > >> 2017-12-12 Rajalakshmi Srinivasaraghavan <raji@linux.vnet.ibm.com> > >> > >> * sysdeps/ieee754/flt-32/s_cosf.c: Move reduced() and > >> constants to s_sincos.h file. > >> * sysdeps/ieee754/flt-32/s_sinf.c: Likewise. > >> * sysdeps/ieee754/flt-32/s_sincosf.c: New implementation. > >> * sysdeps/ieee754/flt-32/s_sincos.h: New file. > > > > OK. > > > > Please put those static constants into a different file and mark them > hidden so that they can be shared when they are compiled twice > for FMA. glibc is build with -fmerge-all-constants, so constants should already be shared at link time; making constants hidden would be relevant only for arrays, not for individual floating-point numbers (and if the individual constants aren't put in appropriate sections with link-time merging, it's probably because not doing so is more efficient on a particular architecture).
On Fri, Dec 15, 2017 at 6:17 AM, Rajalakshmi Srinivasaraghavan <raji@linux.vnet.ibm.com> wrote: > > > On 12/15/2017 07:49 AM, H.J. Lu wrote: >> >> On Thu, Dec 14, 2017 at 4:29 PM, Joseph Myers <joseph@codesourcery.com> >> wrote: >>> >>> On Wed, 13 Dec 2017, Rajalakshmi Srinivasaraghavan wrote: >>> >>>> This implementation is based on generic s_sinf.c and s_cosf.c. >>>> Tested on s390x, powerpc64le and powerpc32. >>>> >>>> 2017-12-12 Rajalakshmi Srinivasaraghavan <raji@linux.vnet.ibm.com> >>>> >>>> * sysdeps/ieee754/flt-32/s_cosf.c: Move reduced() and >>>> constants to s_sincos.h file. >>>> * sysdeps/ieee754/flt-32/s_sinf.c: Likewise. >>>> * sysdeps/ieee754/flt-32/s_sincosf.c: New implementation. >>>> * sysdeps/ieee754/flt-32/s_sincos.h: New file. >>> >>> >>> OK. >>> >> >> Please put those static constants into a different file and mark them >> hidden so that they can be shared when they are compiled twice >> for FMA. >> > > I assume you are asking changes similar to > sysdeps/ieee754/flt-32/e_exp2f_data.c. Could you please explain the Yes. On x86-64, I will compile sinf/cosf/sincosf twice, with and without FMA. Both versions of sinf/cosf/sincosf should share the same set of tables to avoid data duplication. > advantage of doing this or suggest any previous mail thread which I can > refer to? Thanks.
On Fri, Dec 15, 2017 at 6:18 AM, Joseph Myers <joseph@codesourcery.com> wrote: > On Fri, 15 Dec 2017, H.J. Lu wrote: > >> On Thu, Dec 14, 2017 at 4:29 PM, Joseph Myers <joseph@codesourcery.com> wrote: >> > On Wed, 13 Dec 2017, Rajalakshmi Srinivasaraghavan wrote: >> > >> >> This implementation is based on generic s_sinf.c and s_cosf.c. >> >> Tested on s390x, powerpc64le and powerpc32. >> >> >> >> 2017-12-12 Rajalakshmi Srinivasaraghavan <raji@linux.vnet.ibm.com> >> >> >> >> * sysdeps/ieee754/flt-32/s_cosf.c: Move reduced() and >> >> constants to s_sincos.h file. >> >> * sysdeps/ieee754/flt-32/s_sinf.c: Likewise. >> >> * sysdeps/ieee754/flt-32/s_sincosf.c: New implementation. >> >> * sysdeps/ieee754/flt-32/s_sincos.h: New file. >> > >> > OK. >> > >> >> Please put those static constants into a different file and mark them >> hidden so that they can be shared when they are compiled twice >> for FMA. > > glibc is build with -fmerge-all-constants, so constants should already be > shared at link time; making constants hidden would be relevant only for > arrays, not for individual floating-point numbers (and if the individual > constants aren't put in appropriate sections with link-time merging, it's > probably because not doing so is more efficient on a particular > architecture). > What do you suggest for x86-64 to avoid array duplication?
On Fri, 15 Dec 2017, H.J. Lu wrote: > > glibc is build with -fmerge-all-constants, so constants should already be > > shared at link time; making constants hidden would be relevant only for > > arrays, not for individual floating-point numbers (and if the individual > > constants aren't put in appropriate sections with link-time merging, it's > > probably because not doing so is more efficient on a particular > > architecture). > > What do you suggest for x86-64 to avoid array duplication? Well, a followup patch to refactor the arrays into hidden variables would be reasonable once the sincosf patch is in - but it would need benchmarking to make sure it doesn't adversely affect performance.
On Fri, Dec 15, 2017 at 7:48 AM, Joseph Myers <joseph@codesourcery.com> wrote: > On Fri, 15 Dec 2017, H.J. Lu wrote: > >> > glibc is build with -fmerge-all-constants, so constants should already be >> > shared at link time; making constants hidden would be relevant only for >> > arrays, not for individual floating-point numbers (and if the individual >> > constants aren't put in appropriate sections with link-time merging, it's >> > probably because not doing so is more efficient on a particular >> > architecture). >> >> What do you suggest for x86-64 to avoid array duplication? > > Well, a followup patch to refactor the arrays into hidden variables would > be reasonable once the sincosf patch is in - but it would need > benchmarking to make sure it doesn't adversely affect performance. Do you have any suspicions to indicate that hidden array may be slower than local one to access?
On Fri, 15 Dec 2017, H.J. Lu wrote: > On Fri, Dec 15, 2017 at 7:48 AM, Joseph Myers <joseph@codesourcery.com> wrote: > > On Fri, 15 Dec 2017, H.J. Lu wrote: > > > >> > glibc is build with -fmerge-all-constants, so constants should already be > >> > shared at link time; making constants hidden would be relevant only for > >> > arrays, not for individual floating-point numbers (and if the individual > >> > constants aren't put in appropriate sections with link-time merging, it's > >> > probably because not doing so is more efficient on a particular > >> > architecture). > >> > >> What do you suggest for x86-64 to avoid array duplication? > > > > Well, a followup patch to refactor the arrays into hidden variables would > > be reasonable once the sincosf patch is in - but it would need > > benchmarking to make sure it doesn't adversely affect performance. > > Do you have any suspicions to indicate that hidden array may be > slower than local one to access? On some architectures, GCC uses section anchors for access to variables defined in the same source file. In such cases, moving the arrays to being hidden and defined separately means that the address of each array needs to be loaded separately because the offsets between different variables are no longer known when the source files using them are being compiled, whereas if defined in each source file, a function accessing multiple arrays could load the anchor address and them compute all the other addresses with small offsets from it. You should be able to avoid that problem by putting all the data in a single structure, as done in e_exp2f_data.c, rather than separate hidden variables.
On Fri, Dec 15, 2017 at 8:24 AM, Joseph Myers <joseph@codesourcery.com> wrote: > On Fri, 15 Dec 2017, H.J. Lu wrote: > >> On Fri, Dec 15, 2017 at 7:48 AM, Joseph Myers <joseph@codesourcery.com> wrote: >> > On Fri, 15 Dec 2017, H.J. Lu wrote: >> > >> >> > glibc is build with -fmerge-all-constants, so constants should already be >> >> > shared at link time; making constants hidden would be relevant only for >> >> > arrays, not for individual floating-point numbers (and if the individual >> >> > constants aren't put in appropriate sections with link-time merging, it's >> >> > probably because not doing so is more efficient on a particular >> >> > architecture). >> >> >> >> What do you suggest for x86-64 to avoid array duplication? >> > >> > Well, a followup patch to refactor the arrays into hidden variables would >> > be reasonable once the sincosf patch is in - but it would need >> > benchmarking to make sure it doesn't adversely affect performance. >> >> Do you have any suspicions to indicate that hidden array may be >> slower than local one to access? > > On some architectures, GCC uses section anchors for access to variables > defined in the same source file. In such cases, moving the arrays to > being hidden and defined separately means that the address of each array > needs to be loaded separately because the offsets between different > variables are no longer known when the source files using them are being > compiled, whereas if defined in each source file, a function accessing > multiple arrays could load the anchor address and them compute all the > other addresses with small offsets from it. This won't show up on x86. > You should be able to avoid that problem by putting all the data in a > single structure, as done in e_exp2f_data.c, rather than separate hidden > variables. > Good to know. I have no objections now. I will prepare a patch after s_sincosf.c is checked in. Thanks.
On Fri, Dec 15, 2017 at 8:34 AM, H.J. Lu <hjl.tools@gmail.com> wrote: > On Fri, Dec 15, 2017 at 8:24 AM, Joseph Myers <joseph@codesourcery.com> wrote: >> On Fri, 15 Dec 2017, H.J. Lu wrote: >> >>> On Fri, Dec 15, 2017 at 7:48 AM, Joseph Myers <joseph@codesourcery.com> wrote: >>> > On Fri, 15 Dec 2017, H.J. Lu wrote: >>> > >>> >> > glibc is build with -fmerge-all-constants, so constants should already be >>> >> > shared at link time; making constants hidden would be relevant only for >>> >> > arrays, not for individual floating-point numbers (and if the individual >>> >> > constants aren't put in appropriate sections with link-time merging, it's >>> >> > probably because not doing so is more efficient on a particular >>> >> > architecture). >>> >> >>> >> What do you suggest for x86-64 to avoid array duplication? >>> > >>> > Well, a followup patch to refactor the arrays into hidden variables would >>> > be reasonable once the sincosf patch is in - but it would need >>> > benchmarking to make sure it doesn't adversely affect performance. >>> >>> Do you have any suspicions to indicate that hidden array may be >>> slower than local one to access? >> >> On some architectures, GCC uses section anchors for access to variables >> defined in the same source file. In such cases, moving the arrays to >> being hidden and defined separately means that the address of each array >> needs to be loaded separately because the offsets between different >> variables are no longer known when the source files using them are being >> compiled, whereas if defined in each source file, a function accessing >> multiple arrays could load the anchor address and them compute all the >> other addresses with small offsets from it. > > This won't show up on x86. > >> You should be able to avoid that problem by putting all the data in a >> single structure, as done in e_exp2f_data.c, rather than separate hidden >> variables. >> > > Good to know. > > I have no objections now. I will prepare a patch after s_sincosf.c is checked > in. > One comment, should sysdeps/ieee754/flt-32/s_sincos.h be renamed to sysdeps/ieee754/flt-32/s_sincosf.h? I am enclosing a patch to add sysdeps/ieee754/flt-32/s_sincosf_data.c. It saves 352 bytes on x86-64. H.J. From 9ba68929320c3bdb50b3aadc5d024318362fbfba Mon Sep 17 00:00:00 2001 From: "H.J. Lu" <hjl.tools@gmail.com> Date: Fri, 15 Dec 2017 12:08:22 -0800 Subject: [PATCH] Add s_sincosf_data.c Add s_sincosf_data.c to data between sinf, cosf and sincosf. Tested on x86-64. It reduced the read-only data size by 352 bytes. * math/Makefile (type-float-routines): Add s_sincosf_data. * sysdeps/ieee754/flt-32/math_config.h (SINCOSF_PIO2_TABLE): New. (SINCOSF_INVPIO4_TABLE): Likewise. (SINCOSF_DOUBLE_ONES): Likewise. (__sincosf_data): Likewise. * sysdeps/ieee754/flt-32/s_sincos.h: Include "math_config.h". (pio2_table): Replaced with __sincosf_data. (invpio4_table): Likewise. (ones): Likewise. * sysdeps/ieee754/flt-32/s_sincosf_data.c: New file. --- math/Makefile | 2 +- sysdeps/ieee754/flt-32/math_config.h | 11 +++++++++ sysdeps/ieee754/flt-32/s_sincos.h | 26 ++++----------------- sysdeps/ieee754/flt-32/s_sincosf_data.c | 41 +++++++++++++++++++++++++++++++++ 4 files changed, 58 insertions(+), 22 deletions(-) create mode 100644 sysdeps/ieee754/flt-32/s_sincosf_data.c diff --git a/math/Makefile b/math/Makefile index cba9ce9405..8762f20059 100644 --- a/math/Makefile +++ b/math/Makefile @@ -120,7 +120,7 @@ type-double-routines := branred doasin dosincos halfulp mpa mpatan2 \ # float support type-float-suffix := f type-float-routines := k_rem_pio2f math_errf e_exp2f_data e_logf_data \ - e_log2f_data e_powf_log2_data + e_log2f_data e_powf_log2_data s_sincosf_data # _Float128 support type-float128-suffix := f128 diff --git a/sysdeps/ieee754/flt-32/math_config.h b/sysdeps/ieee754/flt-32/math_config.h index e5a830b442..c2bc72805c 100644 --- a/sysdeps/ieee754/flt-32/math_config.h +++ b/sysdeps/ieee754/flt-32/math_config.h @@ -161,4 +161,15 @@ extern const struct powf_log2_data double poly[POWF_LOG2_POLY_ORDER]; } __powf_log2_data attribute_hidden; +/* Shared between sinf, cosf and sincosf. */ +#define SINCOSF_PIO2_TABLE 6 +#define SINCOSF_INVPIO4_TABLE 8 +#define SINCOSF_DOUBLE_ONES 2 +extern const struct sincosf_data +{ + double pio2_table[SINCOSF_PIO2_TABLE]; + double invpio4_table[SINCOSF_INVPIO4_TABLE]; + double ones[SINCOSF_DOUBLE_ONES]; +} __sincosf_data attribute_hidden; + #endif diff --git a/sysdeps/ieee754/flt-32/s_sincos.h b/sysdeps/ieee754/flt-32/s_sincos.h index b0110fc2af..89855c757f 100644 --- a/sysdeps/ieee754/flt-32/s_sincos.h +++ b/sysdeps/ieee754/flt-32/s_sincos.h @@ -16,6 +16,8 @@ License along with the GNU C Library; if not, see <http://www.gnu.org/licenses/>. */ +#include "math_config.h" + /* Chebyshev constants for cos, range -PI/4 - PI/4. */ static const double C0 = -0x1.ffffffffe98aep-2; static const double C1 = 0x1.55555545c50c7p-5; @@ -48,27 +50,9 @@ static const double inv_PI_4 = 0x1.45f306dc9c883p+0; /* 4/PI. */ #define FLOAT_EXPONENT_SHIFT 23 #define FLOAT_EXPONENT_BIAS 127 -static const double pio2_table[] = { - 0 * M_PI_2, - 1 * M_PI_2, - 2 * M_PI_2, - 3 * M_PI_2, - 4 * M_PI_2, - 5 * M_PI_2 -}; - -static const double invpio4_table[] = { - 0x0p+0, - 0x1.45f306cp+0, - 0x1.c9c882ap-28, - 0x1.4fe13a8p-58, - 0x1.f47d4dp-85, - 0x1.bb81b6cp-112, - 0x1.4acc9ep-142, - 0x1.0e4107cp-169 -}; - -static const double ones[] = { 1.0, -1.0 }; +#define pio2_table __sincosf_data.pio2_table +#define invpio4_table __sincosf_data.invpio4_table +#define ones __sincosf_data.ones /* Compute the sine value using Chebyshev polynomials where THETA is the range reduced absolute value of the input diff --git a/sysdeps/ieee754/flt-32/s_sincosf_data.c b/sysdeps/ieee754/flt-32/s_sincosf_data.c new file mode 100644 index 0000000000..d999ac5dfc --- /dev/null +++ b/sysdeps/ieee754/flt-32/s_sincosf_data.c @@ -0,0 +1,41 @@ +/* Shared data between sinf, cosf and sincosf. + Copyright (C) 2017 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 + <http://www.gnu.org/licenses/>. */ + +#include "math_config.h" + +const struct sincosf_data __sincosf_data = { + .pio2_table = { + 0 * M_PI_2, + 1 * M_PI_2, + 2 * M_PI_2, + 3 * M_PI_2, + 4 * M_PI_2, + 5 * M_PI_2 + }, + .invpio4_table = { + 0x0p+0, + 0x1.45f306cp+0, + 0x1.c9c882ap-28, + 0x1.4fe13a8p-58, + 0x1.f47d4dp-85, + 0x1.bb81b6cp-112, + 0x1.4acc9ep-142, + 0x1.0e4107cp-169 + }, + .ones = { 1.0, -1.0 } +};
On Fri, 15 Dec 2017, H.J. Lu wrote: > One comment, should sysdeps/ieee754/flt-32/s_sincos.h be renamed > to sysdeps/ieee754/flt-32/s_sincosf.h? Yes, I think so.
On 12/16/2017 02:20 AM, Joseph Myers wrote: > On Fri, 15 Dec 2017, H.J. Lu wrote: > >> One comment, should sysdeps/ieee754/flt-32/s_sincos.h be renamed >> to sysdeps/ieee754/flt-32/s_sincosf.h? > > Yes, I think so. > Pushed as 984ae9967b49830173490a33ae6130880f3f70d9 with the file name changed as s_sincosf.h. Thanks for the review.
On Sat, Dec 16, 2017 at 12:52 AM, Rajalakshmi Srinivasaraghavan <raji@linux.vnet.ibm.com> wrote: > > > On 12/16/2017 02:20 AM, Joseph Myers wrote: >> >> On Fri, 15 Dec 2017, H.J. Lu wrote: >> >>> One comment, should sysdeps/ieee754/flt-32/s_sincos.h be renamed >>> to sysdeps/ieee754/flt-32/s_sincosf.h? >> >> >> Yes, I think so. >> > > Pushed as 984ae9967b49830173490a33ae6130880f3f70d9 with the file > name changed as s_sincosf.h. Thanks for the review. > I noticed that sysdeps/x86_64/fpu/s_sincosf.S is still faster: Generic: "max": 276.971, "min": 10.813, "mean": 29.3755 SSE2: "max": 138.795, "min": 11.686, "mean": 22.9463 The SSE2 is 28% faster. Do they use the same algorithm?
On 12/16/2017 08:13 PM, H.J. Lu wrote: > On Sat, Dec 16, 2017 at 12:52 AM, Rajalakshmi Srinivasaraghavan > <raji@linux.vnet.ibm.com> wrote: >> >> >> On 12/16/2017 02:20 AM, Joseph Myers wrote: >>> >>> On Fri, 15 Dec 2017, H.J. Lu wrote: >>> >>>> One comment, should sysdeps/ieee754/flt-32/s_sincos.h be renamed >>>> to sysdeps/ieee754/flt-32/s_sincosf.h? >>> >>> >>> Yes, I think so. >>> >> >> Pushed as 984ae9967b49830173490a33ae6130880f3f70d9 with the file >> name changed as s_sincosf.h. Thanks for the review. >> > > I noticed that sysdeps/x86_64/fpu/s_sincosf.S is still faster: > > Generic: > > "max": 276.971, > "min": 10.813, > "mean": 29.3755 > > SSE2: > > "max": 138.795, > "min": 11.686, > "mean": 22.9463 > > The SSE2 is 28% faster. Do they use the same algorithm? Yes, they are same. One small difference is generic version calls reduced_sin() and reduced_cos() whereas asm version handles both in the same branch for reconstruction. > >
On Sun, Dec 17, 2017 at 7:26 PM, Rajalakshmi Srinivasaraghavan <raji@linux.vnet.ibm.com> wrote: > > > On 12/16/2017 08:13 PM, H.J. Lu wrote: >> >> On Sat, Dec 16, 2017 at 12:52 AM, Rajalakshmi Srinivasaraghavan >> <raji@linux.vnet.ibm.com> wrote: >>> >>> >>> >>> On 12/16/2017 02:20 AM, Joseph Myers wrote: >>>> >>>> >>>> On Fri, 15 Dec 2017, H.J. Lu wrote: >>>> >>>>> One comment, should sysdeps/ieee754/flt-32/s_sincos.h be renamed >>>>> to sysdeps/ieee754/flt-32/s_sincosf.h? >>>> >>>> >>>> >>>> Yes, I think so. >>>> >>> >>> Pushed as 984ae9967b49830173490a33ae6130880f3f70d9 with the file >>> name changed as s_sincosf.h. Thanks for the review. >>> >> >> I noticed that sysdeps/x86_64/fpu/s_sincosf.S is still faster: >> >> Generic: >> >> "max": 276.971, >> "min": 10.813, >> "mean": 29.3755 >> >> SSE2: >> >> "max": 138.795, >> "min": 11.686, >> "mean": 22.9463 >> >> The SSE2 is 28% faster. Do they use the same algorithm? > > > Yes, they are same. One small difference is generic version calls > reduced_sin() and reduced_cos() whereas asm version handles both in the same > branch for reconstruction. > Can generic version do the same?
On 12/19/2017 10:05 PM, H.J. Lu wrote: > On Sun, Dec 17, 2017 at 7:26 PM, Rajalakshmi Srinivasaraghavan > <raji@linux.vnet.ibm.com> wrote: >> >> >> On 12/16/2017 08:13 PM, H.J. Lu wrote: >>> >>> On Sat, Dec 16, 2017 at 12:52 AM, Rajalakshmi Srinivasaraghavan >>> <raji@linux.vnet.ibm.com> wrote: >>>> >>>> >>>> >>>> On 12/16/2017 02:20 AM, Joseph Myers wrote: >>>>> >>>>> >>>>> On Fri, 15 Dec 2017, H.J. Lu wrote: >>>>> >>>>>> One comment, should sysdeps/ieee754/flt-32/s_sincos.h be renamed >>>>>> to sysdeps/ieee754/flt-32/s_sincosf.h? >>>>> >>>>> >>>>> >>>>> Yes, I think so. >>>>> >>>> >>>> Pushed as 984ae9967b49830173490a33ae6130880f3f70d9 with the file >>>> name changed as s_sincosf.h. Thanks for the review. >>>> >>> >>> I noticed that sysdeps/x86_64/fpu/s_sincosf.S is still faster: >>> >>> Generic: >>> >>> "max": 276.971, >>> "min": 10.813, >>> "mean": 29.3755 >>> >>> SSE2: >>> >>> "max": 138.795, >>> "min": 11.686, >>> "mean": 22.9463 >>> >>> The SSE2 is 28% faster. Do they use the same algorithm? >> >> >> Yes, they are same. One small difference is generic version calls >> reduced_sin() and reduced_cos() whereas asm version handles both in the same >> branch for reconstruction. >> > > Can generic version do the same? I don't see any difference in mean value with that change. > > >
diff --git a/sysdeps/ieee754/flt-32/s_cosf.c b/sysdeps/ieee754/flt-32/s_cosf.c index ac6d044449..d19397faca 100644 --- a/sysdeps/ieee754/flt-32/s_cosf.c +++ b/sysdeps/ieee754/flt-32/s_cosf.c @@ -20,6 +20,7 @@ #include <math.h> #include <math_private.h> #include <libm-alias-float.h> +#include "s_sincos.h" #ifndef COSF # define COSF_FUNC __cosf @@ -27,95 +28,6 @@ # define COSF_FUNC COSF #endif -/* Chebyshev constants for cos, range -PI/4 - PI/4. */ -static const double C0 = -0x1.ffffffffe98aep-2; -static const double C1 = 0x1.55555545c50c7p-5; -static const double C2 = -0x1.6c16b348b6874p-10; -static const double C3 = 0x1.a00eb9ac43ccp-16; -static const double C4 = -0x1.23c97dd8844d7p-22; - -/* Chebyshev constants for sin, range -PI/4 - PI/4. */ -static const double S0 = -0x1.5555555551cd9p-3; -static const double S1 = 0x1.1111110c2688bp-7; -static const double S2 = -0x1.a019f8b4bd1f9p-13; -static const double S3 = 0x1.71d7264e6b5b4p-19; -static const double S4 = -0x1.a947e1674b58ap-26; - -/* Chebyshev constants for cos, range 2^-27 - 2^-5. */ -static const double CC0 = -0x1.fffffff5cc6fdp-2; -static const double CC1 = 0x1.55514b178dac5p-5; - -/* PI/2 with 98 bits of accuracy. */ -static const double PI_2_hi = 0x1.921fb544p+0; -static const double PI_2_lo = 0x1.0b4611a626332p-34; - -static const double inv_PI_4 = 0x1.45f306dc9c883p+0; /* 4/PI. */ - -#define FLOAT_EXPONENT_SHIFT 23 -#define FLOAT_EXPONENT_BIAS 127 - -static const double pio2_table[] = { - 0 * M_PI_2, - 1 * M_PI_2, - 2 * M_PI_2, - 3 * M_PI_2, - 4 * M_PI_2, - 5 * M_PI_2 -}; - -static const double invpio4_table[] = { - 0x0p+0, - 0x1.45f306cp+0, - 0x1.c9c882ap-28, - 0x1.4fe13a8p-58, - 0x1.f47d4dp-85, - 0x1.bb81b6cp-112, - 0x1.4acc9ep-142, - 0x1.0e4107cp-169 -}; - -static const double ones[] = { 1.0, -1.0 }; - -/* Compute the cosine value using Chebyshev polynomials where - THETA is the range reduced absolute value of the input - and it is less than Pi/4, - N is calculated as trunc(|x|/(Pi/4)) + 1 and it is used to decide - whether a sine or cosine approximation is more accurate and - the sign of the result. */ -static inline float -reduced (double theta, unsigned int n) -{ - double sign, cx; - const double theta2 = theta * theta; - - /* Determine positive or negative primary interval. */ - n += 2; - sign = ones[(n >> 2) & 1]; - - /* Are we in the primary interval of sin or cos? */ - if ((n & 2) == 0) - { - /* Here cosf() is calculated using sin Chebyshev polynomial: - x+x^3*(S0+x^2*(S1+x^2*(S2+x^2*(S3+x^2*S4)))). */ - cx = S3 + theta2 * S4; - cx = S2 + theta2 * cx; - cx = S1 + theta2 * cx; - cx = S0 + theta2 * cx; - cx = theta + theta * theta2 * cx; - } - else - { - /* Here cosf() is calculated using cos Chebyshev polynomial: - 1.0+x^2*(C0+x^2*(C1+x^2*(C2+x^2*(C3+x^2*C4)))). */ - cx = C3 + theta2 * C4; - cx = C2 + theta2 * cx; - cx = C1 + theta2 * cx; - cx = C0 + theta2 * cx; - cx = 1. + theta2 * cx; - } - return sign * cx; -} - float COSF_FUNC (float x) { @@ -161,7 +73,7 @@ COSF_FUNC (float x) pio2_table must go to 5 (9 / 2 + 1). */ unsigned int n = (abstheta * inv_PI_4) + 1; theta = abstheta - pio2_table[n / 2]; - return reduced (theta, n); + return reduced_cos (theta, n); } else if (isless (abstheta, INFINITY)) { @@ -171,7 +83,7 @@ COSF_FUNC (float x) double x = n / 2; theta = (abstheta - x * PI_2_hi) - x * PI_2_lo; /* Argument reduction needed. */ - return reduced (theta, n); + return reduced_cos (theta, n); } else /* |theta| >= 2^23. */ { @@ -199,7 +111,7 @@ COSF_FUNC (float x) e += c; e += d; e *= M_PI_4; - return reduced (e, l + 1); + return reduced_cos (e, l + 1); } else { @@ -209,14 +121,14 @@ COSF_FUNC (float x) if (e <= 1.0) { e *= M_PI_4; - return reduced (e, l + 1); + return reduced_cos (e, l + 1); } else { l++; e -= 2.0; e *= M_PI_4; - return reduced (e, l + 1); + return reduced_cos (e, l + 1); } } } diff --git a/sysdeps/ieee754/flt-32/s_sincos.h b/sysdeps/ieee754/flt-32/s_sincos.h new file mode 100644 index 0000000000..b0110fc2af --- /dev/null +++ b/sysdeps/ieee754/flt-32/s_sincos.h @@ -0,0 +1,155 @@ +/* Used by sinf, cosf and sincosf functions. + Copyright (C) 2017 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 + <http://www.gnu.org/licenses/>. */ + +/* Chebyshev constants for cos, range -PI/4 - PI/4. */ +static const double C0 = -0x1.ffffffffe98aep-2; +static const double C1 = 0x1.55555545c50c7p-5; +static const double C2 = -0x1.6c16b348b6874p-10; +static const double C3 = 0x1.a00eb9ac43ccp-16; +static const double C4 = -0x1.23c97dd8844d7p-22; + +/* Chebyshev constants for sin, range -PI/4 - PI/4. */ +static const double S0 = -0x1.5555555551cd9p-3; +static const double S1 = 0x1.1111110c2688bp-7; +static const double S2 = -0x1.a019f8b4bd1f9p-13; +static const double S3 = 0x1.71d7264e6b5b4p-19; +static const double S4 = -0x1.a947e1674b58ap-26; + +/* Chebyshev constants for sin, range 2^-27 - 2^-5. */ +static const double SS0 = -0x1.555555543d49dp-3; +static const double SS1 = 0x1.110f475cec8c5p-7; + +/* Chebyshev constants for cos, range 2^-27 - 2^-5. */ +static const double CC0 = -0x1.fffffff5cc6fdp-2; +static const double CC1 = 0x1.55514b178dac5p-5; + +/* PI/2 with 98 bits of accuracy. */ +static const double PI_2_hi = 0x1.921fb544p+0; +static const double PI_2_lo = 0x1.0b4611a626332p-34; + +static const double SMALL = 0x1p-50; /* 2^-50. */ +static const double inv_PI_4 = 0x1.45f306dc9c883p+0; /* 4/PI. */ + +#define FLOAT_EXPONENT_SHIFT 23 +#define FLOAT_EXPONENT_BIAS 127 + +static const double pio2_table[] = { + 0 * M_PI_2, + 1 * M_PI_2, + 2 * M_PI_2, + 3 * M_PI_2, + 4 * M_PI_2, + 5 * M_PI_2 +}; + +static const double invpio4_table[] = { + 0x0p+0, + 0x1.45f306cp+0, + 0x1.c9c882ap-28, + 0x1.4fe13a8p-58, + 0x1.f47d4dp-85, + 0x1.bb81b6cp-112, + 0x1.4acc9ep-142, + 0x1.0e4107cp-169 +}; + +static const double ones[] = { 1.0, -1.0 }; + +/* Compute the sine value using Chebyshev polynomials where + THETA is the range reduced absolute value of the input + and it is less than Pi/4, + N is calculated as trunc(|x|/(Pi/4)) + 1 and it is used to decide + whether a sine or cosine approximation is more accurate and + SIGNBIT is used to add the correct sign after the Chebyshev + polynomial is computed. */ +static inline float +reduced_sin (const double theta, const unsigned int n, + const unsigned int signbit) +{ + double sx; + const double theta2 = theta * theta; + /* We are operating on |x|, so we need to add back the original + signbit for sinf. */ + double sign; + /* Determine positive or negative primary interval. */ + sign = ones[((n >> 2) & 1) ^ signbit]; + /* Are we in the primary interval of sin or cos? */ + if ((n & 2) == 0) + { + /* Here sinf() is calculated using sin Chebyshev polynomial: + x+x^3*(S0+x^2*(S1+x^2*(S2+x^2*(S3+x^2*S4)))). */ + sx = S3 + theta2 * S4; /* S3+x^2*S4. */ + sx = S2 + theta2 * sx; /* S2+x^2*(S3+x^2*S4). */ + sx = S1 + theta2 * sx; /* S1+x^2*(S2+x^2*(S3+x^2*S4)). */ + sx = S0 + theta2 * sx; /* S0+x^2*(S1+x^2*(S2+x^2*(S3+x^2*S4))). */ + sx = theta + theta * theta2 * sx; + } + else + { + /* Here sinf() is calculated using cos Chebyshev polynomial: + 1.0+x^2*(C0+x^2*(C1+x^2*(C2+x^2*(C3+x^2*C4)))). */ + sx = C3 + theta2 * C4; /* C3+x^2*C4. */ + sx = C2 + theta2 * sx; /* C2+x^2*(C3+x^2*C4). */ + sx = C1 + theta2 * sx; /* C1+x^2*(C2+x^2*(C3+x^2*C4)). */ + sx = C0 + theta2 * sx; /* C0+x^2*(C1+x^2*(C2+x^2*(C3+x^2*C4))). */ + sx = 1.0 + theta2 * sx; + } + + /* Add in the signbit and assign the result. */ + return sign * sx; +} + +/* Compute the cosine value using Chebyshev polynomials where + THETA is the range reduced absolute value of the input + and it is less than Pi/4, + N is calculated as trunc(|x|/(Pi/4)) + 1 and it is used to decide + whether a sine or cosine approximation is more accurate and + the sign of the result. */ +static inline float +reduced_cos (double theta, unsigned int n) +{ + double sign, cx; + const double theta2 = theta * theta; + + /* Determine positive or negative primary interval. */ + n += 2; + sign = ones[(n >> 2) & 1]; + + /* Are we in the primary interval of sin or cos? */ + if ((n & 2) == 0) + { + /* Here cosf() is calculated using sin Chebyshev polynomial: + x+x^3*(S0+x^2*(S1+x^2*(S2+x^2*(S3+x^2*S4)))). */ + cx = S3 + theta2 * S4; + cx = S2 + theta2 * cx; + cx = S1 + theta2 * cx; + cx = S0 + theta2 * cx; + cx = theta + theta * theta2 * cx; + } + else + { + /* Here cosf() is calculated using cos Chebyshev polynomial: + 1.0+x^2*(C0+x^2*(C1+x^2*(C2+x^2*(C3+x^2*C4)))). */ + cx = C3 + theta2 * C4; + cx = C2 + theta2 * cx; + cx = C1 + theta2 * cx; + cx = C0 + theta2 * cx; + cx = 1. + theta2 * cx; + } + return sign * cx; +} diff --git a/sysdeps/ieee754/flt-32/s_sincosf.c b/sysdeps/ieee754/flt-32/s_sincosf.c index 4946a6eb54..b9f775d45b 100644 --- a/sysdeps/ieee754/flt-32/s_sincosf.c +++ b/sysdeps/ieee754/flt-32/s_sincosf.c @@ -1,7 +1,6 @@ /* Compute sine and cosine of argument. - Copyright (C) 1997-2017 Free Software Foundation, Inc. + Copyright (C) 2017 Free Software Foundation, Inc. This file is part of the GNU C Library. - Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. The GNU C Library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public @@ -19,9 +18,9 @@ #include <errno.h> #include <math.h> - #include <math_private.h> #include <libm-alias-float.h> +#include "s_sincos.h" #ifndef SINCOSF # define SINCOSF_FUNC __sincosf @@ -32,50 +31,137 @@ void SINCOSF_FUNC (float x, float *sinx, float *cosx) { - int32_t ix; - - /* High word of x. */ - GET_FLOAT_WORD (ix, x); - - /* |x| ~< pi/4 */ - ix &= 0x7fffffff; - if (ix <= 0x3f490fd8) - { - *sinx = __kernel_sinf (x, 0.0, 0); - *cosx = __kernel_cosf (x, 0.0); - } - else if (ix>=0x7f800000) + double cx; + double theta = x; + double abstheta = fabs (theta); + /* If |x|< Pi/4. */ + if (isless (abstheta, M_PI_4)) { - /* sin(Inf or NaN) is NaN */ - *sinx = *cosx = x - x; - if (ix == 0x7f800000) - __set_errno (EDOM); + if (abstheta >= 0x1p-5) /* |x| >= 2^-5. */ + { + const double theta2 = theta * theta; + /* Chebyshev polynomial of the form for sin and cos. */ + cx = C3 + theta2 * C4; + cx = C2 + theta2 * cx; + cx = C1 + theta2 * cx; + cx = C0 + theta2 * cx; + cx = 1.0 + theta2 * cx; + *cosx = cx; + cx = S3 + theta2 * S4; + cx = S2 + theta2 * cx; + cx = S1 + theta2 * cx; + cx = S0 + theta2 * cx; + cx = theta + theta * theta2 * cx; + *sinx = cx; + } + else if (abstheta >= 0x1p-27) /* |x| >= 2^-27. */ + { + /* A simpler Chebyshev approximation is close enough for this range: + for sin: x+x^3*(SS0+x^2*SS1) + for cos: 1.0+x^2*(CC0+x^3*CC1). */ + const double theta2 = theta * theta; + cx = CC0 + theta * theta2 * CC1; + cx = 1.0 + theta2 * cx; + *cosx = cx; + cx = SS0 + theta2 * SS1; + cx = theta + theta * theta2 * cx; + *sinx = cx; + } + else + { + /* Handle some special cases. */ + if (theta) + *sinx = theta - (theta * SMALL); + else + *sinx = theta; + *cosx = 1.0 - abstheta; + } } - else + else /* |x| >= Pi/4. */ { - /* Argument reduction needed. */ - float y[2]; - int n; - - n = __ieee754_rem_pio2f (x, y); - switch (n & 3) + unsigned int signbit = isless (x, 0); + if (isless (abstheta, 9 * M_PI_4)) /* |x| < 9*Pi/4. */ + { + /* There are cases where FE_UPWARD rounding mode can + produce a result of abstheta * inv_PI_4 == 9, + where abstheta < 9pi/4, so the domain for + pio2_table must go to 5 (9 / 2 + 1). */ + unsigned int n = (abstheta * inv_PI_4) + 1; + theta = abstheta - pio2_table[n / 2]; + *sinx = reduced_sin (theta, n, signbit); + *cosx = reduced_cos (theta, n); + } + else if (isless (abstheta, INFINITY)) + { + if (abstheta < 0x1p+23) /* |x| < 2^23. */ + { + unsigned int n = ((unsigned int) (abstheta * inv_PI_4)) + 1; + double x = n / 2; + theta = (abstheta - x * PI_2_hi) - x * PI_2_lo; + /* Argument reduction needed. */ + *sinx = reduced_sin (theta, n, signbit); + *cosx = reduced_cos (theta, n); + } + else /* |x| >= 2^23. */ + { + x = fabsf (x); + int exponent; + GET_FLOAT_WORD (exponent, x); + exponent + = (exponent >> FLOAT_EXPONENT_SHIFT) - FLOAT_EXPONENT_BIAS; + exponent += 3; + exponent /= 28; + double a = invpio4_table[exponent] * x; + double b = invpio4_table[exponent + 1] * x; + double c = invpio4_table[exponent + 2] * x; + double d = invpio4_table[exponent + 3] * x; + uint64_t l = a; + l &= ~0x7; + a -= l; + double e = a + b; + l = e; + e = a - l; + if (l & 1) + { + e -= 1.0; + e += b; + e += c; + e += d; + e *= M_PI_4; + *sinx = reduced_sin (e, l + 1, signbit); + *cosx = reduced_cos (e, l + 1); + } + else + { + e += b; + e += c; + e += d; + if (e <= 1.0) + { + e *= M_PI_4; + *sinx = reduced_sin (e, l + 1, signbit); + *cosx = reduced_cos (e, l + 1); + } + else + { + l++; + e -= 2.0; + e *= M_PI_4; + *sinx = reduced_sin (e, l + 1, signbit); + *cosx = reduced_cos (e, l + 1); + } + } + } + } + else { - case 0: - *sinx = __kernel_sinf (y[0], y[1], 1); - *cosx = __kernel_cosf (y[0], y[1]); - break; - case 1: - *sinx = __kernel_cosf (y[0], y[1]); - *cosx = -__kernel_sinf (y[0], y[1], 1); - break; - case 2: - *sinx = -__kernel_sinf (y[0], y[1], 1); - *cosx = -__kernel_cosf (y[0], y[1]); - break; - default: - *sinx = -__kernel_cosf (y[0], y[1]); - *cosx = __kernel_sinf (y[0], y[1], 1); - break; + int32_t ix; + /* High word of x. */ + GET_FLOAT_WORD (ix, abstheta); + /* sin/cos(Inf or NaN) is NaN. */ + *sinx = *cosx = x - x; + if (ix == 0x7f800000) + __set_errno (EDOM); } } } diff --git a/sysdeps/ieee754/flt-32/s_sinf.c b/sysdeps/ieee754/flt-32/s_sinf.c index 418d4487c5..c682bc64fc 100644 --- a/sysdeps/ieee754/flt-32/s_sinf.c +++ b/sysdeps/ieee754/flt-32/s_sinf.c @@ -20,6 +20,7 @@ #include <math.h> #include <math_private.h> #include <libm-alias-float.h> +#include "s_sincos.h" #ifndef SINF # define SINF_FUNC __sinf @@ -27,100 +28,6 @@ # define SINF_FUNC SINF #endif -/* Chebyshev constants for cos, range -PI/4 - PI/4. */ -static const double C0 = -0x1.ffffffffe98aep-2; -static const double C1 = 0x1.55555545c50c7p-5; -static const double C2 = -0x1.6c16b348b6874p-10; -static const double C3 = 0x1.a00eb9ac43ccp-16; -static const double C4 = -0x1.23c97dd8844d7p-22; - -/* Chebyshev constants for sin, range -PI/4 - PI/4. */ -static const double S0 = -0x1.5555555551cd9p-3; -static const double S1 = 0x1.1111110c2688bp-7; -static const double S2 = -0x1.a019f8b4bd1f9p-13; -static const double S3 = 0x1.71d7264e6b5b4p-19; -static const double S4 = -0x1.a947e1674b58ap-26; - -/* Chebyshev constants for sin, range 2^-27 - 2^-5. */ -static const double SS0 = -0x1.555555543d49dp-3; -static const double SS1 = 0x1.110f475cec8c5p-7; - -/* PI/2 with 98 bits of accuracy. */ -static const double PI_2_hi = -0x1.921fb544p+0; -static const double PI_2_lo = -0x1.0b4611a626332p-34; - -static const double SMALL = 0x1p-50; /* 2^-50. */ -static const double inv_PI_4 = 0x1.45f306dc9c883p+0; /* 4/PI. */ - -#define FLOAT_EXPONENT_SHIFT 23 -#define FLOAT_EXPONENT_BIAS 127 - -static const double pio2_table[] = { - 0 * M_PI_2, - 1 * M_PI_2, - 2 * M_PI_2, - 3 * M_PI_2, - 4 * M_PI_2, - 5 * M_PI_2 -}; - -static const double invpio4_table[] = { - 0x0p+0, - 0x1.45f306cp+0, - 0x1.c9c882ap-28, - 0x1.4fe13a8p-58, - 0x1.f47d4dp-85, - 0x1.bb81b6cp-112, - 0x1.4acc9ep-142, - 0x1.0e4107cp-169 -}; - -static const double ones[] = { 1.0, -1.0 }; - -/* Compute the sine value using Chebyshev polynomials where - THETA is the range reduced absolute value of the input - and it is less than Pi/4, - N is calculated as trunc(|x|/(Pi/4)) + 1 and it is used to decide - whether a sine or cosine approximation is more accurate and - SIGNBIT is used to add the correct sign after the Chebyshev - polynomial is computed. */ -static inline float -reduced (const double theta, const unsigned int n, - const unsigned int signbit) -{ - double sx; - const double theta2 = theta * theta; - /* We are operating on |x|, so we need to add back the original - signbit for sinf. */ - double sign; - /* Determine positive or negative primary interval. */ - sign = ones[((n >> 2) & 1) ^ signbit]; - /* Are we in the primary interval of sin or cos? */ - if ((n & 2) == 0) - { - /* Here sinf() is calculated using sin Chebyshev polynomial: - x+x^3*(S0+x^2*(S1+x^2*(S2+x^2*(S3+x^2*S4)))). */ - sx = S3 + theta2 * S4; /* S3+x^2*S4. */ - sx = S2 + theta2 * sx; /* S2+x^2*(S3+x^2*S4). */ - sx = S1 + theta2 * sx; /* S1+x^2*(S2+x^2*(S3+x^2*S4)). */ - sx = S0 + theta2 * sx; /* S0+x^2*(S1+x^2*(S2+x^2*(S3+x^2*S4))). */ - sx = theta + theta * theta2 * sx; - } - else - { - /* Here sinf() is calculated using cos Chebyshev polynomial: - 1.0+x^2*(C0+x^2*(C1+x^2*(C2+x^2*(C3+x^2*C4)))). */ - sx = C3 + theta2 * C4; /* C3+x^2*C4. */ - sx = C2 + theta2 * sx; /* C2+x^2*(C3+x^2*C4). */ - sx = C1 + theta2 * sx; /* C1+x^2*(C2+x^2*(C3+x^2*C4)). */ - sx = C0 + theta2 * sx; /* C0+x^2*(C1+x^2*(C2+x^2*(C3+x^2*C4))). */ - sx = 1.0 + theta2 * sx; - } - - /* Add in the signbit and assign the result. */ - return sign * sx; -} - float SINF_FUNC (float x) { @@ -171,7 +78,7 @@ SINF_FUNC (float x) pio2_table must go to 5 (9 / 2 + 1). */ unsigned int n = (abstheta * inv_PI_4) + 1; theta = abstheta - pio2_table[n / 2]; - return reduced (theta, n, signbit); + return reduced_sin (theta, n, signbit); } else if (isless (abstheta, INFINITY)) { @@ -179,9 +86,9 @@ SINF_FUNC (float x) { unsigned int n = ((unsigned int) (abstheta * inv_PI_4)) + 1; double x = n / 2; - theta = x * PI_2_lo + (x * PI_2_hi + abstheta); + theta = (abstheta - x * PI_2_hi) - x * PI_2_lo; /* Argument reduction needed. */ - return reduced (theta, n, signbit); + return reduced_sin (theta, n, signbit); } else /* |x| >= 2^23. */ { @@ -209,7 +116,7 @@ SINF_FUNC (float x) e += c; e += d; e *= M_PI_4; - return reduced (e, l + 1, signbit); + return reduced_sin (e, l + 1, signbit); } else { @@ -219,14 +126,14 @@ SINF_FUNC (float x) if (e <= 1.0) { e *= M_PI_4; - return reduced (e, l + 1, signbit); + return reduced_sin (e, l + 1, signbit); } else { l++; e -= 2.0; e *= M_PI_4; - return reduced (e, l + 1, signbit); + return reduced_sin (e, l + 1, signbit); } } }