catan.c raw

   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  }
 108