1/* origin: FreeBSD /usr/src/lib/msun/src/e_atan2l.c */
2/*
3 * ====================================================
4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
5 *
6 * Developed at SunSoft, a Sun Microsystems, Inc. business.
7 * Permission to use, copy, modify, and distribute this
8 * software is freely granted, provided that this notice
9 * is preserved.
10 * ====================================================
11 *
12 */
13/*
14 * See comments in atan2.c.
15 * Converted to long double by David Schultz <das@FreeBSD.ORG>.
16 */
17
18#include "libm.h"
19
20#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
21long double atan2l(long double y, long double x)
22{
23 return atan2(y, x);
24}
25#elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
26#include "__invtrigl.h"
27
28long double atan2l(long double y, long double x)
29{
30 union ldshape ux, uy;
31 long double z;
32 int m, ex, ey;
33
34 if (isnan(x) || isnan(y))
35 return x+y;
36 if (x == 1)
37 return atanl(y);
38 ux.f = x;
39 uy.f = y;
40 ex = ux.i.se & 0x7fff;
41 ey = uy.i.se & 0x7fff;
42 m = 2*(ux.i.se>>15) | uy.i.se>>15;
43 if (y == 0) {
44 switch(m) {
45 case 0:
46 case 1: return y; /* atan(+-0,+anything)=+-0 */
47 case 2: return 2*pio2_hi; /* atan(+0,-anything) = pi */
48 case 3: return -2*pio2_hi; /* atan(-0,-anything) =-pi */
49 }
50 }
51 if (x == 0)
52 return m&1 ? -pio2_hi : pio2_hi;
53 if (ex == 0x7fff) {
54 if (ey == 0x7fff) {
55 switch(m) {
56 case 0: return pio2_hi/2; /* atan(+INF,+INF) */
57 case 1: return -pio2_hi/2; /* atan(-INF,+INF) */
58 case 2: return 1.5*pio2_hi; /* atan(+INF,-INF) */
59 case 3: return -1.5*pio2_hi; /* atan(-INF,-INF) */
60 }
61 } else {
62 switch(m) {
63 case 0: return 0.0; /* atan(+...,+INF) */
64 case 1: return -0.0; /* atan(-...,+INF) */
65 case 2: return 2*pio2_hi; /* atan(+...,-INF) */
66 case 3: return -2*pio2_hi; /* atan(-...,-INF) */
67 }
68 }
69 }
70 if (ex+120 < ey || ey == 0x7fff)
71 return m&1 ? -pio2_hi : pio2_hi;
72 /* z = atan(|y/x|) without spurious underflow */
73 if ((m&2) && ey+120 < ex) /* |y/x| < 0x1p-120, x<0 */
74 z = 0.0;
75 else
76 z = atanl(fabsl(y/x));
77 switch (m) {
78 case 0: return z; /* atan(+,+) */
79 case 1: return -z; /* atan(-,+) */
80 case 2: return 2*pio2_hi-(z-2*pio2_lo); /* atan(+,-) */
81 default: /* case 3 */
82 return (z-2*pio2_lo)-2*pio2_hi; /* atan(-,-) */
83 }
84}
85#endif