1181074Sdas/* @(#)s_atan.c 5.1 93/09/24 */ 2181074Sdas/* FreeBSD: head/lib/msun/src/s_atan.c 176451 2008-02-22 02:30:36Z das */ 3181074Sdas/* 4181074Sdas * ==================================================== 5181074Sdas * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 6181074Sdas * 7181074Sdas * Developed at SunPro, a Sun Microsystems, Inc. business. 8181074Sdas * Permission to use, copy, modify, and distribute this 9181074Sdas * software is freely granted, provided that this notice 10181074Sdas * is preserved. 11181074Sdas * ==================================================== 12181074Sdas */ 13181074Sdas 14181074Sdas#include <sys/cdefs.h> 15181074Sdas__FBSDID("$FreeBSD$"); 16181074Sdas 17181074Sdas/* 18181074Sdas * See comments in s_atan.c. 19181074Sdas * Converted to long double by David Schultz <das@FreeBSD.ORG>. 20181074Sdas */ 21181074Sdas 22181074Sdas#include <float.h> 23181074Sdas 24181074Sdas#include "invtrig.h" 25181074Sdas#include "math.h" 26181074Sdas#include "math_private.h" 27181074Sdas 28181074Sdasstatic const long double 29181074Sdasone = 1.0, 30181074Sdashuge = 1.0e300; 31181074Sdas 32181074Sdaslong double 33181074Sdasatanl(long double x) 34181074Sdas{ 35181074Sdas union IEEEl2bits u; 36181074Sdas long double w,s1,s2,z; 37181074Sdas int id; 38181074Sdas int16_t expsign, expt; 39181074Sdas int32_t expman; 40181074Sdas 41181074Sdas u.e = x; 42181074Sdas expsign = u.xbits.expsign; 43181074Sdas expt = expsign & 0x7fff; 44181074Sdas if(expt >= ATAN_CONST) { /* if |x| is large, atan(x)~=pi/2 */ 45181074Sdas if(expt == BIAS + LDBL_MAX_EXP && 46181074Sdas ((u.bits.manh&~LDBL_NBIT)|u.bits.manl)!=0) 47181074Sdas return x+x; /* NaN */ 48181074Sdas if(expsign>0) return atanhi[3]+atanlo[3]; 49181074Sdas else return -atanhi[3]-atanlo[3]; 50181074Sdas } 51181074Sdas /* Extract the exponent and the first few bits of the mantissa. */ 52181074Sdas /* XXX There should be a more convenient way to do this. */ 53181074Sdas expman = (expt << 8) | ((u.bits.manh >> (MANH_SIZE - 9)) & 0xff); 54181074Sdas if (expman < ((BIAS - 2) << 8) + 0xc0) { /* |x| < 0.4375 */ 55181074Sdas if (expt < ATAN_LINEAR) { /* if |x| is small, atanl(x)~=x */ 56181074Sdas if(huge+x>one) return x; /* raise inexact */ 57181074Sdas } 58181074Sdas id = -1; 59181074Sdas } else { 60181074Sdas x = fabsl(x); 61181074Sdas if (expman < (BIAS << 8) + 0x30) { /* |x| < 1.1875 */ 62181074Sdas if (expman < ((BIAS - 1) << 8) + 0x60) { /* 7/16 <=|x|<11/16 */ 63181074Sdas id = 0; x = (2.0*x-one)/(2.0+x); 64181074Sdas } else { /* 11/16<=|x|< 19/16 */ 65181074Sdas id = 1; x = (x-one)/(x+one); 66181074Sdas } 67181074Sdas } else { 68181074Sdas if (expman < ((BIAS + 1) << 8) + 0x38) { /* |x| < 2.4375 */ 69181074Sdas id = 2; x = (x-1.5)/(one+1.5*x); 70181074Sdas } else { /* 2.4375 <= |x| < 2^ATAN_CONST */ 71181074Sdas id = 3; x = -1.0/x; 72181074Sdas } 73181074Sdas }} 74181074Sdas /* end of argument reduction */ 75181074Sdas z = x*x; 76181074Sdas w = z*z; 77181074Sdas /* break sum aT[i]z**(i+1) into odd and even poly */ 78181074Sdas s1 = z*T_even(w); 79181074Sdas s2 = w*T_odd(w); 80181074Sdas if (id<0) return x - x*(s1+s2); 81181074Sdas else { 82181074Sdas z = atanhi[id] - ((x*(s1+s2) - atanlo[id]) - x); 83181074Sdas return (expsign<0)? -z:z; 84181074Sdas } 85181074Sdas} 86