1#ifndef _TGMATH_H
2#define _TGMATH_H
3
4/*
5the return types are only correct with gcc (__GNUC__)
6otherwise they are long double or long double complex
7
8the long double version of a function is never chosen when
9sizeof(double) == sizeof(long double)
10(but the return type is set correctly with gcc)
11*/
12
13#include <math.h>
14#include <complex.h>
15
16#define __IS_FP(x) (sizeof((x)+1ULL) == sizeof((x)+1.0f))
17#define __IS_CX(x) (__IS_FP(x) && sizeof(x) == sizeof((x)+I))
18#define __IS_REAL(x) (__IS_FP(x) && 2*sizeof(x) == sizeof((x)+I))
19
20#define __FLT(x) (__IS_REAL(x) && sizeof(x) == sizeof(float))
21#define __LDBL(x) (__IS_REAL(x) && sizeof(x) == sizeof(long double) && sizeof(long double) != sizeof(double))
22
23#define __FLTCX(x) (__IS_CX(x) && sizeof(x) == sizeof(float complex))
24#define __DBLCX(x) (__IS_CX(x) && sizeof(x) == sizeof(double complex))
25#define __LDBLCX(x) (__IS_CX(x) && sizeof(x) == sizeof(long double complex) && sizeof(long double) != sizeof(double))
26
27/* return type */
28
29#ifdef __GNUC__
30/*
31the result must be casted to the right type
32(otherwise the result type is determined by the conversion
33rules applied to all the function return types so it is long
34double or long double complex except for integral functions)
35
36this cannot be done in c99, so the typeof gcc extension is
37used and that the type of ?: depends on wether an operand is
38a null pointer constant or not
39(in c11 _Generic can be used)
40
41the c arguments below must be integer constant expressions
42so they can be in null pointer constants
43(__IS_FP above was carefully chosen this way)
44*/
45/* if c then t else void */
46#define __type1(c,t) __typeof__(*(0?(t*)0:(void*)!(c)))
47/* if c then t1 else t2 */
48#define __type2(c,t1,t2) __typeof__(*(0?(__type1(c,t1)*)0:(__type1(!(c),t2)*)0))
49/* cast to double when x is integral, otherwise use typeof(x) */
50#define __RETCAST(x) ( \
51 __type2(__IS_FP(x), __typeof__(x), double))
52/* 2 args case, should work for complex types (cpow) */
53#define __RETCAST_2(x, y) ( \
54 __type2(__IS_FP(x) && __IS_FP(y), \
55 __typeof__((x)+(y)), \
56 __typeof__((x)+(y)+1.0)))
57/* 3 args case (fma only) */
58#define __RETCAST_3(x, y, z) ( \
59 __type2(__IS_FP(x) && __IS_FP(y) && __IS_FP(z), \
60 __typeof__((x)+(y)+(z)), \
61 __typeof__((x)+(y)+(z)+1.0)))
62/* drop complex from the type of x */
63/* TODO: wrong when sizeof(long double)==sizeof(double) */
64#define __RETCAST_REAL(x) ( \
65 __type2(__IS_FP(x) && sizeof((x)+I) == sizeof(float complex), float, \
66 __type2(sizeof((x)+1.0+I) == sizeof(double complex), double, \
67 long double)))
68/* add complex to the type of x */
69#define __RETCAST_CX(x) (__typeof__(__RETCAST(x)0+I))
70#else
71#define __RETCAST(x)
72#define __RETCAST_2(x, y)
73#define __RETCAST_3(x, y, z)
74#define __RETCAST_REAL(x)
75#define __RETCAST_CX(x)
76#endif
77
78/* function selection */
79
80#define __tg_real_nocast(fun, x) ( \
81 __FLT(x) ? fun ## f (x) : \
82 __LDBL(x) ? fun ## l (x) : \
83 fun(x) )
84
85#define __tg_real(fun, x) (__RETCAST(x)__tg_real_nocast(fun, x))
86
87#define __tg_real_2_1(fun, x, y) (__RETCAST(x)( \
88 __FLT(x) ? fun ## f (x, y) : \
89 __LDBL(x) ? fun ## l (x, y) : \
90 fun(x, y) ))
91
92#define __tg_real_2(fun, x, y) (__RETCAST_2(x, y)( \
93 __FLT(x) && __FLT(y) ? fun ## f (x, y) : \
94 __LDBL((x)+(y)) ? fun ## l (x, y) : \
95 fun(x, y) ))
96
97#define __tg_complex(fun, x) (__RETCAST_CX(x)( \
98 __FLTCX((x)+I) && __IS_FP(x) ? fun ## f (x) : \
99 __LDBLCX((x)+I) ? fun ## l (x) : \
100 fun(x) ))
101
102#define __tg_complex_retreal(fun, x) (__RETCAST_REAL(x)( \
103 __FLTCX((x)+I) && __IS_FP(x) ? fun ## f (x) : \
104 __LDBLCX((x)+I) ? fun ## l (x) : \
105 fun(x) ))
106
107#define __tg_real_complex(fun, x) (__RETCAST(x)( \
108 __FLTCX(x) ? c ## fun ## f (x) : \
109 __DBLCX(x) ? c ## fun (x) : \
110 __LDBLCX(x) ? c ## fun ## l (x) : \
111 __FLT(x) ? fun ## f (x) : \
112 __LDBL(x) ? fun ## l (x) : \
113 fun(x) ))
114
115/* special cases */
116
117#define __tg_real_remquo(x, y, z) (__RETCAST_2(x, y)( \
118 __FLT(x) && __FLT(y) ? remquof(x, y, z) : \
119 __LDBL((x)+(y)) ? remquol(x, y, z) : \
120 remquo(x, y, z) ))
121
122#define __tg_real_fma(x, y, z) (__RETCAST_3(x, y, z)( \
123 __FLT(x) && __FLT(y) && __FLT(z) ? fmaf(x, y, z) : \
124 __LDBL((x)+(y)+(z)) ? fmal(x, y, z) : \
125 fma(x, y, z) ))
126
127#define __tg_real_complex_pow(x, y) (__RETCAST_2(x, y)( \
128 __FLTCX((x)+(y)) && __IS_FP(x) && __IS_FP(y) ? cpowf(x, y) : \
129 __FLTCX((x)+(y)) ? cpow(x, y) : \
130 __DBLCX((x)+(y)) ? cpow(x, y) : \
131 __LDBLCX((x)+(y)) ? cpowl(x, y) : \
132 __FLT(x) && __FLT(y) ? powf(x, y) : \
133 __LDBL((x)+(y)) ? powl(x, y) : \
134 pow(x, y) ))
135
136#define __tg_real_complex_fabs(x) (__RETCAST_REAL(x)( \
137 __FLTCX(x) ? cabsf(x) : \
138 __DBLCX(x) ? cabs(x) : \
139 __LDBLCX(x) ? cabsl(x) : \
140 __FLT(x) ? fabsf(x) : \
141 __LDBL(x) ? fabsl(x) : \
142 fabs(x) ))
143
144/* suppress any macros in math.h or complex.h */
145
146#undef acos
147#undef acosh
148#undef asin
149#undef asinh
150#undef atan
151#undef atan2
152#undef atanh
153#undef carg
154#undef cbrt
155#undef ceil
156#undef cimag
157#undef conj
158#undef copysign
159#undef cos
160#undef cosh
161#undef cproj
162#undef creal
163#undef erf
164#undef erfc
165#undef exp
166#undef exp2
167#undef expm1
168#undef fabs
169#undef fdim
170#undef floor
171#undef fma
172#undef fmax
173#undef fmin
174#undef fmod
175#undef frexp
176#undef hypot
177#undef ilogb
178#undef ldexp
179#undef lgamma
180#undef llrint
181#undef llround
182#undef log
183#undef log10
184#undef log1p
185#undef log2
186#undef logb
187#undef lrint
188#undef lround
189#undef nearbyint
190#undef nextafter
191#undef nexttoward
192#undef pow
193#undef remainder
194#undef remquo
195#undef rint
196#undef round
197#undef scalbln
198#undef scalbn
199#undef sin
200#undef sinh
201#undef sqrt
202#undef tan
203#undef tanh
204#undef tgamma
205#undef trunc
206
207/* tg functions */
208
209#define acos(x) __tg_real_complex(acos, (x))
210#define acosh(x) __tg_real_complex(acosh, (x))
211#define asin(x) __tg_real_complex(asin, (x))
212#define asinh(x) __tg_real_complex(asinh, (x))
213#define atan(x) __tg_real_complex(atan, (x))
214#define atan2(x,y) __tg_real_2(atan2, (x), (y))
215#define atanh(x) __tg_real_complex(atanh, (x))
216#define carg(x) __tg_complex_retreal(carg, (x))
217#define cbrt(x) __tg_real(cbrt, (x))
218#define ceil(x) __tg_real(ceil, (x))
219#define cimag(x) __tg_complex_retreal(cimag, (x))
220#define conj(x) __tg_complex(conj, (x))
221#define copysign(x,y) __tg_real_2(copysign, (x), (y))
222#define cos(x) __tg_real_complex(cos, (x))
223#define cosh(x) __tg_real_complex(cosh, (x))
224#define cproj(x) __tg_complex(cproj, (x))
225#define creal(x) __tg_complex_retreal(creal, (x))
226#define erf(x) __tg_real(erf, (x))
227#define erfc(x) __tg_real(erfc, (x))
228#define exp(x) __tg_real_complex(exp, (x))
229#define exp2(x) __tg_real(exp2, (x))
230#define expm1(x) __tg_real(expm1, (x))
231#define fabs(x) __tg_real_complex_fabs(x)
232#define fdim(x,y) __tg_real_2(fdim, (x), (y))
233#define floor(x) __tg_real(floor, (x))
234#define fma(x,y,z) __tg_real_fma((x), (y), (z))
235#define fmax(x,y) __tg_real_2(fmax, (x), (y))
236#define fmin(x,y) __tg_real_2(fmin, (x), (y))
237#define fmod(x,y) __tg_real_2(fmod, (x), (y))
238#define frexp(x,y) __tg_real_2_1(frexp, (x), (y))
239#define hypot(x,y) __tg_real_2(hypot, (x), (y))
240#define ilogb(x) __tg_real_nocast(ilogb, (x))
241#define ldexp(x,y) __tg_real_2_1(ldexp, (x), (y))
242#define lgamma(x) __tg_real(lgamma, (x))
243#define llrint(x) __tg_real_nocast(llrint, (x))
244#define llround(x) __tg_real_nocast(llround, (x))
245#define log(x) __tg_real_complex(log, (x))
246#define log10(x) __tg_real(log10, (x))
247#define log1p(x) __tg_real(log1p, (x))
248#define log2(x) __tg_real(log2, (x))
249#define logb(x) __tg_real(logb, (x))
250#define lrint(x) __tg_real_nocast(lrint, (x))
251#define lround(x) __tg_real_nocast(lround, (x))
252#define nearbyint(x) __tg_real(nearbyint, (x))
253#define nextafter(x,y) __tg_real_2(nextafter, (x), (y))
254#define nexttoward(x,y) __tg_real_2(nexttoward, (x), (y))
255#define pow(x,y) __tg_real_complex_pow((x), (y))
256#define remainder(x,y) __tg_real_2(remainder, (x), (y))
257#define remquo(x,y,z) __tg_real_remquo((x), (y), (z))
258#define rint(x) __tg_real(rint, (x))
259#define round(x) __tg_real(round, (x))
260#define scalbln(x,y) __tg_real_2_1(scalbln, (x), (y))
261#define scalbn(x,y) __tg_real_2_1(scalbn, (x), (y))
262#define sin(x) __tg_real_complex(sin, (x))
263#define sinh(x) __tg_real_complex(sinh, (x))
264#define sqrt(x) __tg_real_complex(sqrt, (x))
265#define tan(x) __tg_real_complex(tan, (x))
266#define tanh(x) __tg_real_complex(tanh, (x))
267#define tgamma(x) __tg_real(tgamma, (x))
268#define trunc(x) __tg_real(trunc, (x))
269
270#endif