Branch data Line data Source code
1 : : /* $OpenBSD: s_catan.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 : : /* catan()
19 : : *
20 : : * Complex circular arc tangent
21 : : *
22 : : *
23 : : *
24 : : * SYNOPSIS:
25 : : *
26 : : * double complex catan();
27 : : * double complex z, w;
28 : : *
29 : : * w = catan (z);
30 : : *
31 : : *
32 : : *
33 : : * DESCRIPTION:
34 : : *
35 : : * If
36 : : * z = x + iy,
37 : : *
38 : : * then
39 : : * 1 ( 2x )
40 : : * Re w = - arctan(-----------) + k PI
41 : : * 2 ( 2 2)
42 : : * (1 - x - y )
43 : : *
44 : : * ( 2 2)
45 : : * 1 (x + (y+1) )
46 : : * Im w = - log(------------)
47 : : * 4 ( 2 2)
48 : : * (x + (y-1) )
49 : : *
50 : : * Where k is an arbitrary integer.
51 : : *
52 : : * catan(z) = -i catanh(iz).
53 : : *
54 : : * ACCURACY:
55 : : *
56 : : * Relative error:
57 : : * arithmetic domain # trials peak rms
58 : : * DEC -10,+10 5900 1.3e-16 7.8e-18
59 : : * IEEE -10,+10 30000 2.3e-15 8.5e-17
60 : : * The check catan( ctan(z) ) = z, with |x| and |y| < PI/2,
61 : : * had peak relative error 1.5e-16, rms relative error
62 : : * 2.9e-17. See also clog().
63 : : */
64 : :
65 : : #include <float.h>
66 : : #include <openlibm_complex.h>
67 : : #include <openlibm_math.h>
68 : :
69 : : #include "math_private.h"
70 : :
71 : : #define MAXNUM 1.0e308
72 : :
73 : : static const double DP1 = 3.14159265160560607910E0;
74 : : static const double DP2 = 1.98418714791870343106E-9;
75 : : static const double DP3 = 1.14423774522196636802E-17;
76 : :
77 : : static double
78 : 0 : _redupi(double x)
79 : : {
80 : : double t;
81 : : long i;
82 : :
83 : 0 : t = x/M_PI;
84 [ # # ]: 0 : if(t >= 0.0)
85 : 0 : t += 0.5;
86 : : else
87 : 0 : t -= 0.5;
88 : :
89 : 0 : i = t; /* the multiple */
90 : 0 : t = i;
91 : 0 : t = ((x - t * DP1) - t * DP2) - t * DP3;
92 : 0 : return (t);
93 : : }
94 : :
95 : : double complex
96 : 0 : catan(double complex z)
97 : : {
98 : : double complex w;
99 : : double a, t, x, x2, y;
100 : :
101 : 0 : x = creal (z);
102 : 0 : y = cimag (z);
103 : :
104 [ # # # # ]: 0 : if ((x == 0.0) && (y > 1.0))
105 : 0 : goto ovrf;
106 : :
107 : 0 : x2 = x * x;
108 : 0 : a = 1.0 - x2 - (y * y);
109 [ # # ]: 0 : if (a == 0.0)
110 : 0 : goto ovrf;
111 : :
112 : 0 : t = 0.5 * atan2 (2.0 * x, a);
113 : 0 : w = _redupi (t);
114 : :
115 : 0 : t = y - 1.0;
116 : 0 : a = x2 + (t * t);
117 [ # # ]: 0 : if (a == 0.0)
118 : 0 : goto ovrf;
119 : :
120 : 0 : t = y + 1.0;
121 : 0 : a = (x2 + (t * t))/a;
122 : 0 : w = w + (0.25 * log (a)) * I;
123 : 0 : return (w);
124 : :
125 : 0 : ovrf:
126 : : /*mtherr ("catan", OVERFLOW);*/
127 : 0 : w = MAXNUM + MAXNUM * I;
128 : 0 : return (w);
129 : : }
130 : :
131 : : #if LDBL_MANT_DIG == DBL_MANT_DIG
132 : : openlibm_strong_reference(catan, catanl);
133 : : #endif /* LDBL_MANT_DIG == DBL_MANT_DIG */
|
