1 /* $OpenBSD: s_ctanf.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 /* ctanf() 19 * 20 * Complex circular tangent 21 * 22 * 23 * 24 * SYNOPSIS: 25 * 26 * void ctanf(); 27 * cmplxf z, w; 28 * 29 * ctanf( &z, &w ); 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 * 49 * ACCURACY: 50 * 51 * Relative error: 52 * arithmetic domain # trials peak rms 53 * IEEE -10,+10 30000 3.3e-7 5.1e-8 54 */ 55 56 #include <complex.h> 57 #include <math.h> 58 59 #define MACHEPF 3.0e-8 60 #define MAXNUMF 1.0e38f 61 62 static const double DP1 = 3.140625; 63 static const double DP2 = 9.67502593994140625E-4; 64 static const double DP3 = 1.509957990978376432E-7; 65 66 static float 67 _redupif(float xx) 68 { 69 float x, t; 70 long i; 71 72 x = xx; 73 t = x/(float)M_PI; 74 if(t >= 0.0) 75 t += 0.5; 76 else 77 t -= 0.5; 78 79 i = t; /* the multiple */ 80 t = i; 81 t = ((x - t * DP1) - t * DP2) - t * DP3; 82 return(t); 83 } 84 85 /* Taylor series expansion for cosh(2y) - cos(2x) */ 86 87 static float 88 _ctansf(float complex z) 89 { 90 float f, x, x2, y, y2, rn, t, d; 91 92 x = fabsf(2.0f * creal(z)); 93 y = fabsf(2.0f * cimag(z)); 94 95 x = _redupif(x); 96 97 x = x * x; 98 y = y * y; 99 x2 = 1.0f; 100 y2 = 1.0f; 101 f = 1.0f; 102 rn = 0.0f; 103 d = 0.0f; 104 do { 105 rn += 1.0f; 106 f *= rn; 107 rn += 1.0f; 108 f *= rn; 109 x2 *= x; 110 y2 *= y; 111 t = y2 + x2; 112 t /= f; 113 d += t; 114 115 rn += 1.0f; 116 f *= rn; 117 rn += 1.0f; 118 f *= rn; 119 x2 *= x; 120 y2 *= y; 121 t = y2 - x2; 122 t /= f; 123 d += t; 124 } 125 while (fabsf(t/d) > MACHEPF) 126 ; 127 return(d); 128 } 129 130 float complex 131 ctanf(float complex z) 132 { 133 float complex w; 134 float d; 135 136 d = cosf( 2.0f * creal(z) ) + coshf( 2.0f * cimag(z) ); 137 138 if(fabsf(d) < 0.25f) 139 d = _ctansf(z); 140 141 if (d == 0.0f) { 142 /*mtherr( "ctanf", OVERFLOW );*/ 143 w = MAXNUMF + MAXNUMF * I; 144 return (w); 145 } 146 w = sinf (2.0f * creal(z)) / d + (sinhf (2.0f * cimag(z)) / d) * I; 147 return (w); 148 } 149