nsz Git - libm/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 volatile float
  19 tiny = 1.0e-30;
  20 static const float
  21 zero = 0.0,
  22 pi_o_4 = 7.8539818525e-01, /* 0x3f490fdb */
  23 pi_o_2 = 1.5707963705e+00, /* 0x3fc90fdb */
  24 pi     = 3.1415927410e+00; /* 0x40490fdb */
  25 static volatile float
  26 pi_lo  = -8.7422776573e-08; /* 0xb3bbbd2e */
  27
  28 float atan2f(float y, float x)
  29 {
  30         float z;
  31         int32_t k,m,hx,hy,ix,iy;
  32
  33         GET_FLOAT_WORD(hx, x);
  34         ix = hx & 0x7fffffff;
  35         GET_FLOAT_WORD(hy, y);
  36         iy = hy & 0x7fffffff;
  37         if (ix > 0x7f800000 || iy > 0x7f800000)  /* x or y is NaN */
  38                 return x+y;
  39         if (hx == 0x3f800000)  /* x=1.0 */
  40                 return atanf(y);
  41         m = ((hy>>31)&1) | ((hx>>30)&2);  /* 2*sign(x)+sign(y) */
  42
  43         /* when y = 0 */
  44         if (iy == 0) {
  45                 switch (m) {
  46                 case 0:
  47                 case 1: return y;        /* atan(+-0,+anything)=+-0 */
  48                 case 2: return  pi+tiny; /* atan(+0,-anything) = pi */
  49                 case 3: return -pi-tiny; /* atan(-0,-anything) =-pi */
  50                 }
  51         }
  52         /* when x = 0 */
  53         if (ix == 0)
  54                 return hy < 0 ? -pi_o_2-tiny : pi_o_2+tiny;
  55         /* when x is INF */
  56         if (ix == 0x7f800000) {
  57                 if (iy == 0x7f800000) {
  58                         switch (m) {
  59                         case 0: return  pi_o_4+tiny; /* atan(+INF,+INF) */
  60                         case 1: return -pi_o_4-tiny; /* atan(-INF,+INF) */
  61                         case 2: return (float)3.0*pi_o_4+tiny;  /*atan(+INF,-INF)*/
  62                         case 3: return (float)-3.0*pi_o_4-tiny; /*atan(-INF,-INF)*/
  63                         }
  64                 } else {
  65                         switch (m) {
  66                         case 0: return  zero;    /* atan(+...,+INF) */
  67                         case 1: return -zero;    /* atan(-...,+INF) */
  68                         case 2: return  pi+tiny; /* atan(+...,-INF) */
  69                         case 3: return -pi-tiny; /* atan(-...,-INF) */
  70                         }
  71                 }
  72         }
  73         /* when y is INF */
  74         if (iy == 0x7f800000)
  75                 return hy < 0 ? -pi_o_2-tiny : pi_o_2+tiny;
  76
  77         /* compute y/x */
  78         k = (iy-ix)>>23;
  79         if (k > 26) {                  /* |y/x| >  2**26 */
  80                 z = pi_o_2+(float)0.5*pi_lo;
  81                 m &= 1;
  82         } else if (k < -26 && hx < 0)  /* 0 > |y|/x > -2**-26 */
  83                 z = 0.0;
  84         else                           /* safe to do y/x */
  85                 z = atanf(fabsf(y/x));
  86         switch (m) {
  87         case 0: return z;              /* atan(+,+) */
  88         case 1: return -z;             /* atan(-,+) */
  89         case 2: return pi - (z-pi_lo); /* atan(+,-) */
  90         default: /* case 3 */
  91                 return (z-pi_lo) - pi; /* atan(-,-) */
  92         }
  93 }