1/* origin: FreeBSD /usr/src/lib/msun/src/math_private.h */
2/*
3 * ====================================================
4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
5 *
6 * Developed at SunPro, a Sun Microsystems, Inc. business.
7 * Permission to use, copy, modify, and distribute this
8 * software is freely granted, provided that this notice
9 * is preserved.
10 * ====================================================
11 */
12
13#pragma once
14
15#include <complex.h>
16#include <endian.h>
17#include <float.h>
18#include <math.h>
19#include <stdint.h>
20
21// Clang does not support the C99 rounding mode pragma. Support seems
22// unlikely to be coming soon, but for reference the clang/llvm bug
23// tracking this fact may be found at:
24// https://llvm.org/bugs/show_bug.cgi?id=8100
25#define PRAGMA_STDC_FENV_ACCESS_ON
26
27// GCC complains about the STDC CX_LIMITED_RANGE pragma.
28#define PRAGMA_STDC_CX_LIMITED_RANGE_ON
29
30#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
31#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384 && __BYTE_ORDER == __LITTLE_ENDIAN
32union ldshape {
33    long double f;
34    struct {
35        uint64_t m;
36        uint16_t se;
37    } i;
38};
39#elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384 && __BYTE_ORDER == __LITTLE_ENDIAN
40union ldshape {
41    long double f;
42    struct {
43        uint64_t lo;
44        uint32_t mid;
45        uint16_t top;
46        uint16_t se;
47    } i;
48    struct {
49        uint64_t lo;
50        uint64_t hi;
51    } i2;
52};
53#elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384 && __BYTE_ORDER == __BIG_ENDIAN
54union ldshape {
55    long double f;
56    struct {
57        uint16_t se;
58        uint16_t top;
59        uint32_t mid;
60        uint64_t lo;
61    } i;
62    struct {
63        uint64_t hi;
64        uint64_t lo;
65    } i2;
66};
67#else
68#error Unsupported long double representation
69#endif
70
71#define FORCE_EVAL(x)                             \
72    do {                                          \
73        if (sizeof(x) == sizeof(float)) {         \
74            volatile float __x;                   \
75            __x = (x);                            \
76            (void)__x;                            \
77        } else if (sizeof(x) == sizeof(double)) { \
78            volatile double __x;                  \
79            __x = (x);                            \
80            (void)__x;                            \
81        } else {                                  \
82            volatile long double __x;             \
83            __x = (x);                            \
84            (void)__x;                            \
85        }                                         \
86    } while (0)
87
88/* Get two 32 bit ints from a double.  */
89#define EXTRACT_WORDS(hi, lo, d) \
90    do {                         \
91        union {                  \
92            double f;            \
93            uint64_t i;          \
94        } __u;                   \
95        __u.f = (d);             \
96        (hi) = __u.i >> 32;      \
97        (lo) = (uint32_t)__u.i;  \
98    } while (0)
99
100/* Get the more significant 32 bit int from a double.  */
101#define GET_HIGH_WORD(hi, d) \
102    do {                     \
103        union {              \
104            double f;        \
105            uint64_t i;      \
106        } __u;               \
107        __u.f = (d);         \
108        (hi) = __u.i >> 32;  \
109    } while (0)
110
111/* Get the less significant 32 bit int from a double.  */
112#define GET_LOW_WORD(lo, d)     \
113    do {                        \
114        union {                 \
115            double f;           \
116            uint64_t i;         \
117        } __u;                  \
118        __u.f = (d);            \
119        (lo) = (uint32_t)__u.i; \
120    } while (0)
121
122/* Set a double from two 32 bit ints.  */
123#define INSERT_WORDS(d, hi, lo)                          \
124    do {                                                 \
125        union {                                          \
126            double f;                                    \
127            uint64_t i;                                  \
128        } __u;                                           \
129        __u.i = ((uint64_t)(hi) << 32) | (uint32_t)(lo); \
130        (d) = __u.f;                                     \
131    } while (0)
132
133/* Set the more significant 32 bits of a double from an int.  */
134#define SET_HIGH_WORD(d, hi)           \
135    do {                               \
136        union {                        \
137            double f;                  \
138            uint64_t i;                \
139        } __u;                         \
140        __u.f = (d);                   \
141        __u.i &= 0xffffffff;           \
142        __u.i |= (uint64_t)(hi) << 32; \
143        (d) = __u.f;                   \
144    } while (0)
145
146/* Set the less significant 32 bits of a double from an int.  */
147#define SET_LOW_WORD(d, lo)             \
148    do {                                \
149        union {                         \
150            double f;                   \
151            uint64_t i;                 \
152        } __u;                          \
153        __u.f = (d);                    \
154        __u.i &= 0xffffffff00000000ull; \
155        __u.i |= (uint32_t)(lo);        \
156        (d) = __u.f;                    \
157    } while (0)
158
159/* Get a 32 bit int from a float.  */
160#define GET_FLOAT_WORD(w, d) \
161    do {                     \
162        union {              \
163            float f;         \
164            uint32_t i;      \
165        } __u;               \
166        __u.f = (d);         \
167        (w) = __u.i;         \
168    } while (0)
169
170/* Set a float from a 32 bit int.  */
171#define SET_FLOAT_WORD(d, w) \
172    do {                     \
173        union {              \
174            float f;         \
175            uint32_t i;      \
176        } __u;               \
177        __u.i = (w);         \
178        (d) = __u.f;         \
179    } while (0)
180
181#undef __CMPLX
182#undef CMPLX
183#undef CMPLXF
184#undef CMPLXL
185
186#define __CMPLX(x, y, t)    \
187    ((union {               \
188         _Complex t __z;    \
189         t __xy[2];         \
190     }){.__xy = {(x), (y)}} \
191         .__z)
192
193#define CMPLX(x, y) __CMPLX(x, y, double)
194#define CMPLXF(x, y) __CMPLX(x, y, float)
195#define CMPLXL(x, y) __CMPLX(x, y, long double)
196
197/* fdlibm kernel functions */
198
199int __rem_pio2_large(double*, double*, int, int, int);
200
201int __rem_pio2(double, double*);
202double __sin(double, double, int);
203double __cos(double, double);
204double __tan(double, double, int);
205double __expo2(double);
206double complex __ldexp_cexp(double complex, int);
207
208int __rem_pio2f(float, double*);
209float __sindf(double);
210float __cosdf(double);
211float __tandf(double, int);
212float __expo2f(float);
213float complex __ldexp_cexpf(float complex, int);
214
215int __rem_pio2l(long double, long double*);
216long double __sinl(long double, long double, int);
217long double __cosl(long double, long double);
218long double __tanl(long double, long double, int);
219
220/* polynomial evaluation */
221long double __polevll(long double, const long double*, int);
222long double __p1evll(long double, const long double*, int);
223