github.com/afumu/libc@v0.0.6/musl/src/complex/catanf.c (about)

     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  
    60  static const double DP1 = 3.140625;
    61  static const double DP2 = 9.67502593994140625E-4;
    62  static const double DP3 = 1.509957990978376432E-7;
    63  
    64  static 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  
    82  float 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  	x2 = x * x;
    91  	a = 1.0f - x2 - (y * y);
    92  
    93  	t = 0.5f * atan2f(2.0f * x, a);
    94  	w = _redupif(t);
    95  
    96  	t = y - 1.0f;
    97  	a = x2 + (t * t);
    98  
    99  	t = y + 1.0f;
   100  	a = (x2 + (t * t))/a;
   101  	w = CMPLXF(w, 0.25f * logf(a));
   102  	return w;
   103  }