src/libm/k_cos.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
/*
slouken@2756
    13
 * __kernel_cos( x,  y )
slouken@2756
    14
 * kernel cos function on [-pi/4, pi/4], pi/4 ~ 0.785398164
slouken@2756
    15
 * Input x is assumed to be bounded by ~pi/4 in magnitude.
slouken@2756
    16
 * Input y is the tail of x.
slouken@2756
    17
 *
slouken@2756
    18
 * Algorithm
slouken@2756
    19
 *	1. Since cos(-x) = cos(x), we need only to consider positive x.
slouken@2756
    20
 *	2. if x < 2^-27 (hx<0x3e400000 0), return 1 with inexact if x!=0.
slouken@2756
    21
 *	3. cos(x) is approximated by a polynomial of degree 14 on
slouken@2756
    22
 *	   [0,pi/4]
slouken@2756
    23
 *		  	                 4            14
slouken@2756
    24
 *	   	cos(x) ~ 1 - x*x/2 + C1*x + ... + C6*x
slouken@2756
    25
 *	   where the remez error is
slouken@2756
    26
 *
slouken@2756
    27
 * 	|              2     4     6     8     10    12     14 |     -58
slouken@2756
    28
 * 	|cos(x)-(1-.5*x +C1*x +C2*x +C3*x +C4*x +C5*x  +C6*x  )| <= 2
slouken@2756
    29
 * 	|    					               |
slouken@2756
    30
 *
slouken@2756
    31
 * 	               4     6     8     10    12     14
slouken@2756
    32
 *	4. let r = C1*x +C2*x +C3*x +C4*x +C5*x  +C6*x  , then
slouken@2756
    33
 *	       cos(x) = 1 - x*x/2 + r
slouken@2756
    34
 *	   since cos(x+y) ~ cos(x) - sin(x)*y
slouken@2756
    35
 *			  ~ cos(x) - x*y,
slouken@2756
    36
 *	   a correction term is necessary in cos(x) and hence
slouken@2756
    37
 *		cos(x+y) = 1 - (x*x/2 - (r - x*y))
slouken@2756
    38
 *	   For better accuracy when x > 0.3, let qx = |x|/4 with
slouken@2756
    39
 *	   the last 32 bits mask off, and if x > 0.78125, let qx = 0.28125.
slouken@2756
    40
 *	   Then
slouken@2756
    41
 *		cos(x+y) = (1-qx) - ((x*x/2-qx) - (r-x*y)).
slouken@2756
    42
 *	   Note that 1-qx and (x*x/2-qx) is EXACT here, and the
slouken@2756
    43
 *	   magnitude of the latter is at least a quarter of x*x/2,
slouken@2756
    44
 *	   thus, reducing the rounding error in the subtraction.
slouken@2756
    45
 */
slouken@2756
    46
slouken@6044
    47
#include "math_libm.h"
slouken@2756
    48
#include "math_private.h"
slouken@2756
    49
slouken@2756
    50
static const double
slouken@11683
    51
one =  1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
slouken@11683
    52
C1  =  4.16666666666666019037e-02, /* 0x3FA55555, 0x5555554C */
slouken@11683
    53
C2  = -1.38888888888741095749e-03, /* 0xBF56C16C, 0x16C15177 */
slouken@11683
    54
C3  =  2.48015872894767294178e-05, /* 0x3EFA01A0, 0x19CB1590 */
slouken@11683
    55
C4  = -2.75573143513906633035e-07, /* 0xBE927E4F, 0x809C52AD */
slouken@11683
    56
C5  =  2.08757232129817482790e-09, /* 0x3E21EE9E, 0xBDB4B1C4 */
slouken@11683
    57
C6  = -1.13596475577881948265e-11; /* 0xBDA8FAE9, 0xBE8838D4 */
slouken@2756
    58
slouken@11683
    59
double attribute_hidden __kernel_cos(double x, double y)
slouken@2756
    60
{
slouken@11683
    61
	double a,hz,z,r,qx;
slouken@11683
    62
	int32_t ix;
slouken@11683
    63
	GET_HIGH_WORD(ix,x);
slouken@11683
    64
	ix &= 0x7fffffff;			/* ix = |x|'s high word*/
slouken@11683
    65
	if(ix<0x3e400000) {			/* if x < 2**27 */
slouken@11683
    66
	    if(((int)x)==0) return one;		/* generate inexact */
slouken@11683
    67
	}
slouken@11683
    68
	z  = x*x;
slouken@11683
    69
	r  = z*(C1+z*(C2+z*(C3+z*(C4+z*(C5+z*C6)))));
slouken@11683
    70
	if(ix < 0x3FD33333) 			/* if |x| < 0.3 */
slouken@11683
    71
	    return one - (0.5*z - (z*r - x*y));
slouken@11683
    72
	else {
slouken@11683
    73
	    if(ix > 0x3fe90000) {		/* x > 0.78125 */
slouken@11683
    74
		qx = 0.28125;
slouken@11683
    75
	    } else {
slouken@11683
    76
	        INSERT_WORDS(qx,ix-0x00200000,0);	/* x/4 */
slouken@11683
    77
	    }
slouken@11683
    78
	    hz = 0.5*z-qx;
slouken@11683
    79
	    a  = one-qx;
slouken@11683
    80
	    return a - (hz - (z*r-x*y));
slouken@11683
    81
	}
slouken@2756
    82
}