1  | 
     | 
     | 
    /*	$OpenBSD: s_tanl.c,v 1.1 2008/12/09 20:00:35 martynas Exp $	*/  | 
    
    
    2  | 
     | 
     | 
    /*-  | 
    
    
    3  | 
     | 
     | 
     * Copyright (c) 2007 Steven G. Kargl  | 
    
    
    4  | 
     | 
     | 
     * All rights reserved.  | 
    
    
    5  | 
     | 
     | 
     *  | 
    
    
    6  | 
     | 
     | 
     * Redistribution and use in source and binary forms, with or without  | 
    
    
    7  | 
     | 
     | 
     * modification, are permitted provided that the following conditions  | 
    
    
    8  | 
     | 
     | 
     * are met:  | 
    
    
    9  | 
     | 
     | 
     * 1. Redistributions of source code must retain the above copyright  | 
    
    
    10  | 
     | 
     | 
     *    notice unmodified, this list of conditions, and the following  | 
    
    
    11  | 
     | 
     | 
     *    disclaimer.  | 
    
    
    12  | 
     | 
     | 
     * 2. Redistributions in binary form must reproduce the above copyright  | 
    
    
    13  | 
     | 
     | 
     *    notice, this list of conditions and the following disclaimer in the  | 
    
    
    14  | 
     | 
     | 
     *    documentation and/or other materials provided with the distribution.  | 
    
    
    15  | 
     | 
     | 
     *  | 
    
    
    16  | 
     | 
     | 
     * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR  | 
    
    
    17  | 
     | 
     | 
     * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES  | 
    
    
    18  | 
     | 
     | 
     * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.  | 
    
    
    19  | 
     | 
     | 
     * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,  | 
    
    
    20  | 
     | 
     | 
     * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT  | 
    
    
    21  | 
     | 
     | 
     * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,  | 
    
    
    22  | 
     | 
     | 
     * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY  | 
    
    
    23  | 
     | 
     | 
     * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT  | 
    
    
    24  | 
     | 
     | 
     * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF  | 
    
    
    25  | 
     | 
     | 
     * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.  | 
    
    
    26  | 
     | 
     | 
     */  | 
    
    
    27  | 
     | 
     | 
     | 
    
    
    28  | 
     | 
     | 
    /*  | 
    
    
    29  | 
     | 
     | 
     * Compute tan(x) for x where x is reduced to y = x - k * pi / 2.  | 
    
    
    30  | 
     | 
     | 
     * Limited testing on pseudorandom numbers drawn within [0:4e8] shows  | 
    
    
    31  | 
     | 
     | 
     * an accuracy of <= 1.5 ULP where 247024 values of x out of 40 million  | 
    
    
    32  | 
     | 
     | 
     * possibles resulted in tan(x) that exceeded 0.5 ULP (ie., 0.6%).  | 
    
    
    33  | 
     | 
     | 
     */  | 
    
    
    34  | 
     | 
     | 
     | 
    
    
    35  | 
     | 
     | 
    #include <sys/types.h>  | 
    
    
    36  | 
     | 
     | 
    #include <machine/ieee.h>  | 
    
    
    37  | 
     | 
     | 
    #include <float.h>  | 
    
    
    38  | 
     | 
     | 
    #include <math.h>  | 
    
    
    39  | 
     | 
     | 
     | 
    
    
    40  | 
     | 
     | 
    #include "math_private.h"  | 
    
    
    41  | 
     | 
     | 
     | 
    
    
    42  | 
     | 
     | 
    #if LDBL_MANT_DIG == 64  | 
    
    
    43  | 
     | 
     | 
    #define	NX	3  | 
    
    
    44  | 
     | 
     | 
    #define	PREC	2  | 
    
    
    45  | 
     | 
     | 
    #elif LDBL_MANT_DIG == 113  | 
    
    
    46  | 
     | 
     | 
    #define	NX	5  | 
    
    
    47  | 
     | 
     | 
    #define	PREC	3  | 
    
    
    48  | 
     | 
     | 
    #else  | 
    
    
    49  | 
     | 
     | 
    #error "Unsupported long double format"  | 
    
    
    50  | 
     | 
     | 
    #endif  | 
    
    
    51  | 
     | 
     | 
     | 
    
    
    52  | 
     | 
     | 
    static const long double two24 = 1.67772160000000000000e+07L;  | 
    
    
    53  | 
     | 
     | 
     | 
    
    
    54  | 
     | 
     | 
    long double  | 
    
    
    55  | 
     | 
     | 
    tanl(long double x)  | 
    
    
    56  | 
     | 
     | 
    { | 
    
    
    57  | 
     | 
     | 
    	union { | 
    
    
    58  | 
     | 
     | 
    		long double e;  | 
    
    
    59  | 
     | 
     | 
    		struct ieee_ext bits;  | 
    
    
    60  | 
     | 
     | 
    	} z;  | 
    
    
    61  | 
     | 
     | 
    	int i, e0, s;  | 
    
    
    62  | 
     | 
     | 
    	double xd[NX], yd[PREC];  | 
    
    
    63  | 
     | 
     | 
    	long double hi, lo;  | 
    
    
    64  | 
     | 
     | 
     | 
    
    
    65  | 
     | 
     | 
    	z.e = x;  | 
    
    
    66  | 
     | 
     | 
    	s = z.bits.ext_sign;  | 
    
    
    67  | 
     | 
     | 
    	z.bits.ext_sign = 0;  | 
    
    
    68  | 
     | 
     | 
     | 
    
    
    69  | 
     | 
     | 
    	/* If x = +-0 or x is subnormal, then tan(x) = x. */  | 
    
    
    70  | 
     | 
     | 
    	if (z.bits.ext_exp == 0)  | 
    
    
    71  | 
     | 
     | 
    		return (x);  | 
    
    
    72  | 
     | 
     | 
     | 
    
    
    73  | 
     | 
     | 
    	/* If x = NaN or Inf, then tan(x) = NaN. */  | 
    
    
    74  | 
     | 
     | 
    	if (z.bits.ext_exp == 32767)  | 
    
    
    75  | 
     | 
     | 
    		return ((x - x) / (x - x));  | 
    
    
    76  | 
     | 
     | 
     | 
    
    
    77  | 
     | 
     | 
    	/* Optimize the case where x is already within range. */  | 
    
    
    78  | 
     | 
     | 
    	if (z.e < M_PI_4) { | 
    
    
    79  | 
     | 
     | 
    		hi = __kernel_tanl(z.e, 0, 0);  | 
    
    
    80  | 
     | 
     | 
    		return (s ? -hi : hi);  | 
    
    
    81  | 
     | 
     | 
    	}  | 
    
    
    82  | 
     | 
     | 
     | 
    
    
    83  | 
     | 
     | 
    	/* Split z.e into a 24-bit representation. */  | 
    
    
    84  | 
     | 
     | 
    	e0 = ilogbl(z.e) - 23;  | 
    
    
    85  | 
     | 
     | 
    	z.e = scalbnl(z.e, -e0);  | 
    
    
    86  | 
     | 
     | 
    	for (i = 0; i < NX; i++) { | 
    
    
    87  | 
     | 
     | 
    		xd[i] = (double)((int32_t)z.e);  | 
    
    
    88  | 
     | 
     | 
    		z.e = (z.e - xd[i]) * two24;  | 
    
    
    89  | 
     | 
     | 
    	}  | 
    
    
    90  | 
     | 
     | 
     | 
    
    
    91  | 
     | 
     | 
    	/* yd contains the pieces of xd rem pi/2 such that |yd| < pi/4. */  | 
    
    
    92  | 
     | 
     | 
    	e0 = __kernel_rem_pio2(xd, yd, e0, NX, PREC);  | 
    
    
    93  | 
     | 
     | 
     | 
    
    
    94  | 
     | 
     | 
    #if PREC == 2  | 
    
    
    95  | 
     | 
     | 
    	hi = (long double)yd[0] + yd[1];  | 
    
    
    96  | 
     | 
     | 
    	lo = yd[1] - (hi - yd[0]);  | 
    
    
    97  | 
     | 
     | 
    #else /* PREC == 3 */  | 
    
    
    98  | 
     | 
     | 
    	long double t;  | 
    
    
    99  | 
     | 
     | 
    	t = (long double)yd[2] + yd[1];  | 
    
    
    100  | 
     | 
     | 
    	hi = t + yd[0];  | 
    
    
    101  | 
     | 
     | 
    	lo = yd[0] - (hi - t);  | 
    
    
    102  | 
     | 
     | 
    #endif  | 
    
    
    103  | 
     | 
     | 
     | 
    
    
    104  | 
     | 
     | 
    	switch (e0 & 3) { | 
    
    
    105  | 
     | 
     | 
    	case 0:  | 
    
    
    106  | 
     | 
     | 
    	case 2:  | 
    
    
    107  | 
     | 
     | 
    	    hi = __kernel_tanl(hi, lo, 0);  | 
    
    
    108  | 
     | 
     | 
    	    break;  | 
    
    
    109  | 
     | 
     | 
    	case 1:  | 
    
    
    110  | 
     | 
     | 
    	case 3:  | 
    
    
    111  | 
     | 
     | 
    	    hi = __kernel_tanl(hi, lo, 1);  | 
    
    
    112  | 
     | 
     | 
    	    break;  | 
    
    
    113  | 
     | 
     | 
    	}  | 
    
    
    114  | 
     | 
     | 
     | 
    
    
    115  | 
     | 
     | 
    	return (s ? -hi : hi);  | 
    
    
    116  | 
     | 
     | 
    }  |