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