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.
slouken@2756
     1
/*
slouken@2756
     2
 * ====================================================
slouken@2756
     3
 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
slouken@2756
     4
 *
slouken@2756
     5
 * Developed at SunPro, a Sun Microsystems, Inc. business.
slouken@2756
     6
 * Permission to use, copy, modify, and distribute this
slouken@2756
     7
 * software is freely granted, provided that this notice
slouken@2756
     8
 * is preserved.
slouken@2756
     9
 * ====================================================
slouken@2756
    10
 */
slouken@2756
    11
slouken@2756
    12
/* __kernel_sin( x, y, iy)
slouken@2756
    13
 * kernel sin function on [-pi/4, pi/4], pi/4 ~ 0.7854
slouken@2756
    14
 * Input x is assumed to be bounded by ~pi/4 in magnitude.
slouken@2756
    15
 * Input y is the tail of x.
slouken@2756
    16
 * Input iy indicates whether y is 0. (if iy=0, y assume to be 0).
slouken@2756
    17
 *
slouken@2756
    18
 * Algorithm
slouken@2756
    19
 *	1. Since sin(-x) = -sin(x), we need only to consider positive x.
slouken@2756
    20
 *	2. if x < 2^-27 (hx<0x3e400000 0), return x with inexact if x!=0.
slouken@2756
    21
 *	3. sin(x) is approximated by a polynomial of degree 13 on
slouken@2756
    22
 *	   [0,pi/4]
slouken@2756
    23
 *		  	         3            13
slouken@2756
    24
 *	   	sin(x) ~ x + S1*x + ... + S6*x
slouken@2756
    25
 *	   where
slouken@2756
    26
 *
slouken@2756
    27
 * 	|sin(x)         2     4     6     8     10     12  |     -58
slouken@2756
    28
 * 	|----- - (1+S1*x +S2*x +S3*x +S4*x +S5*x  +S6*x   )| <= 2
slouken@2756
    29
 * 	|  x 					           |
slouken@2756
    30
 *
slouken@2756
    31
 *	4. sin(x+y) = sin(x) + sin'(x')*y
slouken@2756
    32
 *		    ~ sin(x) + (1-x*x/2)*y
slouken@2756
    33
 *	   For better accuracy, let
slouken@2756
    34
 *		     3      2      2      2      2
slouken@2756
    35
 *		r = x *(S2+x *(S3+x *(S4+x *(S5+x *S6))))
slouken@2756
    36
 *	   then                   3    2
slouken@2756
    37
 *		sin(x) = x + (S1*x + (x *(r-y/2)+y))
slouken@2756
    38
 */
slouken@2756
    39
slouken@6044
    40
#include "math_libm.h"
slouken@2756
    41
#include "math_private.h"
slouken@2756
    42
slouken@2756
    43
static const double
slouken@11683
    44
half =  5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
slouken@11683
    45
S1  = -1.66666666666666324348e-01, /* 0xBFC55555, 0x55555549 */
slouken@11683
    46
S2  =  8.33333333332248946124e-03, /* 0x3F811111, 0x1110F8A6 */
slouken@11683
    47
S3  = -1.98412698298579493134e-04, /* 0xBF2A01A0, 0x19C161D5 */
slouken@11683
    48
S4  =  2.75573137070700676789e-06, /* 0x3EC71DE3, 0x57B1FE7D */
slouken@11683
    49
S5  = -2.50507602534068634195e-08, /* 0xBE5AE5E6, 0x8A2B9CEB */
slouken@11683
    50
S6  =  1.58969099521155010221e-10; /* 0x3DE5D93A, 0x5ACFD57C */
slouken@2756
    51
slouken@11683
    52
double attribute_hidden __kernel_sin(double x, double y, int iy)
slouken@2756
    53
{
slouken@11683
    54
	double z,r,v;
slouken@11683
    55
	int32_t ix;
slouken@11683
    56
	GET_HIGH_WORD(ix,x);
slouken@11683
    57
	ix &= 0x7fffffff;			/* high word of x */
slouken@11683
    58
	if(ix<0x3e400000)			/* |x| < 2**-27 */
slouken@11683
    59
	   {if((int)x==0) return x;}		/* generate inexact */
slouken@11683
    60
	z	=  x*x;
slouken@11683
    61
	v	=  z*x;
slouken@11683
    62
	r	=  S2+z*(S3+z*(S4+z*(S5+z*S6)));
slouken@11683
    63
	if(iy==0) return x+v*(S1+z*r);
slouken@11683
    64
	else      return x-((z*(half*y-v*r)-y)-v*S1);
slouken@2756
    65
}