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 /*
  27  * Copyright 2006 Sun Microsystems, Inc.  All rights reserved.
  28  * Use is subject to license terms.
  29  */
  30 
  31 /*
  32  * atan2l(y,x)
  33  *
  34  * Method :
  35  *      1. Reduce y to positive by atan2(y,x)=-atan2(-y,x).
  36  *      2. Reduce x to positive by (if x and y are unexceptional):
  37  *              ARG (x+iy) = arctan(y/x)           ... if x > 0,
  38  *              ARG (x+iy) = pi - arctan[y/(-x)]   ... if x < 0,
  39  *
  40  * Special cases:
  41  *
  42  *      ATAN2((anything), NaN ) is NaN;
  43  *      ATAN2(NAN , (anything) ) is NaN;
  44  *      ATAN2(+-0, +(anything but NaN)) is +-0  ;
  45  *      ATAN2(+-0, -(anything but NaN)) is +-PI ;
  46  *      ATAN2(+-(anything but 0 and NaN), 0) is +-PI/2;
  47  *      ATAN2(+-(anything but INF and NaN), +INF) is +-0 ;
  48  *      ATAN2(+-(anything but INF and NaN), -INF) is +-PI;
  49  *      ATAN2(+-INF,+INF ) is +-PI/4 ;
  50  *      ATAN2(+-INF,-INF ) is +-3PI/4;
  51  *      ATAN2(+-INF, (anything but,0,NaN, and INF)) is +-PI/2;
  52  *
  53  * Constants:
  54  * The hexadecimal values are the intended ones for the following constants.
  55  * The decimal values may be used, provided that the compiler will convert
  56  * from decimal to binary accurately enough to produce the hexadecimal values
  57  * shown.
  58  */
  59 
  60 #pragma weak __atan2l = atan2l
  61 
  62 #include "libm.h"
  63 #include "longdouble.h"
  64 
  65 static const long double 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 
  84         signy = signbitl(y);
  85         signx = signbitl(x);
  86 
  87         if (x == one)
  88                 return (atanl(y));
  89 
  90         /* Ensure sign indicators are boolean */
  91         m = (signy != 0) + (signx != 0) + (signx != 0);
  92 
  93         /* when y = 0 */
  94         if (y == zero) {
  95                 switch (m) {
  96                 case 0:
  97                         return (y);             /* atan(+0,+anything) */
  98                 case 1:
  99                         return (y);             /* atan(-0,+anything) */
 100                 case 2:
 101                         return (PI + tiny);     /* atan(+0,-anything) */
 102                 case 3:
 103                         return (-PI - tiny);    /* atan(-0,-anything) */
 104                 }
 105         }
 106 
 107         /* when x = 0 */
 108         if (x == zero)
 109                 return (signy != 0 ? -PIo2 - tiny : PIo2 + tiny);
 110 
 111         /* when x is INF */
 112         if (!finitel(x)) {
 113                 if (!finitel(y)) {
 114                         switch (m) {
 115                         case 0:
 116                                 return (PIo4 + tiny);   /* atan(+INF,+INF) */
 117                         case 1:
 118                                 return (-PIo4 - tiny);  /* atan(-INF,+INF) */
 119                         case 2:
 120                                 return (PI3o4 + tiny);  /* atan(+INF,-INF) */
 121                         case 3:
 122                                 return (-PI3o4 - tiny); /* atan(-INF,-INF) */
 123                         }
 124                 } else {
 125                         switch (m) {
 126                         case 0:
 127                                 return (zero);          /* atan(+...,+INF) */
 128                         case 1:
 129                                 return (-zero);         /* atan(-...,+INF) */
 130                         case 2:
 131                                 return (PI + tiny);     /* atan(+...,-INF) */
 132                         case 3:
 133                                 return (-PI - tiny);    /* atan(-...,-INF) */
 134                         }
 135                 }
 136         }
 137 
 138         /* when y is INF */
 139         if (!finitel(y))
 140                 return (signy != 0 ? -PIo2 - tiny : PIo2 + tiny);
 141 
 142         /* compute y/x */
 143         x = fabsl(x);
 144         y = fabsl(y);
 145         t = PI_lo;
 146         k = (ilogbl(y) - ilogbl(x));
 147 
 148         if (k > 120)
 149                 z = PIo2 + half * t;
 150         else if (m > 1 && k < -120)
 151                 z = zero;
 152         else
 153                 z = atanl(y / x);
 154 
 155         switch (m) {
 156         case 0:
 157                 return (z);             /* atan(+,+) */
 158         case 1:
 159                 return (-z);            /* atan(-,+) */
 160         case 2:
 161                 return (PI - (z - t));  /* atan(+,-) */
 162         case 3:
 163                 return ((z - t) - PI);  /* atan(-,-) */
 164         }
 165 
 166         /* NOTREACHED */
 167         return (0.0L);
 168 }