1 /* $NetBSD: csqrt.c,v 1.1 2007/08/20 16:01:37 drochner Exp $ */
2
3 /*-
4 * Copyright (c) 2007 The NetBSD Foundation, Inc.
5 * All rights reserved.
6 *
7 * This code is derived from software written by Stephen L. Moshier.
8 * It is redistributed by the NetBSD Foundation by permission of the author.
9 *
10 * Redistribution and use in source and binary forms, with or without
11 * modification, are permitted provided that the following conditions
12 * are met:
13 * 1. Redistributions of source code must retain the above copyright
14 * notice, this list of conditions and the following disclaimer.
15 * 2. Redistributions in binary form must reproduce the above copyright
16 * notice, this list of conditions and the following disclaimer in the
17 * documentation and/or other materials provided with the distribution.
18 *
19 * THIS SOFTWARE IS PROVIDED BY THE NETBSD FOUNDATION, INC. AND CONTRIBUTORS
20 * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
21 * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
22 * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE FOUNDATION OR CONTRIBUTORS
23 * BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
24 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
25 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
26 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
27 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
28 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
29 * POSSIBILITY OF SUCH DAMAGE.
30 *
31 * imported and modified include for newlib 2010/10/03
32 * Marco Atzeri <marco_atzeri@yahoo.it>
33 */
34
35 /*
36 FUNCTION
37 <<csqrt>>, <<csqrtf>>---complex square root
38
39 INDEX
40 csqrt
41 INDEX
42 csqrtf
43
44 SYNOPSIS
45 #include <complex.h>
46 double complex csqrt(double complex <[z]>);
47 float complex csqrtf(float complex <[z]>);
48
49
50 DESCRIPTION
51 These functions compute the complex square root of <[z]>, with
52 a branch cut along the negative real axis.
53
54 <<csqrtf>> is identical to <<csqrt>>, except that it performs
55 its calculations on <<floats complex>>.
56
57 RETURNS
58 The csqrt functions return the complex square root value, in
59 the range of the right halfplane (including the imaginary axis).
60
61 PORTABILITY
62 <<csqrt>> and <<csqrtf>> are ISO C99
63
64 QUICKREF
65 <<csqrt>> and <<csqrtf>> are ISO C99
66
67 */
68
69
70 #include <complex.h>
71 #include <math.h>
72
73 double complex
csqrt(double complex z)74 csqrt(double complex z)
75 {
76 double complex w;
77 double x, y, r, t, scale;
78
79 x = creal (z);
80 y = cimag (z);
81
82 if (y == 0.0) {
83 if (x == 0.0) {
84 w = 0.0 + y * (double complex) I;
85 } else {
86 r = fabs(x);
87 r = sqrt(r);
88 if (x < 0.0) {
89 w = 0.0 + r * (double complex) I;
90 } else {
91 w = r + y * (double complex) I;
92 }
93 }
94 return w;
95 }
96 if (x == 0.0) {
97 r = fabs(y);
98 r = sqrt(0.5 * r);
99 if (y > 0)
100 w = r + r * (double complex) I;
101 else
102 w = r - r * (double complex) I;
103 return w;
104 }
105 /* Rescale to avoid internal overflow or underflow. */
106 if ((fabs(x) > 4.0) || (fabs(y) > 4.0)) {
107 x *= 0.25;
108 y *= 0.25;
109 scale = 2.0;
110 } else {
111 #if 1
112 x *= 1.8014398509481984e16; /* 2^54 */
113 y *= 1.8014398509481984e16;
114 scale = 7.450580596923828125e-9; /* 2^-27 */
115 #else
116 x *= 4.0;
117 y *= 4.0;
118 scale = 0.5;
119 #endif
120 }
121 w = x + y * (double complex) I;
122 r = cabs(w);
123 if (x > 0) {
124 t = sqrt(0.5 * r + 0.5 * x);
125 r = scale * fabs((0.5 * y) / t );
126 t *= scale;
127 } else {
128 r = sqrt(0.5 * r - 0.5 * x);
129 t = scale * fabs((0.5 * y) / r);
130 r *= scale;
131 }
132 if (y < 0)
133 w = t - r * (double complex) I;
134 else
135 w = t + r * (double complex) I;
136 return w;
137 }
138
