libquadmath/math/csinq.c

   1 /* Complex sine function for float types.
   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 csinq (__complex128 x)
  24 {
  25   __complex128 retval;
  26   int negate = signbitq (__real__ x);
  27   int rcls = fpclassifyq (__real__ x);
  28   int icls = fpclassifyq (__imag__ x);
  29
  30   __real__ x = fabsq (__real__ x);
  31
  32   if (__glibc_likely (icls >= QUADFP_ZERO))
  33     {
  34       /* Imaginary part is finite.  */
  35       if (__glibc_likely (rcls >= QUADFP_ZERO))
  36         {
  37           /* Real part is finite.  */
  38           const int t = (int) ((FLT128_MAX_EXP - 1) * M_LN2q);
  39           __float128 sinix, cosix;
  40
  41           if (__glibc_likely (__real__ x > FLT128_MIN))
  42             {
  43               sincosq (__real__ x, &sinix, &cosix);
  44             }
  45           else
  46             {
  47               sinix = __real__ x;
  48               cosix = 1;
  49             }
  50
  51           if (negate)
  52             sinix = -sinix;
  53
  54           if (fabsq (__imag__ x) > t)
  55             {
  56               __float128 exp_t = expq (t);
  57               __float128 ix = fabsq (__imag__ x);
  58               if (signbitq (__imag__ x))
  59                 cosix = -cosix;
  60               ix -= t;
  61               sinix *= exp_t / 2;
  62               cosix *= exp_t / 2;
  63               if (ix > t)
  64                 {
  65                   ix -= t;
  66                   sinix *= exp_t;
  67                   cosix *= exp_t;
  68                 }
  69               if (ix > t)
  70                 {
  71                   /* Overflow (original imaginary part of x > 3t).  */
  72                   __real__ retval = FLT128_MAX * sinix;
  73                   __imag__ retval = FLT128_MAX * cosix;
  74                 }
  75               else
  76                 {
  77                   __float128 exp_val = expq (ix);
  78                   __real__ retval = exp_val * sinix;
  79                   __imag__ retval = exp_val * cosix;
  80                 }
  81             }
  82           else
  83             {
  84               __real__ retval = coshq (__imag__ x) * sinix;
  85               __imag__ retval = sinhq (__imag__ x) * cosix;
  86             }
  87
  88           math_check_force_underflow_complex (retval);
  89         }
  90       else
  91         {
  92           if (icls == QUADFP_ZERO)
  93             {
  94               /* Imaginary part is 0.0.  */
  95               __real__ retval = __real__ x - __real__ x;
  96               __imag__ retval = __imag__ x;
  97             }
  98           else
  99             {
 100               __real__ retval = nanq ("");
 101               __imag__ retval = nanq ("");
 102
 103               feraiseexcept (FE_INVALID);
 104             }
 105         }
 106     }
 107   else if (icls == QUADFP_INFINITE)
 108     {
 109       /* Imaginary part is infinite.  */
 110       if (rcls == QUADFP_ZERO)
 111         {
 112           /* Real part is 0.0.  */
 113           __real__ retval = copysignq (0, negate ? -1 : 1);
 114           __imag__ retval = __imag__ x;
 115         }
 116       else if (rcls > QUADFP_ZERO)
 117         {
 118           /* Real part is finite.  */
 119           __float128 sinix, cosix;
 120
 121           if (__glibc_likely (__real__ x > FLT128_MIN))
 122             {
 123               sincosq (__real__ x, &sinix, &cosix);
 124             }
 125           else
 126             {
 127               sinix = __real__ x;
 128               cosix = 1;
 129             }
 130
 131           __real__ retval = copysignq (HUGE_VALQ, sinix);
 132           __imag__ retval = copysignq (HUGE_VALQ, cosix);
 133
 134           if (negate)
 135             __real__ retval = -__real__ retval;
 136           if (signbitq (__imag__ x))
 137             __imag__ retval = -__imag__ retval;
 138         }
 139       else
 140         {
 141           __real__ retval = __real__ x - __real__ x;
 142           __imag__ retval = HUGE_VALQ;
 143         }
 144     }
 145   else
 146     {
 147       if (rcls == QUADFP_ZERO)
 148         __real__ retval = copysignq (0, negate ? -1 : 1);
 149       else
 150         __real__ retval = nanq ("");
 151       __imag__ retval = nanq ("");
 152     }
 153
 154   return retval;
 155 }