1/* $NetBSD$ */
2
3/*-
4 * Copyright (c) 2007 The NetBSD Foundation, Inc.
5 * All rights reserved.
6 *
7 * This code is derived from software written by Stephen L. Moshier.
8 * It is redistributed by the NetBSD Foundation by permission of the author.
9 *
10 * Redistribution and use in source and binary forms, with or without
11 * modification, are permitted provided that the following conditions
12 * are met:
13 * 1. Redistributions of source code must retain the above copyright
14 *    notice, this list of conditions and the following disclaimer.
15 * 2. Redistributions in binary form must reproduce the above copyright
16 *    notice, this list of conditions and the following disclaimer in the
17 *    documentation and/or other materials provided with the distribution.
18 *
19 * THIS SOFTWARE IS PROVIDED BY THE NETBSD FOUNDATION, INC. AND CONTRIBUTORS
20 * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
21 * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
22 * PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE FOUNDATION OR CONTRIBUTORS
23 * BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
24 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
25 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
26 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
27 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
28 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
29 * POSSIBILITY OF SUCH DAMAGE.
30 */
31
32#include "../src/namespace.h"
33#include <complex.h>
34#include <math.h>
35
36#ifdef __weak_alias
37__weak_alias(casin, _casin)
38#endif
39
40double complex
41casin(double complex z)
42{
43	double complex w;
44	double complex ca, ct, zz, z2;
45	double x, y;
46
47	x = creal(z);
48	y = cimag(z);
49
50#if 0 /* MD: test is incorrect, casin(>1) is defined */
51	if (y == 0.0) {
52		if (fabs(x) > 1.0) {
53			w = M_PI_2 + 0.0 * I;
54#if 0
55			mtherr ("casin", DOMAIN);
56#endif
57		} else {
58			w = asin(x) + 0.0 * I;
59		}
60		return w;
61	}
62#endif
63
64/* Power series expansion */
65/*
66b = cabs(z);
67if( b < 0.125 )
68{
69z2.r = (x - y) * (x + y);
70z2.i = 2.0 * x * y;
71
72cn = 1.0;
73n = 1.0;
74ca.r = x;
75ca.i = y;
76sum.r = x;
77sum.i = y;
78do
79	{
80	ct.r = z2.r * ca.r  -  z2.i * ca.i;
81	ct.i = z2.r * ca.i  +  z2.i * ca.r;
82	ca.r = ct.r;
83	ca.i = ct.i;
84
85	cn *= n;
86	n += 1.0;
87	cn /= n;
88	n += 1.0;
89	b = cn/n;
90
91	ct.r *= b;
92	ct.i *= b;
93	sum.r += ct.r;
94	sum.i += ct.i;
95	b = fabs(ct.r) + fabs(ct.i);
96	}
97while( b > MACHEP );
98w->r = sum.r;
99w->i = sum.i;
100return;
101}
102*/
103
104
105	ca = x + y * I;
106	ct = ca * I;
107	/* sqrt( 1 - z*z) */
108	/* cmul( &ca, &ca, &zz ) */
109	/*x * x  -  y * y */
110	zz = (x - y) * (x + y) + (2.0 * x * y) * I;
111
112	zz = 1.0 - creal(zz) - cimag(zz) * I;
113	z2 = csqrt(zz);
114
115	zz = ct + z2;
116	zz = clog(zz);
117	/* multiply by 1/i = -i */
118	w = zz * (-1.0 * I);
119	return w;
120}
121