ports/sysdeps/m68k/m680x0/fpu/s_cexp.c

   1 /* Complex exponential function.  m68k fpu version
   2    Copyright (C) 1997-2014 Free Software Foundation, Inc.
   3    This file is part of the GNU C Library.
   4    Contributed by Andreas Schwab <schwab@issan.informatik.uni-dortmund.de>
   5
   6    The GNU C Library is free software; you can redistribute it and/or
   7    modify it under the terms of the GNU Lesser General Public
   8    License as published by the Free Software Foundation; either
   9    version 2.1 of the License, or (at your option) any later version.
  10
  11    The GNU C Library is distributed in the hope that it will be useful,
  12    but WITHOUT ANY WARRANTY; without even the implied warranty of
  13    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14    Lesser General Public License for more details.
  15
  16    You should have received a copy of the GNU Lesser General Public
  17    License along with the GNU C Library.  If not, see
  18    <http://www.gnu.org/licenses/>.  */
  19
  20 #include <float.h>
  21 #include <complex.h>
  22 #include <math.h>
  23 #include "mathimpl.h"
  24
  25 #ifndef SUFF
  26 #define SUFF
  27 #endif
  28 #ifndef float_type
  29 #define float_type double
  30 #endif
  31
  32 #define CONCATX(a,b) __CONCAT(a,b)
  33 #define s(name) CONCATX(name,SUFF)
  34 #define m81(func) __m81_u(s(func))
  35
  36 __complex__ float_type
  37 s(__cexp) (__complex__ float_type x)
  38 {
  39   __complex__ float_type retval;
  40   unsigned long ix_cond;
  41
  42   ix_cond = __m81_test (__imag__ x);
  43
  44   if ((ix_cond & (__M81_COND_NAN|__M81_COND_INF)) == 0)
  45     {
  46       /* Imaginary part is finite.  */
  47       unsigned long rx_cond = __m81_test (__real__ x);
  48
  49       if ((rx_cond & (__M81_COND_NAN|__M81_COND_INF)) == 0)
  50         {
  51           const int t = (int) ((LDBL_MAX_EXP - 1) * M_LN2l);
  52           long double sin_ix, cos_ix, exp_val;
  53
  54           __m81_u (__sincosl) (__imag__ x, &sin_ix, &cos_ix);
  55
  56           if (__real__ x > t)
  57             {
  58               long double exp_t = __m81_u(__ieee754_expl) (t);
  59               __real__ x -= t;
  60               sin_ix *= exp_t;
  61               cos_ix *= exp_t;
  62               if (__real__ x > t)
  63                 {
  64                   __real__ x -= t;
  65                   sin_ix *= exp_t;
  66                   cos_ix *= exp_t;
  67                 }
  68             }
  69
  70           exp_val = __m81_u(__ieee754_expl) (__real__ x);
  71           __real__ retval = exp_val * cos_ix;
  72           if (ix_cond & __M81_COND_ZERO)
  73             __imag__ retval = __imag__ x;
  74           else
  75             __imag__ retval = exp_val * sin_ix;
  76         }
  77       else
  78         {
  79           /* Compute the sign of the result.  */
  80           long double remainder, pi_2;
  81           int quadrant;
  82
  83           if ((rx_cond & (__M81_COND_NAN|__M81_COND_NEG)) == __M81_COND_NEG)
  84             __real__ retval = __imag__ retval = 0.0;
  85           else
  86             __real__ retval = __imag__ retval = __real__ x;
  87           __asm ("fmovecr %#0,%0\n\tfscale%.w %#-1,%0" : "=f" (pi_2));
  88           __asm ("fmod%.x %2,%0\n\tfmove%.l %/fpsr,%1"
  89                  : "=f" (remainder), "=dm" (quadrant)
  90                  : "f" (pi_2), "0" (__imag__ x));
  91           quadrant = (quadrant >> 16) & 0x83;
  92           if (quadrant & 0x80)
  93             quadrant ^= 0x83;
  94           switch (quadrant)
  95             {
  96             default:
  97               break;
  98             case 1:
  99               __real__ retval = -__real__ retval;
 100               break;
 101             case 2:
 102               __real__ retval = -__real__ retval;
 103             case 3:
 104               __imag__ retval = -__imag__ retval;
 105               break;
 106             }
 107           if (ix_cond & __M81_COND_ZERO && (rx_cond & __M81_COND_NAN) == 0)
 108             __imag__ retval = __imag__ x;
 109         }
 110     }
 111   else
 112     {
 113       unsigned long rx_cond = __m81_test (__real__ x);
 114
 115       if (rx_cond & __M81_COND_INF)
 116         {
 117           /* Real part is infinite.  */
 118           if (rx_cond & __M81_COND_NEG)
 119             {
 120               __real__ retval = __imag__ retval = 0.0;
 121               if (ix_cond & __M81_COND_NEG)
 122                 __imag__ retval = -__imag__ retval;
 123             }
 124           else
 125             {
 126               __real__ retval = __real__ x;
 127               __imag__ retval = __imag__ x - __imag__ x;
 128             }
 129         }
 130       else
 131         __real__ retval = __imag__ retval = __imag__ x - __imag__ x;
 132     }
 133
 134   return retval;
 135 }
 136 weak_alias (s(__cexp), s(cexp))