1  | 
     | 
     | 
    /*	$OpenBSD: s_ctanl.c,v 1.3 2011/07/20 21:02:51 martynas 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  | 
     | 
     | 
    /*							ctanl()  | 
    
    
    20  | 
     | 
     | 
     *  | 
    
    
    21  | 
     | 
     | 
     *	Complex circular tangent  | 
    
    
    22  | 
     | 
     | 
     *  | 
    
    
    23  | 
     | 
     | 
     *  | 
    
    
    24  | 
     | 
     | 
     *  | 
    
    
    25  | 
     | 
     | 
     * SYNOPSIS:  | 
    
    
    26  | 
     | 
     | 
     *  | 
    
    
    27  | 
     | 
     | 
     * long double complex ctanl();  | 
    
    
    28  | 
     | 
     | 
     * long double complex z, w;  | 
    
    
    29  | 
     | 
     | 
     *  | 
    
    
    30  | 
     | 
     | 
     * w = ctanl( z );  | 
    
    
    31  | 
     | 
     | 
     *  | 
    
    
    32  | 
     | 
     | 
     *  | 
    
    
    33  | 
     | 
     | 
     *  | 
    
    
    34  | 
     | 
     | 
     * DESCRIPTION:  | 
    
    
    35  | 
     | 
     | 
     *  | 
    
    
    36  | 
     | 
     | 
     * If  | 
    
    
    37  | 
     | 
     | 
     *     z = x + iy,  | 
    
    
    38  | 
     | 
     | 
     *  | 
    
    
    39  | 
     | 
     | 
     * then  | 
    
    
    40  | 
     | 
     | 
     *  | 
    
    
    41  | 
     | 
     | 
     *           sin 2x  +  i sinh 2y  | 
    
    
    42  | 
     | 
     | 
     *     w  =  --------------------.  | 
    
    
    43  | 
     | 
     | 
     *            cos 2x  +  cosh 2y  | 
    
    
    44  | 
     | 
     | 
     *  | 
    
    
    45  | 
     | 
     | 
     * On the real axis the denominator is zero at odd multiples  | 
    
    
    46  | 
     | 
     | 
     * of PI/2.  The denominator is evaluated by its Taylor  | 
    
    
    47  | 
     | 
     | 
     * series near these points.  | 
    
    
    48  | 
     | 
     | 
     *  | 
    
    
    49  | 
     | 
     | 
     *  | 
    
    
    50  | 
     | 
     | 
     * ACCURACY:  | 
    
    
    51  | 
     | 
     | 
     *  | 
    
    
    52  | 
     | 
     | 
     *                      Relative error:  | 
    
    
    53  | 
     | 
     | 
     * arithmetic   domain     # trials      peak         rms  | 
    
    
    54  | 
     | 
     | 
     *    DEC       -10,+10      5200       7.1e-17     1.6e-17  | 
    
    
    55  | 
     | 
     | 
     *    IEEE      -10,+10     30000       7.2e-16     1.2e-16  | 
    
    
    56  | 
     | 
     | 
     * Also tested by ctan * ccot = 1 and catan(ctan(z))  =  z.  | 
    
    
    57  | 
     | 
     | 
     */  | 
    
    
    58  | 
     | 
     | 
     | 
    
    
    59  | 
     | 
     | 
    #include <complex.h>  | 
    
    
    60  | 
     | 
     | 
    #include <float.h>  | 
    
    
    61  | 
     | 
     | 
    #include <math.h>  | 
    
    
    62  | 
     | 
     | 
     | 
    
    
    63  | 
     | 
     | 
    #if	LDBL_MANT_DIG == 64  | 
    
    
    64  | 
     | 
     | 
    static const long double MACHEPL= 5.42101086242752217003726400434970855712890625E-20L;  | 
    
    
    65  | 
     | 
     | 
    #elif	LDBL_MANT_DIG == 113  | 
    
    
    66  | 
     | 
     | 
    static const long double MACHEPL = 9.629649721936179265279889712924636592690508e-35L;  | 
    
    
    67  | 
     | 
     | 
    #endif  | 
    
    
    68  | 
     | 
     | 
     | 
    
    
    69  | 
     | 
     | 
    static const long double PIL = 3.141592653589793238462643383279502884197169L;  | 
    
    
    70  | 
     | 
     | 
    static const long double DP1 = 3.14159265358979323829596852490908531763125L;  | 
    
    
    71  | 
     | 
     | 
    static const long double DP2 = 1.6667485837041756656403424829301998703007e-19L;  | 
    
    
    72  | 
     | 
     | 
    static const long double DP3 = 1.8830410776607851167459095484560349402753e-39L;  | 
    
    
    73  | 
     | 
     | 
     | 
    
    
    74  | 
     | 
     | 
    static long double  | 
    
    
    75  | 
     | 
     | 
    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  | 
     | 
     | 
    static long double  | 
    
    
    93  | 
     | 
     | 
    ctansl(long double complex z)  | 
    
    
    94  | 
     | 
     | 
    { | 
    
    
    95  | 
     | 
     | 
    	long double f, x, x2, y, y2, rn, t;  | 
    
    
    96  | 
     | 
     | 
    	long double d;  | 
    
    
    97  | 
     | 
     | 
     | 
    
    
    98  | 
     | 
     | 
    	x = fabsl(2.0L * creall(z));  | 
    
    
    99  | 
     | 
     | 
    	y = fabsl(2.0L * cimagl(z));  | 
    
    
    100  | 
     | 
     | 
     | 
    
    
    101  | 
     | 
     | 
    	x = redupil(x);  | 
    
    
    102  | 
     | 
     | 
     | 
    
    
    103  | 
     | 
     | 
    	x = x * x;  | 
    
    
    104  | 
     | 
     | 
    	y = y * y;  | 
    
    
    105  | 
     | 
     | 
    	x2 = 1.0L;  | 
    
    
    106  | 
     | 
     | 
    	y2 = 1.0L;  | 
    
    
    107  | 
     | 
     | 
    	f = 1.0L;  | 
    
    
    108  | 
     | 
     | 
    	rn = 0.0L;  | 
    
    
    109  | 
     | 
     | 
    	d = 0.0L;  | 
    
    
    110  | 
     | 
     | 
    	do { | 
    
    
    111  | 
     | 
     | 
    		rn += 1.0L;  | 
    
    
    112  | 
     | 
     | 
    		f *= rn;  | 
    
    
    113  | 
     | 
     | 
    		rn += 1.0L;  | 
    
    
    114  | 
     | 
     | 
    		f *= rn;  | 
    
    
    115  | 
     | 
     | 
    		x2 *= x;  | 
    
    
    116  | 
     | 
     | 
    		y2 *= y;  | 
    
    
    117  | 
     | 
     | 
    		t = y2 + x2;  | 
    
    
    118  | 
     | 
     | 
    		t /= f;  | 
    
    
    119  | 
     | 
     | 
    		d += t;  | 
    
    
    120  | 
     | 
     | 
     | 
    
    
    121  | 
     | 
     | 
    		rn += 1.0L;  | 
    
    
    122  | 
     | 
     | 
    		f *= rn;  | 
    
    
    123  | 
     | 
     | 
    		rn += 1.0L;  | 
    
    
    124  | 
     | 
     | 
    		f *= rn;  | 
    
    
    125  | 
     | 
     | 
    		x2 *= x;  | 
    
    
    126  | 
     | 
     | 
    		y2 *= y;  | 
    
    
    127  | 
     | 
     | 
    		t = y2 - x2;  | 
    
    
    128  | 
     | 
     | 
    		t /= f;  | 
    
    
    129  | 
     | 
     | 
    		d += t;  | 
    
    
    130  | 
     | 
     | 
    	}  | 
    
    
    131  | 
     | 
     | 
    	while (fabsl(t/d) > MACHEPL);  | 
    
    
    132  | 
     | 
     | 
    	return(d);  | 
    
    
    133  | 
     | 
     | 
    }  | 
    
    
    134  | 
     | 
     | 
     | 
    
    
    135  | 
     | 
     | 
    long double complex  | 
    
    
    136  | 
     | 
     | 
    ctanl(long double complex z)  | 
    
    
    137  | 
     | 
     | 
    { | 
    
    
    138  | 
     | 
     | 
    	long double complex w;  | 
    
    
    139  | 
     | 
     | 
    	long double d, x, y;  | 
    
    
    140  | 
     | 
     | 
     | 
    
    
    141  | 
     | 
     | 
    	x = creall(z);  | 
    
    
    142  | 
     | 
     | 
    	y = cimagl(z);  | 
    
    
    143  | 
     | 
     | 
    	d = cosl(2.0L * x) + coshl(2.0L * y);  | 
    
    
    144  | 
     | 
     | 
     | 
    
    
    145  | 
     | 
     | 
    	if (fabsl(d) < 0.25L) { | 
    
    
    146  | 
     | 
     | 
    		d = fabsl(d);  | 
    
    
    147  | 
     | 
     | 
    		d = ctansl(z);  | 
    
    
    148  | 
     | 
     | 
    	}  | 
    
    
    149  | 
     | 
     | 
    	if (d == 0.0L) { | 
    
    
    150  | 
     | 
     | 
    		/*mtherr( "ctan", OVERFLOW );*/  | 
    
    
    151  | 
     | 
     | 
    		w = LDBL_MAX + LDBL_MAX * I;  | 
    
    
    152  | 
     | 
     | 
    		return (w);  | 
    
    
    153  | 
     | 
     | 
    	}  | 
    
    
    154  | 
     | 
     | 
     | 
    
    
    155  | 
     | 
     | 
    	w = sinl(2.0L * x) / d + (sinhl(2.0L * y) / d) * I;  | 
    
    
    156  | 
     | 
     | 
    	return (w);  | 
    
    
    157  | 
     | 
     | 
    }  |