1 /* $OpenBSD: s_ctan.c,v 1.1 2008/09/07 20:36:09 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 /* ctan() 19 * 20 * Complex circular tangent 21 * 22 * 23 * 24 * SYNOPSIS: 25 * 26 * double complex ctan(); 27 * double complex z, w; 28 * 29 * w = ctan (z); 30 * 31 * 32 * 33 * DESCRIPTION: 34 * 35 * If 36 * z = x + iy, 37 * 38 * then 39 * 40 * sin 2x + i sinh 2y 41 * w = --------------------. 42 * cos 2x + cosh 2y 43 * 44 * On the real axis the denominator is zero at odd multiples 45 * of PI/2. The denominator is evaluated by its Taylor 46 * series near these points. 47 * 48 * ctan(z) = -i ctanh(iz). 49 * 50 * ACCURACY: 51 * 52 * Relative error: 53 * arithmetic domain # trials peak rms 54 * DEC -10,+10 5200 7.1e-17 1.6e-17 55 * IEEE -10,+10 30000 7.2e-16 1.2e-16 56 * Also tested by ctan * ccot = 1 and catan(ctan(z)) = z. 57 */ 58 59 #include <complex.h> 60 #include <math.h> 61 62 #define MACHEP 1.1e-16 63 #define MAXNUM 1.0e308 64 65 static const double DP1 = 3.14159265160560607910E0; 66 static const double DP2 = 1.98418714791870343106E-9; 67 static const double DP3 = 1.14423774522196636802E-17; 68 69 static double 70 _redupi(double x) 71 { 72 double t; 73 long i; 74 75 t = x/M_PI; 76 if (t >= 0.0) 77 t += 0.5; 78 else 79 t -= 0.5; 80 81 i = t; /* the multiple */ 82 t = i; 83 t = ((x - t * DP1) - t * DP2) - t * DP3; 84 return (t); 85 } 86 87 /* Taylor series expansion for cosh(2y) - cos(2x) */ 88 89 static double 90 _ctans(double complex z) 91 { 92 double f, x, x2, y, y2, rn, t; 93 double d; 94 95 x = fabs (2.0 * creal (z)); 96 y = fabs (2.0 * cimag(z)); 97 98 x = _redupi(x); 99 100 x = x * x; 101 y = y * y; 102 x2 = 1.0; 103 y2 = 1.0; 104 f = 1.0; 105 rn = 0.0; 106 d = 0.0; 107 do { 108 rn += 1.0; 109 f *= rn; 110 rn += 1.0; 111 f *= rn; 112 x2 *= x; 113 y2 *= y; 114 t = y2 + x2; 115 t /= f; 116 d += t; 117 118 rn += 1.0; 119 f *= rn; 120 rn += 1.0; 121 f *= rn; 122 x2 *= x; 123 y2 *= y; 124 t = y2 - x2; 125 t /= f; 126 d += t; 127 } 128 while (fabs(t/d) > MACHEP) 129 ; 130 return (d); 131 } 132 133 double complex 134 ctan(double complex z) 135 { 136 double complex w; 137 double d; 138 139 d = cos (2.0 * creal (z)) + cosh (2.0 * cimag (z)); 140 141 if (fabs(d) < 0.25) 142 d = _ctans (z); 143 144 if (d == 0.0) { 145 /*mtherr ("ctan", OVERFLOW);*/ 146 w = MAXNUM + MAXNUM * I; 147 return (w); 148 } 149 150 w = sin (2.0 * creal(z)) / d + (sinh (2.0 * cimag(z)) / d) * I; 151 return (w); 152 } 153