1  | 
     | 
     | 
    /*	$OpenBSD: s_catanl.c,v 1.4 2016/09/12 19:47:02 guenther Exp $	*/  | 
    
    
    2  | 
     | 
     | 
     | 
    
    
    3  | 
     | 
     | 
    /*  | 
    
    
    4  | 
     | 
     | 
     * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>  | 
    
    
    5  | 
     | 
     | 
     *  | 
    
    
    6  | 
     | 
     | 
     * Permission to use, copy, modify, and distribute this software for any  | 
    
    
    7  | 
     | 
     | 
     * purpose with or without fee is hereby granted, provided that the above  | 
    
    
    8  | 
     | 
     | 
     * copyright notice and this permission notice appear in all copies.  | 
    
    
    9  | 
     | 
     | 
     *  | 
    
    
    10  | 
     | 
     | 
     * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES  | 
    
    
    11  | 
     | 
     | 
     * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF  | 
    
    
    12  | 
     | 
     | 
     * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR  | 
    
    
    13  | 
     | 
     | 
     * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES  | 
    
    
    14  | 
     | 
     | 
     * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN  | 
    
    
    15  | 
     | 
     | 
     * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF  | 
    
    
    16  | 
     | 
     | 
     * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.  | 
    
    
    17  | 
     | 
     | 
     */  | 
    
    
    18  | 
     | 
     | 
     | 
    
    
    19  | 
     | 
     | 
    /*							catanl()  | 
    
    
    20  | 
     | 
     | 
     *  | 
    
    
    21  | 
     | 
     | 
     *	Complex circular arc tangent  | 
    
    
    22  | 
     | 
     | 
     *  | 
    
    
    23  | 
     | 
     | 
     *  | 
    
    
    24  | 
     | 
     | 
     *  | 
    
    
    25  | 
     | 
     | 
     * SYNOPSIS:  | 
    
    
    26  | 
     | 
     | 
     *  | 
    
    
    27  | 
     | 
     | 
     * long double complex catanl();  | 
    
    
    28  | 
     | 
     | 
     * long double complex z, w;  | 
    
    
    29  | 
     | 
     | 
     *  | 
    
    
    30  | 
     | 
     | 
     * w = catanl( z );  | 
    
    
    31  | 
     | 
     | 
     *  | 
    
    
    32  | 
     | 
     | 
     *  | 
    
    
    33  | 
     | 
     | 
     *  | 
    
    
    34  | 
     | 
     | 
     * DESCRIPTION:  | 
    
    
    35  | 
     | 
     | 
     *  | 
    
    
    36  | 
     | 
     | 
     * If  | 
    
    
    37  | 
     | 
     | 
     *     z = x + iy,  | 
    
    
    38  | 
     | 
     | 
     *  | 
    
    
    39  | 
     | 
     | 
     * then  | 
    
    
    40  | 
     | 
     | 
     *          1       (    2x     )  | 
    
    
    41  | 
     | 
     | 
     * Re w  =  - arctan(-----------)  +  k PI  | 
    
    
    42  | 
     | 
     | 
     *          2       (     2    2)  | 
    
    
    43  | 
     | 
     | 
     *                  (1 - x  - y )  | 
    
    
    44  | 
     | 
     | 
     *  | 
    
    
    45  | 
     | 
     | 
     *               ( 2         2)  | 
    
    
    46  | 
     | 
     | 
     *          1    (x  +  (y+1) )  | 
    
    
    47  | 
     | 
     | 
     * Im w  =  - log(------------)  | 
    
    
    48  | 
     | 
     | 
     *          4    ( 2         2)  | 
    
    
    49  | 
     | 
     | 
     *               (x  +  (y-1) )  | 
    
    
    50  | 
     | 
     | 
     *  | 
    
    
    51  | 
     | 
     | 
     * Where k is an arbitrary integer.  | 
    
    
    52  | 
     | 
     | 
     *  | 
    
    
    53  | 
     | 
     | 
     *  | 
    
    
    54  | 
     | 
     | 
     *  | 
    
    
    55  | 
     | 
     | 
     * ACCURACY:  | 
    
    
    56  | 
     | 
     | 
     *  | 
    
    
    57  | 
     | 
     | 
     *                      Relative error:  | 
    
    
    58  | 
     | 
     | 
     * arithmetic   domain     # trials      peak         rms  | 
    
    
    59  | 
     | 
     | 
     *    DEC       -10,+10      5900       1.3e-16     7.8e-18  | 
    
    
    60  | 
     | 
     | 
     *    IEEE      -10,+10     30000       2.3e-15     8.5e-17  | 
    
    
    61  | 
     | 
     | 
     * The check catan( ctan(z) )  =  z, with |x| and |y| < PI/2,  | 
    
    
    62  | 
     | 
     | 
     * had peak relative error 1.5e-16, rms relative error  | 
    
    
    63  | 
     | 
     | 
     * 2.9e-17.  See also clog().  | 
    
    
    64  | 
     | 
     | 
     */  | 
    
    
    65  | 
     | 
     | 
     | 
    
    
    66  | 
     | 
     | 
    #include <complex.h>  | 
    
    
    67  | 
     | 
     | 
    #include <float.h>  | 
    
    
    68  | 
     | 
     | 
    #include <math.h>  | 
    
    
    69  | 
     | 
     | 
     | 
    
    
    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  | 
    
    
    76  | 
     | 
     | 
    redupil(long double x)  | 
    
    
    77  | 
     | 
     | 
    { | 
    
    
    78  | 
     | 
     | 
    	long double t;  | 
    
    
    79  | 
     | 
     | 
    	long i;  | 
    
    
    80  | 
     | 
     | 
     | 
    
    
    81  | 
     | 
     | 
    	t = x / PIL;  | 
    
    
    82  | 
     | 
     | 
    	if (t >= 0.0L)  | 
    
    
    83  | 
     | 
     | 
    		t += 0.5L;  | 
    
    
    84  | 
     | 
     | 
    	else  | 
    
    
    85  | 
     | 
     | 
    		t -= 0.5L;  | 
    
    
    86  | 
     | 
     | 
     | 
    
    
    87  | 
     | 
     | 
    	i = t;	/* the multiple */  | 
    
    
    88  | 
     | 
     | 
    	t = i;  | 
    
    
    89  | 
     | 
     | 
    	t = ((x - t * DP1) - t * DP2) - t * DP3;  | 
    
    
    90  | 
     | 
     | 
    	return (t);  | 
    
    
    91  | 
     | 
     | 
    }  | 
    
    
    92  | 
     | 
     | 
     | 
    
    
    93  | 
     | 
     | 
    long double complex  | 
    
    
    94  | 
     | 
     | 
    catanl(long double complex z)  | 
    
    
    95  | 
     | 
     | 
    { | 
    
    
    96  | 
     | 
     | 
    	long double complex w;  | 
    
    
    97  | 
     | 
     | 
    	long double a, t, x, x2, y;  | 
    
    
    98  | 
     | 
     | 
     | 
    
    
    99  | 
     | 
     | 
    	x = creall(z);  | 
    
    
    100  | 
     | 
     | 
    	y = cimagl(z);  | 
    
    
    101  | 
     | 
     | 
     | 
    
    
    102  | 
     | 
     | 
    	if ((x == 0.0L) && (y > 1.0L))  | 
    
    
    103  | 
     | 
     | 
    		goto ovrf;  | 
    
    
    104  | 
     | 
     | 
     | 
    
    
    105  | 
     | 
     | 
    	x2 = x * x;  | 
    
    
    106  | 
     | 
     | 
    	a = 1.0L - x2 - (y * y);  | 
    
    
    107  | 
     | 
     | 
    	if (a == 0.0L)  | 
    
    
    108  | 
     | 
     | 
    		goto ovrf;  | 
    
    
    109  | 
     | 
     | 
     | 
    
    
    110  | 
     | 
     | 
    	t = atan2l(2.0L * x, a) * 0.5L;  | 
    
    
    111  | 
     | 
     | 
    	w = redupil(t);  | 
    
    
    112  | 
     | 
     | 
     | 
    
    
    113  | 
     | 
     | 
    	t = y - 1.0L;  | 
    
    
    114  | 
     | 
     | 
    	a = x2 + (t * t);  | 
    
    
    115  | 
     | 
     | 
    	if (a == 0.0L)  | 
    
    
    116  | 
     | 
     | 
    		goto ovrf;  | 
    
    
    117  | 
     | 
     | 
     | 
    
    
    118  | 
     | 
     | 
    	t = y + 1.0L;  | 
    
    
    119  | 
     | 
     | 
    	a = (x2 + (t * t)) / a;  | 
    
    
    120  | 
     | 
     | 
    	w = w + (0.25L * logl(a)) * I;  | 
    
    
    121  | 
     | 
     | 
    	return (w);  | 
    
    
    122  | 
     | 
     | 
     | 
    
    
    123  | 
     | 
     | 
    ovrf:  | 
    
    
    124  | 
     | 
     | 
    	/*mtherr( "catanl", OVERFLOW );*/  | 
    
    
    125  | 
     | 
     | 
    	w = LDBL_MAX + LDBL_MAX * I;  | 
    
    
    126  | 
     | 
     | 
    	return (w);  | 
    
    
    127  | 
     | 
     | 
    }  | 
    
    
    128  | 
     | 
     | 
    DEF_STD(catanl);  |