usr/src/lib/libm/common/Q/atan2l.c

   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 }