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 2005 Sun Microsystems, Inc.  All rights reserved.
  28  * Use is subject to license terms.
  29  */
  30 
  31 #pragma weak __atan2f = atan2f
  32 
  33 #include "libm.h"
  34 
  35 #if defined(__i386) && !defined(__amd64)
  36 extern int __swapRP(int);
  37 #endif
  38 
  39 /*
  40  * For i = 0, ..., 192, let x[i] be the double precision number whose
  41  * high order 32 bits are 0x3f900000 + (i << 16) and whose low order
  42  * 32 bits are zero.  Then TBL[i] := atan(x[i]) to double precision.
  43  */
  44 
  45 static const double TBL[] = {
  46         1.56237286204768313e-02,
  47         1.66000375562312640e-02,
  48         1.75763148444955872e-02,
  49         1.85525586258889763e-02,
  50         1.95287670414137082e-02,
  51         2.05049382324763683e-02,
  52         2.14810703409090559e-02,
  53         2.24571615089905717e-02,
  54         2.34332098794675855e-02,
  55         2.44092135955758099e-02,
  56         2.53851708010611396e-02,
  57         2.63610796402007873e-02,
  58         2.73369382578244127e-02,
  59         2.83127447993351995e-02,
  60         2.92884974107309737e-02,
  61         3.02641942386252458e-02,
  62         3.12398334302682774e-02,
  63         3.31909314971115949e-02,
  64         3.51417768027967800e-02,
  65         3.70923545503918164e-02,
  66         3.90426499551669928e-02,
  67         4.09926482452637811e-02,
  68         4.29423346623621707e-02,
  69         4.48916944623464972e-02,
  70         4.68407129159696539e-02,
  71         4.87893753095156174e-02,
  72         5.07376669454602178e-02,
  73         5.26855731431300420e-02,
  74         5.46330792393594777e-02,
  75         5.65801705891457105e-02,
  76         5.85268325663017702e-02,
  77         6.04730505641073168e-02,
  78         6.24188099959573500e-02,
  79         6.63088949198234884e-02,
  80         7.01969710718705203e-02,
  81         7.40829225490337306e-02,
  82         7.79666338315423008e-02,
  83         8.18479898030765457e-02,
  84         8.57268757707448092e-02,
  85         8.96031774848717461e-02,
  86         9.34767811585894698e-02,
  87         9.73475734872236709e-02,
  88         1.01215441667466668e-01,
  89         1.05080273416329528e-01,
  90         1.08941956989865793e-01,
  91         1.12800381201659389e-01,
  92         1.16655435441069349e-01,
  93         1.20507009691224562e-01,
  94         1.24354994546761438e-01,
  95         1.32039761614638762e-01,
  96         1.39708874289163648e-01,
  97         1.47361481088651630e-01,
  98         1.54996741923940973e-01,
  99         1.62613828597948568e-01,
 100         1.70211925285474408e-01,
 101         1.77790228992676075e-01,
 102         1.85347949995694761e-01,
 103         1.92884312257974672e-01,
 104         2.00398553825878512e-01,
 105         2.07889927202262986e-01,
 106         2.15357699697738048e-01,
 107         2.22801153759394521e-01,
 108         2.30219587276843718e-01,
 109         2.37612313865471242e-01,
 110         2.44978663126864143e-01,
 111         2.59629629408257512e-01,
 112         2.74167451119658789e-01,
 113         2.88587361894077410e-01,
 114         3.02884868374971417e-01,
 115         3.17055753209147029e-01,
 116         3.31096076704132103e-01,
 117         3.45002177207105132e-01,
 118         3.58770670270572245e-01,
 119         3.72398446676754202e-01,
 120         3.85882669398073752e-01,
 121         3.99220769575252543e-01,
 122         4.12410441597387323e-01,
 123         4.25449637370042266e-01,
 124         4.38336559857957830e-01,
 125         4.51069655988523499e-01,
 126         4.63647609000806094e-01,
 127         4.88333951056405535e-01,
 128         5.12389460310737732e-01,
 129         5.35811237960463704e-01,
 130         5.58599315343562441e-01,
 131         5.80756353567670414e-01,
 132         6.02287346134964152e-01,
 133         6.23199329934065904e-01,
 134         6.43501108793284371e-01,
 135         6.63202992706093286e-01,
 136         6.82316554874748071e-01,
 137         7.00854407884450192e-01,
 138         7.18829999621624527e-01,
 139         7.36257428981428097e-01,
 140         7.53151280962194414e-01,
 141         7.69526480405658297e-01,
 142         7.85398163397448279e-01,
 143         8.15691923316223422e-01,
 144         8.44153986113171051e-01,
 145         8.70903457075652976e-01,
 146         8.96055384571343927e-01,
 147         9.19719605350416858e-01,
 148         9.42000040379463610e-01,
 149         9.62994330680936206e-01,
 150         9.82793723247329054e-01,
 151         1.00148313569423464e+00,
 152         1.01914134426634972e+00,
 153         1.03584125300880014e+00,
 154         1.05165021254837376e+00,
 155         1.06663036531574362e+00,
 156         1.08083900054116833e+00,
 157         1.09432890732118993e+00,
 158         1.10714871779409041e+00,
 159         1.13095374397916038e+00,
 160         1.15257199721566761e+00,
 161         1.17227388112847630e+00,
 162         1.19028994968253166e+00,
 163         1.20681737028525249e+00,
 164         1.22202532321098967e+00,
 165         1.23605948947808186e+00,
 166         1.24904577239825443e+00,
 167         1.26109338225244039e+00,
 168         1.27229739520871732e+00,
 169         1.28274087974427076e+00,
 170         1.29249666778978534e+00,
 171         1.30162883400919616e+00,
 172         1.31019393504755555e+00,
 173         1.31824205101683711e+00,
 174         1.32581766366803255e+00,
 175         1.33970565959899957e+00,
 176         1.35212738092095464e+00,
 177         1.36330010035969384e+00,
 178         1.37340076694501589e+00,
 179         1.38257482149012589e+00,
 180         1.39094282700241845e+00,
 181         1.39860551227195762e+00,
 182         1.40564764938026987e+00,
 183         1.41214106460849531e+00,
 184         1.41814699839963154e+00,
 185         1.42371797140649403e+00,
 186         1.42889927219073276e+00,
 187         1.43373015248470903e+00,
 188         1.43824479449822262e+00,
 189         1.44247309910910193e+00,
 190         1.44644133224813509e+00,
 191         1.45368758222803240e+00,
 192         1.46013910562100091e+00,
 193         1.46591938806466282e+00,
 194         1.47112767430373470e+00,
 195         1.47584462045214027e+00,
 196         1.48013643959415142e+00,
 197         1.48405798811891154e+00,
 198         1.48765509490645531e+00,
 199         1.49096634108265924e+00,
 200         1.49402443552511865e+00,
 201         1.49685728913695626e+00,
 202         1.49948886200960629e+00,
 203         1.50193983749385196e+00,
 204         1.50422816301907281e+00,
 205         1.50636948736934317e+00,
 206         1.50837751679893928e+00,
 207         1.51204050407917401e+00,
 208         1.51529782154917969e+00,
 209         1.51821326518395483e+00,
 210         1.52083793107295384e+00,
 211         1.52321322351791322e+00,
 212         1.52537304737331958e+00,
 213         1.52734543140336587e+00,
 214         1.52915374769630819e+00,
 215         1.53081763967160667e+00,
 216         1.53235373677370856e+00,
 217         1.53377621092096650e+00,
 218         1.53509721411557254e+00,
 219         1.53632722579538861e+00,
 220         1.53747533091664934e+00,
 221         1.53854944435964280e+00,
 222         1.53955649336462841e+00,
 223         1.54139303859089161e+00,
 224         1.54302569020147562e+00,
 225         1.54448660954197448e+00,
 226         1.54580153317597646e+00,
 227         1.54699130060982659e+00,
 228         1.54807296595325550e+00,
 229         1.54906061995310385e+00,
 230         1.54996600675867957e+00,
 231         1.55079899282174605e+00,
 232         1.55156792769518947e+00,
 233         1.55227992472688747e+00,
 234         1.55294108165534417e+00,
 235         1.55355665560036682e+00,
 236         1.55413120308095598e+00,
 237         1.55466869295126031e+00,
 238         1.55517259817441977e+00,
 239 };
 240 
 241 static const double
 242         pio4 = 7.8539816339744827900e-01,
 243         pio2 = 1.5707963267948965580e+00,
 244         negpi = -3.1415926535897931160e+00,
 245         q1 = -3.3333333333296428046e-01,
 246         q2 = 1.9999999186853752618e-01,
 247         zero = 0.0;
 248 static const float two24 = 16777216.0;
 249 
 250 float
 251 atan2f(float fy, float fx)
 252 {
 253         double a, t, s, dbase;
 254         float x, y, base;
 255         int i, k, hx, hy, ix, iy, sign;
 256 
 257 #if defined(__i386) && !defined(__amd64)
 258         int rp;
 259 #endif
 260 
 261         iy = *(int *)&fy;
 262         ix = *(int *)&fx;
 263         hy = iy & ~0x80000000;
 264         hx = ix & ~0x80000000;
 265 
 266         sign = 0;
 267 
 268         if (hy > hx) {
 269                 x = fy;
 270                 y = fx;
 271                 i = hx;
 272                 hx = hy;
 273                 hy = i;
 274 
 275                 if (iy < 0) {
 276                         x = -x;
 277                         sign = 1;
 278                 }
 279 
 280                 if (ix < 0) {
 281                         y = -y;
 282                         a = pio2;
 283                 } else {
 284                         a = -pio2;
 285                         sign = 1 - sign;
 286                 }
 287         } else {
 288                 y = fy;
 289                 x = fx;
 290 
 291                 if (iy < 0) {
 292                         y = -y;
 293                         sign = 1;
 294                 }
 295 
 296                 if (ix < 0) {
 297                         x = -x;
 298                         a = negpi;
 299                         sign = 1 - sign;
 300                 } else {
 301                         a = zero;
 302                 }
 303         }
 304 
 305         if (hx >= 0x7f800000 || hx - hy >= 0x0c800000) {
 306                 if (hx >= 0x7f800000) {
 307                         if (hx > 0x7f800000) /* nan */
 308                                 return (x * y);
 309                         else if (hy >= 0x7f800000)
 310                                 a += pio4;
 311                 } else if ((int)a == 0) {
 312                         a = (double)y / x;
 313                 }
 314 
 315                 return ((float)((sign) ? -a : a));
 316         }
 317 
 318         if (hy < 0x00800000) {
 319                 if (hy == 0)
 320                         return ((float)((sign) ? -a : a));
 321 
 322                 /* scale subnormal y */
 323                 y *= two24;
 324                 x *= two24;
 325                 hy = *(int *)&y;
 326                 hx = *(int *)&x;
 327         }
 328 
 329 #if defined(__i386) && !defined(__amd64)
 330         rp = __swapRP(fp_extended);
 331 #endif
 332         k = (hy - hx + 0x3f800000) & 0xfff80000;
 333 
 334         if (k >= 0x3c800000) {               /* |y/x| >= 1/64 */
 335                 *(int *)&base = k;
 336                 k = (k - 0x3c800000) >> 19;
 337                 a += TBL[k];
 338         } else {
 339                 /*
 340                  * For some reason this is faster on USIII than just
 341                  * doing t = y/x in this case.
 342                  */
 343                 *(int *)&base = 0;
 344         }
 345 
 346         dbase = (double)base;
 347         t = (y - x * dbase) / (x + y * dbase);
 348         s = t * t;
 349         a = (a + t) + t * s * (q1 + s * q2);
 350 #if defined(__i386) && !defined(__amd64)
 351         if (rp != fp_extended)
 352                 (void) __swapRP(rp);
 353 #endif
 354         return ((float)((sign) ? -a : a));
 355 }