Home | History | Annotate | Line # | Download | only in src 
      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 #include <sys/cdefs.h>
     28 
     29 /*
     30  * See ../src/s_sinpi.c for implementation details.
     31  */
     32 
     33 #define	INLINE_KERNEL_SINDF
     34 #define	INLINE_KERNEL_COSDF
     35 
     36 #include "namespace.h"
     37 __weak_alias(sinpif, _sinpif)
     38 
     39 #include "math.h"
     40 #include "math_private.h"
     41 #include "k_cosdf.c"
     42 #include "k_sindf.c"
     43 
     44 #define	__kernel_cospif(x)	(__kernel_cosdf(M_PI * (x)))
     45 #define	__kernel_sinpif(x)	(__kernel_sindf(M_PI * (x)))
     46 
     47 static const float
     48 pi_hi =  3.14160156e+00F,	/* 0x40491000 */
     49 pi_lo = -8.90890988e-06F;	/* 0xb715777a */
     50 
     51 static volatile const float vzero = 0;
     52 
     53 float
     54 sinpif(float x)
     55 {
     56 	float ax, hi, lo, s;
     57 	uint32_t hx, ix, j0;
     58 
     59 	GET_FLOAT_WORD(hx, x);
     60 	ix = hx & 0x7fffffff;
     61 	SET_FLOAT_WORD(ax, ix);
     62 
     63 	if (ix < 0x3f800000) {			/* |x| < 1 */
     64 		if (ix < 0x3e800000) {		/* |x| < 0.25 */
     65 	 		if (ix < 0x38800000) {	/* |x| < 0x1p-14 */
     66 				if (x == 0)
     67 					return (x);
     68 				SET_FLOAT_WORD(hi, hx & 0xffff0000);
     69 				hi *= 0x1p23F;
     70 				lo = x * 0x1p23F - hi;
     71 				s = (pi_lo + pi_hi) * lo + pi_lo * hi +
     72 				    pi_hi * hi;
     73 				return (s * 0x1p-23F);
     74 			}
     75 
     76 			s = __kernel_sinpif(ax);
     77 			return ((hx & 0x80000000) ? -s : s);
     78 		}
     79 
     80 		if (ix < 0x3f000000)		/* |x| < 0.5 */
     81 			s = __kernel_cospif(0.5F - ax);
     82 		else if (ix < 0x3f400000)	/* |x| < 0.75 */
     83 			s = __kernel_cospif(ax - 0.5F);
     84 		else
     85 			s = __kernel_sinpif(1 - ax);
     86 		return ((hx & 0x80000000) ? -s : s);
     87 	}
     88 
     89 	if (ix < 0x4b000000) {		/* 1 <= |x| < 0x1p23 */
     90 		FFLOORF(x, j0, ix);	/* Integer part of ax. */
     91 		ax -= x;
     92 		GET_FLOAT_WORD(ix, ax);
     93 
     94 		if (ix == 0)
     95 			s = 0;
     96 		else {
     97 			if (ix < 0x3f000000) {		/* |x| < 0.5 */
     98 				if (ix < 0x3e800000)	/* |x| < 0.25 */
     99 					s = __kernel_sinpif(ax);
    100 				else
    101 					s = __kernel_cospif(0.5F - ax);
    102 			} else {
    103 				if (ix < 0x3f400000)	/* |x| < 0.75 */
    104 					s = __kernel_cospif(ax - 0.5F);
    105 				else
    106 					s = __kernel_sinpif(1 - ax);
    107 			}
    108 
    109 			j0 = (uint32_t)x;
    110 			s = (j0 & 1) ? -s : s;
    111 		}
    112 		return ((hx & 0x80000000) ? -s : s);
    113 	}
    114 
    115 	/* x = +-inf or nan. */
    116 	if (ix >= 0x7f800000)
    117 		return (vzero / vzero);
    118 
    119 	/*
    120 	 * |x| >= 0x1p23 is always an integer, so return +-0.
    121 	 */
    122 	return (copysignf(0, x));
    123 }
    124