Home | History | Annotate | Line # | Download | only in ld80 
      1 /*-
      2  * SPDX-License-Identifier: BSD-2-Clause
      3  *
      4  * Copyright (c) 2008 David Schultz <das (at) FreeBSD.ORG>
      5  * All rights reserved.
      6  *
      7  * Redistribution and use in source and binary forms, with or without
      8  * modification, are permitted provided that the following conditions
      9  * are met:
     10  * 1. Redistributions of source code must retain the above copyright
     11  *    notice, this list of conditions and the following 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 AND CONTRIBUTORS ``AS IS'' AND
     17  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
     18  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
     19  * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
     20  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
     21  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
     22  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
     23  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
     24  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
     25  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
     26  * SUCH DAMAGE.
     27  */
     28 
     29 #include <float.h>
     30 
     31 #ifdef __FreeBSD__
     32 #include "fpmath.h"
     33 #endif
     34 
     35 #define	BIAS		(LDBL_MAX_EXP - 1)
     36 #define	MANH_SIZE	EXT_FRACHBITS
     37 
     38 /* Approximation thresholds. */
     39 #define	ASIN_LINEAR	(BIAS - 32)	/* 2**-32 */
     40 #define	ACOS_CONST	(BIAS - 65)	/* 2**-65 */
     41 #define	ATAN_CONST	(BIAS + 65)	/* 2**65 */
     42 #define	ATAN_LINEAR	(BIAS - 32)	/* 2**-32 */
     43 
     44 /* 0.95 */
     45 #define	THRESH	((0xe666666666666666ULL>>(64-(MANH_SIZE-1)))|LDBL_NBIT)
     46 
     47 /* Constants shared by the long double inverse trig functions. */
     48 #define	pS0	_ItL_pS0
     49 #define	pS1	_ItL_pS1
     50 #define	pS2	_ItL_pS2
     51 #define	pS3	_ItL_pS3
     52 #define	pS4	_ItL_pS4
     53 #define	pS5	_ItL_pS5
     54 #define	pS6	_ItL_pS6
     55 #define	qS1	_ItL_qS1
     56 #define	qS2	_ItL_qS2
     57 #define	qS3	_ItL_qS3
     58 #define	qS4	_ItL_qS4
     59 #define	qS5	_ItL_qS5
     60 #define	atanhi	_ItL_atanhi
     61 #define	atanlo	_ItL_atanlo
     62 #define	aT	_ItL_aT
     63 #define	pi_lo	_ItL_pi_lo
     64 
     65 #define	pio2_hi	atanhi[3]
     66 #define	pio2_lo	atanlo[3]
     67 #define	pio4_hi	atanhi[1]
     68 
     69 #ifdef STRUCT_DECLS
     70 typedef struct longdouble {
     71 	uint64_t mant;
     72 	uint16_t expsign;
     73 } LONGDOUBLE;
     74 #else
     75 typedef long double LONGDOUBLE;
     76 #endif
     77 
     78 extern const LONGDOUBLE pS0, pS1, pS2, pS3, pS4, pS5, pS6;
     79 extern const LONGDOUBLE qS1, qS2, qS3, qS4, qS5;
     80 extern const LONGDOUBLE atanhi[], atanlo[], aT[];
     81 extern const LONGDOUBLE pi_lo;
     82 
     83 #ifndef STRUCT_DECLS
     84 
     85 static inline long double
     86 P(long double x)
     87 {
     88 
     89 	return (x * (pS0 + x * (pS1 + x * (pS2 + x * (pS3 + x * \
     90 		(pS4 + x * (pS5 + x * pS6)))))));
     91 }
     92 
     93 static inline long double
     94 Q(long double x)
     95 {
     96 
     97 	return (1.0 + x * (qS1 + x * (qS2 + x * (qS3 + x * (qS4 + x * qS5)))));
     98 }
     99 
    100 static inline long double
    101 T_even(long double x)
    102 {
    103 
    104 	return (aT[0] + x * (aT[2] + x * (aT[4] + x * (aT[6] + x * \
    105 		(aT[8] + x * (aT[10] + x * aT[12]))))));
    106 }
    107 
    108 static inline long double
    109 T_odd(long double x)
    110 {
    111 
    112 	return (aT[1] + x * (aT[3] + x * (aT[5] + x * (aT[7] + x * \
    113 		(aT[9] + x * aT[11])))));
    114 }
    115 
    116 #endif
    117