nsz Git - musl/blob - src/math/atan2f.c

   1 /* origin: FreeBSD /usr/src/lib/msun/src/e_atan2f.c */
   2 /*
   3  * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
   4  */
   5 /*
   6  * ====================================================
   7  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
   8  *
   9  * Developed at SunPro, a Sun Microsystems, Inc. business.
  10  * Permission to use, copy, modify, and distribute this
  11  * software is freely granted, provided that this notice
  12  * is preserved.
  13  * ====================================================
  14  */
  15
  16 #include "libm.h"
  17
  18 static const volatile float
  19 tiny = 1.0e-30;
  20 static const float
  21 pi_o_4 = 7.8539818525e-01, /* 0x3f490fdb */
  22 pi_o_2 = 1.5707963705e+00, /* 0x3fc90fdb */
  23 pi     = 3.1415927410e+00; /* 0x40490fdb */
  24 static const volatile float
  25 pi_lo  = -8.7422776573e-08; /* 0xb3bbbd2e */
  26
  27 float atan2f(float y, float x)
  28 {
  29         float z;
  30         int32_t k,m,hx,hy,ix,iy;
  31
  32         GET_FLOAT_WORD(hx, x);
  33         ix = hx & 0x7fffffff;
  34         GET_FLOAT_WORD(hy, y);
  35         iy = hy & 0x7fffffff;
  36         if (ix > 0x7f800000 || iy > 0x7f800000)  /* x or y is NaN */
  37                 return x+y;
  38         if (hx == 0x3f800000)  /* x=1.0 */
  39                 return atanf(y);
  40         m = ((hy>>31)&1) | ((hx>>30)&2);  /* 2*sign(x)+sign(y) */
  41
  42         /* when y = 0 */
  43         if (iy == 0) {
  44                 switch (m) {
  45                 case 0:
  46                 case 1: return y;        /* atan(+-0,+anything)=+-0 */
  47                 case 2: return  pi+tiny; /* atan(+0,-anything) = pi */
  48                 case 3: return -pi-tiny; /* atan(-0,-anything) =-pi */
  49                 }
  50         }
  51         /* when x = 0 */
  52         if (ix == 0)
  53                 return hy < 0 ? -pi_o_2-tiny : pi_o_2+tiny;
  54         /* when x is INF */
  55         if (ix == 0x7f800000) {
  56                 if (iy == 0x7f800000) {
  57                         switch (m) {
  58                         case 0: return  pi_o_4+tiny; /* atan(+INF,+INF) */
  59                         case 1: return -pi_o_4-tiny; /* atan(-INF,+INF) */
  60                         case 2: return 3.0f*pi_o_4+tiny;  /*atan(+INF,-INF)*/
  61                         case 3: return -3.0f*pi_o_4-tiny; /*atan(-INF,-INF)*/
  62                         }
  63                 } else {
  64                         switch (m) {
  65                         case 0: return  0.0f;    /* atan(+...,+INF) */
  66                         case 1: return -0.0f;    /* atan(-...,+INF) */
  67                         case 2: return  pi+tiny; /* atan(+...,-INF) */
  68                         case 3: return -pi-tiny; /* atan(-...,-INF) */
  69                         }
  70                 }
  71         }
  72         /* when y is INF */
  73         if (iy == 0x7f800000)
  74                 return hy < 0 ? -pi_o_2-tiny : pi_o_2+tiny;
  75
  76         /* compute y/x */
  77         k = (iy-ix)>>23;
  78         if (k > 26) {                  /* |y/x| >  2**26 */
  79                 z = pi_o_2 + 0.5f*pi_lo;
  80                 m &= 1;
  81         } else if (k < -26 && hx < 0)  /* 0 > |y|/x > -2**-26 */
  82                 z = 0.0;
  83         else                           /* safe to do y/x */
  84                 z = atanf(fabsf(y/x));
  85         switch (m) {
  86         case 0: return z;              /* atan(+,+) */
  87         case 1: return -z;             /* atan(-,+) */
  88         case 2: return pi - (z-pi_lo); /* atan(+,-) */
  89         default: /* case 3 */
  90                 return (z-pi_lo) - pi; /* atan(-,-) */
  91         }
  92 }