Home | History | Annotate | Line # | Download | only in ld80 
      1 /*-
      2  * Copyright (c) 2017, 2023 Steven G. Kargl
      3  * All rights reserved.
      4  *
      5  * Redistribution and use in source and binary forms, with or without
      6  * modification, are permitted provided that the following conditions
      7  * are met:
      8  * 1. Redistributions of source code must retain the above copyright
      9  *    notice unmodified, this list of conditions, and the following
     10  *    disclaimer.
     11  * 2. Redistributions in binary form must reproduce the above copyright
     12  *    notice, this list of conditions and the following disclaimer in the
     13  *    documentation and/or other materials provided with the distribution.
     14  *
     15  * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
     16  * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
     17  * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
     18  * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
     19  * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
     20  * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
     21  * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
     22  * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
     23  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
     24  * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
     25  */
     26 
     27 /*
     28  * See ../src/s_sinpi.c for implementation details.
     29  */
     30 
     31 #include <stdint.h>
     32 
     33 #include "math.h"
     34 #include "math_private.h"
     35 
     36 static const union ieee_ext_u
     37 pi_hi_u = LD80C(0xc90fdaa200000000,   1, 3.14159265346825122833e+00L),
     38 pi_lo_u = LD80C(0x85a308d313198a2e, -33, 1.21542010130123852029e-10L);
     39 #define	pi_hi	(pi_hi_u.extu_ld)
     40 #define	pi_lo	(pi_lo_u.extu_ld)
     41 
     42 #include "k_cospil.h"
     43 #include "k_sinpil.h"
     44 
     45 static volatile const double vzero = 0;
     46 
     47 long double
     48 sinpil(long double x)
     49 {
     50 	long double ax, hi, lo, s;
     51 	uint64_t lx;
     52 	uint32_t j0;
     53 	uint16_t hx, ix;
     54 
     55 	EXTRACT_LDBL80_WORDS(hx, lx, x);
     56 	ix = hx & 0x7fff;
     57 	INSERT_LDBL80_WORDS(ax, ix, lx);
     58 
     59 	ENTERI();
     60 
     61 	if (ix < 0x3fff) {			/* |x| < 1 */
     62 		if (ix < 0x3ffd) {		/* |x| < 0.25 */
     63 			if (ix < 0x3fdd) {	/* |x| < 0x1p-34 */
     64 				if (x == 0)
     65 					RETURNI(x);
     66 				INSERT_LDBL80_WORDS(hi, hx,
     67 				    lx & 0xffffffff00000000ull);
     68 				hi *= 0x1p63L;
     69 				lo = x * 0x1p63L - hi;
     70 				s = (pi_lo + pi_hi) * lo + pi_lo * hi +
     71 				    pi_hi * hi;
     72 				RETURNI(s * 0x1p-63L);
     73 			}
     74 			s = __kernel_sinpil(ax);
     75 			RETURNI((hx & 0x8000) ? -s : s);
     76 		}
     77 
     78 		if (ix < 0x3ffe)			/* |x| < 0.5 */
     79 			s = __kernel_cospil(0.5 - ax);
     80 		else if (lx < 0xc000000000000000ull)	/* |x| < 0.75 */
     81 			s = __kernel_cospil(ax - 0.5);
     82 		else
     83 			s = __kernel_sinpil(1 - ax);
     84 		RETURNI((hx & 0x8000) ? -s : s);
     85 	}
     86 
     87 	if (ix < 0x403e) {			/* 1 <= |x| < 0x1p63 */
     88 		FFLOORL80(x, j0, ix, lx);	/* Integer part of ax. */
     89 		ax -= x;
     90 		EXTRACT_LDBL80_WORDS(ix, lx, ax);
     91 
     92 		if (ix == 0) {
     93 			s = 0;
     94 		} else {
     95 			if (ix < 0x3ffe) {		/* |x| < 0.5 */
     96 				if (ix < 0x3ffd)	/* |x| < 0.25 */
     97 					s = __kernel_sinpil(ax);
     98 				else
     99 					s = __kernel_cospil(0.5 - ax);
    100 			} else {
    101 							/* |x| < 0.75 */
    102 				if (lx < 0xc000000000000000ull)
    103 					s = __kernel_cospil(ax - 0.5);
    104 				else
    105 					s = __kernel_sinpil(1 - ax);
    106 			}
    107 
    108 			if (j0 > 40)
    109 				x -= 0x1p40;
    110 			if (j0 > 30)
    111 				x -= 0x1p30;
    112 			j0 = (uint32_t)x;
    113 			if (j0 & 1) s = -s;
    114 		}
    115 		RETURNI((hx & 0x8000) ? -s : s);
    116 	}
    117 
    118 	/* x = +-inf or nan. */
    119 	if (ix >= 0x7fff)
    120 		RETURNI(vzero / vzero);
    121 
    122 	/*
    123 	 * |x| >= 0x1p63 is always an integer, so return +-0.
    124 	 */
    125 	RETURNI(copysignl(0, x));
    126 }
    127