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