2 * IBM Accurate Mathematical Library

3 * written by International Business Machines Corp.

4 * Copyright (C) 2001-2020 Free Software Foundation, Inc.

6 * This program is free software; you can redistribute it and/or modify

7 * it under the terms of the GNU Lesser General Public License as published by

8 * the Free Software Foundation; either version 2.1 of the License, or

9 * (at your option) any later version.

11 * This program 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

14 * GNU Lesser General Public License for more details.

16 * You should have received a copy of the GNU Lesser General Public License

17 * along with this program; if not, see <https://www.gnu.org/licenses/>.

23 /* FUNCTIONS: usin */

24 /* ucos */

25 /* FILES NEEDED: dla.h endian.h mpa.h mydefs.h usncs.h */

28 /* An ultimate sin and cos routine. Given an IEEE double machine number x */

29 /* it computes sin(x) or cos(x) with ~0.55 ULP. */

30 /* Assumption: Machine arithmetic operations are performed in */

31 /* round to nearest mode of IEEE 754 standard. */

49 /* Helper macros to compute sin of the input values. */

54 /* The computed polynomial is a variation of the Taylor series expansion for

57 a - a^3/3! + a^5/5! - a^7/7! + a^9/9! + (1 - a^2) * da / 2

59 The constants s1, s2, s3, etc. are pre-computed values of 1/3!, 1/5! and so

60 on. The result is returned to LHS. */

63 double t = ((POLYNOMIAL (xx) * (a) - 0.5 * (da)) * (xx) + (da)); \

64 double res = (a) + t; \

70 int4 k = u.i[LOW_HALF] << 2; \

71 sn = __sincostab.x[k]; \

72 ssn = __sincostab.x[k + 1]; \

73 cs = __sincostab.x[k + 2]; \

74 ccs = __sincostab.x[k + 3]; \

77 #ifndef SECTION

78 # define SECTION

81 extern const union

[ 880 ]; 83 int4 i

84 double x [ 440 ];

87 static const double

94 int __branred ( double x , double * a , double * aa ); aa

96 /* Given a number partitioned into X and DX, this function computes the cosine

97 of the number by combining the sin and cos of X (as computed by a variation

98 of the Taylor series) with the values looked up from the sin/cos table to

99 get the result. */

100 static __always_inline double __always_inline

101 do_cos ( double x , double dx ) dx

105 if ( x < 0 )

108 u . x = big + fabs ( x ); big

109 x = fabs ( x ) - ( u . x - big ) + dx ; bigdx

111 double xx , s , sn , ssn , c , cs , ccs , cor ; xxsnssncsccscor

= x * x ; 112 xx

113 s = x + x * xx * ( sn3 + xx * sn5 ); xxsn3xxsn5

114 c = xx * ( cs2 + xx * ( cs4 + xx * cs6 )); xxcs2xxcs4xxcs6

115 SINCOS_TABLE_LOOKUP ( u , sn , ssn , cs , ccs ); snssncsccs

= ( ccs - s * ssn - cs * c ) - sn * s ; 116 corccsssncssn

117 return cs + cor ; cscor

120 /* Given a number partitioned into X and DX, this function computes the sine of

121 the number by combining the sin and cos of X (as computed by a variation of

122 the Taylor series) with the values looked up from the sin/cos table to get

123 the result. */

124 static __always_inline double __always_inline

125 do_sin ( double x , double dx ) dx

127 double xold = x ; xold

128 /* Max ULP is 0.501 if |x| < 0.126, otherwise ULP is 0.518. */

129 if ( fabs ( x ) < 0.126 )

130 return TAYLOR_SIN ( x * x , x , dx ); dx

134 if ( x <= 0 )

136 u . x = big + fabs ( x ); big

137 x = fabs ( x ) - ( u . x - big ); big

139 double xx , s , sn , ssn , c , cs , ccs , cor ; xxsnssncsccscor

= x * x ; 140 xx

141 s = x + ( dx + x * xx * ( sn3 + xx * sn5 )); dxxxsn3xxsn5

142 c = x * dx + xx * ( cs2 + xx * ( cs4 + xx * cs6 )); dxxxcs2xxcs4xxcs6

143 SINCOS_TABLE_LOOKUP ( u , sn , ssn , cs , ccs ); snssncsccs

= ( ssn + s * ccs - sn * c ) + cs * s ; 144 corssnccssncs

145 return copysign ( sn + cor , xold ); sncorxold

148 /* Reduce range of x to within PI/2 with abs (x) < 105414350. The high part

149 is written to *a, the low part to *da. Range reduction is accurate to 136

150 bits so that when x is large and *a very close to zero, all 53 bits of *a

151 are correct. */

152 static __always_inline int4 __always_inline int4

153 reduce_sincos ( double x , double * a , double * da ) da

157 double t = ( x * hpinv + toint ); hpinvtoint

