1 /* 2 * CDDL HEADER START 3 * 4 * The contents of this file are subject to the terms of the 5 * Common Development and Distribution License (the "License"). 6 * You may not use this file except in compliance with the License. 7 * 8 * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE 9 * or http://www.opensolaris.org/os/licensing. 10 * See the License for the specific language governing permissions 11 * and limitations under the License. 12 * 13 * When distributing Covered Code, include this CDDL HEADER in each 14 * file and include the License file at usr/src/OPENSOLARIS.LICENSE. 15 * If applicable, add the following below this CDDL HEADER, with the 16 * fields enclosed by brackets "[]" replaced with your own identifying 17 * information: Portions Copyright [yyyy] [name of copyright owner] 18 * 19 * CDDL HEADER END 20 */ 21 22 /* 23 * Copyright 2011 Nexenta Systems, Inc. All rights reserved. 24 */ 25 /* 26 * Copyright 2006 Sun Microsystems, Inc. All rights reserved. 27 * Use is subject to license terms. 28 */ 29 30 /* 31 * atan2l(y,x) 32 * 33 * Method : 34 * 1. Reduce y to positive by atan2(y,x)=-atan2(-y,x). 35 * 2. Reduce x to positive by (if x and y are unexceptional): 36 * ARG (x+iy) = arctan(y/x) ... if x > 0, 37 * ARG (x+iy) = pi - arctan[y/(-x)] ... if x < 0, 38 * 39 * Special cases: 40 * 41 * ATAN2((anything), NaN ) is NaN; 42 * ATAN2(NAN , (anything) ) is NaN; 43 * ATAN2(+-0, +(anything but NaN)) is +-0 ; 44 * ATAN2(+-0, -(anything but NaN)) is +-PI ; 45 * ATAN2(+-(anything but 0 and NaN), 0) is +-PI/2; 46 * ATAN2(+-(anything but INF and NaN), +INF) is +-0 ; 47 * ATAN2(+-(anything but INF and NaN), -INF) is +-PI; 48 * ATAN2(+-INF,+INF ) is +-PI/4 ; 49 * ATAN2(+-INF,-INF ) is +-3PI/4; 50 * ATAN2(+-INF, (anything but,0,NaN, and INF)) is +-PI/2; 51 * 52 * Constants: 53 * The hexadecimal values are the intended ones for the following constants. 54 * The decimal values may be used, provided that the compiler will convert 55 * from decimal to binary accurately enough to produce the hexadecimal values 56 * shown. 57 */ 58 59 #pragma weak __atan2l = atan2l 60 61 #include "libm.h" 62 #include "longdouble.h" 63 64 static const long double 65 zero = 0.0L, 66 tiny = 1.0e-40L, 67 one = 1.0L, 68 half = 0.5L, 69 PI3o4 = 2.356194490192344928846982537459627163148L, 70 PIo4 = 0.785398163397448309615660845819875721049L, 71 PIo2 = 1.570796326794896619231321691639751442099L, 72 PI = 3.141592653589793238462643383279502884197L, 73 PI_lo = 8.671810130123781024797044026043351968762e-35L; 74 75 long double 76 atan2l(long double y, long double x) 77 { 78 long double t, z; 79 int k, m, signy, signx; 80 81 if (x != x || y != y) 82 return (x + y); /* return NaN if x or y is NAN */ 83 signy = signbitl(y); 84 signx = signbitl(x); 85 if (x == one) 86 return (atanl(y)); 87 /* Ensure sign indicators are boolean */ 88 m = (signy != 0) + (signx != 0) + (signx != 0); 89 90 /* when y = 0 */ 91 if (y == zero) 92 switch (m) { 93 case 0: 94 return (y); /* atan(+0,+anything) */ 95 case 1: 96 return (y); /* atan(-0,+anything) */ 97 case 2: 98 return (PI + tiny); /* atan(+0,-anything) */ 99 case 3: 100 return (-PI - tiny); /* atan(-0,-anything) */ 101 } 102 103 /* when x = 0 */ 104 if (x == zero) 105 return (signy != 0 ? -PIo2 - tiny : PIo2 + tiny); 106 107 /* when x is INF */ 108 if (!finitel(x)) { 109 if (!finitel(y)) { 110 switch (m) { 111 case 0: 112 return (PIo4 + tiny); /* atan(+INF,+INF) */ 113 case 1: 114 return (-PIo4 - tiny); /* atan(-INF,+INF) */ 115 case 2: 116 return (PI3o4 + tiny); /* atan(+INF,-INF) */ 117 case 3: 118 return (-PI3o4 - tiny); /* atan(-INF,-INF) */ 119 } 120 } else { 121 switch (m) { 122 case 0: 123 return (zero); /* atan(+...,+INF) */ 124 case 1: 125 return (-zero); /* atan(-...,+INF) */ 126 case 2: 127 return (PI + tiny); /* atan(+...,-INF) */ 128 case 3: 129 return (-PI - tiny); /* atan(-...,-INF) */ 130 } 131 } 132 } 133 /* when y is INF */ 134 if (!finitel(y)) 135 return (signy != 0 ? -PIo2 - tiny : PIo2 + tiny); 136 137 /* compute y/x */ 138 x = fabsl(x); 139 y = fabsl(y); 140 t = PI_lo; 141 k = (ilogbl(y) - ilogbl(x)); 142 143 if (k > 120) 144 z = PIo2 + half * t; 145 else if (m > 1 && k < -120) 146 z = zero; 147 else 148 z = atanl(y / x); 149 150 switch (m) { 151 case 0: 152 return (z); /* atan(+,+) */ 153 case 1: 154 return (-z); /* atan(-,+) */ 155 case 2: 156 return (PI - (z - t)); /* atan(+,-) */ 157 case 3: 158 return ((z - t) - PI); /* atan(-,-) */ 159 } 160 /* NOTREACHED */ 161 return (0.0L); 162 } 163