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