xref: /illumos-gate/usr/src/lib/libm/common/C/atan.c (revision ddc0e0b5)
125c28e83SPiotr Jasiukajtis /*
225c28e83SPiotr Jasiukajtis  * CDDL HEADER START
325c28e83SPiotr Jasiukajtis  *
425c28e83SPiotr Jasiukajtis  * The contents of this file are subject to the terms of the
525c28e83SPiotr Jasiukajtis  * Common Development and Distribution License (the "License").
625c28e83SPiotr Jasiukajtis  * You may not use this file except in compliance with the License.
725c28e83SPiotr Jasiukajtis  *
825c28e83SPiotr Jasiukajtis  * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
925c28e83SPiotr Jasiukajtis  * or http://www.opensolaris.org/os/licensing.
1025c28e83SPiotr Jasiukajtis  * See the License for the specific language governing permissions
1125c28e83SPiotr Jasiukajtis  * and limitations under the License.
1225c28e83SPiotr Jasiukajtis  *
1325c28e83SPiotr Jasiukajtis  * When distributing Covered Code, include this CDDL HEADER in each
1425c28e83SPiotr Jasiukajtis  * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
1525c28e83SPiotr Jasiukajtis  * If applicable, add the following below this CDDL HEADER, with the
1625c28e83SPiotr Jasiukajtis  * fields enclosed by brackets "[]" replaced with your own identifying
1725c28e83SPiotr Jasiukajtis  * information: Portions Copyright [yyyy] [name of copyright owner]
1825c28e83SPiotr Jasiukajtis  *
1925c28e83SPiotr Jasiukajtis  * CDDL HEADER END
2025c28e83SPiotr Jasiukajtis  */
2125c28e83SPiotr Jasiukajtis 
2225c28e83SPiotr Jasiukajtis /*
2325c28e83SPiotr Jasiukajtis  * Copyright 2011 Nexenta Systems, Inc.  All rights reserved.
2425c28e83SPiotr Jasiukajtis  */
2525c28e83SPiotr Jasiukajtis /*
2625c28e83SPiotr Jasiukajtis  * Copyright 2006 Sun Microsystems, Inc.  All rights reserved.
2725c28e83SPiotr Jasiukajtis  * Use is subject to license terms.
2825c28e83SPiotr Jasiukajtis  */
2925c28e83SPiotr Jasiukajtis 
30*ddc0e0b5SRichard Lowe #pragma weak __atan = atan
3125c28e83SPiotr Jasiukajtis 
3225c28e83SPiotr Jasiukajtis /* INDENT OFF */
3325c28e83SPiotr Jasiukajtis /*
3425c28e83SPiotr Jasiukajtis  * atan(x)
3525c28e83SPiotr Jasiukajtis  * Accurate Table look-up algorithm with polynomial approximation in
3625c28e83SPiotr Jasiukajtis  * partially product form.
3725c28e83SPiotr Jasiukajtis  *
3825c28e83SPiotr Jasiukajtis  * -- K.C. Ng, October 17, 2004
3925c28e83SPiotr Jasiukajtis  *
4025c28e83SPiotr Jasiukajtis  * Algorithm
4125c28e83SPiotr Jasiukajtis  *
4225c28e83SPiotr Jasiukajtis  * (1). Purge off Inf and NaN and 0
4325c28e83SPiotr Jasiukajtis  * (2). Reduce x to positive by atan(x) = -atan(-x).
4425c28e83SPiotr Jasiukajtis  * (3). For x <= 1/8 and let z = x*x, return
4525c28e83SPiotr Jasiukajtis  *	(2.1) if x < 2^(-prec/2), atan(x) = x  with inexact flag raised
4625c28e83SPiotr Jasiukajtis  *	(2.2) if x < 2^(-prec/4-1), atan(x) = x+(x/3)(x*x)
4725c28e83SPiotr Jasiukajtis  *	(2.3) if x < 2^(-prec/6-2), atan(x) = x+(z-5/3)(z*x/5)
4825c28e83SPiotr Jasiukajtis  *	(2.4) Otherwise
4925c28e83SPiotr Jasiukajtis  *		atan(x) = poly1(x) = x + A * B,
5025c28e83SPiotr Jasiukajtis  *	where
5125c28e83SPiotr Jasiukajtis  *		A = (p1*x*z) * (p2+z(p3+z))
5225c28e83SPiotr Jasiukajtis  *		B = (p4+z)+z*z) * (p5+z(p6+z))
5325c28e83SPiotr Jasiukajtis  *	Note: (i) domain of poly1 is [0, 1/8], (ii) remez relative
5425c28e83SPiotr Jasiukajtis  *	approximation error of poly1 is bounded by
5525c28e83SPiotr Jasiukajtis  * 		|(atan(x)-poly1(x))/x| <= 2^-57.61
5625c28e83SPiotr Jasiukajtis  * (4). For x >= 8 then
5725c28e83SPiotr Jasiukajtis  *	(3.1) if x >= 2^prec,     atan(x) = atan(inf) - pio2lo
5825c28e83SPiotr Jasiukajtis  *	(3.2) if x >= 2^(prec/3), atan(x) = atan(inf) - 1/x
5925c28e83SPiotr Jasiukajtis  *	(3.3) if x <= 65,	  atan(x) = atan(inf) - poly1(1/x)
6025c28e83SPiotr Jasiukajtis  *	(3.4) otherwise           atan(x) = atan(inf) - poly2(1/x)
6125c28e83SPiotr Jasiukajtis  *	where
6225c28e83SPiotr Jasiukajtis  *		poly2(r) = (q1*r) * (q2+z(q3+z)) * (q4+z),
6325c28e83SPiotr Jasiukajtis  *	its domain is [0, 0.0154]; and its remez absolute
6425c28e83SPiotr Jasiukajtis  *	approximation error is bounded by
6525c28e83SPiotr Jasiukajtis  *		|atan(x)-poly2(x)|<= 2^-59.45
6625c28e83SPiotr Jasiukajtis  *
6725c28e83SPiotr Jasiukajtis  * (5). Now x is in (0.125, 8).
6825c28e83SPiotr Jasiukajtis  *	Recall identity
6925c28e83SPiotr Jasiukajtis  *		atan(x) = atan(y) + atan((x-y)/(1+x*y)).
7025c28e83SPiotr Jasiukajtis  *	Let j = (ix - 0x3fc00000) >> 16, 0 <= j < 96, where ix is the high
7125c28e83SPiotr Jasiukajtis  *	part of x in IEEE double format. Then
7225c28e83SPiotr Jasiukajtis  *		atan(x) = atan(y[j]) + poly2((x-y[j])/(1+x*y[j]))
7325c28e83SPiotr Jasiukajtis  *	where y[j] are carefully chosen so that it matches x to around 4.5
7425c28e83SPiotr Jasiukajtis  *	bits and at the same time atan(y[j]) is very close to an IEEE double
7525c28e83SPiotr Jasiukajtis  *	floating point number. Calculation indicates that
7625c28e83SPiotr Jasiukajtis  *		max|(x-y[j])/(1+x*y[j])| < 0.0154
7725c28e83SPiotr Jasiukajtis  *		j,x
7825c28e83SPiotr Jasiukajtis  *
7925c28e83SPiotr Jasiukajtis  * Accuracy: Maximum error observed is bounded by 0.6 ulp after testing
8025c28e83SPiotr Jasiukajtis  * more than 10 million random arguments
8125c28e83SPiotr Jasiukajtis  */
8225c28e83SPiotr Jasiukajtis /* INDENT ON */
8325c28e83SPiotr Jasiukajtis 
8425c28e83SPiotr Jasiukajtis #include "libm.h"
8525c28e83SPiotr Jasiukajtis #include "libm_protos.h"
8625c28e83SPiotr Jasiukajtis 
8725c28e83SPiotr Jasiukajtis extern const double _TBL_atan[];
8825c28e83SPiotr Jasiukajtis static const double g[] = {
8925c28e83SPiotr Jasiukajtis /* one	= */  1.0,
9025c28e83SPiotr Jasiukajtis /* p1	= */  8.02176624254765935351230154992663301527500152588e-0002,
9125c28e83SPiotr Jasiukajtis /* p2	= */  1.27223421700559402580665846471674740314483642578e+0000,
9225c28e83SPiotr Jasiukajtis /* p3	= */ -1.20606901800503640842521235754247754812240600586e+0000,
9325c28e83SPiotr Jasiukajtis /* p4	= */ -2.36088967922325565496066701598465442657470703125e+0000,
9425c28e83SPiotr Jasiukajtis /* p5	= */  1.38345799501389166152875986881554126739501953125e+0000,
9525c28e83SPiotr Jasiukajtis /* p6	= */  1.06742368078953453469637224770849570631980895996e+0000,
9625c28e83SPiotr Jasiukajtis /* q1   = */ -1.42796626333911796935538518482644576579332351685e-0001,
9725c28e83SPiotr Jasiukajtis /* q2   = */  3.51427110447873227059810477159863497078605962912e+0000,
9825c28e83SPiotr Jasiukajtis /* q3   = */  5.92129112708164262457444237952586263418197631836e-0001,
9925c28e83SPiotr Jasiukajtis /* q4   = */ -1.99272234785683144409063061175402253866195678711e+0000,
10025c28e83SPiotr Jasiukajtis /* pio2hi */  1.570796326794896558e+00,
10125c28e83SPiotr Jasiukajtis /* pio2lo */  6.123233995736765886e-17,
10225c28e83SPiotr Jasiukajtis /* t1   = */ -0.333333333333333333333333333333333,
10325c28e83SPiotr Jasiukajtis /* t2   = */  0.2,
10425c28e83SPiotr Jasiukajtis /* t3   = */ -1.666666666666666666666666666666666,
10525c28e83SPiotr Jasiukajtis };
10625c28e83SPiotr Jasiukajtis 
10725c28e83SPiotr Jasiukajtis #define	one g[0]
10825c28e83SPiotr Jasiukajtis #define	p1 g[1]
10925c28e83SPiotr Jasiukajtis #define	p2 g[2]
11025c28e83SPiotr Jasiukajtis #define	p3 g[3]
11125c28e83SPiotr Jasiukajtis #define	p4 g[4]
11225c28e83SPiotr Jasiukajtis #define	p5 g[5]
11325c28e83SPiotr Jasiukajtis #define	p6 g[6]
11425c28e83SPiotr Jasiukajtis #define	q1 g[7]
11525c28e83SPiotr Jasiukajtis #define	q2 g[8]
11625c28e83SPiotr Jasiukajtis #define	q3 g[9]
11725c28e83SPiotr Jasiukajtis #define	q4 g[10]
11825c28e83SPiotr Jasiukajtis #define	pio2hi g[11]
11925c28e83SPiotr Jasiukajtis #define	pio2lo g[12]
12025c28e83SPiotr Jasiukajtis #define	t1 g[13]
12125c28e83SPiotr Jasiukajtis #define	t2 g[14]
12225c28e83SPiotr Jasiukajtis #define	t3 g[15]
12325c28e83SPiotr Jasiukajtis 
12425c28e83SPiotr Jasiukajtis 
12525c28e83SPiotr Jasiukajtis double
atan(double x)12625c28e83SPiotr Jasiukajtis atan(double x) {
12725c28e83SPiotr Jasiukajtis 	double y, z, r, p, s;
12825c28e83SPiotr Jasiukajtis 	int ix, lx, hx, j;
12925c28e83SPiotr Jasiukajtis 
13025c28e83SPiotr Jasiukajtis 	hx = ((int *) &x)[HIWORD];
13125c28e83SPiotr Jasiukajtis 	lx = ((int *) &x)[LOWORD];
13225c28e83SPiotr Jasiukajtis 	ix = hx & ~0x80000000;
13325c28e83SPiotr Jasiukajtis 	j = ix >> 20;
13425c28e83SPiotr Jasiukajtis 
13525c28e83SPiotr Jasiukajtis 	/* for |x| < 1/8 */
13625c28e83SPiotr Jasiukajtis 	if (j < 0x3fc) {
13725c28e83SPiotr Jasiukajtis 		if (j < 0x3f5) {	/* when |x| < 2**(-prec/6-2) */
13825c28e83SPiotr Jasiukajtis 			if (j < 0x3e3) {	/* if |x| < 2**(-prec/2-2) */
13925c28e83SPiotr Jasiukajtis 				return ((int) x == 0 ? x : one);
14025c28e83SPiotr Jasiukajtis 			}
14125c28e83SPiotr Jasiukajtis 			if (j < 0x3f1) {	/* if |x| < 2**(-prec/4-1) */
14225c28e83SPiotr Jasiukajtis 				return (x + (x * t1) * (x * x));
14325c28e83SPiotr Jasiukajtis 			} else {	/* if |x| < 2**(-prec/6-2) */
14425c28e83SPiotr Jasiukajtis 				z = x * x;
14525c28e83SPiotr Jasiukajtis 				s = t2 * x;
14625c28e83SPiotr Jasiukajtis 				return (x + (t3 + z) * (s * z));
14725c28e83SPiotr Jasiukajtis 			}
14825c28e83SPiotr Jasiukajtis 		}
14925c28e83SPiotr Jasiukajtis 		z = x * x; s = p1 * x;
15025c28e83SPiotr Jasiukajtis 		return (x + ((s * z) * (p2 + z * (p3 + z))) *
15125c28e83SPiotr Jasiukajtis 				(((p4 + z) + z * z) * (p5 + z * (p6 + z))));
15225c28e83SPiotr Jasiukajtis 	}
15325c28e83SPiotr Jasiukajtis 
15425c28e83SPiotr Jasiukajtis 	/* for |x| >= 8.0 */
15525c28e83SPiotr Jasiukajtis 	if (j >= 0x402) {
15625c28e83SPiotr Jasiukajtis 		if (j < 0x436) {
15725c28e83SPiotr Jasiukajtis 			r = one / x;
15825c28e83SPiotr Jasiukajtis 			if (hx >= 0) {
15925c28e83SPiotr Jasiukajtis 				y =  pio2hi; p =  pio2lo;
16025c28e83SPiotr Jasiukajtis 			} else {
16125c28e83SPiotr Jasiukajtis 				y = -pio2hi; p = -pio2lo;
16225c28e83SPiotr Jasiukajtis 			}
16325c28e83SPiotr Jasiukajtis 			if (ix < 0x40504000) {	/* x <  65 */
16425c28e83SPiotr Jasiukajtis 				z = r * r;
16525c28e83SPiotr Jasiukajtis 				s = p1 * r;
16625c28e83SPiotr Jasiukajtis 				return (y + ((p - r) - ((s * z) *
16725c28e83SPiotr Jasiukajtis 					(p2 + z * (p3 + z))) *
16825c28e83SPiotr Jasiukajtis 					(((p4 + z) + z * z) *
16925c28e83SPiotr Jasiukajtis 					(p5 + z * (p6 + z)))));
17025c28e83SPiotr Jasiukajtis 			} else if (j < 0x412) {
17125c28e83SPiotr Jasiukajtis 				z = r * r;
17225c28e83SPiotr Jasiukajtis 				return (y + (p - ((q1 * r) * (q4 + z)) *
17325c28e83SPiotr Jasiukajtis 					(q2 + z * (q3 + z))));
17425c28e83SPiotr Jasiukajtis 			} else
17525c28e83SPiotr Jasiukajtis 				return (y + (p - r));
17625c28e83SPiotr Jasiukajtis 		} else {
17725c28e83SPiotr Jasiukajtis 			if (j >= 0x7ff) /* x is inf or NaN */
17825c28e83SPiotr Jasiukajtis 				if (((ix - 0x7ff00000) | lx) != 0)
17925c28e83SPiotr Jasiukajtis #if defined(FPADD_TRAPS_INCOMPLETE_ON_NAN)
18025c28e83SPiotr Jasiukajtis 					return (ix >= 0x7ff80000 ? x : x - x);
18125c28e83SPiotr Jasiukajtis 					/* assumes sparc-like QNaN */
18225c28e83SPiotr Jasiukajtis #else
18325c28e83SPiotr Jasiukajtis 					return (x - x);
18425c28e83SPiotr Jasiukajtis #endif
18525c28e83SPiotr Jasiukajtis 			y = -pio2lo;
18625c28e83SPiotr Jasiukajtis 			return (hx >= 0 ? pio2hi - y : y - pio2hi);
18725c28e83SPiotr Jasiukajtis 		}
18825c28e83SPiotr Jasiukajtis 	} else {	/* now x is between 1/8 and 8 */
18925c28e83SPiotr Jasiukajtis 		double *w, w0, w1, s, z;
19025c28e83SPiotr Jasiukajtis 		w = (double *) _TBL_atan + (((ix - 0x3fc00000) >> 16) << 1);
19125c28e83SPiotr Jasiukajtis 		w0 = (hx >= 0)? w[0] : -w[0];
19225c28e83SPiotr Jasiukajtis 		s = (x - w0) / (one + x * w0);
19325c28e83SPiotr Jasiukajtis 		w1 = (hx >= 0)? w[1] : -w[1];
19425c28e83SPiotr Jasiukajtis 		z = s * s;
19525c28e83SPiotr Jasiukajtis 		return (((q1 * s) * (q4 + z)) * (q2 + z * (q3 + z)) + w1);
19625c28e83SPiotr Jasiukajtis 	}
19725c28e83SPiotr Jasiukajtis }
198