sysdeps/m68k/m680x0/fpu/s_cexp_template.c

   1 /* Complex exponential function.  m68k fpu version
   2    Copyright (C) 1997-2019 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 #define CONCATX(a,b) __CONCAT(a,b)
  26 #define s(name) M_SUF (name)
  27 #define m81(func) __m81_u(s(func))
  28
  29 CFLOAT
  30 s(__cexp) (CFLOAT x)
  31 {
  32   CFLOAT retval;
  33   unsigned long ix_cond;
  34
  35   ix_cond = __m81_test (__imag__ x);
  36
  37   if ((ix_cond & (__M81_COND_NAN|__M81_COND_INF)) == 0)
  38     {
  39       /* Imaginary part is finite.  */
  40       unsigned long rx_cond = __m81_test (__real__ x);
  41
  42       if ((rx_cond & (__M81_COND_NAN|__M81_COND_INF)) == 0)
  43         {
  44           const int t = (int) ((LDBL_MAX_EXP - 1) * M_LN2l);
  45           long double sin_ix, cos_ix, exp_val;
  46
  47           __m81_u (__sincosl) (__imag__ x, &sin_ix, &cos_ix);
  48
  49           if (__real__ x > t)
  50             {
  51               long double exp_t = __m81_u(__ieee754_expl) (t);
  52               __real__ x -= t;
  53               sin_ix *= exp_t;
  54               cos_ix *= exp_t;
  55               if (__real__ x > t)
  56                 {
  57                   __real__ x -= t;
  58                   sin_ix *= exp_t;
  59                   cos_ix *= exp_t;
  60                 }
  61             }
  62
  63           exp_val = __m81_u(__ieee754_expl) (__real__ x);
  64           __real__ retval = exp_val * cos_ix;
  65           if (ix_cond & __M81_COND_ZERO)
  66             __imag__ retval = __imag__ x;
  67           else
  68             __imag__ retval = exp_val * sin_ix;
  69         }
  70       else
  71         {
  72           /* Compute the sign of the result.  */
  73           long double remainder, pi_2;
  74           int quadrant;
  75
  76           if ((rx_cond & (__M81_COND_NAN|__M81_COND_NEG)) == __M81_COND_NEG)
  77             __real__ retval = __imag__ retval = 0.0;
  78           else
  79             __real__ retval = __imag__ retval = __real__ x;
  80           __asm ("fmovecr %#0,%0\n\tfscale%.w %#-1,%0" : "=f" (pi_2));
  81           __asm ("fmod%.x %2,%0\n\tfmove%.l %/fpsr,%1"
  82                  : "=f" (remainder), "=dm" (quadrant)
  83                  : "f" (pi_2), "0" (__imag__ x));
  84           quadrant = (quadrant >> 16) & 0x83;
  85           if (quadrant & 0x80)
  86             quadrant ^= 0x83;
  87           switch (quadrant)
  88             {
  89             default:
  90               break;
  91             case 1:
  92               __real__ retval = -__real__ retval;
  93               break;
  94             case 2:
  95               __real__ retval = -__real__ retval;
  96             case 3:
  97               __imag__ retval = -__imag__ retval;
  98               break;
  99             }
 100           if (ix_cond & __M81_COND_ZERO && (rx_cond & __M81_COND_NAN) == 0)
 101             __imag__ retval = __imag__ x;
 102         }
 103     }
 104   else
 105     {
 106       unsigned long rx_cond = __m81_test (__real__ x);
 107
 108       if (rx_cond & __M81_COND_INF)
 109         {
 110           /* Real part is infinite.  */
 111           if (rx_cond & __M81_COND_NEG)
 112             {
 113               __real__ retval = __imag__ retval = 0.0;
 114               if (ix_cond & __M81_COND_NEG)
 115                 __imag__ retval = -__imag__ retval;
 116             }
 117           else
 118             {
 119               __real__ retval = __real__ x;
 120               __imag__ retval = __imag__ x - __imag__ x;
 121             }
 122         }
 123       else
 124         __real__ retval = __imag__ retval = __imag__ x - __imag__ x;
 125     }
 126
 127   return retval;
 128 }
 129 declare_mgen_alias (__cexp, cexp)