158 double xn = t - toint ; xntoint

159 v . x = t ;

160 double y = ( x - xn * mp1 ) - xn * mp2 ; xnmp1xnmp2

= v . i [ LOW_HALF ] & 3 ; 161 int4 nLOW_HALF

163 double b , db , t1 , t2 ; dbt1t2

= xn * pp3 ; 164 t1xnpp3

= y - t1 ; 165 t2t1

= ( y - t2 ) - t1 ; 166 dbt2t1

= xn * pp4 ; 168 t1xnpp4

169 b = t2 - t1 ; t2t1

+= ( t2 - b ) - t1 ; 170 dbt2t1

172 * a = b ;

173 * da = db ; dadb

174 return n ;

177 /* Compute sin or cos (A + DA) for the given quadrant N. */

178 static __always_inline double __always_inline

179 do_sincos ( double a , double da , int4 n ) daint4 n

181 double retval ; retval

183 if ( n & 1 )

184 /* Max ULP is 0.513. */

= do_cos ( a , da ); 185 retvalda

186 else

187 /* Max ULP is 0.501 if xx < 0.01588, otherwise ULP is 0.518. */

= do_sin ( a , da ); 188 retvalda

190 return ( n & 2 ) ? - retval : retval ; retvalretval

195 /* An ultimate sin routine. Given an IEEE double machine number x */

196 /* it computes the rounded value of sin(x). */

199 double

201 __sin ( double x )

203 double t , a , da ; da

, m , n ; 205 int4 k

206 double retval = 0 ; retval

210 u . x = x ;

211 m = u . i [ HIGH_HALF ]; HIGH_HALF

212 k = 0x7fffffff & m ; /* no sign */

213 if ( k < 0x3e500000 ) /* if x->0 =>sin(x)=x */

215 math_check_force_underflow ( x );

219 else if ( k < 0x3feb6000 )

221 /* Max ULP is 0.548. */

= do_sin ( x , 0 ); 222 retval

223 } /* else if (k < 0x3feb6000) */

226 else if ( k < 0x400368fd )

228 t = hp0 - fabs ( x ); hp0

229 /* Max ULP is 0.51. */

= copysign ( do_cos ( t , hp1 ), x ); 230 retvalhp1

231 } /* else if (k < 0x400368fd) */

233 /*-------------------------- 2.426265<|x|< 105414350 ----------------------*/

234 else if ( k < 0x419921FB )

236 n = reduce_sincos ( x , & a , & da ); da

= do_sincos ( a , da , n ); 237 retvalda

238 } /* else if (k < 0x419921FB ) */

241 else if ( k < 0x7ff00000 )

243 n = __branred ( x , & a , & da ); da

= do_sincos ( a , da , n ); 244 retvalda

247 else

249 if ( k == 0x7ff00000 && u . i [ LOW_HALF ] == 0 ) LOW_HALF

250 __set_errno ( EDOM ); EDOM

= x / x ; 251 retval

254 return retval ; retval

259 /* An ultimate cos routine. Given an IEEE double machine number x */

260 /* it computes the rounded value of cos(x). */

263 double

265 __cos ( double x )

267 double y , a , da ; da

, m , n ; 269 int4 k

271 double retval = 0 ; retval

275 u . x = x ;

276 m = u . i [ HIGH_HALF ]; HIGH_HALF

277 k = 0x7fffffff & m ;

280 if ( k < 0x3e400000 )

283 else if ( k < 0x3feb6000 )

285 /* Max ULP is 0.51. */

= do_cos ( x , 0 ); 286 retval

287 } /* else if (k < 0x3feb6000) */

289 else if ( k < 0x400368fd )

291 y = hp0 - fabs ( x ); hp0

292 a = y + hp1 ; hp1

= ( y - a ) + hp1 ; 293 dahp1

294 /* Max ULP is 0.501 if xx < 0.01588 or 0.518 otherwise.

295 Range reduction uses 106 bits here which is sufficient. */

= do_sin ( a , da ); 296 retvalda

297 } /* else if (k < 0x400368fd) */

299 else if ( k < 0x419921FB )

300 { /* 2.426265<|x|< 105414350 */

301 n = reduce_sincos ( x , & a , & da ); da

= do_sincos ( a , da , n + 1 ); 302 retvalda

303 } /* else if (k < 0x419921FB ) */

305 /* 105414350 <|x| <2^1024 */

306 else if ( k < 0x7ff00000 )

308 n = __branred ( x , & a , & da ); da

= do_sincos ( a , da , n + 1 ); 309 retvalda

312 else

314 if ( k == 0x7ff00000 && u . i [ LOW_HALF ] == 0 ) LOW_HALF

315 __set_errno ( EDOM ); EDOM

= x / x ; /* |x| > 2^1024 */ 316 retval

319 return retval ; retval

323 libm_alias_double ( __cos , cos ) __coscos