Line data Source code
1 : /*
2 : * SPDX-License-Identifier: MPL-2.0
3 : *
4 : * This Source Code Form is subject to the terms of the Mozilla Public
5 : * License, v. 2.0. If a copy of the MPL was not distributed with this
6 : * file, You can obtain one at http://mozilla.org/MPL/2.0/.
7 : *
8 : * Copyright 2024 MonetDB Foundation;
9 : * Copyright August 2008 - 2023 MonetDB B.V.;
10 : * Copyright 1997 - July 2008 CWI.
11 : */
12 :
13 : /*
14 : * N.J. Nes, M. Kersten
15 : * 07/01/1996
16 : * The math module
17 : * This module contains the math commands. The implementation is very simply,
18 : * the c math library functions are called. See for documentation the
19 : * ANSI-C/POSIX manuals of the equally named functions.
20 : *
21 : * NOTE: the operand itself is being modified, rather than that we produce
22 : * a new BAT. This to save the expensive copying.
23 : */
24 : #include "monetdb_config.h"
25 : #include "mal.h"
26 : #include "mal_exception.h"
27 : #include <fenv.h>
28 : #include "mmath_private.h"
29 :
30 : double
31 10 : cot(double x)
32 : {
33 10 : return 1.0 / tan(x);
34 : }
35 :
36 : float
37 0 : cotf(float x)
38 : {
39 0 : return 1.0f / tanf(x);
40 : }
41 :
42 : double
43 44 : radians(double x)
44 : {
45 44 : return x * (M_PI / 180.0);
46 : }
47 :
48 : float
49 0 : radiansf(float x)
50 : {
51 0 : return x * (M_PIF / 180.0f);
52 : }
53 :
54 : double
55 26 : degrees(double x)
56 : {
57 26 : return x * (180.0 / M_PI);
58 : }
59 :
60 : float
61 0 : degreesf(float x)
62 : {
63 0 : return x * (180.0f / M_PIF);
64 : }
65 :
66 : double
67 8 : logbs(double base, double x)
68 : {
69 8 : if (base == 1) {
70 2 : feraiseexcept(FE_DIVBYZERO);
71 2 : return INFINITY;
72 : }
73 6 : return log(x) / log(base);
74 : }
75 :
76 : float
77 0 : logbsf(float base, float x)
78 : {
79 0 : if (base == 1) {
80 0 : feraiseexcept(FE_DIVBYZERO);
81 0 : return INFINITY;
82 : }
83 0 : return logf(x) / logf(base);
84 : }
85 :
86 : #define unopbaseM5(NAME, FUNC, TYPE) \
87 : static str \
88 : MATHunary##NAME##TYPE(TYPE *res, const TYPE *a) \
89 : { \
90 : if (is_##TYPE##_nil(*a)) { \
91 : *res = TYPE##_nil; \
92 : } else { \
93 : int e = 0, ex = 0; \
94 : errno = 0; \
95 : feclearexcept(FE_ALL_EXCEPT); \
96 : *res = FUNC(*a); \
97 : if ((e = errno) != 0 || \
98 : (ex = fetestexcept(FE_INVALID | FE_DIVBYZERO | \
99 : FE_OVERFLOW)) != 0) { \
100 : const char *err; \
101 : char buf[128]; \
102 : if (e) { \
103 : err = GDKstrerror(e, buf, 128); \
104 : } else if (ex & FE_DIVBYZERO) \
105 : err = "Divide by zero"; \
106 : else if (ex & FE_OVERFLOW) \
107 : err = "Overflow"; \
108 : else \
109 : err = "Invalid result"; \
110 : throw(MAL, "mmath." #FUNC, "Math exception: %s", err); \
111 : } \
112 : } \
113 : return MAL_SUCCEED; \
114 : }
115 :
116 : #define unopM5(NAME, FUNC) \
117 : unopbaseM5(NAME, FUNC, dbl) \
118 : unopbaseM5(NAME, FUNC##f, flt)
119 :
120 : #define binopbaseM5(NAME, FUNC, TYPE) \
121 : static str \
122 : MATHbinary##NAME##TYPE(TYPE *res, const TYPE *a, const TYPE *b) \
123 : { \
124 : if (is_##TYPE##_nil(*a) || is_##TYPE##_nil(*b)) { \
125 : *res = TYPE##_nil; \
126 : } else { \
127 : int e = 0, ex = 0; \
128 : errno = 0; \
129 : feclearexcept(FE_ALL_EXCEPT); \
130 : *res = FUNC(*a, *b); \
131 : if ((e = errno) != 0 || \
132 : (ex = fetestexcept(FE_INVALID | FE_DIVBYZERO | \
133 : FE_OVERFLOW)) != 0) { \
134 : const char *err; \
135 : char buf[128]; \
136 : if (e) { \
137 : err = GDKstrerror(e, buf, 128); \
138 : } else if (ex & FE_DIVBYZERO) \
139 : err = "Divide by zero"; \
140 : else if (ex & FE_OVERFLOW) \
141 : err = "Overflow"; \
142 : else \
143 : err = "Invalid result"; \
144 : throw(MAL, "mmath." #FUNC, "Math exception: %s", err); \
145 : } \
146 : } \
147 : return MAL_SUCCEED; \
148 : }
149 :
150 : #define binopM5(NAME, FUNC) \
151 : binopbaseM5(NAME, FUNC, dbl) \
152 : binopbaseM5(NAME, FUNC##f, flt)
153 :
154 : #define roundM5(TYPE) \
155 : static str \
156 : MATHbinary_ROUND##TYPE(TYPE *res, const TYPE *x, const int *y) \
157 : { \
158 : if (is_##TYPE##_nil(*x) || is_int_nil(*y)) { \
159 : *res = TYPE##_nil; \
160 : } else { \
161 : dbl factor = pow(10,*y), integral; \
162 : dbl tmp = *y > 0 ? modf(*x, &integral) : *x; \
163 : \
164 : tmp *= factor; \
165 : if (tmp >= 0) \
166 : tmp = floor(tmp + 0.5); \
167 : else \
168 : tmp = ceil(tmp - 0.5); \
169 : tmp /= factor; \
170 : \
171 : if (*y > 0) \
172 : tmp += integral; \
173 : \
174 : *res = (TYPE) tmp; \
175 : } \
176 : return MAL_SUCCEED; \
177 : }
178 :
179 :
180 3 : unopM5(_ACOS, acos)
181 4 : unopM5(_ASIN, asin)
182 3 : unopM5(_ATAN, atan)
183 23 : unopM5(_COS, cos)
184 10015 : unopM5(_SIN, sin)
185 2 : unopM5(_TAN, tan)
186 8 : unopM5(_COT, cot)
187 31 : unopM5(_RADIANS, radians)
188 12 : unopM5(_DEGREES, degrees)
189 2 : unopM5(_COSH, cosh)
190 2 : unopM5(_SINH, sinh)
191 2 : unopM5(_TANH, tanh)
192 6 : unopM5(_EXP, exp)
193 4 : unopM5(_LOG, log)
194 2 : unopM5(_LOG10, log10)
195 1 : unopM5(_LOG2, log2)
196 7 : unopM5(_SQRT, sqrt)
197 4 : unopM5(_CBRT, cbrt)
198 9 : unopM5(_CEIL, ceil)
199 14 : unopM5(_FLOOR, floor)
200 3 : binopM5(_ATAN2, atan2)
201 15 : binopM5(_POW, pow)
202 8 : binopM5(_LOG, logbs)
203 20 : binopM5(_NEXTAFTER, nextafter)
204 :
205 : static str
206 0 : MATHunary_FABSdbl(dbl *res, const dbl *a)
207 : {
208 0 : *res = is_dbl_nil(*a) ? dbl_nil : fabs(*a);
209 0 : return MAL_SUCCEED;
210 : }
211 :
212 0 : roundM5(dbl)
213 0 : roundM5(flt)
214 : static str
215 0 : MATHunary_ISNAN(bit *res, const dbl *a)
216 : {
217 0 : if (is_dbl_nil(*a)) {
218 0 : *res = bit_nil;
219 : } else {
220 0 : *res = isnan(*a) != 0;
221 : }
222 0 : return MAL_SUCCEED;
223 : }
224 :
225 : static str
226 0 : MATHunary_ISINF(int *res, const dbl *a)
227 : {
228 0 : if (is_dbl_nil(*a)) {
229 0 : *res = int_nil;
230 : } else {
231 0 : if (isinf(*a)) {
232 0 : *res = (*a < 0.0) ? -1 : 1;
233 : } else {
234 0 : *res = 0;
235 : }
236 : }
237 0 : return MAL_SUCCEED;
238 : }
239 :
240 : static str
241 0 : MATHunary_FINITE(bit *res, const dbl *a)
242 : {
243 0 : if (is_dbl_nil(*a)) {
244 0 : *res = bit_nil;
245 : } else {
246 0 : *res = isfinite(*a) != 0;
247 : }
248 0 : return MAL_SUCCEED;
249 : }
250 :
251 : /* global pseudo random generator state */
252 : random_state_engine mmath_rse;
253 : /* serialize access to state */
254 : MT_Lock mmath_rse_lock = MT_LOCK_INITIALIZER(mmath_rse_lock);
255 :
256 : static str
257 330 : MATHprelude(void)
258 : {
259 330 : init_random_state_engine(mmath_rse, (uint64_t) GDKusec());
260 330 : return MAL_SUCCEED;
261 : }
262 :
263 : static str
264 10004 : MATHrandint(int *res)
265 : {
266 : #ifdef __COVERITY__
267 : *res = 0;
268 : #else
269 10004 : MT_lock_set(&mmath_rse_lock);
270 10004 : *res = (int) (next(mmath_rse) >> 33);
271 10004 : MT_lock_unset(&mmath_rse_lock);
272 : #endif
273 10004 : return MAL_SUCCEED;
274 : }
275 :
276 : static str
277 0 : MATHrandintarg(int *res, const int *dummy)
278 : {
279 0 : (void) dummy;
280 : #ifdef __COVERITY__
281 : *res = 0;
282 : #else
283 0 : MT_lock_set(&mmath_rse_lock);
284 0 : *res = (int) (next(mmath_rse) >> 33);
285 0 : MT_lock_unset(&mmath_rse_lock);
286 : #endif
287 0 : return MAL_SUCCEED;
288 : }
289 :
290 : static str
291 2 : MATHsrandint(void *ret, const int *seed)
292 : {
293 2 : (void) ret;
294 2 : MT_lock_set(&mmath_rse_lock);
295 2 : init_random_state_engine(mmath_rse, (uint64_t) * seed);
296 2 : MT_lock_unset(&mmath_rse_lock);
297 2 : return MAL_SUCCEED;
298 : }
299 :
300 : static str
301 2 : MATHsqlrandint(int *res, const int *seed)
302 : {
303 : #ifdef __COVERITY__
304 : (void) seed;
305 : *res = 0;
306 : #else
307 2 : MT_lock_set(&mmath_rse_lock);
308 2 : init_random_state_engine(mmath_rse, (uint64_t) * seed);
309 2 : *res = (int) (next(mmath_rse) >> 33);
310 2 : MT_lock_unset(&mmath_rse_lock);
311 : #endif
312 2 : return MAL_SUCCEED;
313 : }
314 :
315 : static str
316 84 : MATHpi(dbl *pi)
317 : {
318 84 : *pi = M_PI;
319 84 : return MAL_SUCCEED;
320 : }
321 :
322 : #include "mel.h"
323 : mel_func mmath_init_funcs[] = {
324 : command("mmath", "acos", MATHunary_ACOSflt, false, "", args(1,2, arg("",flt),arg("x",flt))),
325 : command("mmath", "acos", MATHunary_ACOSdbl, false, "The acos(x) function calculates the arc cosine of x, that is the \nvalue whose cosine is x. The value is returned in radians and is \nmathematically defined to be between 0 and PI (inclusive).", args(1,2, arg("",dbl),arg("x",dbl))),
326 : command("mmath", "asin", MATHunary_ASINflt, false, "", args(1,2, arg("",flt),arg("x",flt))),
327 : command("mmath", "asin", MATHunary_ASINdbl, false, "The asin(x) function calculates the arc sine of x, that is the value \nwhose sine is x. The value is returned in radians and is mathematically \ndefined to be between -PI/20 and -PI/2 (inclusive).", args(1,2, arg("",dbl),arg("x",dbl))),
328 : command("mmath", "atan", MATHunary_ATANflt, false, "", args(1,2, arg("",flt),arg("x",flt))),
329 : command("mmath", "atan", MATHunary_ATANdbl, false, "The atan(x) function calculates the arc tangent of x, that is the value \nwhose tangent is x. The value is returned in radians and is mathematically \ndefined to be between -PI/2 and PI/2 (inclusive).", args(1,2, arg("",dbl),arg("x",dbl))),
330 : command("mmath", "atan2", MATHbinary_ATAN2flt, false, "", args(1,3, arg("",flt),arg("x",flt),arg("y",flt))),
331 : command("mmath", "atan2", MATHbinary_ATAN2dbl, false, "The atan2(x,y) function calculates the arc tangent of the two \nvariables x and y. It is similar to calculating the arc\ntangent of y / x, except that the signs of both arguments are \nused to determine the quadrant of the result. The value is \nreturned in radians and is mathematically defined to be between \n-PI/2 and PI/2 (inclusive).", args(1,3, arg("",dbl),arg("x",dbl),arg("y",dbl))),
332 : command("mmath", "cos", MATHunary_COSflt, false, "", args(1,2, arg("",flt),arg("x",flt))),
333 : command("mmath", "cos", MATHunary_COSdbl, false, "The cos(x) function returns the cosine of x, where x is given in \nradians. The return value is between -1 and 1.", args(1,2, arg("",dbl),arg("x",dbl))),
334 : command("mmath", "sin", MATHunary_SINflt, false, "", args(1,2, arg("",flt),arg("x",flt))),
335 : command("mmath", "sin", MATHunary_SINdbl, false, "The sin(x) function returns the cosine of x, where x is given in \nradians. The return value is between -1 and 1.", args(1,2, arg("",dbl),arg("x",dbl))),
336 : command("mmath", "tan", MATHunary_TANflt, false, "", args(1,2, arg("",flt),arg("x",flt))),
337 : command("mmath", "tan", MATHunary_TANdbl, false, "The tan(x) function returns the tangent of x,\nwhere x is given in radians", args(1,2, arg("",dbl),arg("x",dbl))),
338 : command("mmath", "cot", MATHunary_COTflt, false, "", args(1,2, arg("",flt),arg("x",flt))),
339 : command("mmath", "cot", MATHunary_COTdbl, false, "The cot(x) function returns the Cotangent of x,\nwhere x is given in radians", args(1,2, arg("",dbl),arg("x",dbl))),
340 : command("mmath", "cosh", MATHunary_COSHflt, false, "", args(1,2, arg("",flt),arg("x",flt))),
341 : command("mmath", "cosh", MATHunary_COSHdbl, false, "The cosh() function returns the hyperbolic cosine of x, which is \ndefined mathematically as (exp(x) + exp(-x)) / 2.", args(1,2, arg("",dbl),arg("x",dbl))),
342 : command("mmath", "sinh", MATHunary_SINHflt, false, "", args(1,2, arg("",flt),arg("x",flt))),
343 : command("mmath", "sinh", MATHunary_SINHdbl, false, "The sinh() function returns the hyperbolic sine of x, which is \ndefined mathematically as (exp(x) - exp(-x)) / 2.", args(1,2, arg("",dbl),arg("x",dbl))),
344 : command("mmath", "tanh", MATHunary_TANHflt, false, "", args(1,2, arg("",flt),arg("x",flt))),
345 : command("mmath", "tanh", MATHunary_TANHdbl, false, "The tanh() function returns the hyperbolic tangent of x, which is \ndefined mathematically as sinh(x) / cosh(x).", args(1,2, arg("",dbl),arg("x",dbl))),
346 : command("mmath", "radians", MATHunary_RADIANSflt, false, "", args(1,2, arg("",flt),arg("x",flt))),
347 : command("mmath", "radians", MATHunary_RADIANSdbl, false, "The radians() function converts degrees into radians", args(1,2, arg("",dbl),arg("x",dbl))),
348 : command("mmath", "degrees", MATHunary_DEGREESflt, false, "", args(1,2, arg("",flt),arg("x",flt))),
349 : command("mmath", "degrees", MATHunary_DEGREESdbl, false, "The degrees() function converts radians into degrees", args(1,2, arg("",dbl),arg("x",dbl))),
350 : command("mmath", "exp", MATHunary_EXPflt, false, "", args(1,2, arg("",flt),arg("x",flt))),
351 : command("mmath", "exp", MATHunary_EXPdbl, false, "The exp(x) function returns the value of e (the base of \nnatural logarithms) raised to the power of x.", args(1,2, arg("",dbl),arg("x",dbl))),
352 : command("mmath", "log", MATHunary_LOGflt, false, "", args(1,2, arg("",flt),arg("x",flt))),
353 : command("mmath", "log", MATHunary_LOGdbl, false, "The log(x) function returns the natural logarithm of x.", args(1,2, arg("",dbl),arg("x",dbl))),
354 : command("mmath", "log2arg", MATHbinary_LOGflt, false, "The log(x) function returns the logarithm of x in the given base.", args(1,3, arg("",flt),arg("x",flt),arg("base",flt))),
355 : command("mmath", "log2arg", MATHbinary_LOGdbl, false, "The log(x) function returns the logarithm of x in the given base.", args(1,3, arg("",dbl),arg("x",dbl),arg("base",dbl))),
356 : command("mmath", "log10", MATHunary_LOG10flt, false, "", args(1,2, arg("",flt),arg("x",flt))),
357 : command("mmath", "log10", MATHunary_LOG10dbl, false, "The log10(x) function returns the base-10 logarithm of x.", args(1,2, arg("",dbl),arg("x",dbl))),
358 : command("mmath", "log2", MATHunary_LOG2flt, false, "", args(1,2, arg("",flt),arg("x",flt))),
359 : command("mmath", "log2", MATHunary_LOG2dbl, false, "The log2(x) function returns the base-2 logarithm of x.", args(1,2, arg("",dbl),arg("x",dbl))),
360 : command("mmath", "pow", MATHbinary_POWflt, false, "", args(1,3, arg("",flt),arg("x",flt),arg("y",flt))),
361 : command("mmath", "pow", MATHbinary_POWdbl, false, "The pow(x,y) function returns the value of x raised to the power of y.", args(1,3, arg("",dbl),arg("x",dbl),arg("y",dbl))),
362 : command("mmath", "nextafter", MATHbinary_NEXTAFTERflt, false, "The nextafter(x,y) function returns the next representable floating-point value of x in the direction of y.", args(1,3, arg("",flt),arg("x",flt),arg("y",flt))),
363 : command("mmath", "nextafter", MATHbinary_NEXTAFTERdbl, false, "The nextafter(x,y) function returns the next representable floating-point value of x in the direction of y.", args(1,3, arg("",dbl),arg("x",dbl),arg("y",dbl))),
364 : command("mmath", "sqrt", MATHunary_SQRTflt, false, "", args(1,2, arg("",flt),arg("y",flt))),
365 : command("mmath", "sqrt", MATHunary_SQRTdbl, false, "The sqrt(x) function returns the non-negative root of x.", args(1,2, arg("",dbl),arg("y",dbl))),
366 : command("mmath", "cbrt", MATHunary_CBRTflt, false, "", args(1,2, arg("",flt),arg("y",flt))),
367 : command("mmath", "cbrt", MATHunary_CBRTdbl, false, "The cbrt(x) function returns the cube root of x.", args(1,2, arg("",dbl),arg("y",dbl))),
368 : command("mmath", "ceil", MATHunary_CEILflt, false, "", args(1,2, arg("",flt),arg("y",flt))),
369 : command("mmath", "ceil", MATHunary_CEILdbl, false, "The ceil(x) function rounds x upwards to the nearest integer.", args(1,2, arg("",dbl),arg("y",dbl))),
370 : command("mmath", "fabs", MATHunary_FABSdbl, false, "The fabs(x) function returns the absolute value of the floating-point number x.", args(1,2, arg("",dbl),arg("y",dbl))),
371 : command("mmath", "floor", MATHunary_FLOORflt, false, "", args(1,2, arg("",flt),arg("y",flt))),
372 : command("mmath", "floor", MATHunary_FLOORdbl, false, "The floor(x) function rounds x downwards to the nearest integer.", args(1,2, arg("",dbl),arg("y",dbl))),
373 : command("mmath", "round", MATHbinary_ROUNDflt, false, "", args(1,3, arg("",flt),arg("x",flt),arg("y",int))),
374 : command("mmath", "round", MATHbinary_ROUNDdbl, false, "The round(n, m) returns n rounded to m places to the right \nof the decimal point; if m is omitted, to 0 places. m can be \nnegative to round off digits left of the decimal point. \nm must be an integer.", args(1,3, arg("",dbl),arg("x",dbl),arg("y",int))),
375 : command("mmath", "isnan", MATHunary_ISNAN, false, "The isnan(x) function returns true if x is 'not-a-number' \n(NaN), and false otherwise.", args(1,2, arg("",bit),arg("d",dbl))),
376 : command("mmath", "isinf", MATHunary_ISINF, false, "The isinf(x) function returns -1 if x represents negative \ninfinity, 1 if x represents positive infinity, and 0 otherwise.", args(1,2, arg("",int),arg("d",dbl))),
377 : command("mmath", "finite", MATHunary_FINITE, false, "The finite(x) function returns true if x is neither infinite \nnor a 'not-a-number' (NaN) value, and false otherwise.", args(1,2, arg("",bit),arg("d",dbl))),
378 : command("mmath", "rand", MATHrandint, true, "return a random number", args(1,1, arg("",int))),
379 : command("mmath", "rand", MATHrandintarg, true, "return a random number", args(1,2, arg("",int),arg("v",int))),
380 : command("mmath", "srand", MATHsrandint, false, "initialize the rand() function with a seed", args(1,2, arg("",void),arg("seed",int))),
381 : command("mmath", "sqlrand", MATHsqlrandint, false, "initialize the rand() function with a seed and call rand()", args(1,2, arg("",int),arg("seed",int))),
382 : command("mmath", "pi", MATHpi, false, "return an important mathematical value", args(1,1, arg("",dbl))),
383 : { .imp=NULL }
384 : };
385 : #include "mal_import.h"
386 : #ifdef _MSC_VER
387 : #undef read
388 : #pragma section(".CRT$XCU",read)
389 : #endif
390 323 : LIB_STARTUP_FUNC(init_mmath_mal)
391 323 : { mal_module2("mmath", NULL, mmath_init_funcs, MATHprelude, NULL); }
|
