sysdeps/m68k/fpu/e_atan2.c

   1 /* Copyright (C) 1997, 1999 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, write to the Free
  16    Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
  17    02111-1307 USA.  */
  18
  19 #include <math.h>
  20 #include "math_private.h"
  21 #include "mathimpl.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 }