Elliott Hughes | a0ee078 | 2013-01-30 19:06:37 -0800 | [diff] [blame] | 1 | /* @(#)s_atan.c 5.1 93/09/24 */ |
| 2 | /* FreeBSD: head/lib/msun/src/s_atan.c 176451 2008-02-22 02:30:36Z das */ |
| 3 | /* |
| 4 | * ==================================================== |
| 5 | * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. |
| 6 | * |
| 7 | * Developed at SunPro, a Sun Microsystems, Inc. business. |
| 8 | * Permission to use, copy, modify, and distribute this |
| 9 | * software is freely granted, provided that this notice |
| 10 | * is preserved. |
| 11 | * ==================================================== |
| 12 | */ |
| 13 | |
| 14 | #include <sys/cdefs.h> |
| 15 | __FBSDID("$FreeBSD$"); |
| 16 | |
| 17 | /* |
| 18 | * See comments in s_atan.c. |
| 19 | * Converted to long double by David Schultz <das@FreeBSD.ORG>. |
| 20 | */ |
| 21 | |
| 22 | #include <float.h> |
| 23 | |
| 24 | #include "invtrig.h" |
| 25 | #include "math.h" |
| 26 | #include "math_private.h" |
| 27 | |
| 28 | static const long double |
| 29 | one = 1.0, |
| 30 | huge = 1.0e300; |
| 31 | |
| 32 | long double |
| 33 | atanl(long double x) |
| 34 | { |
| 35 | union IEEEl2bits u; |
| 36 | long double w,s1,s2,z; |
| 37 | int id; |
| 38 | int16_t expsign, expt; |
| 39 | int32_t expman; |
| 40 | |
| 41 | u.e = x; |
| 42 | expsign = u.xbits.expsign; |
| 43 | expt = expsign & 0x7fff; |
| 44 | if(expt >= ATAN_CONST) { /* if |x| is large, atan(x)~=pi/2 */ |
| 45 | if(expt == BIAS + LDBL_MAX_EXP && |
| 46 | ((u.bits.manh&~LDBL_NBIT)|u.bits.manl)!=0) |
| 47 | return x+x; /* NaN */ |
| 48 | if(expsign>0) return atanhi[3]+atanlo[3]; |
| 49 | else return -atanhi[3]-atanlo[3]; |
| 50 | } |
| 51 | /* Extract the exponent and the first few bits of the mantissa. */ |
| 52 | /* XXX There should be a more convenient way to do this. */ |
| 53 | expman = (expt << 8) | ((u.bits.manh >> (MANH_SIZE - 9)) & 0xff); |
| 54 | if (expman < ((BIAS - 2) << 8) + 0xc0) { /* |x| < 0.4375 */ |
| 55 | if (expt < ATAN_LINEAR) { /* if |x| is small, atanl(x)~=x */ |
| 56 | if(huge+x>one) return x; /* raise inexact */ |
| 57 | } |
| 58 | id = -1; |
| 59 | } else { |
| 60 | x = fabsl(x); |
| 61 | if (expman < (BIAS << 8) + 0x30) { /* |x| < 1.1875 */ |
| 62 | if (expman < ((BIAS - 1) << 8) + 0x60) { /* 7/16 <=|x|<11/16 */ |
| 63 | id = 0; x = (2.0*x-one)/(2.0+x); |
| 64 | } else { /* 11/16<=|x|< 19/16 */ |
| 65 | id = 1; x = (x-one)/(x+one); |
| 66 | } |
| 67 | } else { |
| 68 | if (expman < ((BIAS + 1) << 8) + 0x38) { /* |x| < 2.4375 */ |
| 69 | id = 2; x = (x-1.5)/(one+1.5*x); |
| 70 | } else { /* 2.4375 <= |x| < 2^ATAN_CONST */ |
| 71 | id = 3; x = -1.0/x; |
| 72 | } |
| 73 | }} |
| 74 | /* end of argument reduction */ |
| 75 | z = x*x; |
| 76 | w = z*z; |
| 77 | /* break sum aT[i]z**(i+1) into odd and even poly */ |
| 78 | s1 = z*T_even(w); |
| 79 | s2 = w*T_odd(w); |
| 80 | if (id<0) return x - x*(s1+s2); |
| 81 | else { |
| 82 | z = atanhi[id] - ((x*(s1+s2) - atanlo[id]) - x); |
| 83 | return (expsign<0)? -z:z; |
| 84 | } |
| 85 | } |