GCC Code Coverage Report
Directory: ./ Exec Total Coverage
File: lib/libm/src/s_csqrtf.c Lines: 0 35 0.0 %
Date: 2017-11-13 Branches: 0 14 0.0 %

Line Branch Exec Source
1
/*	$OpenBSD: s_csqrtf.c,v 1.4 2016/09/12 19:47:02 guenther Exp $	*/
2
/*
3
 * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
4
 *
5
 * Permission to use, copy, modify, and distribute this software for any
6
 * purpose with or without fee is hereby granted, provided that the above
7
 * copyright notice and this permission notice appear in all copies.
8
 *
9
 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
10
 * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
11
 * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
12
 * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
13
 * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
14
 * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
15
 * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
16
 */
17
18
/*							csqrtf()
19
 *
20
 *	Complex square root
21
 *
22
 *
23
 *
24
 * SYNOPSIS:
25
 *
26
 * float complex csqrtf();
27
 * float complex z, w;
28
 *
29
 * w = csqrtf( z );
30
 *
31
 *
32
 *
33
 * DESCRIPTION:
34
 *
35
 *
36
 * If z = x + iy,  r = |z|, then
37
 *
38
 *                       1/2
39
 * Re w  =  [ (r + x)/2 ]   ,
40
 *
41
 *                       1/2
42
 * Im w  =  [ (r - x)/2 ]   .
43
 *
44
 * Cancellation error in r-x or r+x is avoided by using the
45
 * identity  2 Re w Im w  =  y.
46
 *
47
 * Note that -w is also a square root of z.  The root chosen
48
 * is always in the right half plane and Im w has the same sign as y.
49
 *
50
 *
51
 *
52
 * ACCURACY:
53
 *
54
 *
55
 *                      Relative error:
56
 * arithmetic   domain     # trials      peak         rms
57
 *    IEEE      -10,+10    1,000,000    1.8e-7       3.5e-8
58
 *
59
 */
60
61
#include <complex.h>
62
#include <math.h>
63
64
float complex
65
csqrtf(float complex z)
66
{
67
	float complex w;
68
	float x, y, r, t, scale;
69
70
	x = crealf(z);
71
	y = cimagf(z);
72
73
	if(y == 0.0f) {
74
		if (x < 0.0f) {
75
			w = 0.0f + copysign(sqrtf(-x), y) * I;
76
			return (w);
77
		}
78
		else if (x == 0.0f) {
79
			return (0.0f + y * I);
80
		}
81
		else {
82
			w = sqrtf(x) + y * I;
83
			return (w);
84
		}
85
	}
86
87
	if (x == 0.0f) {
88
		r = fabsf(y);
89
		r = sqrtf(0.5f*r);
90
		if(y > 0)
91
			w = r + r * I;
92
		else
93
			w = r - r * I;
94
		return (w);
95
	}
96
97
	/* Rescale to avoid internal overflow or underflow.  */
98
	if ((fabsf(x) > 4.0f) || (fabsf(y) > 4.0f)) {
99
		x *= 0.25f;
100
		y *= 0.25f;
101
		scale = 2.0f;
102
	}
103
	else {
104
		x *= 6.7108864e7f; /* 2^26 */
105
		y *= 6.7108864e7f;
106
		scale = 1.220703125e-4f; /* 2^-13 */
107
#if 0
108
		x *= 4.0f;
109
		y *= 4.0f;
110
		scale = 0.5f;
111
#endif
112
	}
113
	w = x + y * I;
114
	r = cabsf(w);
115
	if (x > 0) {
116
		t = sqrtf( 0.5f * r + 0.5f * x );
117
		r = scale * fabsf((0.5f * y) / t);
118
		t *= scale;
119
	}
120
	else {
121
		r = sqrtf(0.5f * r - 0.5f * x);
122
		t = scale * fabsf((0.5f * y) / r);
123
		r *= scale;
124
	}
125
126
	if (y < 0)
127
		w = t - r * I;
128
	else
129
		w = t + r * I;
130
	return (w);
131
}
132
DEF_STD(csqrtf);