src/libm/k_sin.c
 author Ryan C. Gordon Thu, 21 Apr 2016 03:16:44 -0400 changeset 11729 d1ce8396c356 parent 11683 48bcba563d9c permissions -rw-r--r--
Initial shot at a renderer target for Apple's Metal API.

This isn't complete, but is enough to run testsprite2. It's currently
Mac-only; with a little work to figure out how to properly glue in a Metal
layer to a UIView, this will likely work on iOS, too.

This is only wired up to the configure script right now, and disabled by
default. CMake and Xcode still need their bits filled in as appropriate.
1 /*
2  * ====================================================
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  */
12 /* __kernel_sin( x, y, iy)
13  * kernel sin function on [-pi/4, pi/4], pi/4 ~ 0.7854
14  * Input x is assumed to be bounded by ~pi/4 in magnitude.
15  * Input y is the tail of x.
16  * Input iy indicates whether y is 0. (if iy=0, y assume to be 0).
17  *
18  * Algorithm
19  *	1. Since sin(-x) = -sin(x), we need only to consider positive x.
20  *	2. if x < 2^-27 (hx<0x3e400000 0), return x with inexact if x!=0.
21  *	3. sin(x) is approximated by a polynomial of degree 13 on
22  *	   [0,pi/4]
23  *		  	         3            13
24  *	   	sin(x) ~ x + S1*x + ... + S6*x
25  *	   where
26  *
27  * 	|sin(x)         2     4     6     8     10     12  |     -58
28  * 	|----- - (1+S1*x +S2*x +S3*x +S4*x +S5*x  +S6*x   )| <= 2
29  * 	|  x 					           |
30  *
31  *	4. sin(x+y) = sin(x) + sin'(x')*y
32  *		    ~ sin(x) + (1-x*x/2)*y
33  *	   For better accuracy, let
34  *		     3      2      2      2      2
35  *		r = x *(S2+x *(S3+x *(S4+x *(S5+x *S6))))
36  *	   then                   3    2
37  *		sin(x) = x + (S1*x + (x *(r-y/2)+y))
38  */
40 #include "math_libm.h"
41 #include "math_private.h"
43 static const double
44 half =  5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
45 S1  = -1.66666666666666324348e-01, /* 0xBFC55555, 0x55555549 */
46 S2  =  8.33333333332248946124e-03, /* 0x3F811111, 0x1110F8A6 */
47 S3  = -1.98412698298579493134e-04, /* 0xBF2A01A0, 0x19C161D5 */
48 S4  =  2.75573137070700676789e-06, /* 0x3EC71DE3, 0x57B1FE7D */
49 S5  = -2.50507602534068634195e-08, /* 0xBE5AE5E6, 0x8A2B9CEB */
50 S6  =  1.58969099521155010221e-10; /* 0x3DE5D93A, 0x5ACFD57C */
52 double attribute_hidden __kernel_sin(double x, double y, int iy)
53 {
54 	double z,r,v;
55 	int32_t ix;
56 	GET_HIGH_WORD(ix,x);
57 	ix &= 0x7fffffff;			/* high word of x */
58 	if(ix<0x3e400000)			/* |x| < 2**-27 */
59 	   {if((int)x==0) return x;}		/* generate inexact */
60 	z	=  x*x;
61 	v	=  z*x;
62 	r	=  S2+z*(S3+z*(S4+z*(S5+z*S6)));
63 	if(iy==0) return x+v*(S1+z*r);
64 	else      return x-((z*(half*y-v*r)-y)-v*S1);
65 }