Home | History | Annotate | Line # | Download | only in src 
      1 /*-
      2  * ====================================================
      3  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
      4  *
      5  * Developed at SunPro, a Sun Microsystems, Inc. business.
      6  * Permission to use, copy, modify, and distribute this
      7  * software is freely granted, provided that this notice
      8  * is preserved.
      9  * ====================================================
     10  *
     11  * s_sin.c and s_cos.c merged by Steven G. Kargl.  Descriptions of the
     12  * algorithms are contained in the original files.
     13  */
     14 
     15 #include <sys/cdefs.h>
     16 #if defined(LIBM_SCCS) && !defined(lint)
     17 __RCSID("$NetBSD: s_sincos.c,v 1.5 2022/08/29 01:48:34 riastradh Exp $");
     18 #endif
     19 #if 0
     20 __FBSDID("$FreeBSD: head/lib/msun/src/s_sincos.c 319047 2017-05-28 06:13:38Z mmel $");
     21 #endif
     22 
     23 #include "namespace.h"
     24 #include <float.h>
     25 
     26 #include "math.h"
     27 // #define INLINE_REM_PIO2
     28 #include "math_private.h"
     29 // #include "e_rem_pio2.c"
     30 #include "k_sincos.h"
     31 
     32 #ifdef __weak_alias
     33 __weak_alias(sincos,_sincos)
     34 #endif
     35 
     36 void
     37 sincos(double x, double *sn, double *cs)
     38 {
     39 	double y[2];
     40 	int32_t n, ix;
     41 
     42 	/* High word of x. */
     43 	GET_HIGH_WORD(ix, x);
     44 
     45 	/* |x| ~< pi/4 */
     46 	ix &= 0x7fffffff;
     47 	if (ix <= 0x3fe921fb) {
     48 		if (ix < 0x3e400000) {		/* |x| < 2**-27 */
     49 			if ((int)x == 0) {	/* Generate inexact. */
     50 				*sn = x;
     51 				*cs = 1;
     52 				return;
     53 			}
     54 		}
     55 		__kernel_sincos(x, 0, 0, sn, cs);
     56 		return;
     57 	}
     58 
     59 	/* If x = Inf or NaN, then sin(x) = NaN and cos(x) = NaN. */
     60 	if (ix >= 0x7ff00000) {
     61 		*sn = x - x;
     62 		*cs = x - x;
     63 		return;
     64 	}
     65 
     66 	/* Argument reduction. */
     67 	n = __ieee754_rem_pio2(x, y);
     68 
     69 	switch(n & 3) {
     70 	case 0:
     71 		__kernel_sincos(y[0], y[1], 1, sn, cs);
     72 		break;
     73 	case 1:
     74 		__kernel_sincos(y[0], y[1], 1, cs, sn);
     75 		*cs = -*cs;
     76 		break;
     77 	case 2:
     78 		__kernel_sincos(y[0], y[1], 1, sn, cs);
     79 		*sn = -*sn;
     80 		*cs = -*cs;
     81 		break;
     82 	default:
     83 		__kernel_sincos(y[0], y[1], 1, cs, sn);
     84 		*sn = -*sn;
     85 	}
     86 }
     87 
     88 #if !defined(__HAVE_LONG_DOUBLE)
     89 __weak_alias(sincosl, sincos);
     90 #endif
     91