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 45 : radians(double x)
44 : {
45 45 : 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 32 : 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 : static str
204 0 : MATHunary_FABSdbl(dbl *res, const dbl *a)
205 : {
206 0 : *res = is_dbl_nil(*a) ? dbl_nil : fabs(*a);
207 0 : return MAL_SUCCEED;
208 : }
209 :
210 0 : roundM5(dbl)
211 0 : roundM5(flt)
212 : static str
213 0 : MATHunary_ISNAN(bit *res, const dbl *a)
214 : {
215 0 : if (is_dbl_nil(*a)) {
216 0 : *res = bit_nil;
217 : } else {
218 0 : *res = isnan(*a) != 0;
219 : }
220 0 : return MAL_SUCCEED;
221 : }
222 :
223 : static str
224 0 : MATHunary_ISINF(int *res, const dbl *a)
225 : {
226 0 : if (is_dbl_nil(*a)) {
227 0 : *res = int_nil;
228 : } else {
229 0 : if (isinf(*a)) {
230 0 : *res = (*a < 0.0) ? -1 : 1;
231 : } else {
232 0 : *res = 0;
233 : }
234 : }
235 0 : return MAL_SUCCEED;
236 : }
237 :
238 : static str
239 0 : MATHunary_FINITE(bit *res, const dbl *a)
240 : {
241 0 : if (is_dbl_nil(*a)) {
242 0 : *res = bit_nil;
243 : } else {
244 0 : *res = isfinite(*a) != 0;
245 : }
246 0 : return MAL_SUCCEED;
247 : }
248 :
249 : /* global pseudo random generator state */
250 : random_state_engine mmath_rse;
251 : /* serialize access to state */
252 : MT_Lock mmath_rse_lock = MT_LOCK_INITIALIZER(mmath_rse_lock);
253 :
254 : static str
255 332 : MATHprelude(void)
256 : {
257 332 : init_random_state_engine(mmath_rse, (uint64_t) GDKusec());
258 332 : return MAL_SUCCEED;
259 : }
260 :
261 : static str
262 1010004 : MATHrandint(int *res)
263 : {
264 : #ifdef __COVERITY__
265 : *res = 0;
266 : #else
267 1010004 : MT_lock_set(&mmath_rse_lock);
268 1010004 : *res = (int) (next(mmath_rse) >> 33);
269 1010004 : MT_lock_unset(&mmath_rse_lock);
270 : #endif
271 1010004 : return MAL_SUCCEED;
272 : }
273 :
274 : static str
275 0 : MATHrandintarg(int *res, const int *dummy)
276 : {
277 0 : (void) dummy;
278 : #ifdef __COVERITY__
279 : *res = 0;
280 : #else
281 0 : MT_lock_set(&mmath_rse_lock);
282 0 : *res = (int) (next(mmath_rse) >> 33);
283 0 : MT_lock_unset(&mmath_rse_lock);
284 : #endif
285 0 : return MAL_SUCCEED;
286 : }
287 :
288 : static str
289 3 : MATHsrandint(void *ret, const int *seed)
290 : {
291 3 : (void) ret;
292 3 : MT_lock_set(&mmath_rse_lock);
293 3 : init_random_state_engine(mmath_rse, (uint64_t) * seed);
294 3 : MT_lock_unset(&mmath_rse_lock);
295 3 : return MAL_SUCCEED;
296 : }
297 :
298 : static str
299 2 : MATHsqlrandint(int *res, const int *seed)
300 : {
301 : #ifdef __COVERITY__
302 : (void) seed;
303 : *res = 0;
304 : #else
305 2 : MT_lock_set(&mmath_rse_lock);
306 2 : init_random_state_engine(mmath_rse, (uint64_t) * seed);
307 2 : *res = (int) (next(mmath_rse) >> 33);
308 2 : MT_lock_unset(&mmath_rse_lock);
309 : #endif
310 2 : return MAL_SUCCEED;
311 : }
312 :
313 : static str
314 84 : MATHpi(dbl *pi)
315 : {
316 84 : *pi = M_PI;
317 84 : return MAL_SUCCEED;
318 : }
319 :
320 : #include "mel.h"
321 : mel_func mmath_init_funcs[] = {
322 : command("mmath", "acos", MATHunary_ACOSflt, false, "", args(1,2, arg("",flt),arg("x",flt))),
323 : 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))),
324 : command("mmath", "asin", MATHunary_ASINflt, false, "", args(1,2, arg("",flt),arg("x",flt))),
325 : 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))),
326 : command("mmath", "atan", MATHunary_ATANflt, false, "", args(1,2, arg("",flt),arg("x",flt))),
327 : 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))),
328 : command("mmath", "atan2", MATHbinary_ATAN2flt, false, "", args(1,3, arg("",flt),arg("x",flt),arg("y",flt))),
329 : 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))),
330 : command("mmath", "cos", MATHunary_COSflt, false, "", args(1,2, arg("",flt),arg("x",flt))),
331 : 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))),
332 : command("mmath", "sin", MATHunary_SINflt, false, "", args(1,2, arg("",flt),arg("x",flt))),
333 : 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))),
334 : command("mmath", "tan", MATHunary_TANflt, false, "", args(1,2, arg("",flt),arg("x",flt))),
335 : 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))),
336 : command("mmath", "cot", MATHunary_COTflt, false, "", args(1,2, arg("",flt),arg("x",flt))),
337 : 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))),
338 : command("mmath", "cosh", MATHunary_COSHflt, false, "", args(1,2, arg("",flt),arg("x",flt))),
339 : 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))),
340 : command("mmath", "sinh", MATHunary_SINHflt, false, "", args(1,2, arg("",flt),arg("x",flt))),
341 : command("mmath", "sinh", MATHunary_SINHdbl, false, "The sinh() function returns the hyperbolic sine of x, which \nis defined mathematically as (exp(x) - exp(-x)) / 2.", args(1,2, arg("",dbl),arg("x",dbl))),
342 : command("mmath", "tanh", MATHunary_TANHflt, false, "", args(1,2, arg("",flt),arg("x",flt))),
343 : 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))),
344 : command("mmath", "radians", MATHunary_RADIANSflt, false, "", args(1,2, arg("",flt),arg("x",flt))),
345 : command("mmath", "radians", MATHunary_RADIANSdbl, false, "The radians() function converts degrees into radians", args(1,2, arg("",dbl),arg("x",dbl))),
346 : command("mmath", "degrees", MATHunary_DEGREESflt, false, "", args(1,2, arg("",flt),arg("x",flt))),
347 : command("mmath", "degrees", MATHunary_DEGREESdbl, false, "The degrees() function converts radians into degrees", args(1,2, arg("",dbl),arg("x",dbl))),
348 : command("mmath", "exp", MATHunary_EXPflt, false, "", args(1,2, arg("",flt),arg("x",flt))),
349 : 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))),
350 : command("mmath", "log", MATHunary_LOGflt, false, "", args(1,2, arg("",flt),arg("x",flt))),
351 : command("mmath", "log", MATHunary_LOGdbl, false, "The log(x) function returns the natural logarithm of x.", args(1,2, arg("",dbl),arg("x",dbl))),
352 : 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))),
353 : 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))),
354 : command("mmath", "log10", MATHunary_LOG10flt, false, "", args(1,2, arg("",flt),arg("x",flt))),
355 : 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))),
356 : command("mmath", "log2", MATHunary_LOG2flt, false, "", args(1,2, arg("",flt),arg("x",flt))),
357 : 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))),
358 : command("mmath", "pow", MATHbinary_POWflt, false, "", args(1,3, arg("",flt),arg("x",flt),arg("y",flt))),
359 : 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))),
360 : command("mmath", "sqrt", MATHunary_SQRTflt, false, "", args(1,2, arg("",flt),arg("y",flt))),
361 : 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))),
362 : command("mmath", "cbrt", MATHunary_CBRTflt, false, "", args(1,2, arg("",flt),arg("y",flt))),
363 : command("mmath", "cbrt", MATHunary_CBRTdbl, false, "The cbrt(x) function returns the cube root of x.", args(1,2, arg("",dbl),arg("y",dbl))),
364 : command("mmath", "ceil", MATHunary_CEILflt, false, "", args(1,2, arg("",flt),arg("y",flt))),
365 : 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))),
366 : command("mmath", "fabs", MATHunary_FABSdbl, false, "The fabs(x) function returns the absolute value of the \nfloating-point number x.", args(1,2, arg("",dbl),arg("y",dbl))),
367 : command("mmath", "floor", MATHunary_FLOORflt, false, "", args(1,2, arg("",flt),arg("y",flt))),
368 : 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))),
369 : command("mmath", "round", MATHbinary_ROUNDflt, false, "", args(1,3, arg("",flt),arg("x",flt),arg("y",int))),
370 : 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))),
371 : 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))),
372 : 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))),
373 : 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))),
374 : command("mmath", "rand", MATHrandint, true, "return a random number", args(1,1, arg("",int))),
375 : command("mmath", "rand", MATHrandintarg, true, "return a random number", args(1,2, arg("",int),arg("v",int))),
376 : command("mmath", "srand", MATHsrandint, false, "initialize the rand() function with a seed", args(1,2, arg("",void),arg("seed",int))),
377 : command("mmath", "sqlrand", MATHsqlrandint, false, "initialize the rand() function with a seed and call rand()", args(1,2, arg("",int),arg("seed",int))),
378 : command("mmath", "pi", MATHpi, false, "return an important mathematical value", args(1,1, arg("",dbl))),
379 : { .imp=NULL }
380 : };
381 : #include "mal_import.h"
382 : #ifdef _MSC_VER
383 : #undef read
384 : #pragma section(".CRT$XCU",read)
385 : #endif
386 325 : LIB_STARTUP_FUNC(init_mmath_mal)
387 325 : { mal_module2("mmath", NULL, mmath_init_funcs, MATHprelude, NULL); }
|
