Branch data Line data Source code
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 <float.h>
53 : : #include <openlibm_complex.h>
54 : : #include <openlibm_math.h>
55 : :
56 : : #include "math_private.h"
57 : :
58 : : double complex
59 : 0 : casin(double complex z)
60 : : {
61 : : double complex w;
62 : : static double complex ca, ct, zz, z2;
63 : : double x, y;
64 : :
65 : 0 : x = creal (z);
66 : 0 : y = cimag (z);
67 : :
68 [ # # ]: 0 : if (y == 0.0) {
69 [ # # ]: 0 : if (fabs(x) > 1.0) {
70 : 0 : w = M_PI_2 + 0.0 * I;
71 : : /*mtherr ("casin", DOMAIN);*/
72 : : }
73 : : else {
74 : 0 : w = asin (x) + 0.0 * I;
75 : : }
76 : 0 : return (w);
77 : : }
78 : :
79 : : /* Power series expansion */
80 : : /*
81 : : b = cabs(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 = fabs(ct.r) + fabs(ct.i);
109 : : }
110 : : while( b > MACHEP );
111 : : w->r = sum.r;
112 : : w->i = sum.i;
113 : : return;
114 : : }
115 : : */
116 : :
117 : 0 : ca = x + y * I;
118 : 0 : ct = ca * I;
119 : : /* sqrt( 1 - z*z) */
120 : : /* cmul( &ca, &ca, &zz ) */
121 : : /*x * x - y * y */
122 : 0 : zz = (x - y) * (x + y) + (2.0 * x * y) * I;
123 : :
124 : 0 : zz = 1.0 - creal(zz) - cimag(zz) * I;
125 : 0 : z2 = csqrt (zz);
126 : :
127 : 0 : zz = ct + z2;
128 : 0 : zz = clog (zz);
129 : : /* multiply by 1/i = -i */
130 : 0 : w = zz * (-1.0 * I);
131 : 0 : return (w);
132 : : }
133 : :
134 : : #if LDBL_MANT_DIG == DBL_MANT_DIG
135 : : openlibm_strong_reference(casin, casinl);
136 : : #endif /* LDBL_MANT_DIG == DBL_MANT_DIG */
|
