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

Line Branch Exec Source
1
/*	$OpenBSD: s_casinl.c,v 1.5 2016/09/12 19:47:02 guenther Exp $	*/
2
3
/*
4
 * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
5
 *
6
 * Permission to use, copy, modify, and distribute this software for any
7
 * purpose with or without fee is hereby granted, provided that the above
8
 * copyright notice and this permission notice appear in all copies.
9
 *
10
 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
11
 * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
12
 * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
13
 * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
14
 * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
15
 * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
16
 * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
17
 */
18
19
/*							casinl()
20
 *
21
 *	Complex circular arc sine
22
 *
23
 *
24
 *
25
 * SYNOPSIS:
26
 *
27
 * long double complex casinl();
28
 * long double complex z, w;
29
 *
30
 * w = casinl( z );
31
 *
32
 *
33
 *
34
 * DESCRIPTION:
35
 *
36
 * Inverse complex sine:
37
 *
38
 *                               2
39
 * w = -i clog( iz + csqrt( 1 - z ) ).
40
 *
41
 *
42
 * ACCURACY:
43
 *
44
 *                      Relative error:
45
 * arithmetic   domain     # trials      peak         rms
46
 *    DEC       -10,+10     10100       2.1e-15     3.4e-16
47
 *    IEEE      -10,+10     30000       2.2e-14     2.7e-15
48
 * Larger relative error can be observed for z near zero.
49
 * Also tested by csin(casin(z)) = z.
50
 */
51
52
#include <complex.h>
53
#include <float.h>
54
#include <math.h>
55
56
#if	LDBL_MANT_DIG == 64
57
static const long double MACHEPL= 5.42101086242752217003726400434970855712890625E-20L;
58
#elif	LDBL_MANT_DIG == 113
59
static const long double MACHEPL = 9.629649721936179265279889712924636592690508e-35L;
60
#endif
61
62
static const long double PIO2L = 1.570796326794896619231321691639751442098585L;
63
64
long double complex
65
casinl(long double complex z)
66
{
67
	long double complex w;
68
	long double x, y, b;
69
	static long double complex ca, ct, zz, z2;
70
71
	x = creall(z);
72
	y = cimagl(z);
73
74
#if 0
75
	if (y == 0.0L) {
76
		if (fabsl(x) > 1.0L) {
77
			w = PIO2L + 0.0L * I;
78
			/*mtherr( "casinl", DOMAIN );*/
79
		}
80
		else {
81
			w = asinl(x) + 0.0L * I;
82
		}
83
		return (w);
84
	}
85
#endif
86
87
	/* Power series expansion */
88
	b = cabsl(z);
89
	if (b < 0.125L) {
90
		long double complex sum;
91
		long double n, cn;
92
93
		z2 = (x - y) * (x + y) + (2.0L * x * y) * I;
94
		cn = 1.0L;
95
		n = 1.0L;
96
		ca = x + y * I;
97
		sum = x + y * I;
98
		do {
99
			ct = z2 * ca;
100
			ca = ct;
101
102
			cn *= n;
103
			n += 1.0L;
104
			cn /= n;
105
			n += 1.0L;
106
			b = cn/n;
107
108
			ct *= b;
109
			sum += ct;
110
			b = cabsl(ct);
111
		}
112
113
		while (b > MACHEPL);
114
		w = sum;
115
		return w;
116
	}
117
118
	ca = x + y * I;
119
	ct = ca * I;	/* iz */
120
	/* sqrt(1 - z*z) */
121
	/* cmul(&ca, &ca, &zz) */
122
	/* x * x  -  y * y */
123
	zz = (x - y) * (x + y) + (2.0L * x * y) * I;
124
	zz = 1.0L - creall(zz) - cimagl(zz) * I;
125
	z2 = csqrtl(zz);
126
127
	zz = ct + z2;
128
	zz = clogl(zz);
129
	/* multiply by 1/i = -i */
130
	w = zz * (-1.0L * I);
131
	return (w);
132
}
133
DEF_STD(casinl);