1 /*	$OpenBSD: s_casinf.c,v 1.3 2011/07/20 19:28:33 martynas 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 /*							casinf()
19  *
20  *	Complex circular arc sine
21  *
22  *
23  *
24  * SYNOPSIS:
25  *
26  * void casinf();
27  * cmplxf z, w;
28  *
29  * casinf( &z, &w );
30  *
31  *
32  *
33  * DESCRIPTION:
34  *
35  * Inverse complex sine:
36  *
37  *                               2
38  * w = -i clog( iz + csqrt( 1 - z ) ).
39  *
40  *
41  * ACCURACY:
42  *
43  *                      Relative error:
44  * arithmetic   domain     # trials      peak         rms
45  *    IEEE      -10,+10     30000       1.1e-5      1.5e-6
46  * Larger relative error can be observed for z near zero.
47  *
48  */
49 
50 #include <complex.h>
51 #include <math.h>
52 
53 float complex
54 casinf(float complex z)
55 {
56 	float complex w;
57 	float x, y;
58 	static float complex ca, ct, zz, z2;
59 	/*
60 	float cn, n;
61 	static float a, b, s, t, u, v, y2;
62 	static cmplxf sum;
63 	*/
64 
65 	x = crealf(z);
66 	y = cimagf(z);
67 
68 	if(y == 0.0f) {
69 		if(fabsf(x) > 1.0f) {
70 			w = (float)M_PI_2 + 0.0f * I;
71 			/*mtherr( "casinf", DOMAIN );*/
72 		}
73 		else {
74 			w = asinf (x) + 0.0f * I;
75 		}
76 		return (w);
77 	}
78 
79 	/* Power series expansion */
80 	/*
81 	b = cabsf(z);
82 	if(b < 0.125) {
83 		z2.r = (x - y) * (x + y);
84 		z2.i = 2.0 * x * y;
85 
86 		cn = 1.0;
87 		n = 1.0;
88 		ca.r = x;
89 		ca.i = y;
90 		sum.r = x;
91 		sum.i = y;
92 		do {
93 			ct.r = z2.r * ca.r  -  z2.i * ca.i;
94 			ct.i = z2.r * ca.i  +  z2.i * ca.r;
95 			ca.r = ct.r;
96 			ca.i = ct.i;
97 
98 			cn *= n;
99 			n += 1.0;
100 			cn /= n;
101 			n += 1.0;
102 			b = cn/n;
103 
104 			ct.r *= b;
105 			ct.i *= b;
106 			sum.r += ct.r;
107 			sum.i += ct.i;
108 			b = fabsf(ct.r) + fabsf(ct.i);
109 		}
110 		while(b > MACHEPF);
111 		w->r = sum.r;
112 		w->i = sum.i;
113 		return;
114 	}
115 	*/
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.0f * x * y) * I;
124 	zz = 1.0f - crealf(zz) - cimagf(zz) * I;
125 	z2 = csqrtf (zz);
126 
127 	zz = ct + z2;
128 	zz = clogf (zz);
129 	/* multiply by 1/i = -i */
130 	w = zz * (-1.0f * I);
131 	return (w);
132 }
133