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 "complex_impl.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 const float float_pi = M_PI;
 65
 66static float _redupif(float xx)
 67{
 68	float x, t;
 69	long i;
 70
 71	x = xx;
 72	t = x/float_pi;
 73	if (t >= 0.0f)
 74		t += 0.5f;
 75	else
 76		t -= 0.5f;
 77
 78	i = t;  /* the multiple */
 79	t = i;
 80	t = ((x - t * DP1) - t * DP2) - t * DP3;
 81	return t;
 82}
 83
 84float complex catanf(float complex z)
 85{
 86	float complex w;
 87	float a, t, x, x2, y;
 88
 89	x = crealf(z);
 90	y = cimagf(z);
 91
 92	x2 = x * x;
 93	a = 1.0f - x2 - (y * y);
 94
 95	t = 0.5f * atan2f(2.0f * x, a);
 96	w = _redupif(t);
 97
 98	t = y - 1.0f;
 99	a = x2 + (t * t);
100
101	t = y + 1.0f;
102	a = (x2 + (t * t))/a;
103	w = CMPLXF(w, 0.25f * logf(a));
104	return w;
105}