k_tan.c revision 97413
118334Speter/* @(#)k_tan.c 5.1 93/09/24 */ 218334Speter/* 318334Speter * ==================================================== 490075Sobrien * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 590075Sobrien * 618334Speter * Developed at SunPro, a Sun Microsystems, Inc. business. 718334Speter * Permission to use, copy, modify, and distribute this 890075Sobrien * software is freely granted, provided that this notice 918334Speter * is preserved. 1090075Sobrien * ==================================================== 1190075Sobrien */ 1290075Sobrien 1390075Sobrien#ifndef lint 1418334Speterstatic char rcsid[] = "$FreeBSD: head/lib/msun/src/k_tan.c 97413 2002-05-28 18:15:04Z alfred $"; 1590075Sobrien#endif 1690075Sobrien 1790075Sobrien/* __kernel_tan( x, y, k ) 1890075Sobrien * kernel tan function on [-pi/4, pi/4], pi/4 ~ 0.7854 1918334Speter * Input x is assumed to be bounded by ~pi/4 in magnitude. 2018334Speter * Input y is the tail of x. 2190075Sobrien * Input k indicates whether tan (if k=1) or 2290075Sobrien * -1/tan (if k= -1) is returned. 2390075Sobrien * 2418334Speter * Algorithm 2518334Speter * 1. Since tan(-x) = -tan(x), we need only to consider positive x. 2650397Sobrien * 2. if x < 2^-28 (hx<0x3e300000 0), return x with inexact if x!=0. 2752284Sobrien * 3. tan(x) is approximated by a odd polynomial of degree 27 on 2818334Speter * [0,0.67434] 2918334Speter * 3 27 3018334Speter * tan(x) ~ x + T1*x + ... + T13*x 3118334Speter * where 3218334Speter * 3318334Speter * |tan(x) 2 4 26 | -59.2 3418334Speter * |----- - (1+T1*x +T2*x +.... +T13*x )| <= 2 3518334Speter * | x | 3618334Speter * 3718334Speter * Note: tan(x+y) = tan(x) + tan'(x)*y 3818334Speter * ~ tan(x) + (1+x*x)*y 3918334Speter * Therefore, for better accuracy in computing tan(x+y), let 4052284Sobrien * 3 2 2 2 2 4118334Speter * r = x *(T2+x *(T3+x *(...+x *(T12+x *T13)))) 4290075Sobrien * then 4390075Sobrien * 3 2 4490075Sobrien * tan(x+y) = x + (T1*x + (x *(r+y)+y)) 4590075Sobrien * 4690075Sobrien * 4. For x in [0.67434,pi/4], let y = pi/4 - x, then 4790075Sobrien * tan(x) = tan(pi/4-y) = (1-tan(y))/(1+tan(y)) 4818334Speter * = 1 - 2*(tan(y) - (tan(y)^2)/(1+tan(y))) 4918334Speter */ 5018334Speter 5118334Speter#include "math.h" 5218334Speter#include "math_private.h" 5318334Speterstatic const double 5418334Speterone = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */ 5518334Speterpio4 = 7.85398163397448278999e-01, /* 0x3FE921FB, 0x54442D18 */ 5618334Speterpio4lo= 3.06161699786838301793e-17, /* 0x3C81A626, 0x33145C07 */ 5718334SpeterT[] = { 5818334Speter 3.33333333333334091986e-01, /* 0x3FD55555, 0x55555563 */ 5918334Speter 1.33333333333201242699e-01, /* 0x3FC11111, 0x1110FE7A */ 6018334Speter 5.39682539762260521377e-02, /* 0x3FABA1BA, 0x1BB341FE */ 6118334Speter 2.18694882948595424599e-02, /* 0x3F9664F4, 0x8406D637 */ 6218334Speter 8.86323982359930005737e-03, /* 0x3F8226E3, 0xE96E8493 */ 6350397Sobrien 3.59207910759131235356e-03, /* 0x3F6D6D22, 0xC9560328 */ 6452284Sobrien 1.45620945432529025516e-03, /* 0x3F57DBC8, 0xFEE08315 */ 6552284Sobrien 5.88041240820264096874e-04, /* 0x3F4344D8, 0xF2F26501 */ 6618334Speter 2.46463134818469906812e-04, /* 0x3F3026F7, 0x1A8D1068 */ 6790075Sobrien 7.81794442939557092300e-05, /* 0x3F147E88, 0xA03792A6 */ 6818334Speter 7.14072491382608190305e-05, /* 0x3F12B80F, 0x32F0A7E9 */ 6918334Speter -1.85586374855275456654e-05, /* 0xBEF375CB, 0xDB605373 */ 7018334Speter 2.59073051863633712884e-05, /* 0x3EFB2A70, 0x74BF7AD4 */ 7118334Speter}; 7218334Speter 7318334Speterdouble 7418334Speter__kernel_tan(double x, double y, int iy) 7518334Speter{ 7618334Speter double z,r,v,w,s; 7718334Speter int32_t ix,hx; 7818334Speter GET_HIGH_WORD(hx,x); 7918334Speter ix = hx&0x7fffffff; /* high word of |x| */ 8018334Speter if(ix<0x3e300000) /* x < 2**-28 */ 8118334Speter {if((int)x==0) { /* generate inexact */ 8218334Speter u_int32_t low; 8318334Speter GET_LOW_WORD(low,x); 8418334Speter if(((ix|low)|(iy+1))==0) return one/fabs(x); 8518334Speter else return (iy==1)? x: -one/x; 8618334Speter } 8718334Speter } 8818334Speter if(ix>=0x3FE59428) { /* |x|>=0.6744 */ 8918334Speter if(hx<0) {x = -x; y = -y;} 9018334Speter z = pio4-x; 9118334Speter w = pio4lo-y; 9218334Speter x = z+w; y = 0.0; 9390075Sobrien } 9418334Speter z = x*x; 9518334Speter w = z*z; 9618334Speter /* Break x^5*(T[1]+x^2*T[2]+...) into 9790075Sobrien * x^5(T[1]+x^4*T[3]+...+x^20*T[11]) + 9818334Speter * x^5(x^2*(T[2]+x^4*T[4]+...+x^22*[T12])) 9990075Sobrien */ 10090075Sobrien r = T[1]+w*(T[3]+w*(T[5]+w*(T[7]+w*(T[9]+w*T[11])))); 10118334Speter v = z*(T[2]+w*(T[4]+w*(T[6]+w*(T[8]+w*(T[10]+w*T[12]))))); 10218334Speter s = z*x; 10318334Speter r = y + z*(s*(r+v)+y); 10418334Speter r += T[0]*s; 10518334Speter w = x+r; 10618334Speter if(ix>=0x3FE59428) { 10718334Speter v = (double)iy; 10850397Sobrien return (double)(1-((hx>>30)&2))*(v-2.0*(x-(w*w/(w+v)-r))); 10918334Speter } 11052284Sobrien if(iy==1) return w; 11118334Speter else { /* if allow error up to 2 ulp, 11218334Speter simply return -1.0/(x+r) here */ 11318334Speter /* compute -1.0/(x+r) accurately */ 11418334Speter double a,t; 11552284Sobrien z = w; 11618334Speter SET_LOW_WORD(z,0); 11718334Speter v = r-(z - x); /* z+v = r+x */ 11818334Speter t = a = -1.0/w; /* a = -1.0/w */ 11918334Speter SET_LOW_WORD(t,0); 12018334Speter s = 1.0+t*z; 12118334Speter return t+a*(s+t*v); 12218334Speter } 12318334Speter} 12452284Sobrien