libquadmath/math/cexpq.c

   1 /* Return value of complex exponential function for a float type.
   2    Copyright (C) 1997-2018 Free Software Foundation, Inc.
   3    This file is part of the GNU C Library.
   4    Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
   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 "quadmath-imp.h"
  21
  22 __complex128
  23 cexpq (__complex128 x)
  24 {
  25   __complex128 retval;
  26   int rcls = fpclassifyq (__real__ x);
  27   int icls = fpclassifyq (__imag__ x);
  28
  29   if (__glibc_likely (rcls >= QUADFP_ZERO))
  30     {
  31       /* Real part is finite.  */
  32       if (__glibc_likely (icls >= QUADFP_ZERO))
  33         {
  34           /* Imaginary part is finite.  */
  35           const int t = (int) ((FLT128_MAX_EXP - 1) * M_LN2q);
  36           __float128 sinix, cosix;
  37
  38           if (__glibc_likely (fabsq (__imag__ x) > FLT128_MIN))
  39             {
  40               sincosq (__imag__ x, &sinix, &cosix);
  41             }
  42           else
  43             {
  44               sinix = __imag__ x;
  45               cosix = 1;
  46             }
  47
  48           if (__real__ x > t)
  49             {
  50               __float128 exp_t = expq (t);
  51               __real__ x -= t;
  52               sinix *= exp_t;
  53               cosix *= exp_t;
  54               if (__real__ x > t)
  55                 {
  56                   __real__ x -= t;
  57                   sinix *= exp_t;
  58                   cosix *= exp_t;
  59                 }
  60             }
  61           if (__real__ x > t)
  62             {
  63               /* Overflow (original real part of x > 3t).  */
  64               __real__ retval = FLT128_MAX * cosix;
  65               __imag__ retval = FLT128_MAX * sinix;
  66             }
  67           else
  68             {
  69               __float128 exp_val = expq (__real__ x);
  70               __real__ retval = exp_val * cosix;
  71               __imag__ retval = exp_val * sinix;
  72             }
  73           math_check_force_underflow_complex (retval);
  74         }
  75       else
  76         {
  77           /* If the imaginary part is +-inf or NaN and the real part
  78              is not +-inf the result is NaN + iNaN.  */
  79           __real__ retval = nanq ("");
  80           __imag__ retval = nanq ("");
  81
  82           feraiseexcept (FE_INVALID);
  83         }
  84     }
  85   else if (__glibc_likely (rcls == QUADFP_INFINITE))
  86     {
  87       /* Real part is infinite.  */
  88       if (__glibc_likely (icls >= QUADFP_ZERO))
  89         {
  90           /* Imaginary part is finite.  */
  91           __float128 value = signbitq (__real__ x) ? 0 : HUGE_VALQ;
  92
  93           if (icls == QUADFP_ZERO)
  94             {
  95               /* Imaginary part is 0.0.  */
  96               __real__ retval = value;
  97               __imag__ retval = __imag__ x;
  98             }
  99           else
 100             {
 101               __float128 sinix, cosix;
 102
 103               if (__glibc_likely (fabsq (__imag__ x) > FLT128_MIN))
 104                 {
 105                   sincosq (__imag__ x, &sinix, &cosix);
 106                 }
 107               else
 108                 {
 109                   sinix = __imag__ x;
 110                   cosix = 1;
 111                 }
 112
 113               __real__ retval = copysignq (value, cosix);
 114               __imag__ retval = copysignq (value, sinix);
 115             }
 116         }
 117       else if (signbitq (__real__ x) == 0)
 118         {
 119           __real__ retval = HUGE_VALQ;
 120           __imag__ retval = __imag__ x - __imag__ x;
 121         }
 122       else
 123         {
 124           __real__ retval = 0;
 125           __imag__ retval = copysignq (0, __imag__ x);
 126         }
 127     }
 128   else
 129     {
 130       /* If the real part is NaN the result is NaN + iNaN unless the
 131          imaginary part is zero.  */
 132       __real__ retval = nanq ("");
 133       if (icls == QUADFP_ZERO)
 134         __imag__ retval = __imag__ x;
 135       else
 136         {
 137           __imag__ retval = nanq ("");
 138
 139           if (rcls != QUADFP_NAN || icls != QUADFP_NAN)
 140             feraiseexcept (FE_INVALID);
 141         }
 142     }
 143
 144   return retval;
 145 }