Update.
authorUlrich Drepper <drepper@redhat.com>
Thu, 15 Jul 1999 18:26:25 +0000 (18:26 +0000)
committerUlrich Drepper <drepper@redhat.com>
Thu, 15 Jul 1999 18:26:25 +0000 (18:26 +0000)
1999-07-15  Jakub Jelinek  <jj@ultra.linux.cz>

* math/Makefile: Add t_sincosl and k_sincosl support routines.
* math/math_private.h (__kernel_sincosl): New declaration.
* sysdeps/generic/t_sincosl.c: New file.
* sysdeps/generic/k_sincosl.c: New file.
* sysdeps/ieee754/ldbl-128/k_cosl.c: New file.
* sysdeps/ieee754/ldbl-128/k_sinl.c: New file.
* sysdeps/ieee754/ldbl-128/k_sincosl.c: New file.
* sysdeps/ieee754/ldbl-128/t_sincosl.c: New file.
* sysdeps/ieee754/ldbl-128/e_rem_pio2l.c: New file.
* sysdeps/ieee754/ldbl-128/s_sincosl.c (__sincosl): Use
__kernel_sincosl.
* sysdeps/ieee754/ldbl-128/math_ldbl.h (GET_LDOUBLE_LSW64): New
definition.

ChangeLog
math/Makefile
math/math_private.h
sysdeps/generic/k_sincosl.c [new file with mode: 0644]
sysdeps/generic/t_sincosl.c [new file with mode: 0644]
sysdeps/ieee754/ldbl-128/e_rem_pio2l.c [new file with mode: 0644]
sysdeps/ieee754/ldbl-128/k_cosl.c [new file with mode: 0644]
sysdeps/ieee754/ldbl-128/k_sincosl.c [new file with mode: 0644]
sysdeps/ieee754/ldbl-128/k_sinl.c [new file with mode: 0644]
sysdeps/ieee754/ldbl-128/math_ldbl.h
sysdeps/ieee754/ldbl-128/s_sincosl.c

index 8c764df..5dd0ab4 100644 (file)
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,19 @@
+1999-07-15  Jakub Jelinek  <jj@ultra.linux.cz>
+
+       * math/Makefile: Add t_sincosl and k_sincosl support routines.
+       * math/math_private.h (__kernel_sincosl): New declaration.
+       * sysdeps/generic/t_sincosl.c: New file.
+       * sysdeps/generic/k_sincosl.c: New file.
+       * sysdeps/ieee754/ldbl-128/k_cosl.c: New file.
+       * sysdeps/ieee754/ldbl-128/k_sinl.c: New file.
+       * sysdeps/ieee754/ldbl-128/k_sincosl.c: New file.
+       * sysdeps/ieee754/ldbl-128/t_sincosl.c: New file.
+       * sysdeps/ieee754/ldbl-128/e_rem_pio2l.c: New file.
+       * sysdeps/ieee754/ldbl-128/s_sincosl.c (__sincosl): Use
+       __kernel_sincosl.
+       * sysdeps/ieee754/ldbl-128/math_ldbl.h (GET_LDOUBLE_LSW64): New
+       definition.
+
 1999-07-15  Ulrich Drepper  <drepper@cygnus.com>
 
        * posix/unistd.h: Use __PMT for exit.
index e1b3431..d5434f7 100644 (file)
@@ -61,7 +61,8 @@ libm-routines = $(strip $(libm-support) $(libm-calls) \
                        $(patsubst %_rf,%f_r,$(libm-calls:=f))  \
                        $(long-m-$(long-double-fcts)))
 long-m-routines = $(patsubst %_rl,%l_r,$(libm-calls:=l))
-long-m-yes = $(long-m-routines)
+long-m-support = t_sincosl k_sincosl
+long-m-yes = $(long-m-routines) $(long-m-support)
 distribute += $(long-m-yes:=.c)
 
 # These functions are in libc instead of libm because __printf_fp
index 35e3e4e..29dbb46 100644 (file)
@@ -264,6 +264,8 @@ extern long double __ieee754_scalbl __P((long double,long double));
 extern long double __kernel_sinl __P((long double,long double,int));
 extern long double __kernel_cosl __P((long double,long double));
 extern long double __kernel_tanl __P((long double,long double,int));
+extern void __kernel_sincosl __P((long double,long double,
+                                 long double *,long double *, int));
 extern int   __kernel_rem_pio2l __P((long double*,long double*,int,int,
                                     int,const int*));
 
