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

     1  /* origin: OpenBSD /usr/src/lib/libm/src/s_catan.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   * double complex catan();
    24   * double complex z, w;
    25   *
    26   * w = catan (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   * catan(z) = -i catanh(iz).
    49   *
    50   * ACCURACY:
    51   *
    52   *                      Relative error:
    53   * arithmetic   domain     # trials      peak         rms
    54   *    DEC       -10,+10      5900       1.3e-16     7.8e-18
    55   *    IEEE      -10,+10     30000       2.3e-15     8.5e-17
    56   * The check catan( ctan(z) )  =  z, with |x| and |y| < PI/2,
    57   * had peak relative error 1.5e-16, rms relative error
    58   * 2.9e-17.  See also clog().
    59   */
    60  
    61  #include "complex_impl.h"
    62  
    63  #define MAXNUM 1.0e308
    64  
    65  static const double DP1 = 3.14159265160560607910E0;
    66  static const double DP2 = 1.98418714791870343106E-9;
    67  static const double DP3 = 1.14423774522196636802E-17;
    68  
    69  static double _redupi(double x)
    70  {
    71  	double t;
    72  	long i;
    73  
    74  	t = x/M_PI;
    75  	if (t >= 0.0)
    76  		t += 0.5;
    77  	else
    78  		t -= 0.5;
    79  
    80  	i = t;  /* the multiple */
    81  	t = i;
    82  	t = ((x - t * DP1) - t * DP2) - t * DP3;
    83  	return t;
    84  }
    85  
    86  double complex catan(double complex z)
    87  {
    88  	double complex w;
    89  	double a, t, x, x2, y;
    90  
    91  	x = creal(z);
    92  	y = cimag(z);
    93  
    94  	x2 = x * x;
    95  	a = 1.0 - x2 - (y * y);
    96  
    97  	t = 0.5 * atan2(2.0 * x, a);
    98  	w = _redupi(t);
    99  
   100  	t = y - 1.0;
   101  	a = x2 + (t * t);
   102  
   103  	t = y + 1.0;
   104  	a = (x2 + t * t)/a;
   105  	w = CMPLX(w, 0.25 * log(a));
   106  	return w;
   107  }