src/libm/k_sin.c
author Ryan C. Gordon <icculus@icculus.org>
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  * ====================================================
     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 
    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  */
    39 
    40 #include "math_libm.h"
    41 #include "math_private.h"
    42 
    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 */
    51 
    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 }