1/* origin: OpenBSD /usr/src/lib/libm/src/s_catanf.c */ 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 * Complex circular arc tangent 19 * 20 * 21 * SYNOPSIS: 22 * 23 * float complex catanf(); 24 * float complex z, w; 25 * 26 * w = catanf( z ); 27 * 28 * 29 * DESCRIPTION: 30 * 31 * If 32 * z = x + iy, 33 * 34 * then 35 * 1 ( 2x ) 36 * Re w = - arctan(-----------) + k PI 37 * 2 ( 2 2) 38 * (1 - x - y ) 39 * 40 * ( 2 2) 41 * 1 (x + (y+1) ) 42 * Im w = - log(------------) 43 * 4 ( 2 2) 44 * (x + (y-1) ) 45 * 46 * Where k is an arbitrary integer. 47 * 48 * 49 * ACCURACY: 50 * 51 * Relative error: 52 * arithmetic domain # trials peak rms 53 * IEEE -10,+10 30000 2.3e-6 5.2e-8 54 */ 55 56#include "libm.h" 57 58#define MAXNUMF 1.0e38F 59 60static const double DP1 = 3.140625; 61static const double DP2 = 9.67502593994140625E-4; 62static const double DP3 = 1.509957990978376432E-7; 63 64static float _redupif(float xx) 65{ 66 float x, t; 67 long i; 68 69 x = xx; 70 t = x/(float)M_PI; 71 if (t >= 0.0f) 72 t += 0.5f; 73 else 74 t -= 0.5f; 75 76 i = t; /* the multiple */ 77 t = i; 78 t = ((x - t * DP1) - t * DP2) - t * DP3; 79 return t; 80} 81 82float complex catanf(float complex z) 83{ 84 float complex w; 85 float a, t, x, x2, y; 86 87 x = crealf(z); 88 y = cimagf(z); 89 90 if ((x == 0.0f) && (y > 1.0f)) 91 goto ovrf; 92 93 x2 = x * x; 94 a = 1.0f - x2 - (y * y); 95 if (a == 0.0f) 96 goto ovrf; 97 98 t = 0.5f * atan2f(2.0f * x, a); 99 w = _redupif(t); 100 101 t = y - 1.0f; 102 a = x2 + (t * t); 103 if (a == 0.0f) 104 goto ovrf; 105 106 t = y + 1.0f; 107 a = (x2 + (t * t))/a; 108 w = w + (0.25f * logf (a)) * I; 109 return w; 110 111ovrf: 112 // FIXME 113 w = MAXNUMF + MAXNUMF * I; 114 return w; 115} 116