src/libm/k_tan.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_tan( x, y, k )
    13  * kernel tan 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 k indicates whether tan (if k=1) or
    17  * -1/tan (if k= -1) is returned.
    18  *
    19  * Algorithm
    20  *	1. Since tan(-x) = -tan(x), we need only to consider positive x.
    21  *	2. if x < 2^-28 (hx<0x3e300000 0), return x with inexact if x!=0.
    22  *	3. tan(x) is approximated by a odd polynomial of degree 27 on
    23  *	   [0,0.67434]
    24  *		  	         3             27
    25  *	   	tan(x) ~ x + T1*x + ... + T13*x
    26  *	   where
    27  *
    28  * 	        |tan(x)         2     4            26   |     -59.2
    29  * 	        |----- - (1+T1*x +T2*x +.... +T13*x    )| <= 2
    30  * 	        |  x 					|
    31  *
    32  *	   Note: tan(x+y) = tan(x) + tan'(x)*y
    33  *		          ~ tan(x) + (1+x*x)*y
    34  *	   Therefore, for better accuracy in computing tan(x+y), let
    35  *		     3      2      2       2       2
    36  *		r = x *(T2+x *(T3+x *(...+x *(T12+x *T13))))
    37  *	   then
    38  *		 		    3    2
    39  *		tan(x+y) = x + (T1*x + (x *(r+y)+y))
    40  *
    41  *      4. For x in [0.67434,pi/4],  let y = pi/4 - x, then
    42  *		tan(x) = tan(pi/4-y) = (1-tan(y))/(1+tan(y))
    43  *		       = 1 - 2*(tan(y) - (tan(y)^2)/(1+tan(y)))
    44  */
    45 
    46 #include "math_libm.h"
    47 #include "math_private.h"
    48 
    49 static const double
    50 one   =  1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
    51 pio4  =  7.85398163397448278999e-01, /* 0x3FE921FB, 0x54442D18 */
    52 pio4lo=  3.06161699786838301793e-17, /* 0x3C81A626, 0x33145C07 */
    53 T[] =  {
    54   3.33333333333334091986e-01, /* 0x3FD55555, 0x55555563 */
    55   1.33333333333201242699e-01, /* 0x3FC11111, 0x1110FE7A */
    56   5.39682539762260521377e-02, /* 0x3FABA1BA, 0x1BB341FE */
    57   2.18694882948595424599e-02, /* 0x3F9664F4, 0x8406D637 */
    58   8.86323982359930005737e-03, /* 0x3F8226E3, 0xE96E8493 */
    59   3.59207910759131235356e-03, /* 0x3F6D6D22, 0xC9560328 */
    60   1.45620945432529025516e-03, /* 0x3F57DBC8, 0xFEE08315 */
    61   5.88041240820264096874e-04, /* 0x3F4344D8, 0xF2F26501 */
    62   2.46463134818469906812e-04, /* 0x3F3026F7, 0x1A8D1068 */
    63   7.81794442939557092300e-05, /* 0x3F147E88, 0xA03792A6 */
    64   7.14072491382608190305e-05, /* 0x3F12B80F, 0x32F0A7E9 */
    65  -1.85586374855275456654e-05, /* 0xBEF375CB, 0xDB605373 */
    66   2.59073051863633712884e-05, /* 0x3EFB2A70, 0x74BF7AD4 */
    67 };
    68 
    69 double attribute_hidden __kernel_tan(double x, double y, int iy)
    70 {
    71 	double z,r,v,w,s;
    72 	int32_t ix,hx;
    73 	GET_HIGH_WORD(hx,x);
    74 	ix = hx&0x7fffffff;	/* high word of |x| */
    75 	if(ix<0x3e300000)			/* x < 2**-28 */
    76 	    {if((int)x==0) {			/* generate inexact */
    77 	        u_int32_t low;
    78 		GET_LOW_WORD(low,x);
    79 		if(((ix|low)|(iy+1))==0) return one/fabs(x);
    80 		else return (iy==1)? x: -one/x;
    81 	    }
    82 	    }
    83 	if(ix>=0x3FE59428) { 			/* |x|>=0.6744 */
    84 	    if(hx<0) {x = -x; y = -y;}
    85 	    z = pio4-x;
    86 	    w = pio4lo-y;
    87 	    x = z+w; y = 0.0;
    88 	}
    89 	z	=  x*x;
    90 	w 	=  z*z;
    91     /* Break x^5*(T[1]+x^2*T[2]+...) into
    92      *	  x^5(T[1]+x^4*T[3]+...+x^20*T[11]) +
    93      *	  x^5(x^2*(T[2]+x^4*T[4]+...+x^22*[T12]))
    94      */
    95 	r = T[1]+w*(T[3]+w*(T[5]+w*(T[7]+w*(T[9]+w*T[11]))));
    96 	v = z*(T[2]+w*(T[4]+w*(T[6]+w*(T[8]+w*(T[10]+w*T[12])))));
    97 	s = z*x;
    98 	r = y + z*(s*(r+v)+y);
    99 	r += T[0]*s;
   100 	w = x+r;
   101 	if(ix>=0x3FE59428) {
   102 	    v = (double)iy;
   103 	    return (double)(1-((hx>>30)&2))*(v-2.0*(x-(w*w/(w+v)-r)));
   104 	}
   105 	if(iy==1) return w;
   106 	else {		/* if allow error up to 2 ulp,
   107 			   simply return -1.0/(x+r) here */
   108      /*  compute -1.0/(x+r) accurately */
   109 	    double a,t;
   110 	    z  = w;
   111 	    SET_LOW_WORD(z,0);
   112 	    v  = r-(z - x); 	/* z+v = r+x */
   113 	    t = a  = -1.0/w;	/* a = -1.0/w */
   114 	    SET_LOW_WORD(t,0);
   115 	    s  = 1.0+t*z;
   116 	    return t+a*(s+t*v);
   117 	}
   118 }