117303c62SDavid Schultz /*-
2*4d846d26SWarner Losh * SPDX-License-Identifier: BSD-2-Clause
35e53a4f9SPedro F. Giffuni *
417303c62SDavid Schultz * Copyright (c) 2008 David Schultz <das@FreeBSD.ORG>
517303c62SDavid Schultz * All rights reserved.
617303c62SDavid Schultz *
717303c62SDavid Schultz * Redistribution and use in source and binary forms, with or without
817303c62SDavid Schultz * modification, are permitted provided that the following conditions
917303c62SDavid Schultz * are met:
1017303c62SDavid Schultz * 1. Redistributions of source code must retain the above copyright
1117303c62SDavid Schultz * notice, this list of conditions and the following disclaimer.
1217303c62SDavid Schultz * 2. Redistributions in binary form must reproduce the above copyright
1317303c62SDavid Schultz * notice, this list of conditions and the following disclaimer in the
1417303c62SDavid Schultz * documentation and/or other materials provided with the distribution.
1517303c62SDavid Schultz *
1617303c62SDavid Schultz * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
1717303c62SDavid Schultz * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
1817303c62SDavid Schultz * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
1917303c62SDavid Schultz * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
2017303c62SDavid Schultz * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
2117303c62SDavid Schultz * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
2217303c62SDavid Schultz * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
2317303c62SDavid Schultz * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
2417303c62SDavid Schultz * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
2517303c62SDavid Schultz * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
2617303c62SDavid Schultz * SUCH DAMAGE.
2717303c62SDavid Schultz */
2817303c62SDavid Schultz
2917303c62SDavid Schultz #include <float.h>
3017303c62SDavid Schultz
3117303c62SDavid Schultz #include "fpmath.h"
3217303c62SDavid Schultz
3317303c62SDavid Schultz #define BIAS (LDBL_MAX_EXP - 1)
3417303c62SDavid Schultz #define MANH_SIZE LDBL_MANH_SIZE
3517303c62SDavid Schultz
3617303c62SDavid Schultz /* Approximation thresholds. */
3717303c62SDavid Schultz #define ASIN_LINEAR (BIAS - 32) /* 2**-32 */
3817303c62SDavid Schultz #define ACOS_CONST (BIAS - 65) /* 2**-65 */
3917303c62SDavid Schultz #define ATAN_CONST (BIAS + 65) /* 2**65 */
4017303c62SDavid Schultz #define ATAN_LINEAR (BIAS - 32) /* 2**-32 */
4117303c62SDavid Schultz
4217303c62SDavid Schultz /* 0.95 */
4317303c62SDavid Schultz #define THRESH ((0xe666666666666666ULL>>(64-(MANH_SIZE-1)))|LDBL_NBIT)
4417303c62SDavid Schultz
4517303c62SDavid Schultz /* Constants shared by the long double inverse trig functions. */
4617303c62SDavid Schultz #define pS0 _ItL_pS0
4717303c62SDavid Schultz #define pS1 _ItL_pS1
4817303c62SDavid Schultz #define pS2 _ItL_pS2
4917303c62SDavid Schultz #define pS3 _ItL_pS3
5017303c62SDavid Schultz #define pS4 _ItL_pS4
5117303c62SDavid Schultz #define pS5 _ItL_pS5
5217303c62SDavid Schultz #define pS6 _ItL_pS6
5317303c62SDavid Schultz #define qS1 _ItL_qS1
5417303c62SDavid Schultz #define qS2 _ItL_qS2
5517303c62SDavid Schultz #define qS3 _ItL_qS3
5617303c62SDavid Schultz #define qS4 _ItL_qS4
5717303c62SDavid Schultz #define qS5 _ItL_qS5
5817303c62SDavid Schultz #define atanhi _ItL_atanhi
5917303c62SDavid Schultz #define atanlo _ItL_atanlo
6017303c62SDavid Schultz #define aT _ItL_aT
6117303c62SDavid Schultz #define pi_lo _ItL_pi_lo
6217303c62SDavid Schultz
6317303c62SDavid Schultz #define pio2_hi atanhi[3]
6417303c62SDavid Schultz #define pio2_lo atanlo[3]
6517303c62SDavid Schultz #define pio4_hi atanhi[1]
6617303c62SDavid Schultz
671192a80eSDavid Schultz #ifdef STRUCT_DECLS
681192a80eSDavid Schultz typedef struct longdouble {
691192a80eSDavid Schultz uint64_t mant;
701192a80eSDavid Schultz uint16_t expsign;
711192a80eSDavid Schultz } LONGDOUBLE;
721192a80eSDavid Schultz #else
731192a80eSDavid Schultz typedef long double LONGDOUBLE;
741192a80eSDavid Schultz #endif
751192a80eSDavid Schultz
761192a80eSDavid Schultz extern const LONGDOUBLE pS0, pS1, pS2, pS3, pS4, pS5, pS6;
771192a80eSDavid Schultz extern const LONGDOUBLE qS1, qS2, qS3, qS4, qS5;
781192a80eSDavid Schultz extern const LONGDOUBLE atanhi[], atanlo[], aT[];
791192a80eSDavid Schultz extern const LONGDOUBLE pi_lo;
801192a80eSDavid Schultz
811192a80eSDavid Schultz #ifndef STRUCT_DECLS
8217303c62SDavid Schultz
8317303c62SDavid Schultz static inline long double
P(long double x)8417303c62SDavid Schultz P(long double x)
8517303c62SDavid Schultz {
8617303c62SDavid Schultz
8717303c62SDavid Schultz return (x * (pS0 + x * (pS1 + x * (pS2 + x * (pS3 + x * \
8817303c62SDavid Schultz (pS4 + x * (pS5 + x * pS6)))))));
8917303c62SDavid Schultz }
9017303c62SDavid Schultz
9117303c62SDavid Schultz static inline long double
Q(long double x)9217303c62SDavid Schultz Q(long double x)
9317303c62SDavid Schultz {
9417303c62SDavid Schultz
9517303c62SDavid Schultz return (1.0 + x * (qS1 + x * (qS2 + x * (qS3 + x * (qS4 + x * qS5)))));
9617303c62SDavid Schultz }
9717303c62SDavid Schultz
9817303c62SDavid Schultz static inline long double
T_even(long double x)9917303c62SDavid Schultz T_even(long double x)
10017303c62SDavid Schultz {
10117303c62SDavid Schultz
10217303c62SDavid Schultz return (aT[0] + x * (aT[2] + x * (aT[4] + x * (aT[6] + x * \
10317303c62SDavid Schultz (aT[8] + x * (aT[10] + x * aT[12]))))));
10417303c62SDavid Schultz }
10517303c62SDavid Schultz
10617303c62SDavid Schultz static inline long double
T_odd(long double x)10717303c62SDavid Schultz T_odd(long double x)
10817303c62SDavid Schultz {
10917303c62SDavid Schultz
11017303c62SDavid Schultz return (aT[1] + x * (aT[3] + x * (aT[5] + x * (aT[7] + x * \
11117303c62SDavid Schultz (aT[9] + x * aT[11])))));
11217303c62SDavid Schultz }
1131192a80eSDavid Schultz
1141192a80eSDavid Schultz #endif
115