1 /* $OpenBSD: s_ctanl.c,v 1.3 2011/07/20 21:02:51 martynas Exp $ */ 2 3 /* 4 * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net> 5 * 6 * Permission to use, copy, modify, and distribute this software for any 7 * purpose with or without fee is hereby granted, provided that the above 8 * copyright notice and this permission notice appear in all copies. 9 * 10 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES 11 * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF 12 * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR 13 * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES 14 * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN 15 * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF 16 * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. 17 */ 18 19 /* ctanl() 20 * 21 * Complex circular tangent 22 * 23 * 24 * 25 * SYNOPSIS: 26 * 27 * long double complex ctanl(); 28 * long double complex z, w; 29 * 30 * w = ctanl( z ); 31 * 32 * 33 * 34 * DESCRIPTION: 35 * 36 * If 37 * z = x + iy, 38 * 39 * then 40 * 41 * sin 2x + i sinh 2y 42 * w = --------------------. 43 * cos 2x + cosh 2y 44 * 45 * On the real axis the denominator is zero at odd multiples 46 * of PI/2. The denominator is evaluated by its Taylor 47 * series near these points. 48 * 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 #if LDBL_MANT_DIG == 64 64 static const long double MACHEPL= 5.42101086242752217003726400434970855712890625E-20L; 65 #elif LDBL_MANT_DIG == 113 66 static const long double MACHEPL = 9.629649721936179265279889712924636592690508e-35L; 67 #endif 68 69 static const long double PIL = 3.141592653589793238462643383279502884197169L; 70 static const long double DP1 = 3.14159265358979323829596852490908531763125L; 71 static const long double DP2 = 1.6667485837041756656403424829301998703007e-19L; 72 static const long double DP3 = 1.8830410776607851167459095484560349402753e-39L; 73 74 static long double 75 redupil(long double x) 76 { 77 long double t; 78 long i; 79 80 t = x / PIL; 81 if (t >= 0.0L) 82 t += 0.5L; 83 else 84 t -= 0.5L; 85 86 i = t; /* the multiple */ 87 t = i; 88 t = ((x - t * DP1) - t * DP2) - t * DP3; 89 return (t); 90 } 91 92 static long double 93 ctansl(long double complex z) 94 { 95 long double f, x, x2, y, y2, rn, t; 96 long double d; 97 98 x = fabsl(2.0L * creall(z)); 99 y = fabsl(2.0L * cimagl(z)); 100 101 x = redupil(x); 102 103 x = x * x; 104 y = y * y; 105 x2 = 1.0L; 106 y2 = 1.0L; 107 f = 1.0L; 108 rn = 0.0L; 109 d = 0.0L; 110 do { 111 rn += 1.0L; 112 f *= rn; 113 rn += 1.0L; 114 f *= rn; 115 x2 *= x; 116 y2 *= y; 117 t = y2 + x2; 118 t /= f; 119 d += t; 120 121 rn += 1.0L; 122 f *= rn; 123 rn += 1.0L; 124 f *= rn; 125 x2 *= x; 126 y2 *= y; 127 t = y2 - x2; 128 t /= f; 129 d += t; 130 } 131 while (fabsl(t/d) > MACHEPL); 132 return(d); 133 } 134 135 long double complex 136 ctanl(long double complex z) 137 { 138 long double complex w; 139 long double d, x, y; 140 141 x = creall(z); 142 y = cimagl(z); 143 d = cosl(2.0L * x) + coshl(2.0L * y); 144 145 if (fabsl(d) < 0.25L) { 146 d = fabsl(d); 147 d = ctansl(z); 148 } 149 if (d == 0.0L) { 150 /*mtherr( "ctan", OVERFLOW );*/ 151 w = LDBL_MAX + LDBL_MAX * I; 152 return (w); 153 } 154 155 w = sinl(2.0L * x) / d + (sinhl(2.0L * y) / d) * I; 156 return (w); 157 } 158