1  /* Complex square root of a float type.
       2     Copyright (C) 1997-2018 Free Software Foundation, Inc.
       3     This file is part of the GNU C Library.
       4     Based on an algorithm by Stephen L. Moshier <moshier@world.std.com>.
       5     Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
       6  
       7     The GNU C Library is free software; you can redistribute it and/or
       8     modify it under the terms of the GNU Lesser General Public
       9     License as published by the Free Software Foundation; either
      10     version 2.1 of the License, or (at your option) any later version.
      11  
      12     The GNU C Library is distributed in the hope that it will be useful,
      13     but WITHOUT ANY WARRANTY; without even the implied warranty of
      14     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
      15     Lesser General Public License for more details.
      16  
      17     You should have received a copy of the GNU Lesser General Public
      18     License along with the GNU C Library; if not, see
      19     <http://www.gnu.org/licenses/>.  */
      20  
      21  #include "quadmath-imp.h"
      22  
      23  __complex128
      24  csqrtq (__complex128 x)
      25  {
      26    __complex128 res;
      27    int rcls = fpclassifyq (__real__ x);
      28    int icls = fpclassifyq (__imag__ x);
      29  
      30    if (__glibc_unlikely (rcls <= QUADFP_INFINITE || icls <= QUADFP_INFINITE))
      31      {
      32        if (icls == QUADFP_INFINITE)
      33  	{
      34  	  __real__ res = HUGE_VALQ;
      35  	  __imag__ res = __imag__ x;
      36  	}
      37        else if (rcls == QUADFP_INFINITE)
      38  	{
      39  	  if (__real__ x < 0)
      40  	    {
      41  	      __real__ res = icls == QUADFP_NAN ? nanq ("") : 0;
      42  	      __imag__ res = copysignq (HUGE_VALQ, __imag__ x);
      43  	    }
      44  	  else
      45  	    {
      46  	      __real__ res = __real__ x;
      47  	      __imag__ res = (icls == QUADFP_NAN
      48  			      ? nanq ("") : copysignq (0, __imag__ x));
      49  	    }
      50  	}
      51        else
      52  	{
      53  	  __real__ res = nanq ("");
      54  	  __imag__ res = nanq ("");
      55  	}
      56      }
      57    else
      58      {
      59        if (__glibc_unlikely (icls == QUADFP_ZERO))
      60  	{
      61  	  if (__real__ x < 0)
      62  	    {
      63  	      __real__ res = 0;
      64  	      __imag__ res = copysignq (sqrtq (-__real__ x), __imag__ x);
      65  	    }
      66  	  else
      67  	    {
      68  	      __real__ res = fabsq (sqrtq (__real__ x));
      69  	      __imag__ res = copysignq (0, __imag__ x);
      70  	    }
      71  	}
      72        else if (__glibc_unlikely (rcls == QUADFP_ZERO))
      73  	{
      74  	  __float128 r;
      75  	  if (fabsq (__imag__ x) >= 2 * FLT128_MIN)
      76  	    r = sqrtq (0.5Q * fabsq (__imag__ x));
      77  	  else
      78  	    r = 0.5Q * sqrtq (2 * fabsq (__imag__ x));
      79  
      80  	  __real__ res = r;
      81  	  __imag__ res = copysignq (r, __imag__ x);
      82  	}
      83        else
      84  	{
      85  	  __float128 d, r, s;
      86  	  int scale = 0;
      87  
      88  	  if (fabsq (__real__ x) > FLT128_MAX / 4)
      89  	    {
      90  	      scale = 1;
      91  	      __real__ x = scalbnq (__real__ x, -2 * scale);
      92  	      __imag__ x = scalbnq (__imag__ x, -2 * scale);
      93  	    }
      94  	  else if (fabsq (__imag__ x) > FLT128_MAX / 4)
      95  	    {
      96  	      scale = 1;
      97  	      if (fabsq (__real__ x) >= 4 * FLT128_MIN)
      98  		__real__ x = scalbnq (__real__ x, -2 * scale);
      99  	      else
     100  		__real__ x = 0;
     101  	      __imag__ x = scalbnq (__imag__ x, -2 * scale);
     102  	    }
     103  	  else if (fabsq (__real__ x) < 2 * FLT128_MIN
     104  		   && fabsq (__imag__ x) < 2 * FLT128_MIN)
     105  	    {
     106  	      scale = -((FLT128_MANT_DIG + 1) / 2);
     107  	      __real__ x = scalbnq (__real__ x, -2 * scale);
     108  	      __imag__ x = scalbnq (__imag__ x, -2 * scale);
     109  	    }
     110  
     111  	  d = hypotq (__real__ x, __imag__ x);
     112  	  /* Use the identity   2  Re res  Im res = Im x
     113  	     to avoid cancellation error in  d +/- Re x.  */
     114  	  if (__real__ x > 0)
     115  	    {
     116  	      r = sqrtq (0.5Q * (d + __real__ x));
     117  	      if (scale == 1 && fabsq (__imag__ x) < 1)
     118  		{
     119  		  /* Avoid possible intermediate underflow.  */
     120  		  s = __imag__ x / r;
     121  		  r = scalbnq (r, scale);
     122  		  scale = 0;
     123  		}
     124  	      else
     125  		s = 0.5Q * (__imag__ x / r);
     126  	    }
     127  	  else
     128  	    {
     129  	      s = sqrtq (0.5Q * (d - __real__ x));
     130  	      if (scale == 1 && fabsq (__imag__ x) < 1)
     131  		{
     132  		  /* Avoid possible intermediate underflow.  */
     133  		  r = fabsq (__imag__ x / s);
     134  		  s = scalbnq (s, scale);
     135  		  scale = 0;
     136  		}
     137  	      else
     138  		r = fabsq (0.5Q * (__imag__ x / s));
     139  	    }
     140  
     141  	  if (scale)
     142  	    {
     143  	      r = scalbnq (r, scale);
     144  	      s = scalbnq (s, scale);
     145  	    }
     146  
     147  	  math_check_force_underflow (r);
     148  	  math_check_force_underflow (s);
     149  
     150  	  __real__ res = r;
     151  	  __imag__ res = copysignq (s, __imag__ x);
     152  	}
     153      }
     154  
     155    return res;
     156  }