diff --git a/sysdeps/generic/k_sincosl.c b/sysdeps/generic/k_sincosl.c
new file mode 100644 (file)
index 0000000..aa038c2
--- /dev/null
@@ -0,0 +1 @@
+/* Empty.  Not needed.  */
diff --git a/sysdeps/generic/t_sincosl.c b/sysdeps/generic/t_sincosl.c
new file mode 100644 (file)
index 0000000..6b271e6
--- /dev/null
@@ -0,0 +1 @@
+/* Empty.  Not needed unless ldbl __kernel_* functions use it. */
diff --git a/sysdeps/ieee754/ldbl-128/e_rem_pio2l.c b/sysdeps/ieee754/ldbl-128/e_rem_pio2l.c
new file mode 100644 (file)
index 0000000..0da0560
--- /dev/null
@@ -0,0 +1,274 @@
+/* Quad-precision floating point argument reduction.
+   Copyright (C) 1999 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+   Contributed by Jakub Jelinek <jj@ultra.linux.cz>
+
+   The GNU C Library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Library General Public License as
+   published by the Free Software Foundation; either version 2 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
+   Library General Public License for more details.
+
+   You should have received a copy of the GNU Library General Public
+   License along with the GNU C Library; see the file COPYING.LIB.  If not,
+   write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
+   Boston, MA 02111-1307, USA.  */
+
+#include "math.h"
+#include "math_private.h"
+
+/*
+ * Table of constants for 2/pi, 5628 hexadecimal digits of 2/pi 
+ */
+static const int32_t two_over_pi[] = {
+0xa2f983, 0x6e4e44, 0x1529fc, 0x2757d1, 0xf534dd, 0xc0db62, 
+0x95993c, 0x439041, 0xfe5163, 0xabdebb, 0xc561b7, 0x246e3a, 
+0x424dd2, 0xe00649, 0x2eea09, 0xd1921c, 0xfe1deb, 0x1cb129, 
+0xa73ee8, 0x8235f5, 0x2ebb44, 0x84e99c, 0x7026b4, 0x5f7e41, 
+0x3991d6, 0x398353, 0x39f49c, 0x845f8b, 0xbdf928, 0x3b1ff8, 
+0x97ffde, 0x05980f, 0xef2f11, 0x8b5a0a, 0x6d1f6d, 0x367ecf, 
+0x27cb09, 0xb74f46, 0x3f669e, 0x5fea2d, 0x7527ba, 0xc7ebe5, 
+0xf17b3d, 0x0739f7, 0x8a5292, 0xea6bfb, 0x5fb11f, 0x8d5d08, 
+0x560330, 0x46fc7b, 0x6babf0, 0xcfbc20, 0x9af436, 0x1da9e3, 
+0x91615e, 0xe61b08, 0x659985, 0x5f14a0, 0x68408d, 0xffd880, 
+0x4d7327, 0x310606, 0x1556ca, 0x73a8c9, 0x60e27b, 0xc08c6b, 
+0x47c419, 0xc367cd, 0xdce809, 0x2a8359, 0xc4768b, 0x961ca6, 
+0xddaf44, 0xd15719, 0x053ea5, 0xff0705, 0x3f7e33, 0xe832c2, 
+0xde4f98, 0x327dbb, 0xc33d26, 0xef6b1e, 0x5ef89f, 0x3a1f35, 
+0xcaf27f, 0x1d87f1, 0x21907c, 0x7c246a, 0xfa6ed5, 0x772d30, 
+0x433b15, 0xc614b5, 0x9d19c3, 0xc2c4ad, 0x414d2c, 0x5d000c, 
+0x467d86, 0x2d71e3, 0x9ac69b, 0x006233, 0x7cd2b4, 0x97a7b4, 
+0xd55537, 0xf63ed7, 0x1810a3, 0xfc764d, 0x2a9d64, 0xabd770, 
+0xf87c63, 0x57b07a, 0xe71517, 0x5649c0, 0xd9d63b, 0x3884a7, 
+0xcb2324, 0x778ad6, 0x23545a, 0xb91f00, 0x1b0af1, 0xdfce19, 
+0xff319f, 0x6a1e66, 0x615799, 0x47fbac, 0xd87f7e, 0xb76522, 
+0x89e832, 0x60bfe6, 0xcdc4ef, 0x09366c, 0xd43f5d, 0xd7de16, 
+0xde3b58, 0x929bde, 0x2822d2, 0xe88628, 0x4d58e2, 0x32cac6, 
+0x16e308, 0xcb7de0, 0x50c017, 0xa71df3, 0x5be018, 0x34132e, 
+0x621283, 0x014883, 0x5b8ef5, 0x7fb0ad, 0xf2e91e, 0x434a48, 
+0xd36710, 0xd8ddaa, 0x425fae, 0xce616a, 0xa4280a, 0xb499d3, 
+0xf2a606, 0x7f775c, 0x83c2a3, 0x883c61, 0x78738a, 0x5a8caf, 
+0xbdd76f, 0x63a62d, 0xcbbff4, 0xef818d, 0x67c126, 0x45ca55, 
+0x36d9ca, 0xd2a828, 0x8d61c2, 0x77c912, 0x142604, 0x9b4612, 
+0xc459c4, 0x44c5c8, 0x91b24d, 0xf31700, 0xad43d4, 0xe54929, 
+0x10d5fd, 0xfcbe00, 0xcc941e, 0xeece70, 0xf53e13, 0x80f1ec, 
+0xc3e7b3, 0x28f8c7, 0x940593, 0x3e71c1, 0xb3092e, 0xf3450b, 
+0x9c1288, 0x7b20ab, 0x9fb52e, 0xc29247, 0x2f327b, 0x6d550c, 
+0x90a772, 0x1fe76b, 0x96cb31, 0x4a1679, 0xe27941, 0x89dff4, 
+0x9794e8, 0x84e6e2, 0x973199, 0x6bed88, 0x365f5f, 0x0efdbb, 
+0xb49a48, 0x6ca467, 0x427271, 0x325d8d, 0xb8159f, 0x09e5bc, 
+0x25318d, 0x3974f7, 0x1c0530, 0x010c0d, 0x68084b, 0x58ee2c, 
+0x90aa47, 0x02e774, 0x24d6bd, 0xa67df7, 0x72486e, 0xef169f, 
+0xa6948e, 0xf691b4, 0x5153d1, 0xf20acf, 0x339820, 0x7e4bf5, 
+0x6863b2, 0x5f3edd, 0x035d40, 0x7f8985, 0x295255, 0xc06437, 
+0x10d86d, 0x324832, 0x754c5b, 0xd4714e, 0x6e5445, 0xc1090b, 
+0x69f52a, 0xd56614, 0x9d0727, 0x50045d, 0xdb3bb4, 0xc576ea, 
+0x17f987, 0x7d6b49, 0xba271d, 0x296996, 0xacccc6, 0x5414ad, 
+0x6ae290, 0x89d988, 0x50722c, 0xbea404, 0x940777, 0x7030f3, 
+0x27fc00, 0xa871ea, 0x49c266, 0x3de064, 0x83dd97, 0x973fa3, 
+0xfd9443, 0x8c860d, 0xde4131, 0x9d3992, 0x8c70dd, 0xe7b717, 
+0x3bdf08, 0x2b3715, 0xa0805c, 0x93805a, 0x921110, 0xd8e80f, 
+0xaf806c, 0x4bffdb, 0x0f9038, 0x761859, 0x15a562, 0xbbcb61, 
+0xb989c7, 0xbd4010, 0x04f2d2, 0x277549, 0xf6b6eb, 0xbb22db, 
+0xaa140a, 0x2f2689, 0x768364, 0x333b09, 0x1a940e, 0xaa3a51, 
+0xc2a31d, 0xaeedaf, 0x12265c, 0x4dc26d, 0x9c7a2d, 0x9756c0, 
+0x833f03, 0xf6f009, 0x8c402b, 0x99316d, 0x07b439, 0x15200c, 
+0x5bc3d8, 0xc492f5, 0x4badc6, 0xa5ca4e, 0xcd37a7, 0x36a9e6, 
+0x9492ab, 0x6842dd, 0xde6319, 0xef8c76, 0x528b68, 0x37dbfc, 
+0xaba1ae, 0x3115df, 0xa1ae00, 0xdafb0c, 0x664d64, 0xb705ed, 
+0x306529, 0xbf5657, 0x3aff47, 0xb9f96a, 0xf3be75, 0xdf9328, 
+0x3080ab, 0xf68c66, 0x15cb04, 0x0622fa, 0x1de4d9, 0xa4b33d, 
+0x8f1b57, 0x09cd36, 0xe9424e, 0xa4be13, 0xb52333, 0x1aaaf0, 
+0xa8654f, 0xa5c1d2, 0x0f3f0b, 0xcd785b, 0x76f923, 0x048b7b, 
+0x721789, 0x53a6c6, 0xe26e6f, 0x00ebef, 0x584a9b, 0xb7dac4, 
+0xba66aa, 0xcfcf76, 0x1d02d1, 0x2df1b1, 0xc1998c, 0x77adc3, 
+0xda4886, 0xa05df7, 0xf480c6, 0x2ff0ac, 0x9aecdd, 0xbc5c3f, 
+0x6dded0, 0x1fc790, 0xb6db2a, 0x3a25a3, 0x9aaf00, 0x9353ad, 
+0x0457b6, 0xb42d29, 0x7e804b, 0xa707da, 0x0eaa76, 0xa1597b, 
+0x2a1216, 0x2db7dc, 0xfde5fa, 0xfedb89, 0xfdbe89, 0x6c76e4, 
+0xfca906, 0x70803e, 0x156e85, 0xff87fd, 0x073e28, 0x336761, 
+0x86182a, 0xeabd4d, 0xafe7b3, 0x6e6d8f, 0x396795, 0x5bbf31, 
+0x48d784, 0x16df30, 0x432dc7, 0x356125, 0xce70c9, 0xb8cb30, 
+0xfd6cbf, 0xa200a4, 0xe46c05, 0xa0dd5a, 0x476f21, 0xd21262, 
+0x845cb9, 0x496170, 0xe0566b, 0x015299, 0x375550, 0xb7d51e, 
+0xc4f133, 0x5f6e13, 0xe4305d, 0xa92e85, 0xc3b21d, 0x3632a1, 
+0xa4b708, 0xd4b1ea, 0x21f716, 0xe4698f, 0x77ff27, 0x80030c, 
+0x2d408d, 0xa0cd4f, 0x99a520, 0xd3a2b3, 0x0a5d2f, 0x42f9b4, 
+0xcbda11, 0xd0be7d, 0xc1db9b, 0xbd17ab, 0x81a2ca, 0x5c6a08, 
+0x17552e, 0x550027, 0xf0147f, 0x8607e1, 0x640b14, 0x8d4196, 
+0xdebe87, 0x2afdda, 0xb6256b, 0x34897b, 0xfef305, 0x9ebfb9, 
+0x4f6a68, 0xa82a4a, 0x5ac44f, 0xbcf82d, 0x985ad7, 0x95c7f4, 
+0x8d4d0d, 0xa63a20, 0x5f57a4, 0xb13f14, 0x953880, 0x0120cc, 
+0x86dd71, 0xb6dec9, 0xf560bf, 0x11654d, 0x6b0701, 0xacb08c, 
+0xd0c0b2, 0x485551, 0x0efb1e, 0xc37295, 0x3b06a3, 0x3540c0, 
+0x7bdc06, 0xcc45e0, 0xfa294e, 0xc8cad6, 0x41f3e8, 0xde647c, 
+0xd8649b, 0x31bed9, 0xc397a4, 0xd45877, 0xc5e369, 0x13daf0, 
+0x3c3aba, 0x461846, 0x5f7555, 0xf5bdd2, 0xc6926e, 0x5d2eac, 
+0xed440e, 0x423e1c, 0x87c461, 0xe9fd29, 0xf3d6e7, 0xca7c22, 
+0x35916f, 0xc5e008, 0x8dd7ff, 0xe26a6e, 0xc6fdb0, 0xc10893, 
+0x745d7c, 0xb2ad6b, 0x9d6ecd, 0x7b723e, 0x6a11c6, 0xa9cff7, 
+0xdf7329, 0xbac9b5, 0x5100b7, 0x0db2e2, 0x24ba74, 0x607de5, 
+0x8ad874, 0x2c150d, 0x0c1881, 0x94667e, 0x162901, 0x767a9f, 
+0xbefdfd, 0xef4556, 0x367ed9, 0x13d9ec, 0xb9ba8b, 0xfc97c4, 
+0x27a831, 0xc36ef1, 0x36c594, 0x56a8d8, 0xb5a8b4, 0x0ecccf, 
+0x2d8912, 0x34576f, 0x89562c, 0xe3ce99, 0xb920d6, 0xaa5e6b, 
+0x9c2a3e, 0xcc5f11, 0x4a0bfd, 0xfbf4e1, 0x6d3b8e, 0x2c86e2, 
+0x84d4e9, 0xa9b4fc, 0xd1eeef, 0xc9352e, 0x61392f, 0x442138, 
+0xc8d91b, 0x0afc81, 0x6a4afb, 0xd81c2f, 0x84b453, 0x8c994e, 
+0xcc2254, 0xdc552a, 0xd6c6c0, 0x96190b, 0xb8701a, 0x649569, 
+0x605a26, 0xee523f, 0x0f117f, 0x11b5f4, 0xf5cbfc, 0x2dbc34, 
+0xeebc34, 0xcc5de8, 0x605edd, 0x9b8e67, 0xef3392, 0xb817c9, 
+0x9b5861, 0xbc57e1, 0xc68351, 0x103ed8, 0x4871dd, 0xdd1c2d, 
+0xa118af, 0x462c21, 0xd7f359, 0x987ad9, 0xc0549e, 0xfa864f, 
+0xfc0656, 0xae79e5, 0x362289, 0x22ad38, 0xdc9367, 0xaae855, 
+0x382682, 0x9be7ca, 0xa40d51, 0xb13399, 0x0ed7a9, 0x480569, 
+0xf0b265, 0xa7887f, 0x974c88, 0x36d1f9, 0xb39221, 0x4a827b, 
+0x21cf98, 0xdc9f40, 0x5547dc, 0x3a74e1, 0x42eb67, 0xdf9dfe, 
+0x5fd45e, 0xa4677b, 0x7aacba, 0xa2f655, 0x23882b, 0x55ba41, 
+0x086e59, 0x862a21, 0x834739, 0xe6e389, 0xd49ee5, 0x40fb49, 
+0xe956ff, 0xca0f1c, 0x8a59c5, 0x2bfa94, 0xc5c1d3, 0xcfc50f, 
+0xae5adb, 0x86c547, 0x624385, 0x3b8621, 0x94792c, 0x876110, 
+0x7b4c2a, 0x1a2c80, 0x12bf43, 0x902688, 0x893c78, 0xe4c4a8, 
+0x7bdbe5, 0xc23ac4, 0xeaf426, 0x8a67f7, 0xbf920d, 0x2ba365, 
+0xb1933d, 0x0b7cbd, 0xdc51a4, 0x63dd27, 0xdde169, 0x19949a, 
+0x9529a8, 0x28ce68, 0xb4ed09, 0x209f44, 0xca984e, 0x638270, 
+0x237c7e, 0x32b90f, 0x8ef5a7, 0xe75614, 0x08f121, 0x2a9db5, 
+0x4d7e6f, 0x5119a5, 0xabf9b5, 0xd6df82, 0x61dd96, 0x023616, 
+0x9f3ac4, 0xa1a283, 0x6ded72, 0x7a8d39, 0xa9b882, 0x5c326b, 
+0x5b2746, 0xed3400, 0x7700d2, 0x55f4fc, 0x4d5901, 0x8071e0, 
+0xe13f89, 0xb295f3, 0x64a8f1, 0xaea74b, 0x38fc4c, 0xeab2bb, 
+0x47270b, 0xabc3a7, 0x34ba60, 0x52dd34, 0xf8563a, 0xeb7e8a, 
+0x31bb36, 0x5895b7, 0x47f7a9, 0x94c3aa, 0xd39225, 0x1e7f3e, 
+0xd8974e, 0xbba94f, 0xd8ae01, 0xe661b4, 0x393d8e, 0xa523aa, 
+0x33068e, 0x1633b5, 0x3bb188, 0x1d3a9d, 0x4013d0, 0xcc1be5, 
+0xf862e7, 0x3bf28f, 0x39b5bf, 0x0bc235, 0x22747e, 0xa247c0, 
+0xd52d1f, 0x19add3, 0x9094df, 0x9311d0, 0xb42b25, 0x496db2, 
+0xe264b2, 0x5ef135, 0x3bc6a4, 0x1a4ad0, 0xaac92e, 0x64e886, 
+0x573091, 0x982cfb, 0x311b1a, 0x08728b, 0xbdcee1, 0x60e142, 
+0xeb641d, 0xd0bba3, 0xe559d4, 0x597b8c, 0x2a4483, 0xf332ba, 
+0xf84867, 0x2c8d1b, 0x2fa9b0, 0x50f3dd, 0xf9f573, 0xdb61b4, 
+0xfe233e, 0x6c41a6, 0xeea318, 0x775a26, 0xbc5e5c, 0xcea708, 
+0x94dc57, 0xe20196, 0xf1e839, 0xbe4851, 0x5d2d2f, 0x4e9555, 
+0xd96ec2, 0xe7d755, 0x6304e0, 0xc02e0e, 0xfc40a0, 0xbbf9b3, 
+0x7125a7, 0x222dfb, 0xf619d8, 0x838c1c, 0x6619e6, 0xb20d55, 
+0xbb5137, 0x79e809, 0xaf9149, 0x0d73de, 0x0b0da5, 0xce7f58, 
+0xac1934, 0x724667, 0x7a1a13, 0x9e26bc, 0x4555e7, 0x585cb5, 
+0x711d14, 0x486991, 0x480d60, 0x56adab, 0xd62f64, 0x96ee0c, 
+0x212ff3, 0x5d6d88, 0xa67684, 0x95651e, 0xab9e0a, 0x4ddefe, 
+0x571010, 0x836a39, 0xf8ea31, 0x9e381d, 0xeac8b1, 0xcac96b, 
+0x37f21e, 0xd505e9, 0x984743, 0x9fc56c, 0x0331b7, 0x3b8bf8, 
+0x86e56a, 0x8dc343, 0x6230e7, 0x93cfd5, 0x6a8f2d, 0x733005, 
+0x1af021, 0xa09fcb, 0x7415a1, 0xd56b23, 0x6ff725, 0x2f4bc7, 
+0xb8a591, 0x7fac59, 0x5c55de, 0x212c38, 0xb13296, 0x5cff50, 
+0x366262, 0xfa7b16, 0xf4d9a6, 0x2acfe7, 0xf07403, 0xd4d604, 
+0x6fd916, 0x31b1bf, 0xcbb450, 0x5bd7c8, 0x0ce194, 0x6bd643, 
+0x4fd91c, 0xdf4543, 0x5f3453, 0xe2b5aa, 0xc9aec8, 0x131485, 
+0xf9d2bf, 0xbadb9e, 0x76f5b9, 0xaf15cf, 0xca3182, 0x14b56d, 
+0xe9fe4d, 0x50fc35, 0xf5aed5, 0xa2d0c1, 0xc96057, 0x192eb6, 
+0xe91d92, 0x07d144, 0xaea3c6, 0x343566, 0x26d5b4, 0x3161e2, 
+0x37f1a2, 0x209eff, 0x958e23, 0x493798, 0x35f4a6, 0x4bdc02, 
+0xc2be13, 0xbe80a0, 0x0b72a3, 0x115c5f, 0x1e1bd1, 0x0db4d3, 
+0x869e85, 0x96976b, 0x2ac91f, 0x8a26c2, 0x3070f0, 0x041412, 
+0xfc9fa5, 0xf72a38, 0x9c6878, 0xe2aa76, 0x50cfe1, 0x559274, 
+0x934e38, 0x0a92f7, 0x5533f0, 0xa63db4, 0x399971, 0xe2b755, 
+0xa98a7c, 0x008f19, 0xac54d2, 0x2ea0b4, 0xf5f3e0, 0x60c849, 
+0xffd269, 0xae52ce, 0x7a5fdd, 0xe9ce06, 0xfb0ae8, 0xa50cce, 
+0xea9d3e, 0x3766dd, 0xb834f5, 0x0da090, 0x846f88, 0x4ae3d5, 
+0x099a03, 0x2eae2d, 0xfcb40a, 0xfb9b33, 0xe281dd, 0x1b16ba, 
+0xd8c0af, 0xd96b97, 0xb52dc9, 0x9c277f, 0x5951d5, 0x21ccd6, 
+0xb6496b, 0x584562, 0xb3baf2, 0xa1a5c4, 0x7ca2cf, 0xa9b93d, 
+0x7b7b89, 0x483d38, 
+};
+
+static const long double c[] = {
+/* 93 bits of pi/2 */
+#define PI_2_1 c[0]
+ 1.57079632679489661923132169155131424e+00L, /* 3fff921fb54442d18469898cc5100000 */
+
+/* pi/2 - PI_2_1 */
+#define PI_2_1t c[1]
+ 8.84372056613570112025531863263659260e-29L, /* 3fa1c06e0e68948127044533e63a0106 */
+};
+
+int32_t __ieee754_rem_pio2l(long double x, long double *y)
+{
+  long double z, w, t;
+  double tx[8];
+  int64_t exp, n, ix, hx;
+  u_int64_t lx;
+
+  GET_LDOUBLE_WORDS64 (hx, lx, x);
+  ix = hx & 0x7fffffffffffffffLL;
+  if (ix <= 0x3ffe921fb54442d1LL)      /* x in <-pi/4, pi/4> */
+    {
+      y[0] = x;
+      y[1] = 0;
+      return 0;
+    }
+
+  if (ix < 0x40002d97c7f3321dLL)       /* |x| in <pi/4, 3pi/4) */
+    {
+      if (hx > 0)
+       { 
+         /* 113 + 93 bit PI is ok */
+         z = x - PI_2_1;
+         y[0] = z - PI_2_1t;
+         y[1] = (z - y[0]) - PI_2_1t;
+         return 1;
+       }
+      else
+        {
+         /* 113 + 93 bit PI is ok */
+         z = x + PI_2_1;
+         y[0] = z + PI_2_1t;
+         y[1] = (z - y[0]) + PI_2_1t;
+         return -1;
+       }
+    }
+
+  if (ix >= 0x7fff000000000000LL)      /* x is +=oo or NaN */
+    {
+      y[0] = x - x;
+      y[1] = y[0];
+      return 0;
+    }
+
+  /* Handle large arguments.
+     We split the 113 bits of the mantissa into 5 24bit integers
+     stored in a double array.  */
+  exp = (ix >> 48) - 16383 - 23;
+
+  /* This is faster than doing this in floating point, because we
+     have to convert it to integers anyway and like this we can keep
+     both integer and floating point units busy.  */
+  tx [0] = (double)(((ix >> 25) & 0x7fffff) | 0x800000);
+  tx [1] = (double)((ix >> 1) & 0xffffff);
+  tx [2] = (double)(((ix << 23) | (lx >> 41)) & 0xffffff);
+  tx [3] = (double)((lx >> 17) & 0xffffff);
+  tx [4] = (double)((lx << 7) & 0xffffff);
+
+  n = __kernel_rem_pio2 (tx, tx + 5, exp, ((lx << 7) & 0xffffff) ? 5 : 4,
+                        3, two_over_pi);
+
+  /* The result is now stored in 3 double values, we need to convert it into
+     two long double values.  */
+  t = (long double) tx [6] + (long double) tx [7];
+  w = (long double) tx [5];
+
+  if (hx >= 0)
+    {
+      y[0] = w + t;
+      y[1] = t - (y[0] - w);
+      return n;
+    }
+  else
+    {
+      y[0] = -(w + t);
+      y[1] = -t - (y[0] + w);
+      return -n;
+    }
+}
diff --git a/sysdeps/ieee754/ldbl-128/k_cosl.c b/sysdeps/ieee754/ldbl-128/k_cosl.c
new file mode 100644 (file)
index 0000000..9a74710
--- /dev/null
@@ -0,0 +1,128 @@
+/* Quad-precision floating point cosine on <-pi/4,pi/4>.
+   Copyright (C) 1999 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+   Contributed by Jakub Jelinek <jj@ultra.linux.cz>
+
+   The GNU C Library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Library General Public License as
+   published by the Free Software Foundation; either version 2 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
+   Library General Public License for more details.
+
+   You should have received a copy of the GNU Library General Public
+   License along with the GNU C Library; see the file COPYING.LIB.  If not,
+   write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
+   Boston, MA 02111-1307, USA.  */
+
+#include "math.h"
+#include "math_private.h"
+
+static const long double c[] = {
+#define ONE c[0]
+ 1.00000000000000000000000000000000000E+00L, /* 3fff0000000000000000000000000000 */
+
+/* cos x ~ ONE + x^2 ( SCOS1 + SCOS2 * x^2 + ... + SCOS4 * x^6 + SCOS5 * x^8 )
+   x in <0,1/256>  */
+#define SCOS1 c[1]
+#define SCOS2 c[2]
+#define SCOS3 c[3]
+#define SCOS4 c[4]
+#define SCOS5 c[5]
+-5.00000000000000000000000000000000000E-01L, /* bffe0000000000000000000000000000 */
+ 4.16666666666666666666666666556146073E-02L, /* 3ffa5555555555555555555555395023 */
+-1.38888888888888888888309442601939728E-03L, /* bff56c16c16c16c16c16a566e42c0375 */
+ 2.48015873015862382987049502531095061E-05L, /* 3fefa01a01a019ee02dcf7da2d6d5444 */
+-2.75573112601362126593516899592158083E-07L, /* bfe927e4f5dce637cb0b54908754bde0 */
+
+/* cos x ~ ONE + x^2 ( COS1 + COS2 * x^2 + ... + COS7 * x^12 + COS8 * x^14 )
+   x in <0,0.1484375>  */
+#define COS1 c[6]
+#define COS2 c[7]
+#define COS3 c[8]
+#define COS4 c[9]
+#define COS5 c[10]
+#define COS6 c[11]
+#define COS7 c[12]
+#define COS8 c[13]
+-4.99999999999999999999999999999999759E-01L, /* bffdfffffffffffffffffffffffffffb */
+ 4.16666666666666666666666666651287795E-02L, /* 3ffa5555555555555555555555516f30 */
+-1.38888888888888888888888742314300284E-03L, /* bff56c16c16c16c16c16c16a463dfd0d */
+ 2.48015873015873015867694002851118210E-05L, /* 3fefa01a01a01a01a0195cebe6f3d3a5 */
+-2.75573192239858811636614709689300351E-07L, /* bfe927e4fb7789f5aa8142a22044b51f */
+ 2.08767569877762248667431926878073669E-09L, /* 3fe21eed8eff881d1e9262d7adff4373 */
+-1.14707451049343817400420280514614892E-11L, /* bfda9397496922a9601ed3d4ca48944b */
+ 4.77810092804389587579843296923533297E-14L, /* 3fd2ae5f8197cbcdcaf7c3fb4523414c */
+
+/* sin x ~ ONE * x + x^3 ( SSIN1 + SSIN2 * x^2 + ... + SSIN4 * x^6 + SSIN5 * x^8 )
+   x in <0,1/256>  */
+#define SSIN1 c[14]
+#define SSIN2 c[15]
+#define SSIN3 c[16]
+#define SSIN4 c[17]
+#define SSIN5 c[18]
+-1.66666666666666666666666666666666659E-01L, /* bffc5555555555555555555555555555 */
+ 8.33333333333333333333333333146298442E-03L, /* 3ff81111111111111111111110fe195d */
+-1.98412698412698412697726277416810661E-04L, /* bff2a01a01a01a01a019e7121e080d88 */
+ 2.75573192239848624174178393552189149E-06L, /* 3fec71de3a556c640c6aaa51aa02ab41 */
+-2.50521016467996193495359189395805639E-08L, /* bfe5ae644ee90c47dc71839de75b2787 */
+};
+
+#define SINCOSL_COS_HI 0
+#define SINCOSL_COS_LO 1
+#define SINCOSL_SIN_HI 2
+#define SINCOSL_SIN_LO 3
+extern const long double __sincosl_table[];
+
+long double
+__kernel_cosl(long double x, long double y)
+{
+  long double h, l, z, sin_l, cos_l_m1;
+  int64_t ix;
+  u_int32_t tix, hix, index;
+  GET_LDOUBLE_MSW64 (ix, x);
+  tix = ((u_int64_t)ix) >> 32;
+  tix &= ~0x80000000;                  /* tix = |x|'s high 32 bits */
+  if (tix < 0x3ffc3000)                        /* |x| < 0.1484375 */
+    {
+      /* Argument is small enough to approximate it by a Chebyshev
+        polynomial of degree 16.  */
+      if (tix < 0x3fc60000)            /* |x| < 2^-57 */
+       if (!((int)x)) return ONE;      /* generate inexact */
+      z = x * x;
+      return ONE + (z*(COS1+z*(COS2+z*(COS3+z*(COS4+
+                   z*(COS5+z*(COS6+z*(COS7+z*COS8))))))));
+    }
+  else
+    {
+      /* So that we don't have to use too large polynomial,  we find
+        l and h such that x = l + h,  where fabsl(l) <= 1.0/256 with 83
+        possible values for h.  We look up cosl(h) and sinl(h) in
+        pre-computed tables,  compute cosl(l) and sinl(l) using a
+        Chebyshev polynomial of degree 10(11) and compute
+        cosl(h+l) = cosl(h)cosl(l) - sinl(h)sinl(l).  */
+      index = 0x3ffe - (tix >> 16);
+      hix = (tix + (0x200 << index)) & (0xfffffc00 << index);
+      x = fabsl (x);
+      switch (index)
+       {
+       case 0: index = ((45 << 10) + hix - 0x3ffe0000) >> 8; break;
+       case 1: index = ((13 << 11) + hix - 0x3ffd0000) >> 9; break;
+       default:
+       case 2: index = (hix - 0x3ffc3000) >> 10; break;
+       }
+
+      SET_LDOUBLE_WORDS64(h, ((u_int64_t)hix) << 32, 0);
+      l = y - (h - x);
+      z = l * l;
+      sin_l = l*(ONE+z*(SSIN1+z*(SSIN2+z*(SSIN3+z*(SSIN4+z*SSIN5)))));
+      cos_l_m1 = z*(SCOS1+z*(SCOS2+z*(SCOS3+z*(SCOS4+z*SCOS5))));
+      return __sincosl_table [index + SINCOSL_COS_HI]
+            + (__sincosl_table [index + SINCOSL_COS_LO]
+               - (__sincosl_table [index + SINCOSL_SIN_HI] * sin_l
+                  - __sincosl_table [index + SINCOSL_COS_HI] * cos_l_m1));
+    }
+}
diff --git a/sysdeps/ieee754/ldbl-128/k_sincosl.c b/sysdeps/ieee754/ldbl-128/k_sincosl.c
new file mode 100644 (file)
index 0000000..2c8dcb8
--- /dev/null
@@ -0,0 +1,163 @@
+/* Quad-precision floating point sine and cosine on <-pi/4,pi/4>.
+   Copyright (C) 1999 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+   Contributed by Jakub Jelinek <jj@ultra.linux.cz>
+
+   The GNU C Library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Library General Public License as
+   published by the Free Software Foundation; either version 2 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
+   Library General Public License for more details.
+
+   You should have received a copy of the GNU Library General Public
+   License along with the GNU C Library; see the file COPYING.LIB.  If not,
+   write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
+   Boston, MA 02111-1307, USA.  */
+
+#include "math.h"
+#include "math_private.h"
+
+static const long double c[] = {
+#define ONE c[0]
+ 1.00000000000000000000000000000000000E+00L, /* 3fff0000000000000000000000000000 */
+
+/* cos x ~ ONE + x^2 ( SCOS1 + SCOS2 * x^2 + ... + SCOS4 * x^6 + SCOS5 * x^8 )
+   x in <0,1/256>  */
+#define SCOS1 c[1]
+#define SCOS2 c[2]
+#define SCOS3 c[3]
+#define SCOS4 c[4]
+#define SCOS5 c[5]
+-5.00000000000000000000000000000000000E-01L, /* bffe0000000000000000000000000000 */
+ 4.16666666666666666666666666556146073E-02L, /* 3ffa5555555555555555555555395023 */
+-1.38888888888888888888309442601939728E-03L, /* bff56c16c16c16c16c16a566e42c0375 */
+ 2.48015873015862382987049502531095061E-05L, /* 3fefa01a01a019ee02dcf7da2d6d5444 */
+-2.75573112601362126593516899592158083E-07L, /* bfe927e4f5dce637cb0b54908754bde0 */
+
+/* cos x ~ ONE + x^2 ( COS1 + COS2 * x^2 + ... + COS7 * x^12 + COS8 * x^14 )
+   x in <0,0.1484375>  */
+#define COS1 c[6]
+#define COS2 c[7]
+#define COS3 c[8]
+#define COS4 c[9]
+#define COS5 c[10]
+#define COS6 c[11]
+#define COS7 c[12]
+#define COS8 c[13]
+-4.99999999999999999999999999999999759E-01L, /* bffdfffffffffffffffffffffffffffb */
+ 4.16666666666666666666666666651287795E-02L, /* 3ffa5555555555555555555555516f30 */
+-1.38888888888888888888888742314300284E-03L, /* bff56c16c16c16c16c16c16a463dfd0d */
+ 2.48015873015873015867694002851118210E-05L, /* 3fefa01a01a01a01a0195cebe6f3d3a5 */
+-2.75573192239858811636614709689300351E-07L, /* bfe927e4fb7789f5aa8142a22044b51f */
+ 2.08767569877762248667431926878073669E-09L, /* 3fe21eed8eff881d1e9262d7adff4373 */
+-1.14707451049343817400420280514614892E-11L, /* bfda9397496922a9601ed3d4ca48944b */
+ 4.77810092804389587579843296923533297E-14L, /* 3fd2ae5f8197cbcdcaf7c3fb4523414c */
+
+/* sin x ~ ONE * x + x^3 ( SSIN1 + SSIN2 * x^2 + ... + SSIN4 * x^6 + SSIN5 * x^8 )
+   x in <0,1/256>  */
+#define SSIN1 c[14]
+#define SSIN2 c[15]
+#define SSIN3 c[16]
+#define SSIN4 c[17]
+#define SSIN5 c[18]
+-1.66666666666666666666666666666666659E-01L, /* bffc5555555555555555555555555555 */
+ 8.33333333333333333333333333146298442E-03L, /* 3ff81111111111111111111110fe195d */
+-1.98412698412698412697726277416810661E-04L, /* bff2a01a01a01a01a019e7121e080d88 */
+ 2.75573192239848624174178393552189149E-06L, /* 3fec71de3a556c640c6aaa51aa02ab41 */
+-2.50521016467996193495359189395805639E-08L, /* bfe5ae644ee90c47dc71839de75b2787 */
+
+/* sin x ~ ONE * x + x^3 ( SIN1 + SIN2 * x^2 + ... + SIN7 * x^12 + SIN8 * x^14 )
+   x in <0,0.1484375>  */
+#define SIN1 c[19]
+#define SIN2 c[20]
+#define SIN3 c[21]
+#define SIN4 c[22]
+#define SIN5 c[23]
+#define SIN6 c[24]
+#define SIN7 c[25]
+#define SIN8 c[26]
+-1.66666666666666666666666666666666538e-01L, /* bffc5555555555555555555555555550 */
+ 8.33333333333333333333333333307532934e-03L, /* 3ff811111111111111111111110e7340 */
+-1.98412698412698412698412534478712057e-04L, /* bff2a01a01a01a01a01a019e7a626296 */
+ 2.75573192239858906520896496653095890e-06L, /* 3fec71de3a556c7338fa38527474b8f5 */
+-2.50521083854417116999224301266655662e-08L, /* bfe5ae64567f544e16c7de65c2ea551f */
+ 1.60590438367608957516841576404938118e-10L, /* 3fde6124613a811480538a9a41957115 */
+-7.64716343504264506714019494041582610e-13L, /* bfd6ae7f3d5aef30c7bc660b060ef365 */
+ 2.81068754939739570236322404393398135e-15L, /* 3fce9510115aabf87aceb2022a9a9180 */
+};
+
+#define SINCOSL_COS_HI 0
+#define SINCOSL_COS_LO 1
+#define SINCOSL_SIN_HI 2
+#define SINCOSL_SIN_LO 3
+extern const long double __sincosl_table[];
+
+void
+__kernel_sincosl(long double x, long double y, long double *sinx, long double *cosx, int iy)
+{
+  long double h, l, z, sin_l, cos_l_m1;
+  int64_t ix;
+  u_int32_t tix, hix, index;
+  GET_LDOUBLE_MSW64 (ix, x);
+  tix = ((u_int64_t)ix) >> 32;
+  tix &= ~0x80000000;                  /* tix = |x|'s high 32 bits */
+  if (tix < 0x3ffc3000)                        /* |x| < 0.1484375 */
+    {
+      /* Argument is small enough to approximate it by a Chebyshev
+        polynomial of degree 16(17).  */
+      if (tix < 0x3fc60000)            /* |x| < 2^-57 */
+       if (!((int)x))                  /* generate inexact */
+         {
+           *sinx = x;
+           *cosx = ONE;
+           return;
+         }
+      z = x * x;
+      *sinx = x + (x * (z*(SIN1+z*(SIN2+z*(SIN3+z*(SIN4+
+                       z*(SIN5+z*(SIN6+z*(SIN7+z*SIN8)))))))));
+      *cosx = ONE + (z*(COS1+z*(COS2+z*(COS3+z*(COS4+
+                    z*(COS5+z*(COS6+z*(COS7+z*COS8))))))));
+    }
+  else
+    {
+      /* So that we don't have to use too large polynomial,  we find
+        l and h such that x = l + h,  where fabsl(l) <= 1.0/256 with 83
+        possible values for h.  We look up cosl(h) and sinl(h) in
+        pre-computed tables,  compute cosl(l) and sinl(l) using a
+        Chebyshev polynomial of degree 10(11) and compute
+        sinl(h+l) = sinl(h)cosl(l) + cosl(h)sinl(l) and
+        cosl(h+l) = cosl(h)cosl(l) - sinl(h)sinl(l).  */
+      index = 0x3ffe - (tix >> 16);
+      hix = (tix + (0x200 << index)) & (0xfffffc00 << index);
+      x = fabsl (x);
+      switch (index)
+       {
+       case 0: index = ((45 << 10) + hix - 0x3ffe0000) >> 8; break;
+       case 1: index = ((13 << 11) + hix - 0x3ffd0000) >> 9; break;
+       default:
+       case 2: index = (hix - 0x3ffc3000) >> 10; break;
+       }
+
+      SET_LDOUBLE_WORDS64(h, ((u_int64_t)hix) << 32, 0);
+      if (iy)
+       l = y - (h - x);
+      else
+       l = x - h;
+      z = l * l;
+      sin_l = l*(ONE+z*(SSIN1+z*(SSIN2+z*(SSIN3+z*(SSIN4+z*SSIN5)))));
+      cos_l_m1 = z*(SCOS1+z*(SCOS2+z*(SCOS3+z*(SCOS4+z*SCOS5))));
+      z = __sincosl_table [index + SINCOSL_SIN_HI]
+         + (__sincosl_table [index + SINCOSL_SIN_LO]
+            + (__sincosl_table [index + SINCOSL_SIN_HI] * cos_l_m1)
+            + (__sincosl_table [index + SINCOSL_COS_HI] * sin_l));
+      *sinx = (ix < 0) ? -z : z;
+      *cosx = __sincosl_table [index + SINCOSL_COS_HI]
+             + (__sincosl_table [index + SINCOSL_COS_LO]
+                - (__sincosl_table [index + SINCOSL_SIN_HI] * sin_l
+                   - __sincosl_table [index + SINCOSL_COS_HI] * cos_l_m1));
+    }
+}
diff --git a/sysdeps/ieee754/ldbl-128/k_sinl.c b/sysdeps/ieee754/ldbl-128/k_sinl.c
new file mode 100644 (file)
index 0000000..2247e06
--- /dev/null
@@ -0,0 +1,132 @@
+/* Quad-precision floating point sine on <-pi/4,pi/4>.
+   Copyright (C) 1999 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+   Contributed by Jakub Jelinek <jj@ultra.linux.cz>
+
+   The GNU C Library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Library General Public License as
+   published by the Free Software Foundation; either version 2 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
+   Library General Public License for more details.
+
+   You should have received a copy of the GNU Library General Public
+   License along with the GNU C Library; see the file COPYING.LIB.  If not,
+   write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
+   Boston, MA 02111-1307, USA.  */
+
+#include "math.h"
+#include "math_private.h"
+
+static const long double c[] = {
+#define ONE c[0]
+ 1.00000000000000000000000000000000000E+00L, /* 3fff0000000000000000000000000000 */
+
+/* cos x ~ ONE + x^2 ( SCOS1 + SCOS2 * x^2 + ... + SCOS4 * x^6 + SCOS5 * x^8 )
+   x in <0,1/256>  */
+#define SCOS1 c[1]
+#define SCOS2 c[2]
+#define SCOS3 c[3]
+#define SCOS4 c[4]
+#define SCOS5 c[5]
+-5.00000000000000000000000000000000000E-01L, /* bffe0000000000000000000000000000 */
+ 4.16666666666666666666666666556146073E-02L, /* 3ffa5555555555555555555555395023 */
+-1.38888888888888888888309442601939728E-03L, /* bff56c16c16c16c16c16a566e42c0375 */
+ 2.48015873015862382987049502531095061E-05L, /* 3fefa01a01a019ee02dcf7da2d6d5444 */
+-2.75573112601362126593516899592158083E-07L, /* bfe927e4f5dce637cb0b54908754bde0 */
+
+/* sin x ~ ONE * x + x^3 ( SIN1 + SIN2 * x^2 + ... + SIN7 * x^12 + SIN8 * x^14 )
+   x in <0,0.1484375>  */
+#define SIN1 c[6]
+#define SIN2 c[7]
+#define SIN3 c[8]
+#define SIN4 c[9]
+#define SIN5 c[10]
+#define SIN6 c[11]
+#define SIN7 c[12]
+#define SIN8 c[13]
+-1.66666666666666666666666666666666538e-01L, /* bffc5555555555555555555555555550 */
+ 8.33333333333333333333333333307532934e-03L, /* 3ff811111111111111111111110e7340 */
+-1.98412698412698412698412534478712057e-04L, /* bff2a01a01a01a01a01a019e7a626296 */
+ 2.75573192239858906520896496653095890e-06L, /* 3fec71de3a556c7338fa38527474b8f5 */
+-2.50521083854417116999224301266655662e-08L, /* bfe5ae64567f544e16c7de65c2ea551f */
+ 1.60590438367608957516841576404938118e-10L, /* 3fde6124613a811480538a9a41957115 */
+-7.64716343504264506714019494041582610e-13L, /* bfd6ae7f3d5aef30c7bc660b060ef365 */
+ 2.81068754939739570236322404393398135e-15L, /* 3fce9510115aabf87aceb2022a9a9180 */
+
+/* sin x ~ ONE * x + x^3 ( SSIN1 + SSIN2 * x^2 + ... + SSIN4 * x^6 + SSIN5 * x^8 )
+   x in <0,1/256>  */
+#define SSIN1 c[14]
+#define SSIN2 c[15]
+#define SSIN3 c[16]
+#define SSIN4 c[17]
+#define SSIN5 c[18]
+-1.66666666666666666666666666666666659E-01L, /* bffc5555555555555555555555555555 */
+ 8.33333333333333333333333333146298442E-03L, /* 3ff81111111111111111111110fe195d */
+-1.98412698412698412697726277416810661E-04L, /* bff2a01a01a01a01a019e7121e080d88 */
+ 2.75573192239848624174178393552189149E-06L, /* 3fec71de3a556c640c6aaa51aa02ab41 */
+-2.50521016467996193495359189395805639E-08L, /* bfe5ae644ee90c47dc71839de75b2787 */
+};
+
+#define SINCOSL_COS_HI 0
+#define SINCOSL_COS_LO 1
+#define SINCOSL_SIN_HI 2
+#define SINCOSL_SIN_LO 3
+extern const long double __sincosl_table[];
+
+long double
+__kernel_sinl(long double x, long double y, int iy)
+{
+  long double h, l, z, sin_l, cos_l_m1;
+  int64_t ix;
+  u_int32_t tix, hix, index;
+  GET_LDOUBLE_MSW64 (ix, x);
+  tix = ((u_int64_t)ix) >> 32;
+  tix &= ~0x80000000;                  /* tix = |x|'s high 32 bits */
+  if (tix < 0x3ffc3000)                        /* |x| < 0.1484375 */
+    {
+      /* Argument is small enough to approximate it by a Chebyshev
+        polynomial of degree 17.  */
+      if (tix < 0x3fc60000)            /* |x| < 2^-57 */
+       if (!((int)x)) return x;        /* generate inexact */
+      z = x * x;
+      return x + (x * (z*(SIN1+z*(SIN2+z*(SIN3+z*(SIN4+
+                      z*(SIN5+z*(SIN6+z*(SIN7+z*SIN8)))))))));
+    }
+  else
+    {
+      /* So that we don't have to use too large polynomial,  we find
+        l and h such that x = l + h,  where fabsl(l) <= 1.0/256 with 83
+        possible values for h.  We look up cosl(h) and sinl(h) in
+        pre-computed tables,  compute cosl(l) and sinl(l) using a
+        Chebyshev polynomial of degree 10(11) and compute
+        sinl(h+l) = sinl(h)cosl(l) + cosl(h)sinl(l).  */
+      index = 0x3ffe - (tix >> 16);
+      hix = (tix + (0x200 << index)) & (0xfffffc00 << index);
+      x = fabsl (x);
+      switch (index)
+       {
+       case 0: index = ((45 << 10) + hix - 0x3ffe0000) >> 8; break;
+       case 1: index = ((13 << 11) + hix - 0x3ffd0000) >> 9; break;
+       default:
+       case 2: index = (hix - 0x3ffc3000) >> 10; break;
+       }
+
+      SET_LDOUBLE_WORDS64(h, ((u_int64_t)hix) << 32, 0);
+      if (iy)
+       l = y - (h - x);
+      else
+       l = x - h;
+      z = l * l;
+      sin_l = l*(ONE+z*(SSIN1+z*(SSIN2+z*(SSIN3+z*(SSIN4+z*SSIN5)))));
+      cos_l_m1 = z*(SCOS1+z*(SCOS2+z*(SCOS3+z*(SCOS4+z*SCOS5))));
+      z = __sincosl_table [index + SINCOSL_SIN_HI]
+         + (__sincosl_table [index + SINCOSL_SIN_LO]
+            + (__sincosl_table [index + SINCOSL_SIN_HI] * cos_l_m1)
+            + (__sincosl_table [index + SINCOSL_COS_HI] * sin_l));
+      return (ix < 0) ? -z : z;
+    }
+}
index aecb20a..b3faa04 100644 (file)
@@ -80,3 +80,11 @@ do {                                                         \
   (d) = sh_u.value;                                            \
 } while (0)
 
