1 /* $OpenBSD: e_atan2l.c,v 1.1 2008/12/09 20:00:35 martynas Exp $ */ 2 /* @(#)e_atan2.c 1.3 95/01/18 */ 3 /* FreeBSD: head/lib/msun/src/e_atan2.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 SunSoft, 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 /* 17 * See comments in e_atan2.c. 18 * Converted to long double by David Schultz <das@FreeBSD.ORG>. 19 * Adapted for OpenBSD by Martynas Venckus <martynas@openbsd.org>. 20 */ 21 22 #include <float.h> 23 #include <math.h> 24 25 #include "invtrig.h" 26 #include "math.h" 27 #include "math_private.h" 28 29 #ifdef EXT_IMPLICIT_NBIT 30 #define LDBL_NBIT 0 31 #else /* EXT_IMPLICIT_NBIT */ 32 #define LDBL_NBIT 0x80000000 33 #endif /* EXT_IMPLICIT_NBIT */ 34 35 static volatile long double 36 tiny = 1.0e-300; 37 static const long double 38 zero = 0.0; 39 40 #if defined(__amd64__) || defined(__i386__) 41 /* XXX Work around the fact that gcc truncates long double constants on i386 */ 42 static volatile double 43 pi1 = 3.14159265358979311600e+00, /* 0x1.921fb54442d18p+1 */ 44 pi2 = 1.22514845490862001043e-16; /* 0x1.1a80000000000p-53 */ 45 #define pi ((long double)pi1 + pi2) 46 #else 47 static const long double 48 pi = 3.14159265358979323846264338327950280e+00L; 49 #endif 50 51 long double 52 atan2l(long double y, long double x) 53 { 54 union { 55 long double e; 56 struct ieee_ext bits; 57 } ux, uy; 58 long double z; 59 int32_t k,m; 60 int16_t exptx, expsignx, expty, expsigny; 61 62 uy.e = y; 63 expsigny = (uy.bits.ext_sign << 15) | uy.bits.ext_exp; 64 expty = expsigny & 0x7fff; 65 ux.e = x; 66 expsignx = (ux.bits.ext_sign << 15) | ux.bits.ext_exp; 67 exptx = expsignx & 0x7fff; 68 69 if ((exptx==BIAS+LDBL_MAX_EXP && 70 ((ux.bits.ext_frach&~LDBL_NBIT) 71 #ifdef EXT_FRACHMBITS 72 | ux.bits.ext_frachm 73 #endif /* EXT_FRACHMBITS */ 74 #ifdef EXT_FRACLMBITS 75 | ux.bits.ext_fraclm 76 #endif /* EXT_FRACLMBITS */ 77 | ux.bits.ext_fracl)!=0) || /* x is NaN */ 78 (expty==BIAS+LDBL_MAX_EXP && 79 ((uy.bits.ext_frach&~LDBL_NBIT) 80 #ifdef EXT_FRACHMBITS 81 | uy.bits.ext_frachm 82 #endif /* EXT_FRACHMBITS */ 83 #ifdef EXT_FRACLMBITS 84 | uy.bits.ext_fraclm 85 #endif /* EXT_FRACLMBITS */ 86 | uy.bits.ext_fracl)!=0)) /* y is NaN */ 87 return x+y; 88 if (expsignx==BIAS && ((ux.bits.ext_frach&~LDBL_NBIT) 89 #ifdef EXT_FRACHMBITS 90 | ux.bits.ext_frachm 91 #endif /* EXT_FRACHMBITS */ 92 #ifdef EXT_FRACLMBITS 93 | ux.bits.ext_fraclm 94 #endif /* EXT_FRACLMBITS */ 95 | ux.bits.ext_fracl)==0) 96 return atanl(y); /* x=1.0 */ 97 m = ((expsigny>>15)&1)|((expsignx>>14)&2); /* 2*sign(x)+sign(y) */ 98 99 /* when y = 0 */ 100 if(expty==0 && ((uy.bits.ext_frach&~LDBL_NBIT) 101 #ifdef EXT_FRACHMBITS 102 | uy.bits.ext_frachm 103 #endif /* EXT_FRACHMBITS */ 104 #ifdef EXT_FRACLMBITS 105 | uy.bits.ext_fraclm 106 #endif /* EXT_FRACLMBITS */ 107 | uy.bits.ext_fracl)==0) { 108 switch(m) { 109 case 0: 110 case 1: return y; /* atan(+-0,+anything)=+-0 */ 111 case 2: return pi+tiny;/* atan(+0,-anything) = pi */ 112 case 3: return -pi-tiny;/* atan(-0,-anything) =-pi */ 113 } 114 } 115 /* when x = 0 */ 116 if(exptx==0 && ((ux.bits.ext_frach&~LDBL_NBIT) 117 #ifdef EXT_FRACHMBITS 118 | ux.bits.ext_frachm 119 #endif /* EXT_FRACHMBITS */ 120 #ifdef EXT_FRACLMBITS 121 | ux.bits.ext_fraclm 122 #endif /* EXT_FRACLMBITS */ 123 | ux.bits.ext_fracl)==0) 124 return (expsigny<0)? -pio2_hi-tiny: pio2_hi+tiny; 125 126 /* when x is INF */ 127 if(exptx==BIAS+LDBL_MAX_EXP) { 128 if(expty==BIAS+LDBL_MAX_EXP) { 129 switch(m) { 130 case 0: return pio2_hi*0.5+tiny;/* atan(+INF,+INF) */ 131 case 1: return -pio2_hi*0.5-tiny;/* atan(-INF,+INF) */ 132 case 2: return 1.5*pio2_hi+tiny;/*atan(+INF,-INF)*/ 133 case 3: return -1.5*pio2_hi-tiny;/*atan(-INF,-INF)*/ 134 } 135 } else { 136 switch(m) { 137 case 0: return zero ; /* atan(+...,+INF) */ 138 case 1: return -zero ; /* atan(-...,+INF) */ 139 case 2: return pi+tiny ; /* atan(+...,-INF) */ 140 case 3: return -pi-tiny ; /* atan(-...,-INF) */ 141 } 142 } 143 } 144 /* when y is INF */ 145 if(expty==BIAS+LDBL_MAX_EXP) 146 return (expsigny<0)? -pio2_hi-tiny: pio2_hi+tiny; 147 148 /* compute y/x */ 149 k = expty-exptx; 150 if(k > LDBL_MANT_DIG+2) { /* |y/x| huge */ 151 z=pio2_hi+pio2_lo; 152 m&=1; 153 } 154 else if(expsignx<0&&k<-LDBL_MANT_DIG-2) z=0.0; /* |y/x| tiny, x<0 */ 155 else z=atanl(fabsl(y/x)); /* safe to do y/x */ 156 switch (m) { 157 case 0: return z ; /* atan(+,+) */ 158 case 1: return -z ; /* atan(-,+) */ 159 case 2: return pi-(z-pi_lo);/* atan(+,-) */ 160 default: /* case 3 */ 161 return (z-pi_lo)-pi;/* atan(-,-) */ 162 } 163 } 164