1 /*	$OpenBSD: s_casin.c,v 1.6 2013/07/03 04:46:36 espie 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 /*							casin()
19  *
20  *	Complex circular arc sine
21  *
22  *
23  *
24  * SYNOPSIS:
25  *
26  * double complex casin();
27  * double complex z, w;
28  *
29  * w = casin (z);
30  *
31  *
32  *
33  * DESCRIPTION:
34  *
35  * Inverse complex sine:
36  *
37  *                               2
38  * w = -i clog( iz + csqrt( 1 - z ) ).
39  *
40  * casin(z) = -i casinh(iz)
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 double complex
57 casin(double complex z)
58 {
59 	double complex w;
60 	static double complex ca, ct, zz, z2;
61 	double x, y;
62 
63 	x = creal (z);
64 	y = cimag (z);
65 
66 	if (y == 0.0) {
67 		if (fabs(x) > 1.0) {
68 			w = M_PI_2 + 0.0 * I;
69 			/*mtherr ("casin", DOMAIN);*/
70 		}
71 		else {
72 			w = asin (x) + 0.0 * I;
73 		}
74 		return (w);
75 	}
76 
77 	/* Power series expansion */
78 	/*
79 	b = cabs(z);
80 	if( b < 0.125 ) {
81 		z2.r = (x - y) * (x + y);
82 		z2.i = 2.0 * x * y;
83 
84 		cn = 1.0;
85 		n = 1.0;
86 		ca.r = x;
87 		ca.i = y;
88 		sum.r = x;
89 		sum.i = y;
90 		do {
91 			ct.r = z2.r * ca.r  -  z2.i * ca.i;
92 			ct.i = z2.r * ca.i  +  z2.i * ca.r;
93 			ca.r = ct.r;
94 			ca.i = ct.i;
95 
96 			cn *= n;
97 			n += 1.0;
98 			cn /= n;
99 			n += 1.0;
100 			b = cn/n;
101 
102 			ct.r *= b;
103 			ct.i *= b;
104 			sum.r += ct.r;
105 			sum.i += ct.i;
106 			b = fabs(ct.r) + fabs(ct.i);
107 		}
108 		while( b > MACHEP );
109 		w->r = sum.r;
110 		w->i = sum.i;
111 		return;
112 	}
113 	*/
114 
115 	ca = x + y * I;
116 	ct = ca * I;
117 	/* sqrt( 1 - z*z) */
118 	/* cmul( &ca, &ca, &zz ) */
119 	/*x * x  -  y * y */
120 	zz = (x - y) * (x + y) + (2.0 * x * y) * I;
121 
122 	zz = 1.0 - creal(zz) - cimag(zz) * I;
123 	z2 = csqrt (zz);
124 
125 	zz = ct + z2;
126 	zz = clog (zz);
127 	/* multiply by 1/i = -i */
128 	w = zz * (-1.0 * I);
129 	return (w);
130 }
131 
132 #if	LDBL_MANT_DIG == DBL_MANT_DIG
133 __strong_alias(casinl, casin);
134 #endif	/* LDBL_MANT_DIG == DBL_MANT_DIG */
135