e_atan2l.c revision 181204
136108Sjb 243818Swes/* @(#)e_atan2.c 1.3 95/01/18 */ 336108Sjb/* FreeBSD: head/lib/msun/src/e_atan2.c 176451 2008-02-22 02:30:36Z das */ 436108Sjb/* 543818Swes * ==================================================== 643818Swes * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 743818Swes * 843818Swes * Developed at SunSoft, a Sun Microsystems, Inc. business. 943818Swes * Permission to use, copy, modify, and distribute this 1043818Swes * software is freely granted, provided that this notice 1143818Swes * is preserved. 1243818Swes * ==================================================== 1343818Swes * 1443818Swes */ 1543818Swes 1643818Swes#include <sys/cdefs.h> 1743818Swes__FBSDID("$FreeBSD: head/lib/msun/src/e_atan2l.c 181204 2008-08-02 19:17:00Z das $"); 1843818Swes 1943818Swes/* 2043818Swes * See comments in e_atan2.c. 2143818Swes * Converted to long double by David Schultz <das@FreeBSD.ORG>. 2243818Swes */ 2343818Swes 2443818Swes#include <float.h> 2543818Swes 2643818Swes#include "invtrig.h" 2736108Sjb#include "math.h" 2836108Sjb#include "math_private.h" 2936108Sjb 3036108Sjbstatic volatile long double 3136108Sjbtiny = 1.0e-300; 3236108Sjbstatic const long double 3336108Sjbzero = 0.0; 3436108Sjb 3536108Sjb#ifdef __i386__ 3636108Sjb/* XXX Work around the fact that gcc truncates long double constants on i386 */ 3736108Sjbstatic volatile double 3836108Sjbpi1 = 3.14159265358979311600e+00, /* 0x1.921fb54442d18p+1 */ 3936108Sjbpi2 = 1.22514845490862001043e-16; /* 0x1.1a80000000000p-53 */ 4036108Sjb#define pi ((long double)pi1 + pi2) 4136108Sjb#else 4236108Sjbstatic const long double 4336108Sjbpi = 3.14159265358979323846264338327950280e+00L; 4436108Sjb#endif 4536108Sjb 4636108Sjblong double 4736108Sjbatan2l(long double y, long double x) 4836108Sjb{ 4936108Sjb union IEEEl2bits ux, uy; 5036108Sjb long double z; 5136108Sjb int32_t k,m; 5236108Sjb int16_t exptx, expsignx, expty, expsigny; 5336108Sjb 5436108Sjb uy.e = y; 5536108Sjb expsigny = uy.xbits.expsign; 5636108Sjb expty = expsigny & 0x7fff; 5736108Sjb ux.e = x; 5836108Sjb expsignx = ux.xbits.expsign; 5936108Sjb exptx = expsignx & 0x7fff; 6036108Sjb 6136108Sjb if ((exptx==BIAS+LDBL_MAX_EXP && 6236108Sjb ((ux.bits.manh&~LDBL_NBIT)|ux.bits.manl)!=0) || /* x is NaN */ 6336108Sjb (expty==BIAS+LDBL_MAX_EXP && 6436108Sjb ((uy.bits.manh&~LDBL_NBIT)|uy.bits.manl)!=0)) /* y is NaN */ 6536108Sjb return x+y; 6636108Sjb if (expsignx==BIAS && ((ux.bits.manh&~LDBL_NBIT)|ux.bits.manl)==0) 6736108Sjb return atanl(y); /* x=1.0 */ 6836108Sjb m = ((expsigny>>15)&1)|((expsignx>>14)&2); /* 2*sign(x)+sign(y) */ 6936108Sjb 7036108Sjb /* when y = 0 */ 7136108Sjb if(expty==0 && ((uy.bits.manh&~LDBL_NBIT)|uy.bits.manl)==0) { 7236108Sjb switch(m) { 7336108Sjb case 0: 7436108Sjb case 1: return y; /* atan(+-0,+anything)=+-0 */ 7536108Sjb case 2: return pi+tiny;/* atan(+0,-anything) = pi */ 7636108Sjb case 3: return -pi-tiny;/* atan(-0,-anything) =-pi */ 7736108Sjb } 7836108Sjb } 7936108Sjb /* when x = 0 */ 8036108Sjb if(exptx==0 && ((ux.bits.manh&~LDBL_NBIT)|ux.bits.manl)==0) 8136108Sjb return (expsigny<0)? -pio2_hi-tiny: pio2_hi+tiny; 8236108Sjb 8336108Sjb /* when x is INF */ 8436108Sjb if(exptx==BIAS+LDBL_MAX_EXP) { 8536108Sjb if(expty==BIAS+LDBL_MAX_EXP) { 8636108Sjb switch(m) { 8736108Sjb case 0: return pio2_hi*0.5+tiny;/* atan(+INF,+INF) */ 8836108Sjb case 1: return -pio2_hi*0.5-tiny;/* atan(-INF,+INF) */ 8936108Sjb case 2: return 1.5*pio2_hi+tiny;/*atan(+INF,-INF)*/ 9036108Sjb case 3: return -1.5*pio2_hi-tiny;/*atan(-INF,-INF)*/ 9136108Sjb } 9236108Sjb } else { 9336108Sjb switch(m) { 9436108Sjb case 0: return zero ; /* atan(+...,+INF) */ 9536108Sjb case 1: return -zero ; /* atan(-...,+INF) */ 9636108Sjb case 2: return pi+tiny ; /* atan(+...,-INF) */ 9736108Sjb case 3: return -pi-tiny ; /* atan(-...,-INF) */ 9836108Sjb } 9936108Sjb } 10036108Sjb } 10136108Sjb /* when y is INF */ 10236108Sjb if(expty==BIAS+LDBL_MAX_EXP) 10336108Sjb return (expsigny<0)? -pio2_hi-tiny: pio2_hi+tiny; 10436108Sjb 10536108Sjb /* compute y/x */ 10636108Sjb k = expty-exptx; 10736108Sjb if(k > LDBL_MANT_DIG+2) { /* |y/x| huge */ 10836108Sjb z=pio2_hi+pio2_lo; 10936108Sjb m&=1; 11036108Sjb } 11136108Sjb else if(expsignx<0&&k<-LDBL_MANT_DIG-2) z=0.0; /* |y/x| tiny, x<0 */ 11236108Sjb else z=atanl(fabsl(y/x)); /* safe to do y/x */ 11336108Sjb switch (m) { 11436108Sjb case 0: return z ; /* atan(+,+) */ 11536108Sjb case 1: return -z ; /* atan(-,+) */ 11636108Sjb case 2: return pi-(z-pi_lo);/* atan(+,-) */ 11736108Sjb default: /* case 3 */ 11836108Sjb return (z-pi_lo)-pi;/* atan(-,-) */ 11936108Sjb } 12036108Sjb} 12136108Sjb