Back to index

glibc  2.9
e_atan2.c
Go to the documentation of this file.
00001 /* Copyright (C) 1997, 1999 Free Software Foundation, Inc.
00002    This file is part of the GNU C Library.
00003 
00004    The GNU C Library is free software; you can redistribute it and/or
00005    modify it under the terms of the GNU Lesser General Public
00006    License as published by the Free Software Foundation; either
00007    version 2.1 of the License, or (at your option) any later version.
00008 
00009    The GNU C Library is distributed in the hope that it will be useful,
00010    but WITHOUT ANY WARRANTY; without even the implied warranty of
00011    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00012    Lesser General Public License for more details.
00013 
00014    You should have received a copy of the GNU Lesser General Public
00015    License along with the GNU C Library; if not, write to the Free
00016    Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
00017    02111-1307 USA.  */
00018 
00019 #include <math.h>
00020 #include "math_private.h"
00021 #include "mathimpl.h"
00022 
00023 #ifndef SUFF
00024 #define SUFF
00025 #endif
00026 #ifndef float_type
00027 #define float_type double
00028 #endif
00029 
00030 #define CONCATX(a,b) __CONCAT(a,b)
00031 #define s(name) CONCATX(name,SUFF)
00032 #define m81(func) __m81_u(s(func))
00033 
00034 float_type
00035 s(__ieee754_atan2) (float_type y, float_type x)
00036 {
00037   float_type pi, pi_2, z;
00038   unsigned long y_cond, x_cond;
00039 
00040   __asm ("fmovecr%.x %#0, %0" : "=f" (pi));
00041   __asm ("fscale%.w %#-1, %0" : "=f" (pi_2) : "0" (pi));
00042   y_cond = __m81_test (y);
00043   x_cond = __m81_test (x);
00044 
00045   if ((x_cond | y_cond) & __M81_COND_NAN)
00046     z = x + y;
00047   else if (y_cond & __M81_COND_ZERO)
00048     {
00049       if (x_cond & __M81_COND_NEG)
00050        z = y_cond & __M81_COND_NEG ? -pi : pi;
00051       else
00052        z = y;
00053     }
00054   else if (x_cond & __M81_COND_INF)
00055     {
00056       if (y_cond & __M81_COND_INF)
00057        {
00058          float_type pi_4;
00059          __asm ("fscale%.w %#-2, %0" : "=f" (pi_4) : "0" (pi));
00060          z = x_cond & __M81_COND_NEG ? 3 * pi_4 : pi_4;
00061        }
00062       else
00063        z = x_cond & __M81_COND_NEG ? pi : 0;
00064       if (y_cond & __M81_COND_NEG)
00065        z = -z;
00066     }
00067   else if (y_cond & __M81_COND_INF)
00068     z = y_cond & __M81_COND_NEG ? -pi_2 : pi_2;
00069   else if (x_cond & __M81_COND_NEG)
00070     {
00071       if (y_cond & __M81_COND_NEG)
00072        {
00073          if (-x > -y)
00074            z = -pi + m81(__atan) (y / x);
00075          else
00076            z = -pi_2 - m81(__atan) (x / y);
00077        }
00078       else
00079        {
00080          if (-x > y)
00081            z = pi + m81(__atan) (y / x);
00082          else
00083            z = pi_2 - m81(__atan) (x / y);
00084        }
00085     }
00086   else
00087     {
00088       if (y_cond & __M81_COND_NEG)
00089        {
00090          if (x > -y)
00091            z = m81(__atan) (y / x);
00092          else
00093            z = -pi_2 - m81(__atan) (x / y);
00094        }
00095       else
00096        {
00097          if (x > y)
00098            z = m81(__atan) (y / x);
00099          else
00100            z = pi_2 - m81(__atan) (x / y);
00101        }
00102     }
00103   return z;
00104 }