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