1 /* $OpenBSD: s_ctan.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 /* 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 <float.h> 61 #include <math.h> 62 63 #define MACHEP 1.1e-16 64 #define MAXNUM 1.0e308 65 66 static const double DP1 = 3.14159265160560607910E0; 67 static const double DP2 = 1.98418714791870343106E-9; 68 static const double DP3 = 1.14423774522196636802E-17; 69 70 static double 71 _redupi(double x) 72 { 73 double t; 74 long i; 75 76 t = x/M_PI; 77 if (t >= 0.0) 78 t += 0.5; 79 else 80 t -= 0.5; 81 82 i = t; /* the multiple */ 83 t = i; 84 t = ((x - t * DP1) - t * DP2) - t * DP3; 85 return (t); 86 } 87 88 /* Taylor series expansion for cosh(2y) - cos(2x) */ 89 90 static double 91 _ctans(double complex z) 92 { 93 double f, x, x2, y, y2, rn, t; 94 double d; 95 96 x = fabs (2.0 * creal (z)); 97 y = fabs (2.0 * cimag(z)); 98 99 x = _redupi(x); 100 101 x = x * x; 102 y = y * y; 103 x2 = 1.0; 104 y2 = 1.0; 105 f = 1.0; 106 rn = 0.0; 107 d = 0.0; 108 do { 109 rn += 1.0; 110 f *= rn; 111 rn += 1.0; 112 f *= rn; 113 x2 *= x; 114 y2 *= y; 115 t = y2 + x2; 116 t /= f; 117 d += t; 118 119 rn += 1.0; 120 f *= rn; 121 rn += 1.0; 122 f *= rn; 123 x2 *= x; 124 y2 *= y; 125 t = y2 - x2; 126 t /= f; 127 d += t; 128 } 129 while (fabs(t/d) > MACHEP) 130 ; 131 return (d); 132 } 133 134 double complex 135 ctan(double complex z) 136 { 137 double complex w; 138 double d; 139 140 d = cos (2.0 * creal (z)) + cosh (2.0 * cimag (z)); 141 142 if (fabs(d) < 0.25) 143 d = _ctans (z); 144 145 if (d == 0.0) { 146 /*mtherr ("ctan", OVERFLOW);*/ 147 w = MAXNUM + MAXNUM * I; 148 return (w); 149 } 150 151 w = sin (2.0 * creal(z)) / d + (sinh (2.0 * cimag(z)) / d) * I; 152 return (w); 153 } 154 155 #if LDBL_MANT_DIG == DBL_MANT_DIG 156 __strong_alias(ctanl, ctan); 157 #endif /* LDBL_MANT_DIG == DBL_MANT_DIG */ 158