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 long double t, z; 78 int k, m, signy, signx; 79 80 if (x != x || y != y) 81 return (x + y); /* return NaN if x or y is NAN */ 82 signy = signbitl(y); 83 signx = signbitl(x); 84 if (x == one) 85 return (atanl(y)); 86 m = signy + signx + signx; 87 88 /* when y = 0 */ 89 if (y == zero) 90 switch (m) { 91 case 0: 92 return (y); /* atan(+0,+anything) */ 93 case 1: 94 return (y); /* atan(-0,+anything) */ 95 case 2: 96 return (PI + tiny); /* atan(+0,-anything) */ 97 case 3: 98 return (-PI - tiny); /* atan(-0,-anything) */ 99 } 100 101 /* when x = 0 */ 102 if (x == zero) 103 return (signy == 1 ? -PIo2 - tiny : PIo2 + tiny); 104 105 /* when x is INF */ 106 if (!finitel(x)) { 107 if (!finitel(y)) { 108 switch (m) { 109 case 0: 110 return (PIo4 + tiny); /* atan(+INF,+INF) */ 111 case 1: 112 return (-PIo4 - tiny); /* atan(-INF,+INF) */ 113 case 2: 114 return (PI3o4 + tiny); /* atan(+INF,-INF) */ 115 case 3: 116 return (-PI3o4 - tiny); /* atan(-INF,-INF) */ 117 } 118 } else { 119 switch (m) { 120 case 0: 121 return (zero); /* atan(+...,+INF) */ 122 case 1: 123 return (-zero); /* atan(-...,+INF) */ 124 case 2: 125 return (PI + tiny); /* atan(+...,-INF) */ 126 case 3: 127 return (-PI - tiny); /* atan(-...,-INF) */ 128 } 129 } 130 } 131 /* when y is INF */ 132 if (!finitel(y)) 133 return (signy == 1 ? -PIo2 - tiny : PIo2 + tiny); 134 135 /* compute y/x */ 136 x = fabsl(x); 137 y = fabsl(y); 138 t = PI_lo; 139 k = (ilogbl(y) - ilogbl(x)); 140 141 if (k > 120) 142 z = PIo2 + half * t; 143 else if (m > 1 && k < -120) 144 z = zero; 145 else 146 z = atanl(y / x); 147 148 switch (m) { 149 case 0: 150 return (z); /* atan(+,+) */ 151 case 1: 152 return (-z); /* atan(-,+) */ 153 case 2: 154 return (PI - (z - t)); /* atan(+,-) */ 155 case 3: 156 return ((z - t) - PI); /* atan(-,-) */ 157 } 158 /* NOTREACHED */ 159 return 0.0L; 160 } 161