1 /* Quad-precision floating point sine and cosine on <-pi/4,pi/4>.
2 Copyright (C) 1999,2004,2006 Free Software Foundation, Inc.
3 This file is part of the GNU C Library.
4 Contributed by Jakub Jelinek <jj@ultra.linux.cz>
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, write to the Free
18 Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
22 #include "math_private.h"
24 static const long double c
[] = {
26 1.00000000000000000000000000000000000E+00L, /* 3fff0000000000000000000000000000 */
28 /* cos x ~ ONE + x^2 ( SCOS1 + SCOS2 * x^2 + ... + SCOS4 * x^6 + SCOS5 * x^8 )
35 -5.00000000000000000000000000000000000E-01L, /* bffe0000000000000000000000000000 */
36 4.16666666666666666666666666556146073E-02L, /* 3ffa5555555555555555555555395023 */
37 -1.38888888888888888888309442601939728E-03L, /* bff56c16c16c16c16c16a566e42c0375 */
38 2.48015873015862382987049502531095061E-05L, /* 3fefa01a01a019ee02dcf7da2d6d5444 */
39 -2.75573112601362126593516899592158083E-07L, /* bfe927e4f5dce637cb0b54908754bde0 */
41 /* cos x ~ ONE + x^2 ( COS1 + COS2 * x^2 + ... + COS7 * x^12 + COS8 * x^14 )
51 -4.99999999999999999999999999999999759E-01L, /* bffdfffffffffffffffffffffffffffb */
52 4.16666666666666666666666666651287795E-02L, /* 3ffa5555555555555555555555516f30 */
53 -1.38888888888888888888888742314300284E-03L, /* bff56c16c16c16c16c16c16a463dfd0d */
54 2.48015873015873015867694002851118210E-05L, /* 3fefa01a01a01a01a0195cebe6f3d3a5 */
55 -2.75573192239858811636614709689300351E-07L, /* bfe927e4fb7789f5aa8142a22044b51f */
56 2.08767569877762248667431926878073669E-09L, /* 3fe21eed8eff881d1e9262d7adff4373 */
57 -1.14707451049343817400420280514614892E-11L, /* bfda9397496922a9601ed3d4ca48944b */
58 4.77810092804389587579843296923533297E-14L, /* 3fd2ae5f8197cbcdcaf7c3fb4523414c */
60 /* sin x ~ ONE * x + x^3 ( SSIN1 + SSIN2 * x^2 + ... + SSIN4 * x^6 + SSIN5 * x^8 )
67 -1.66666666666666666666666666666666659E-01L, /* bffc5555555555555555555555555555 */
68 8.33333333333333333333333333146298442E-03L, /* 3ff81111111111111111111110fe195d */
69 -1.98412698412698412697726277416810661E-04L, /* bff2a01a01a01a01a019e7121e080d88 */
70 2.75573192239848624174178393552189149E-06L, /* 3fec71de3a556c640c6aaa51aa02ab41 */
71 -2.50521016467996193495359189395805639E-08L, /* bfe5ae644ee90c47dc71839de75b2787 */
73 /* sin x ~ ONE * x + x^3 ( SIN1 + SIN2 * x^2 + ... + SIN7 * x^12 + SIN8 * x^14 )
83 -1.66666666666666666666666666666666538e-01L, /* bffc5555555555555555555555555550 */
84 8.33333333333333333333333333307532934e-03L, /* 3ff811111111111111111111110e7340 */
85 -1.98412698412698412698412534478712057e-04L, /* bff2a01a01a01a01a01a019e7a626296 */
86 2.75573192239858906520896496653095890e-06L, /* 3fec71de3a556c7338fa38527474b8f5 */
87 -2.50521083854417116999224301266655662e-08L, /* bfe5ae64567f544e16c7de65c2ea551f */
88 1.60590438367608957516841576404938118e-10L, /* 3fde6124613a811480538a9a41957115 */
89 -7.64716343504264506714019494041582610e-13L, /* bfd6ae7f3d5aef30c7bc660b060ef365 */
90 2.81068754939739570236322404393398135e-15L, /* 3fce9510115aabf87aceb2022a9a9180 */
93 #define SINCOSL_COS_HI 0
94 #define SINCOSL_COS_LO 1
95 #define SINCOSL_SIN_HI 2
96 #define SINCOSL_SIN_LO 3
97 extern const long double __sincosl_table
[];
100 __kernel_sincosl(long double x
, long double y
, long double *sinx
, long double *cosx
, int iy
)
102 long double h
, l
, z
, sin_l
, cos_l_m1
;
104 u_int32_t tix
, hix
, index
;
105 GET_LDOUBLE_MSW64 (ix
, x
);
106 tix
= ((u_int64_t
)ix
) >> 32;
107 tix
&= ~0x80000000; /* tix = |x|'s high 32 bits */
108 if (tix
< 0x3fc30000) /* |x| < 0.1484375 */
110 /* Argument is small enough to approximate it by a Chebyshev
111 polynomial of degree 16(17). */
112 if (tix
< 0x3c600000) /* |x| < 2^-57 */
113 if (!((int)x
)) /* generate inexact */
120 *sinx
= x
+ (x
* (z
*(SIN1
+z
*(SIN2
+z
*(SIN3
+z
*(SIN4
+
121 z
*(SIN5
+z
*(SIN6
+z
*(SIN7
+z
*SIN8
)))))))));
122 *cosx
= ONE
+ (z
*(COS1
+z
*(COS2
+z
*(COS3
+z
*(COS4
+
123 z
*(COS5
+z
*(COS6
+z
*(COS7
+z
*COS8
))))))));
127 /* So that we don't have to use too large polynomial, we find
128 l and h such that x = l + h, where fabsl(l) <= 1.0/256 with 83
129 possible values for h. We look up cosl(h) and sinl(h) in
130 pre-computed tables, compute cosl(l) and sinl(l) using a
131 Chebyshev polynomial of degree 10(11) and compute
132 sinl(h+l) = sinl(h)cosl(l) + cosl(h)sinl(l) and
133 cosl(h+l) = cosl(h)cosl(l) - sinl(h)sinl(l). */
135 tix
= ((six
- 0x3ff00000) >> 4) + 0x3fff0000;
136 index
= 0x3ffe - (tix
>> 16);
137 hix
= (tix
+ (0x200 << index
)) & (0xfffffc00 << index
);
141 case 0: index
= ((45 << 10) + hix
- 0x3ffe0000) >> 8; break;
142 case 1: index
= ((13 << 11) + hix
- 0x3ffd0000) >> 9; break;
144 case 2: index
= (hix
- 0x3ffc3000) >> 10; break;
146 hix
= (hix
<< 4) & 0x3fffffff;
148 The following should work for double but generates the wrong index.
149 For now the code above converts double to ieee extended to compute
150 the index back to double for the h value.
153 index = 0x3fe - (tix >> 20);
154 hix = (tix + (0x2000 << index)) & (0xffffc000 << index);
158 case 0: index = ((45 << 14) + hix - 0x3fe00000) >> 12; break;
159 case 1: index = ((13 << 15) + hix - 0x3fd00000) >> 13; break;
161 case 2: index = (hix - 0x3fc30000) >> 14; break;
164 SET_LDOUBLE_WORDS64(h
, ((u_int64_t
)hix
) << 32, 0);
170 sin_l
= l
*(ONE
+z
*(SSIN1
+z
*(SSIN2
+z
*(SSIN3
+z
*(SSIN4
+z
*SSIN5
)))));
171 cos_l_m1
= z
*(SCOS1
+z
*(SCOS2
+z
*(SCOS3
+z
*(SCOS4
+z
*SCOS5
))));
172 z
= __sincosl_table
[index
+ SINCOSL_SIN_HI
]
173 + (__sincosl_table
[index
+ SINCOSL_SIN_LO
]
174 + (__sincosl_table
[index
+ SINCOSL_SIN_HI
] * cos_l_m1
)
175 + (__sincosl_table
[index
+ SINCOSL_COS_HI
] * sin_l
));
176 *sinx
= (ix
< 0) ? -z
: z
;
177 *cosx
= __sincosl_table
[index
+ SINCOSL_COS_HI
]
178 + (__sincosl_table
[index
+ SINCOSL_COS_LO
]
179 - (__sincosl_table
[index
+ SINCOSL_SIN_HI
] * sin_l
180 - __sincosl_table
[index
+ SINCOSL_COS_HI
] * cos_l_m1
));