1  /* Complex exponential function.  m68k fpu version
       2     Copyright (C) 1997-2023 Free Software Foundation, Inc.
       3     This file is part of the GNU C Library.
       4  
       5     The GNU C Library is free software; you can redistribute it and/or
       6     modify it under the terms of the GNU Lesser General Public
       7     License as published by the Free Software Foundation; either
       8     version 2.1 of the License, or (at your option) any later version.
       9  
      10     The GNU C Library is distributed in the hope that it will be useful,
      11     but WITHOUT ANY WARRANTY; without even the implied warranty of
      12     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
      13     Lesser General Public License for more details.
      14  
      15     You should have received a copy of the GNU Lesser General Public
      16     License along with the GNU C Library.  If not, see
      17     <https://www.gnu.org/licenses/>.  */
      18  
      19  #include <float.h>
      20  #include <complex.h>
      21  #include <math.h>
      22  #include "mathimpl.h"
      23  
      24  #define CONCATX(a,b) __CONCAT(a,b)
      25  #define s(name) M_SUF (name)
      26  #define m81(func) __m81_u(s(func))
      27  
      28  CFLOAT
      29  s(__cexp) (CFLOAT x)
      30  {
      31    CFLOAT retval;
      32    unsigned long ix_cond;
      33  
      34    ix_cond = __m81_test (__imag__ x);
      35  
      36    if ((ix_cond & (__M81_COND_NAN|__M81_COND_INF)) == 0)
      37      {
      38        /* Imaginary part is finite.  */
      39        unsigned long rx_cond = __m81_test (__real__ x);
      40  
      41        if ((rx_cond & (__M81_COND_NAN|__M81_COND_INF)) == 0)
      42  	{
      43  	  const int t = (int) ((LDBL_MAX_EXP - 1) * M_LN2l);
      44  	  long double sin_ix, cos_ix, exp_val;
      45  
      46  	  __m81_u (__sincosl) (__imag__ x, &sin_ix, &cos_ix);
      47  
      48  	  if (__real__ x > t)
      49  	    {
      50  	      long double exp_t = __m81_u(__ieee754_expl) (t);
      51  	      __real__ x -= t;
      52  	      sin_ix *= exp_t;
      53  	      cos_ix *= exp_t;
      54  	      if (__real__ x > t)
      55  		{
      56  		  __real__ x -= t;
      57  		  sin_ix *= exp_t;
      58  		  cos_ix *= exp_t;
      59  		}
      60  	    }
      61  
      62  	  exp_val = __m81_u(__ieee754_expl) (__real__ x);
      63  	  __real__ retval = exp_val * cos_ix;
      64  	  if (ix_cond & __M81_COND_ZERO)
      65  	    __imag__ retval = __imag__ x;
      66  	  else
      67  	    __imag__ retval = exp_val * sin_ix;
      68  	}
      69        else
      70  	{
      71  	  /* Compute the sign of the result.  */
      72  	  long double remainder, pi_2;
      73  	  int quadrant;
      74  
      75  	  if ((rx_cond & (__M81_COND_NAN|__M81_COND_NEG)) == __M81_COND_NEG)
      76  	    __real__ retval = __imag__ retval = 0.0;
      77  	  else
      78  	    __real__ retval = __imag__ retval = __real__ x;
      79  	  __asm ("fmovecr %#0,%0\n\tfscale%.w %#-1,%0" : "=f" (pi_2));
      80  	  __asm ("fmod%.x %2,%0\n\tfmove%.l %/fpsr,%1"
      81  		 : "=f" (remainder), "=dm" (quadrant)
      82  		 : "f" (pi_2), "0" (__imag__ x));
      83  	  quadrant = (quadrant >> 16) & 0x83;
      84  	  if (quadrant & 0x80)
      85  	    quadrant ^= 0x83;
      86  	  switch (quadrant)
      87  	    {
      88  	    default:
      89  	      break;
      90  	    case 1:
      91  	      __real__ retval = -__real__ retval;
      92  	      break;
      93  	    case 2:
      94  	      __real__ retval = -__real__ retval;
      95  	      /* Fall through.  */
      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)