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