atan2f.c 2.5 KB

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