amath  1.8.5
Simple command line calculator
atan2.c
Go to the documentation of this file.
1 /*-
2  * Copyright (c) 2014-2018 Carsten Sonne Larsen <cs@innolan.net>
3  * All rights reserved.
4  *
5  * Redistribution and use in source and binary forms, with or without
6  * modification, are permitted provided that the following conditions
7  * are met:
8  * 1. Redistributions of source code must retain the above copyright
9  * notice, this list of conditions and the following disclaimer.
10  * 2. Redistributions in binary form must reproduce the above copyright
11  * notice, this list of conditions and the following disclaimer in the
12  * documentation and/or other materials provided with the distribution.
13  *
14  * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
15  * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
16  * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
17  * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
18  * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
19  * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
20  * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
21  * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
22  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
23  * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
24  *
25  * Project homepage:
26  * https://amath.innolan.net
27  *
28  * The original source code can be obtained from:
29  * http://www.netlib.org/fdlibm/e_atan2.c
30  *
31  * =================================================================
32  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
33  *
34  * Developed at SunSoft, a Sun Microsystems, Inc. business.
35  * Permission to use, copy, modify, and distribute this
36  * software is freely granted, provided that this notice
37  * is preserved.
38  * =================================================================
39  */
40 
41 /**
42  * @file atan2.c
43  * @brief Inverse tangent function
44  */
45 
46 #include "prim.h"
47 
48 static const double
49  tiny = 1.0e-300,
50  zero = 0.0,
51  pi_o_4 = 7.8539816339744827900E-01, /* 0x3FE921FB, 0x54442D18 */
52  pi_o_2 = 1.5707963267948965580E+00, /* 0x3FF921FB, 0x54442D18 */
53  pi = 3.1415926535897931160E+00, /* 0x400921FB, 0x54442D18 */
54  pi_lo = 1.2246467991473531772E-16; /* 0x3CA1A626, 0x33145C07 */
55 
56 
57 /**
58  * @brief Inverse tangent function
59  * @param y,x
60  * @details
61  * <pre>
62  * Method
63  * 1. Reduce y to positive by atan2(y,x)=-atan2(-y,x).
64  * 2. Reduce x to positive by (if x and y are unexceptional):
65  * ARG (x+iy) = arctan(y/x) ... if x > 0,
66  * ARG (x+iy) = pi - arctan[y/(-x)] ... if x < 0,
67  *
68  * Special cases
69  * ATAN2((anything), NaN ) is NaN;
70  * ATAN2(NAN , (anything) ) is NaN;
71  * ATAN2(+-0, +(anything but NaN)) is +-0 ;
72  * ATAN2(+-0, -(anything but NaN)) is +-pi ;
73  * ATAN2(+-(anything but 0 and NaN), 0) is +-pi/2;
74  * ATAN2(+-(anything but INF and NaN), +INF) is +-0 ;
75  * ATAN2(+-(anything but INF and NaN), -INF) is +-pi;
76  * ATAN2(+-INF,+INF ) is +-pi/4 ;
77  * ATAN2(+-INF,-INF ) is +-3pi/4;
78  * ATAN2(+-INF, (anything but,0,NaN, and INF)) is +-pi/2;
79  *
80  * Constants
81  * The hexadecimal values are the intended ones for the following
82  * constants. The decimal values may be used, provided that the
83  * compiler will convert from decimal to binary accurately enough
84  * to produce the hexadecimal values shown.
85  * </pre>
86  */
87 double atan2(double y, double x)
88 {
89  double z;
90  int32_t k, m, hx, hy, ix, iy;
91  uint32_t lx, ly;
92 
93  EXTRACT_WORDS(hx, lx, x);
94  ix = hx & 0x7fffffff;
95  EXTRACT_WORDS(hy, ly, y);
96  iy = hy & 0x7fffffff;
97  if (((ix | ((lx | -lx) >> 31)) > 0x7ff00000) ||
98  ((iy | ((ly | -ly) >> 31)) > 0x7ff00000)) /* x or y is NaN */
99  return x + y;
100  if (((hx - 0x3ff00000) | lx) == 0)
101  return atan(y); /* x=1.0 */
102  m = ((hy >> 31) & 1) | ((hx >> 30) & 2); /* 2*sign(x)+sign(y) */
103 
104  /* when y = 0 */
105  if ((iy | ly) == 0)
106  {
107  switch (m)
108  {
109  case 0:
110  case 1:
111  return y; /* atan(+-0,+anything)=+-0 */
112  case 2:
113  return pi + tiny; /* atan(+0,-anything) = pi */
114  case 3:
115  return -pi - tiny; /* atan(-0,-anything) =-pi */
116  }
117  }
118  /* when x = 0 */
119  if ((ix | lx) == 0)
120  return (hy < 0) ? -pi_o_2 - tiny : pi_o_2 + tiny;
121 
122  /* when x is INF */
123  if (ix == 0x7ff00000)
124  {
125  if (iy == 0x7ff00000)
126  {
127  switch (m)
128  {
129  case 0:
130  return pi_o_4 + tiny; /* atan(+INF,+INF) */
131  case 1:
132  return -pi_o_4 - tiny; /* atan(-INF,+INF) */
133  case 2:
134  return 3.0 * pi_o_4 + tiny; /*atan(+INF,-INF)*/
135  case 3:
136  return -3.0 * pi_o_4 - tiny; /*atan(-INF,-INF)*/
137  }
138  }
139  else
140  {
141  switch (m)
142  {
143  case 0:
144  return zero; /* atan(+...,+INF) */
145  case 1:
146  return -zero; /* atan(-...,+INF) */
147  case 2:
148  return pi + tiny; /* atan(+...,-INF) */
149  case 3:
150  return -pi - tiny; /* atan(-...,-INF) */
151  }
152  }
153  }
154  /* when y is INF */
155  if (iy == 0x7ff00000)
156  return (hy < 0) ? -pi_o_2 - tiny : pi_o_2 + tiny;
157 
158  /* compute y/x */
159  k = (iy - ix) >> 20;
160  if (k > 60)
161  z = pi_o_2 + 0.5 * pi_lo; /* |y/x| > 2**60 */
162  else if (hx < 0 && k < -60)
163  z = 0.0; /* |y|/x < -2**60 */
164  else
165  z = atan(fabs(y / x)); /* safe to do y/x */
166  switch (m)
167  {
168  case 0:
169  return z; /* atan(+,+) */
170  case 1:
171  {
172  uint32_t zh;
173  GET_HIGH_WORD(zh, z);
174  SET_HIGH_WORD(z, zh ^ 0x80000000);
175  }
176  return z; /* atan(-,+) */
177  case 2:
178  return pi - (z - pi_lo); /* atan(+,-) */
179  default: /* case 3 */
180  return (z - pi_lo) - pi; /* atan(-,-) */
181  }
182 }
static const double pi
Definition: atan2.c:53
#define GET_HIGH_WORD(i, d)
Get the more significant 32 bit int from a double.
Definition: prim.h:167
#define EXTRACT_WORDS(ix0, ix1, d)
Get two 32 bit ints from a double.
Definition: prim.h:156
static const double pi_lo
Definition: atan2.c:54
static const double pi_o_2
Definition: atan2.c:52
double fabs(double x)
Returns the absolute value of x.
Definition: fabs.c:51
static const double tiny
Definition: atan2.c:49
double atan(double x)
Inverse tangent function.
Definition: atan.c:103
double atan2(double y, double x)
Inverse tangent function.
Definition: atan2.c:87
static const double zero
Definition: atan2.c:50
static const double pi_o_4
Definition: atan2.c:51
#define SET_HIGH_WORD(d, v)
Set the more significant 32 bits of a double from an int.
Definition: prim.h:198