Back to index

glibc  2.9
e_atan2l.c
Go to the documentation of this file.
00001 /* e_atan2l.c -- long double version of e_atan2.c.
00002  * Conversion to long double by Ulrich Drepper,
00003  * Cygnus Support, drepper@cygnus.com.
00004  */
00005 
00006 /*
00007  * ====================================================
00008  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
00009  *
00010  * Developed at SunPro, a Sun Microsystems, Inc. business.
00011  * Permission to use, copy, modify, and distribute this
00012  * software is freely granted, provided that this notice
00013  * is preserved.
00014  * ====================================================
00015  */
00016 
00017 #if defined(LIBM_SCCS) && !defined(lint)
00018 static char rcsid[] = "$NetBSD: $";
00019 #endif
00020 
00021 /* __ieee754_atan2l(y,x)
00022  * Method :
00023  *     1. Reduce y to positive by atan2(y,x)=-atan2(-y,x).
00024  *     2. Reduce x to positive by (if x and y are unexceptional):
00025  *            ARG (x+iy) = arctan(y/x)       ... if x > 0,
00026  *            ARG (x+iy) = pi - arctan[y/(-x)]   ... if x < 0,
00027  *
00028  * Special cases:
00029  *
00030  *     ATAN2((anything), NaN ) is NaN;
00031  *     ATAN2(NAN , (anything) ) is NaN;
00032  *     ATAN2(+-0, +(anything but NaN)) is +-0  ;
00033  *     ATAN2(+-0, -(anything but NaN)) is +-pi ;
00034  *     ATAN2(+-(anything but 0 and NaN), 0) is +-pi/2;
00035  *     ATAN2(+-(anything but INF and NaN), +INF) is +-0 ;
00036  *     ATAN2(+-(anything but INF and NaN), -INF) is +-pi;
00037  *     ATAN2(+-INF,+INF ) is +-pi/4 ;
00038  *     ATAN2(+-INF,-INF ) is +-3pi/4;
00039  *     ATAN2(+-INF, (anything but,0,NaN, and INF)) is +-pi/2;
00040  *
00041  * Constants:
00042  * The hexadecimal values are the intended ones for the following
00043  * constants. The decimal values may be used, provided that the
00044  * compiler will convert from decimal to binary accurately enough
00045  * to produce the hexadecimal values shown.
00046  */
00047 
00048 #include "math.h"
00049 #include "math_private.h"
00050 
00051 #ifdef __STDC__
00052 static const long double
00053 #else
00054 static long double
00055 #endif
00056 tiny  = 1.0e-4900L,
00057 zero  = 0.0,
00058 pi_o_4  = 7.85398163397448309628202E-01L, /* 0x3FFE, 0xC90FDAA2, 0x2168C235 */
00059 pi_o_2  = 1.5707963267948966192564E+00L,  /* 0x3FFF, 0xC90FDAA2, 0x2168C235 */
00060 pi      = 3.14159265358979323851281E+00L, /* 0x4000, 0xC90FDAA2, 0x2168C235 */
00061 pi_lo   = -5.01655761266833202345176e-20L;/* 0xBFBE, 0xECE675D1, 0xFC8F8CBB */
00062 
00063 #ifdef __STDC__
00064        long double __ieee754_atan2l(long double y, long double x)
00065 #else
00066        long double __ieee754_atan2l(y,x)
00067        long double  y,x;
00068 #endif
00069 {
00070        long double z;
00071        int32_t k,m,hx,hy,ix,iy;
00072        u_int32_t sx,sy,lx,ly;
00073 
00074        GET_LDOUBLE_WORDS(sx,hx,lx,x);
00075        ix = sx&0x7fff;
00076        lx |= hx & 0x7fffffff;
00077        GET_LDOUBLE_WORDS(sy,hy,ly,y);
00078        iy = sy&0x7fff;
00079        ly |= hy & 0x7fffffff;
00080        if(((2*ix|((lx|-lx)>>31))>0xfffe)||
00081           ((2*iy|((ly|-ly)>>31))>0xfffe)) /* x or y is NaN */
00082           return x+y;
00083        if(((sx-0x3fff)|lx)==0) return __atanl(y);   /* x=1.0 */
00084        m = ((sy>>15)&1)|((sx>>14)&2);     /* 2*sign(x)+sign(y) */
00085 
00086     /* when y = 0 */
00087        if((iy|ly)==0) {
00088            switch(m) {
00089               case 0:
00090               case 1: return y;    /* atan(+-0,+anything)=+-0 */
00091               case 2: return  pi+tiny;/* atan(+0,-anything) = pi */
00092               case 3: return -pi-tiny;/* atan(-0,-anything) =-pi */
00093            }
00094        }
00095     /* when x = 0 */
00096        if((ix|lx)==0) return (sy>=0x8000)?  -pi_o_2-tiny: pi_o_2+tiny;
00097 
00098     /* when x is INF */
00099        if(ix==0x7fff) {
00100            if(iy==0x7fff) {
00101               switch(m) {
00102                   case 0: return  pi_o_4+tiny;/* atan(+INF,+INF) */
00103                   case 1: return -pi_o_4-tiny;/* atan(-INF,+INF) */
00104                   case 2: return  3.0*pi_o_4+tiny;/*atan(+INF,-INF)*/
00105                   case 3: return -3.0*pi_o_4-tiny;/*atan(-INF,-INF)*/
00106               }
00107            } else {
00108               switch(m) {
00109                   case 0: return  zero  ; /* atan(+...,+INF) */
00110                   case 1: return -zero  ; /* atan(-...,+INF) */
00111                   case 2: return  pi+tiny  ;     /* atan(+...,-INF) */
00112                   case 3: return -pi-tiny  ;     /* atan(-...,-INF) */
00113               }
00114            }
00115        }
00116     /* when y is INF */
00117        if(iy==0x7fff) return (sy>=0x8000)? -pi_o_2-tiny: pi_o_2+tiny;
00118 
00119     /* compute y/x */
00120        k = sy-sx;
00121        if(k > 70) z=pi_o_2+0.5*pi_lo;     /* |y/x| >  2**70 */
00122        else if(sx>=0x8000&&k<-70) z=0.0;  /* |y|/x < -2**70 */
00123        else z=__atanl(fabsl(y/x)); /* safe to do y/x */
00124        switch (m) {
00125            case 0: return       z  ;      /* atan(+,+) */
00126            case 1: {
00127                     u_int32_t sz;
00128                     GET_LDOUBLE_EXP(sz,z);
00129                     SET_LDOUBLE_EXP(z,sz ^ 0x8000);
00130                   }
00131                   return       z  ;       /* atan(-,+) */
00132            case 2: return  pi-(z-pi_lo);/* atan(+,-) */
00133            default: /* case 3 */
00134                   return  (z-pi_lo)-pi;/* atan(-,-) */
00135        }
00136 }