1  /* Copyright (C) 1997-2023 Free Software Foundation, Inc.
       2     This file is part of the GNU C Library.
       3  
       4     The GNU C Library is free software; you can redistribute it and/or
       5     modify it under the terms of the GNU Lesser General Public
       6     License as published by the Free Software Foundation; either
       7     version 2.1 of the License, or (at your option) any later version.
       8  
       9     The GNU C Library is distributed in the hope that it will be useful,
      10     but WITHOUT ANY WARRANTY; without even the implied warranty of
      11     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
      12     Lesser General Public License for more details.
      13  
      14     You should have received a copy of the GNU Lesser General Public
      15     License along with the GNU C Library.  If not, see
      16     <https://www.gnu.org/licenses/>.  */
      17  
      18  #include <math.h>
      19  #include <math_private.h>
      20  #include "mathimpl.h"
      21  #include <libm-alias-finite.h>
      22  
      23  #ifndef SUFF
      24  #define SUFF
      25  #endif
      26  #ifndef float_type
      27  #define float_type double
      28  #endif
      29  
      30  #define CONCATX(a,b) __CONCAT(a,b)
      31  #define s(name) CONCATX(name,SUFF)
      32  #define m81(func) __m81_u(s(func))
      33  
      34  float_type
      35  s(__ieee754_atan2) (float_type y, float_type x)
      36  {
      37    float_type pi, pi_2, z;
      38    unsigned long y_cond, x_cond;
      39  
      40    __asm ("fmovecr%.x %#0, %0" : "=f" (pi));
      41    __asm ("fscale%.w %#-1, %0" : "=f" (pi_2) : "0" (pi));
      42    y_cond = __m81_test (y);
      43    x_cond = __m81_test (x);
      44  
      45    if ((x_cond | y_cond) & __M81_COND_NAN)
      46      z = x + y;
      47    else if (y_cond & __M81_COND_ZERO)
      48      {
      49        if (x_cond & __M81_COND_NEG)
      50  	z = y_cond & __M81_COND_NEG ? -pi : pi;
      51        else
      52  	z = y;
      53      }
      54    else if (x_cond & __M81_COND_INF)
      55      {
      56        if (y_cond & __M81_COND_INF)
      57  	{
      58  	  float_type pi_4;
      59  	  __asm ("fscale%.w %#-2, %0" : "=f" (pi_4) : "0" (pi));
      60  	  z = x_cond & __M81_COND_NEG ? 3 * pi_4 : pi_4;
      61  	}
      62        else
      63  	z = x_cond & __M81_COND_NEG ? pi : 0;
      64        if (y_cond & __M81_COND_NEG)
      65  	z = -z;
      66      }
      67    else if (y_cond & __M81_COND_INF)
      68      z = y_cond & __M81_COND_NEG ? -pi_2 : pi_2;
      69    else if (x_cond & __M81_COND_NEG)
      70      {
      71        if (y_cond & __M81_COND_NEG)
      72  	{
      73  	  if (-x > -y)
      74  	    z = -pi + m81(__atan) (y / x);
      75  	  else
      76  	    z = -pi_2 - m81(__atan) (x / y);
      77  	}
      78        else
      79  	{
      80  	  if (-x > y)
      81  	    z = pi + m81(__atan) (y / x);
      82  	  else
      83  	    z = pi_2 - m81(__atan) (x / y);
      84  	}
      85      }
      86    else
      87      {
      88        if (y_cond & __M81_COND_NEG)
      89  	{
      90  	  if (x > -y)
      91  	    z = m81(__atan) (y / x);
      92  	  else
      93  	    z = -pi_2 - m81(__atan) (x / y);
      94  	}
      95        else
      96  	{
      97  	  if (x > y)
      98  	    z = m81(__atan) (y / x);
      99  	  else
     100  	    z = pi_2 - m81(__atan) (x / y);
     101  	}
     102      }
     103    return z;
     104  }
     105  libm_alias_finite (s(__ieee754_atan2), s (__atan2))