]>
Commit | Line | Data |
---|---|---|
5ce24954 BK |
1 | /* Complex square root of double value. */ |
2 | ||
3 | /* Copyright (C) 1997-1999 Free Software Foundation, Inc. | |
4 | ||
5 | This file is part of the GNU ISO C++ Library. This library is free | |
6 | software; you can redistribute it and/or modify it under the | |
7 | terms of the GNU General Public License as published by the | |
8 | Free Software Foundation; either version 2, or (at your option) | |
9 | any later version. | |
10 | ||
11 | This 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 | |
14 | GNU General Public License for more details. | |
15 | ||
16 | You should have received a copy of the GNU General Public License along | |
17 | with this library; see the file COPYING. If not, write to the Free | |
18 | Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, | |
19 | USA. | |
20 | ||
21 | As a special exception, you may use this file as part of a free software | |
22 | library without restriction. Specifically, if other files instantiate | |
23 | templates or use macros or inline functions from this file, or you compile | |
24 | this file and link it with other files to produce an executable, this | |
25 | file does not by itself cause the resulting executable to be covered by | |
26 | the GNU General Public License. This exception does not however | |
27 | invalidate any other reasons why the executable file might be covered by | |
28 | the GNU General Public License. */ | |
29 | ||
30 | ||
31 | #include <math.h> | |
32 | #include "mathconf.h" | |
33 | ||
34 | ||
35 | __complex__ double | |
36 | csqrt (__complex__ double x) | |
37 | { | |
38 | __complex__ double res; | |
39 | ||
40 | if (!FINITE_P (__real__ x) || !FINITE_P (__imag__ x)) | |
41 | { | |
42 | if (INFINITE_P (__imag__ x)) | |
43 | { | |
44 | __real__ res = HUGE_VAL; | |
45 | __imag__ res = __imag__ x; | |
46 | } | |
47 | else if (INFINITE_P (__real__ x)) | |
48 | { | |
49 | if (__real__ x < 0.0) | |
50 | { | |
51 | __real__ res = __imag__ x != __imag__ x ? NAN : 0; | |
52 | __imag__ res = copysign (HUGE_VAL, __imag__ x); | |
53 | } | |
54 | else | |
55 | { | |
56 | __real__ res = __real__ x; | |
57 | __imag__ res = (__imag__ x != __imag__ x | |
58 | ? NAN : copysign (0.0, __imag__ x)); | |
59 | } | |
60 | } | |
61 | else | |
62 | { | |
63 | __real__ res = NAN; | |
64 | __imag__ res = NAN; | |
65 | } | |
66 | } | |
67 | else | |
68 | { | |
69 | if (__imag__ x == 0.0) | |
70 | { | |
71 | if (__real__ x < 0.0) | |
72 | { | |
73 | __real__ res = 0.0; | |
74 | __imag__ res = copysign (sqrt (-__real__ x), __imag__ x); | |
75 | } | |
76 | else | |
77 | { | |
78 | __real__ res = fabs (sqrt (__real__ x)); | |
79 | __imag__ res = copysign (0.0, __imag__ x); | |
80 | } | |
81 | } | |
82 | else if (__real__ x == 0.0) | |
83 | { | |
84 | double r = sqrt (0.5 * fabs (__imag__ x)); | |
85 | ||
86 | __real__ res = copysign (r, __imag__ x); | |
87 | __imag__ res = r; | |
88 | } | |
89 | else | |
90 | { | |
91 | __complex__ double q; | |
92 | double t, r; | |
93 | ||
94 | if (fabs (__imag__ x) < 2.0e-4 * fabs (__real__ x)) | |
95 | t = 0.25 * __imag__ x * (__imag__ x / __real__ x); | |
96 | else | |
97 | t = 0.5 * (hypot (__real__ x, __imag__ x) - __real__ x); | |
98 | ||
99 | r = sqrt (t); | |
100 | ||
101 | __real__ q = __imag__ x / (2.0 * r); | |
102 | __imag__ q = r; | |
103 | ||
104 | /* Heron iteration in complex arithmetic. */ | |
105 | res = 0.5 * (q + q / x); | |
106 | } | |
107 | } | |
108 | ||
109 | return res; | |
110 | } |