1/*
2 * This Source Code Form is subject to the terms of the Mozilla Public
3 * License, v. 2.0. If a copy of the MPL was not distributed with this
4 * file, You can obtain one at http://mozilla.org/MPL/2.0/.
5 *
6 * Copyright 1997 - July 2008 CWI, August 2008 - 2019 MonetDB B.V.
7 */
8
9/*
10 * N.J. Nes, M. Kersten
11 * 07/01/1996
12 * The math module
13 * This module contains the math commands. The implementation is very simply,
14 * the c math library functions are called. See for documentation the
15 * ANSI-C/POSIX manuals of the equaly named functions.
16 *
17 * NOTE: the operand itself is being modified, rather than that we produce
18 * a new BAT. This to save the expensive copying.
19 */
20#include "monetdb_config.h"
21#include "mmath.h"
22#include <fenv.h>
23#include "mmath_private.h"
24#ifndef FE_INVALID
25#define FE_INVALID 0
26#endif
27#ifndef FE_DIVBYZERO
28#define FE_DIVBYZERO 0
29#endif
30#ifndef FE_OVERFLOW
31#define FE_OVERFLOW 0
32#endif
33
34#define cot(x) (1 / tan(x))
35#define radians(x) ((x) * 3.14159265358979323846 / 180.0)
36#define degrees(x) ((x) * 180.0 / 3.14159265358979323846)
37
38double
39logbs(double x, double base)
40{
41 return log(x) / log(base);
42}
43
44float
45logbsf(float x, float base)
46{
47 return logf(x) / logf(base);
48}
49
50#define unopbaseM5(NAME, FUNC, TYPE) \
51str \
52MATHunary##NAME##TYPE(TYPE *res , const TYPE *a) \
53{ \
54 if (is_##TYPE##_nil(*a)) { \
55 *res = TYPE##_nil; \
56 } else { \
57 double a1 = *a, r; \
58 int e = 0, ex = 0; \
59 errno = 0; \
60 feclearexcept(FE_ALL_EXCEPT); \
61 r = FUNC(a1); \
62 if ((e = errno) != 0 || \
63 (ex = fetestexcept(FE_INVALID | FE_DIVBYZERO | \
64 FE_OVERFLOW)) != 0) { \
65 const char *err; \
66 if (e) { \
67 err = strerror(e); \
68 } else if (ex & FE_DIVBYZERO) \
69 err = "Divide by zero"; \
70 else if (ex & FE_OVERFLOW) \
71 err = "Overflow"; \
72 else \
73 err = "Invalid result"; \
74 throw(MAL, "mmath." #FUNC, "Math exception: %s", err); \
75 } \
76 *res = (TYPE) r; \
77 } \
78 return MAL_SUCCEED; \
79}
80
81#define unopM5(NAME, FUNC) \
82 unopbaseM5(NAME, FUNC, dbl) \
83 unopbaseM5(NAME, FUNC, flt)
84
85#define binopbaseM5(NAME, FUNC, TYPE) \
86str \
87MATHbinary##NAME##TYPE(TYPE *res, const TYPE *a, const TYPE *b) \
88{ \
89 if (is_##TYPE##_nil(*a) || is_##TYPE##_nil(*b)) { \
90 *res = TYPE##_nil; \
91 } else { \
92 double r1, a1 = *a, b1 = *b; \
93 int e = 0, ex = 0; \
94 errno = 0; \
95 feclearexcept(FE_ALL_EXCEPT); \
96 r1 = FUNC(a1, b1); \
97 if ((e = errno) != 0 || \
98 (ex = fetestexcept(FE_INVALID | FE_DIVBYZERO | \
99 FE_OVERFLOW)) != 0) { \
100 const char *err; \
101 if (e) { \
102 err = strerror(e); \
103 } else if (ex & FE_DIVBYZERO) \
104 err = "Divide by zero"; \
105 else if (ex & FE_OVERFLOW) \
106 err = "Overflow"; \
107 else \
108 err = "Invalid result"; \
109 throw(MAL, "mmath." #FUNC, "Math exception: %s", err); \
110 } \
111 *res= (TYPE) r1; \
112 } \
113 return MAL_SUCCEED; \
114}
115
116#define unopM5NOT(NAME, FUNC) \
117str \
118MATHunary##NAME##dbl(dbl *res , const dbl *a) \
119{ \
120 (void)res; \
121 (void)a; \
122 throw(MAL, "mmath." #FUNC, SQLSTATE(0A000) PROGRAM_NYI); \
123} \
124str \
125MATHunary##NAME##flt(flt *res , const flt *a) \
126{ \
127 (void)res; \
128 (void)a; \
129 throw(MAL, "mmath." #FUNC, SQLSTATE(0A000) PROGRAM_NYI); \
130}
131
132#define binopM5(NAME, FUNC) \
133 binopbaseM5(NAME, FUNC, dbl) \
134 binopbaseM5(NAME, FUNC, flt)
135
136#define roundM5(TYPE) \
137str \
138MATHbinary_ROUND##TYPE(TYPE *res, const TYPE *x, const int *y) \
139{ \
140 if (is_##TYPE##_nil(*x) || is_int_nil(*y)) { \
141 *res = TYPE##_nil; \
142 } else { \
143 dbl factor = pow(10,*y), integral; \
144 dbl tmp = *y > 0 ? modf(*x, &integral) : *x; \
145 \
146 tmp *= factor; \
147 if (tmp >= 0) \
148 tmp = floor(tmp + 0.5); \
149 else \
150 tmp = ceil(tmp - 0.5); \
151 tmp /= factor; \
152 \
153 if (*y > 0) \
154 tmp += integral; \
155 \
156 *res = (TYPE) tmp; \
157 } \
158 return MAL_SUCCEED; \
159}
160
161
162unopM5(_ACOS,acos)
163unopM5(_ASIN,asin)
164unopM5(_ATAN,atan)
165binopM5(_ATAN2,atan2)
166unopM5(_COS,cos)
167unopM5(_SIN,sin)
168unopM5(_TAN,tan)
169unopM5(_COT,cot)
170unopM5(_RADIANS,radians)
171unopM5(_DEGREES,degrees)
172
173unopM5(_COSH,cosh)
174unopM5(_SINH,sinh)
175unopM5(_TANH,tanh)
176
177unopM5(_EXP,exp)
178unopM5(_LOG,log)
179unopM5(_LOG10,log10)
180unopM5(_LOG2,log2)
181
182binopM5(_LOG,logbs)
183
184binopM5(_POW,pow)
185unopM5(_SQRT,sqrt)
186#ifdef HAVE_CBRT
187unopM5(_CBRT,cbrt)
188#else
189unopM5NOT(_CBRT,cbrt)
190#endif
191
192unopM5(_CEIL,ceil)
193unopM5(_FLOOR,floor)
194
195str
196MATHunary_FABSdbl(dbl *res , const dbl *a)
197{
198 *res = is_dbl_nil(*a) ? dbl_nil : fabs(*a);
199 return MAL_SUCCEED;
200}
201
202roundM5(dbl)
203roundM5(flt)
204
205str
206MATHunary_ISNAN(bit *res, const dbl *a)
207{
208 if (is_dbl_nil(*a)) {
209 *res = bit_nil;
210 } else {
211 *res = isnan(*a) != 0;
212 }
213 return MAL_SUCCEED;
214}
215
216str
217MATHunary_ISINF(int *res, const dbl *a)
218{
219 if (is_dbl_nil(*a)) {
220 *res = int_nil;
221 } else {
222 if (isinf(*a)) {
223 *res = (*a < 0.0) ? -1 : 1;
224 } else {
225 *res = 0;
226 }
227 }
228 return MAL_SUCCEED;
229}
230
231str
232MATHunary_FINITE(bit *res, const dbl *a)
233{
234 if (is_dbl_nil(*a)) {
235 *res = bit_nil;
236 } else {
237 *res = isfinite(*a) != 0;
238 }
239 return MAL_SUCCEED;
240}
241
242#include "xoshiro256starstar.h"
243
244/* global pseudo random generator state */
245static random_state_engine mmath_rse;
246/* serialize access to state */
247static MT_Lock mmath_rse_lock = MT_LOCK_INITIALIZER("mmath_rse_lock");
248
249str
250MATHprelude(void *ret)
251{
252 (void) ret;
253 init_random_state_engine(mmath_rse, (uint64_t) GDKusec());
254 return MAL_SUCCEED;
255}
256
257str
258MATHrandint(int *res)
259{
260#ifdef STATIC_CODE_ANALYSIS
261 *res = 0;
262#else
263 MT_lock_set(&mmath_rse_lock);
264 *res = (int) (next(mmath_rse) >> 33);
265 MT_lock_unset(&mmath_rse_lock);
266#endif
267 return MAL_SUCCEED;
268}
269
270str
271MATHrandintarg(int *res, const int *dummy)
272{
273 (void) dummy;
274#ifdef STATIC_CODE_ANALYSIS
275 *res = 0;
276#else
277 MT_lock_set(&mmath_rse_lock);
278 *res = (int) (next(mmath_rse) >> 33);
279 MT_lock_unset(&mmath_rse_lock);
280#endif
281 return MAL_SUCCEED;
282}
283
284str
285MATHsrandint(void *ret, const int *seed)
286{
287 (void) ret;
288 MT_lock_set(&mmath_rse_lock);
289 init_random_state_engine(mmath_rse, (uint64_t) *seed);
290 MT_lock_unset(&mmath_rse_lock);
291 return MAL_SUCCEED;
292}
293
294str
295MATHsqlrandint(int *res, const int *seed)
296{
297#ifdef STATIC_CODE_ANALYSIS
298 *res = 0;
299#else
300 MT_lock_set(&mmath_rse_lock);
301 init_random_state_engine(mmath_rse, (uint64_t) *seed);
302 *res = (int) (next(mmath_rse) >> 33);
303 MT_lock_unset(&mmath_rse_lock);
304#endif
305 return MAL_SUCCEED;
306}
307
308str
309MATHpi(dbl *pi)
310{
311 *pi = 3.14159265358979323846;
312 return MAL_SUCCEED;
313}
314
315