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
57static const long double MACHEPL= 5.42101086242752217003726400434970855712890625E-20L;
58#elif	LDBL_MANT_DIG == 113
59static const long double MACHEPL = 9.629649721936179265279889712924636592690508e-35L;
60#endif
61
62static const long double PIO2L = 1.570796326794896619231321691639751442098585L;
63
64long double complex
65casinl(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}
133DEF_STD(casinl);
134