1/* Complex square root of long double value.
2   Copyright (C) 1997, 1998 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, write to the Free
19   Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
20   02111-1307 USA.  */
21
22#include <complex.h>
23#include <math.h>
24
25#include "math_private.h"
26
27
28__complex__ long double
29__csqrtl (__complex__ long double x)
30{
31  __complex__ long double res;
32  int rcls = fpclassify (__real__ x);
33  int icls = fpclassify (__imag__ x);
34
35  if (rcls <= FP_INFINITE || icls <= FP_INFINITE)
36    {
37      if (icls == FP_INFINITE)
38	{
39	  __real__ res = HUGE_VALL;
40	  __imag__ res = __imag__ x;
41	}
42      else if (rcls == FP_INFINITE)
43	{
44	  if (__real__ x < 0.0)
45	    {
46	      __real__ res = icls == FP_NAN ? nanl ("") : 0;
47	      __imag__ res = copysignl (HUGE_VALL, __imag__ x);
48	    }
49	  else
50	    {
51	      __real__ res = __real__ x;
52	      __imag__ res = (icls == FP_NAN
53			      ? nanl ("") : copysignl (0.0, __imag__ x));
54	    }
55	}
56      else
57	{
58	  __real__ res = nanl ("");
59	  __imag__ res = nanl ("");
60	}
61    }
62  else
63    {
64      if (icls == FP_ZERO)
65	{
66	  if (__real__ x < 0.0)
67	    {
68	      __real__ res = 0.0;
69	      __imag__ res = copysignl (sqrtl (-__real__ x),
70					  __imag__ x);
71	    }
72	  else
73	    {
74	      __real__ res = fabsl (sqrtl (__real__ x));
75	      __imag__ res = copysignl (0.0, __imag__ x);
76	    }
77	}
78      else if (rcls == FP_ZERO)
79	{
80	  long double r = sqrtl (0.5 * fabsl (__imag__ x));
81
82	  __real__ res = copysignl (r, __imag__ x);
83	  __imag__ res = r;
84	}
85      else
86	{
87	  long double d, r, s;
88
89	  d = hypotl (__real__ x, __imag__ x);
90	  /* Use the identity   2  Re res  Im res = Im x
91	     to avoid cancellation error in  d +/- Re x.  */
92	  if (__real__ x > 0)
93	    {
94	      r = sqrtl (0.5L * d + 0.5L * __real__ x);
95	      s = (0.5L * __imag__ x) / r;
96	    }
97	  else
98	    {
99	      s = sqrtl (0.5L * d - 0.5L * __real__ x);
100	      r = fabsl ((0.5L * __imag__ x) / s);
101	    }
102
103	  __real__ res = r;
104	  __imag__ res = copysignl (s, __imag__ x);
105	}
106    }
107
108  return res;
109}
110weak_alias (__csqrtl, csqrtl)
111