+/* Get the least significant 64 bits of a long double mantissa.  */
+
+#define GET_LDOUBLE_LSW64(v,d)                                 \
+do {                                                           \
+  ieee854_long_double_shape_type sh_u;                         \
+  sh_u.value = (d);                                            \
+  (v) = sh_u.parts64.lsw;                                      \
+} while (0)
index 2d74e72..623b3b5 100644 (file)
@@ -23,9 +23,6 @@
 
 #include "math_private.h"
 
-/* Note: We should probably introduce __kernel_sincosl to speed things up,
-   because __kernel_{cos,sin}l sometimes compute both sine and cosine.  */
-
 void
 __sincosl (long double x, long double *sinx, long double *cosx)
 {
@@ -37,10 +34,7 @@ __sincosl (long double x, long double *sinx, long double *cosx)
   /* |x| ~< pi/4 */
   ix &= 0x7fffffffffffffffLL;
   if (ix <= 0x3ffe921fb54442d1LL)
-    {
-      *sinx = __kernel_sinl (x, 0.0, 0);
-      *cosx = __kernel_cosl (x, 0.0);
-    }
+    __kernel_sincosl (x, 0.0L, sinx, cosx, 0);
   else if (ix >= 0x7fff000000000000LL)
     {
       /* sin(Inf or NaN) is NaN */
@@ -56,20 +50,20 @@ __sincosl (long double x, long double *sinx, long double *cosx)
       switch (n & 3)
        {
        case 0:
-         *sinx = __kernel_sinl (y[0], y[1], 1);
-         *cosx = __kernel_cosl (y[0], y[1]);
+         __kernel_sincosl (y[0], y[1], sinx, cosx, 1);
          break;
        case 1:
-         *sinx = __kernel_cosl (y[0], y[1]);
-         *cosx = -__kernel_sinl (y[0], y[1], 1);
+         __kernel_sincosl (y[0], y[1], cosx, sinx, 1);
+         *cosx = -*cosx;
          break;
        case 2:
-         *sinx = -__kernel_sinl (y[0], y[1], 1);
-         *cosx = -__kernel_cosl (y[0], y[1]);
+         __kernel_sincosl (y[0], y[1], sinx, cosx, 1);
+         *sinx = -*sinx;
+         *cosx = -*cosx;
          break;
        default:
-         *sinx = -__kernel_cosl (y[0], y[1]);
-         *cosx = __kernel_sinl (y[0], y[1], 1);
+         __kernel_sincosl (y[0], y[1], cosx, sinx, 1);
+         *sinx = -*sinx;
          break;
        }
     }