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

     1  /* origin: OpenBSD /usr/src/lib/libm/src/s_catanl.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   * long double complex catanl();
    24   * long double complex z, w;
    25   *
    26   * w = catanl( 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   *    DEC       -10,+10      5900       1.3e-16     7.8e-18
    54   *    IEEE      -10,+10     30000       2.3e-15     8.5e-17
    55   * The check catan( ctan(z) )  =  z, with |x| and |y| < PI/2,
    56   * had peak relative error 1.5e-16, rms relative error
    57   * 2.9e-17.  See also clog().
    58   */
    59  
    60  #include <complex.h>
    61  #include <float.h>
    62  #include "complex_impl.h"
    63  
    64  #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
    65  long double complex catanl(long double complex z)
    66  {
    67  	return catan(z);
    68  }
    69  #else
    70  static const long double PIL = 3.141592653589793238462643383279502884197169L;
    71  static const long double DP1 = 3.14159265358979323829596852490908531763125L;
    72  static const long double DP2 = 1.6667485837041756656403424829301998703007e-19L;
    73  static const long double DP3 = 1.8830410776607851167459095484560349402753e-39L;
    74  
    75  static long double redupil(long double x)
    76  {
    77  	long double t;
    78  	long i;
    79  
    80  	t = x / PIL;
    81  	if (t >= 0.0L)
    82  		t += 0.5L;
    83  	else
    84  		t -= 0.5L;
    85  
    86  	i = t;  /* the multiple */
    87  	t = i;
    88  	t = ((x - t * DP1) - t * DP2) - t * DP3;
    89  	return t;
    90  }
    91  
    92  long double complex catanl(long double complex z)
    93  {
    94  	long double complex w;
    95  	long double a, t, x, x2, y;
    96  
    97  	x = creall(z);
    98  	y = cimagl(z);
    99  
   100  	x2 = x * x;
   101  	a = 1.0L - x2 - (y * y);
   102  
   103  	t = atan2l(2.0L * x, a) * 0.5L;
   104  	w = redupil(t);
   105  
   106  	t = y - 1.0L;
   107  	a = x2 + (t * t);
   108  
   109  	t = y + 1.0L;
   110  	a = (x2 + (t * t)) / a;
   111  	w = CMPLXF(w, 0.25L * logl(a));
   112  	return w;
   113  }
   114  #endif