source: branches/dev-cw-2658/RuntimeLibrary/src/jmi/jmi_math.c @ 14003

Last change on this file since 14003 was 14003, checked in by Christian Andersson, 4 weeks ago

Fixed segfault and minor refactoring of math.h. Related to ticket:5868

File size: 12.9 KB
Line 
1/*
2    Copyright (C) 2016 Modelon AB
3
4    This program is free software: you can redistribute it and/or modify
5    it under the terms of the GNU General Public License version 3 as published
6    by the Free Software Foundation, or optionally, under the terms of the
7    Common Public License version 1.0 as published by IBM.
8
9    This program is distributed in the hope that it will be useful,
10    but WITHOUT ANY WARRANTY; without even the implied warranty of
11    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12    GNU General Public License, or the Common Public License, for more details.
13
14    You should have received copies of the GNU General Public License
15    and the Common Public License along with this program.  If not,
16    see <http://www.gnu.org/licenses/> or
17    <http://www.ibm.com/developerworks/library/os-cpl.html/> respectively.
18*/
19
20#include "jmi.h"
21#include "jmi_math.h"
22
23
24void jmi_log_func_or_eq(jmi_t *jmi, const char category_name[], const char func_name[], const char msg[], const char val[]) {
25    if (jmi == NULL) jmi = jmi_get_current();
26   
27    if (func_name != NULL) {
28        char buf[64];
29        sprintf(buf, "%s%s", category_name, "InFunc");
30        jmi_log_node(jmi->log, logWarning, buf, "<func: %s, exp: %s, val:%s>", func_name, msg, val);
31    } else {
32        jmi_log_node(jmi->log, logWarning, category_name, "<exp:%s, val: %s>", msg, val);
33    }
34}
35
36int jmi_check_nan(jmi_t *jmi, jmi_real_t* val, size_t n_val, jmi_int_t* index_of_nan) {
37    size_t i = 0;
38    for (i = 0; i < n_val; i++) {
39        if ( val[i] - val[i] != 0) {
40            *index_of_nan = i;
41            return JMI_ERROR;
42        }
43    }
44    return JMI_OK;
45}
46
47void jmi_inf_log(jmi_t *jmi, const char func_name[], const char msg[], jmi_real_t res, jmi_real_t x) {
48    if (((res - res) != 0)) {
49       
50        if (res > 0) {
51            /* res is +inf */
52            char val[64];
53            sprintf(val, "%.14E", x);
54            jmi_log_func_or_eq(jmi, "RangeError", func_name, msg, val);
55        } else if (res < 0){
56            /* res is -inf */
57            char val[64];
58            sprintf(val, "%.14E", x);
59            jmi_log_func_or_eq(jmi, "RangeError", func_name, msg, val);
60        }
61    }
62}
63
64/*Some of these functions return types are a temporary remnant of CppAD*/
65jmi_real_t jmi_divide(jmi_t *jmi, const char func_name[], jmi_real_t num, jmi_real_t den, const char msg[]) {
66    if (den == 0) {
67        char val[64];
68        sprintf(val, "%.14E, %.14E", num, den);
69       
70        jmi_log_func_or_eq(jmi, "DivideByZero", func_name, msg, val);
71    }
72   
73    return num/den;
74}
75
76jmi_real_t jmi_divide_function(const char func_name[], jmi_real_t num, jmi_real_t den, const char msg[]) {
77    return jmi_divide(NULL, func_name, num, den, msg);
78}
79
80jmi_real_t jmi_divide_equation(jmi_t *jmi, jmi_real_t num, jmi_real_t den, const char msg[]) {
81    return jmi_divide(jmi, NULL, num, den, msg);
82}
83
84jmi_real_t jmi_sqrt(jmi_t *jmi, const char func_name[], jmi_real_t x, const char msg[]) {
85
86    jmi_real_t to_return = sqrt(x);
87
88    if (x < 0.0) {
89        /* Range problem, will return NAN */
90        char val[64];
91        sprintf(val, "%.14E", x);
92        jmi_log_func_or_eq(jmi, "RangeError", func_name, msg, val);
93    }
94
95    return to_return;
96}
97
98jmi_real_t jmi_sqrt_function(const char func_name[], jmi_real_t x, const char msg[]) {
99    return jmi_sqrt(NULL, func_name, x, msg);
100}
101
102jmi_real_t jmi_sqrt_equation(jmi_t *jmi, jmi_real_t x, const char msg[]) {
103    return jmi_sqrt(jmi, NULL, x, msg);
104}
105
106jmi_real_t jmi_asin(jmi_t *jmi, const char func_name[], jmi_real_t x, const char msg[]) {
107
108    jmi_real_t to_return = asin(x);
109
110    if ((to_return - to_return) != 0) {
111        /* The returned value is not a number */
112        char val[64];
113        sprintf(val, "%.14E", x);
114        jmi_log_func_or_eq(jmi, "RangeError", func_name, msg, val);
115    }
116    return to_return;
117}
118
119jmi_real_t jmi_asin_function(const char func_name[], jmi_real_t x, const char msg[]) {
120    return jmi_asin(NULL, func_name, x, msg);
121}
122
123jmi_real_t jmi_asin_equation(jmi_t *jmi, jmi_real_t x, const char msg[]) {
124    return jmi_asin(jmi, NULL, x, msg);
125}
126
127jmi_real_t jmi_acos(jmi_t *jmi, const char func_name[], jmi_real_t x, const char msg[]) {
128
129    jmi_real_t to_return = acos(x);
130
131    if ((to_return - to_return) != 0) {
132        /* The returned value is not a number */
133        char val[64];
134        sprintf(val, "%.14E", x);
135        jmi_log_func_or_eq(jmi, "RangeError", func_name, msg, val);
136    }
137    return to_return;
138}
139
140jmi_real_t jmi_acos_function(const char func_name[], jmi_real_t x, const char msg[]) {
141    return jmi_acos(NULL, func_name, x, msg);
142}
143
144jmi_real_t jmi_acos_equation(jmi_t *jmi, jmi_real_t x, const char msg[]) {
145    return jmi_acos(jmi, NULL, x, msg);
146}
147
148jmi_real_t jmi_atan2(jmi_t *jmi, const char func_name[], jmi_real_t x, jmi_real_t y, const char msg[]) {
149    jmi_real_t to_return = atan2(x, y);
150    if (x == 0 && y == 0) {
151        char val[64];
152        sprintf(val, "%.14E, %.14E", x, y);
153       
154        jmi_log_func_or_eq(jmi, "IllegalAtan2Input", func_name, msg, val);
155    }
156   
157    return to_return;
158}
159
160jmi_real_t jmi_atan2_function(const char func_name[], jmi_real_t x, jmi_real_t y, const char msg[]) {
161    return jmi_atan2(NULL, func_name, x, y, msg);
162}
163
164jmi_real_t jmi_atan2_equation(jmi_t *jmi, jmi_real_t x, jmi_real_t y, const char msg[]) {
165    return jmi_atan2(jmi, NULL, x, y, msg);
166}
167
168jmi_real_t jmi_pow(jmi_t *jmi, const char func_name[], jmi_real_t x, jmi_real_t y, const char msg[]) {
169    jmi_real_t to_return = pow(x, y);
170
171    /* The returned value is not a number */
172    if ((to_return - to_return) != 0) {
173
174        /* Check that the inputs are in the domain of the function*/
175        if (x > 0 || (x == 0 && y > 0) || (x < 0 && (int) y == y)) {
176            /* Range problem, will return JMI_INF or -JMI_INF */
177            char val[64];
178            sprintf(val, "%.14E, %.14E", x, y);
179            jmi_log_func_or_eq(jmi, "RangeError", func_name, msg, val);
180        } else if (x == 0 && y < 0) {
181            /* Pole error */
182            char val[64];
183            sprintf(val, "%.14E, %.14E", x, y);
184            jmi_log_func_or_eq(jmi, "DivideByZero", func_name, msg, val);
185        }
186    }
187    /* jmi_inf_log(jmi, func_name, msg, to_return); */
188    return to_return;
189}
190
191jmi_real_t jmi_pow_function(const char func_name[], jmi_real_t x, jmi_real_t y, const char msg[]) {
192    return jmi_pow(NULL, func_name, x, y, msg);
193}
194
195jmi_real_t jmi_pow_equation(jmi_t *jmi, jmi_real_t x, jmi_real_t y, const char msg[]) {
196    return jmi_pow(jmi, NULL, x, y, msg);
197}
198
199jmi_real_t jmi_exp(jmi_t *jmi, const char func_name[], jmi_real_t x, const char msg[]) {
200
201    jmi_real_t to_return = exp(x);
202    jmi_inf_log(jmi, func_name, msg, to_return, x);
203    return to_return;
204}
205
206jmi_real_t jmi_exp_function(const char func_name[], jmi_real_t x, const char msg[]) {
207    return jmi_exp(NULL, func_name, x, msg);
208}
209
210jmi_real_t jmi_exp_equation(jmi_t *jmi, jmi_real_t x, const char msg[]) {
211    return jmi_exp(jmi, NULL, x, msg);
212}
213
214jmi_real_t jmi_log(jmi_t *jmi, const char func_name[], jmi_real_t x, const char msg[]) {
215
216    jmi_real_t to_return = log(x);
217   
218    /* The returned value is not a number */
219    if ((to_return - to_return) != 0) {
220       
221        if (x == 0) {
222            /* Pole problem, will return -JMI_INF */
223            char val[64];
224            sprintf(val, "%.14E", x);
225            jmi_log_func_or_eq(jmi, "LogarithmOfZero", func_name, msg, val);
226        } else if (x > 0) {
227            /* Range problem, will return JMI_INF */
228            char val[64];
229            sprintf(val, "%.14E", x);
230            jmi_log_func_or_eq(jmi, "LogarithmOfInf", func_name, msg, val);
231        }
232    }
233    return to_return;
234}
235
236jmi_real_t jmi_log_function(const char func_name[], jmi_real_t x, const char msg[]) {
237    return jmi_log(NULL, func_name, x, msg);
238}
239
240jmi_real_t jmi_log_equation(jmi_t *jmi, jmi_real_t x, const char msg[]) {
241    return jmi_log(jmi, NULL, x, msg);
242}
243
244jmi_real_t jmi_log10(jmi_t *jmi, const char func_name[], jmi_real_t x, const char msg[]) {
245
246    jmi_real_t to_return = log10(x);
247   
248    /* The returned value is not a number */
249    if ((to_return - to_return) != 0) {
250
251        if (x == 0) {
252            /* Pole problem, will return -JMI_INF */
253            char val[64];
254            sprintf(val, "%.14E", x);
255            jmi_log_func_or_eq(jmi, "LogarithmOfZero", func_name, msg, val);
256        } else if (x > 0) {
257            /* Infinity problem, will return JMI_INF */
258            char val[64];
259            sprintf(val, "%.14E", x);
260            jmi_log_func_or_eq(jmi, "LogarithmOfInf", func_name, msg, val);
261        }
262    }
263    return to_return;
264}
265
266jmi_real_t jmi_log10_function(const char func_name[], jmi_real_t x, const char msg[]) {
267    return jmi_log10(NULL, func_name, x, msg);
268}
269
270jmi_real_t jmi_log10_equation(jmi_t *jmi, jmi_real_t x, const char msg[]) {
271    return jmi_log10(jmi, NULL, x, msg);
272}
273
274jmi_real_t jmi_sinh(jmi_t *jmi, const char func_name[], jmi_real_t x, const char msg[]) {
275
276    jmi_real_t to_return = sinh(x);
277    jmi_inf_log(jmi, func_name, msg, to_return, x);
278    return to_return;
279}
280
281jmi_real_t jmi_sinh_function(const char func_name[], jmi_real_t x, const char msg[]) {
282    return jmi_sinh(NULL, func_name, x, msg);
283}
284
285jmi_real_t jmi_sinh_equation(jmi_t *jmi, jmi_real_t x, const char msg[]) {
286    return jmi_sinh(jmi, NULL, x, msg);
287}
288
289jmi_real_t jmi_cosh(jmi_t *jmi, const char func_name[], jmi_real_t x, const char msg[]) {
290
291    jmi_real_t to_return = cosh(x);
292    jmi_inf_log(jmi, func_name, msg, to_return, x);
293    return to_return;
294}
295
296jmi_real_t jmi_cosh_function(const char func_name[], jmi_real_t x, const char msg[]) {
297    return jmi_cosh(NULL, func_name, x, msg);
298}
299
300jmi_real_t jmi_cosh_equation(jmi_t *jmi, jmi_real_t x, const char msg[]) {
301    return jmi_cosh(jmi, NULL, x, msg);
302}
303
304jmi_real_t jmi_tan(jmi_t *jmi, const char func_name[], jmi_real_t x, const char msg[]) {
305
306    jmi_real_t to_return = tan(x);
307    jmi_inf_log(jmi, func_name, msg, to_return, x);
308    return to_return;
309}
310
311jmi_real_t jmi_tan_function(const char func_name[], jmi_real_t x, const char msg[]) {
312    return jmi_tan(NULL, func_name, x, msg);
313}
314
315jmi_real_t jmi_tan_equation(jmi_t *jmi, jmi_real_t x, const char msg[]) {
316    return jmi_tan(jmi, NULL, x, msg);
317}
318
319jmi_real_t jmi_sin(jmi_t *jmi, const char func_name[], jmi_real_t x, const char msg[]) {
320
321    jmi_real_t to_return = sin(x);
322    jmi_inf_log(jmi, func_name, msg, to_return, x);
323    return to_return;
324}
325
326jmi_real_t jmi_sin_function(const char func_name[], jmi_real_t x, const char msg[]) {
327    return jmi_sin(NULL, func_name, x, msg);
328}
329
330jmi_real_t jmi_sin_equation(jmi_t *jmi, jmi_real_t x, const char msg[]) {
331    return jmi_sin(jmi, NULL, x, msg);
332}
333
334jmi_real_t jmi_cos(jmi_t *jmi, const char func_name[], jmi_real_t x, const char msg[]) {
335
336    jmi_real_t to_return = cos(x);
337    jmi_inf_log(jmi, func_name, msg, to_return, x);
338    return to_return;
339}
340
341jmi_real_t jmi_cos_function(const char func_name[], jmi_real_t x, const char msg[]) {
342    return jmi_cos(NULL, func_name, x, msg);
343}
344
345jmi_real_t jmi_cos_equation(jmi_t *jmi, jmi_real_t x, const char msg[]) {
346    return jmi_cos(jmi, NULL, x, msg);
347}
348
349jmi_real_t jmi_atan(jmi_t *jmi, const char func_name[], jmi_real_t x, const char msg[]) {
350
351    jmi_real_t to_return = atan(x);
352    jmi_inf_log(jmi, func_name, msg, to_return, x);
353    return to_return;
354}
355
356jmi_real_t jmi_atan_function(const char func_name[], jmi_real_t x, const char msg[]) {
357    return jmi_atan(NULL, func_name, x, msg);
358}
359
360jmi_real_t jmi_atan_equation(jmi_t *jmi, jmi_real_t x, const char msg[]) {
361    return jmi_atan(jmi, NULL, x, msg);
362}
363
364jmi_real_t jmi_tanh(jmi_t *jmi, const char func_name[], jmi_real_t x, const char msg[]) {
365
366    jmi_real_t to_return = tanh(x);
367    jmi_inf_log(jmi, func_name, msg, to_return, x);
368    return to_return;
369}
370
371jmi_real_t jmi_tanh_function(const char func_name[], jmi_real_t x, const char msg[]) {
372    return jmi_tanh(NULL, func_name, x, msg);
373}
374
375jmi_real_t jmi_tanh_equation(jmi_t *jmi, jmi_real_t x, const char msg[]) {
376    return jmi_tanh(jmi, NULL, x, msg);
377}
378
379jmi_real_t jmi_abs(jmi_real_t v) {
380    return COND_EXP_GE(v, 0.0, v, -v);
381}
382
383jmi_real_t jmi_sign(jmi_real_t v) {
384    return COND_EXP_GT(v, 0.0, 1.0, COND_EXP_LT(v, 0.0, -1.0, 0.0));
385}
386
387jmi_real_t jmi_min(jmi_real_t x, jmi_real_t y) {
388    return COND_EXP_LT(x, y, x ,y);
389}
390
391jmi_real_t jmi_max(jmi_real_t x, jmi_real_t y) {
392    return COND_EXP_GT(x, y, x ,y);
393}
394
395jmi_real_t jmi_dround(jmi_real_t x) {
396    return (x >= 0)? floor(x + 0.5) : floor(x - 0.5);
397}
398
399jmi_real_t jmi_dremainder(jmi_t* jmi, jmi_real_t x, jmi_real_t y) {
400    jmi_real_t res = fmod(x,y);
401    jmi_real_t scaling = jmi_max(1.0, jmi_max(x,y));
402    return ((jmi_abs(res-y)/scaling)<jmi->time_events_epsilon)? (res-y)/scaling : res/scaling;
403}
404
405jmi_real_t jmi_sample(jmi_t* jmi, jmi_real_t offset, jmi_real_t h) {
406    jmi_real_t t = jmi_get_t(jmi)[0];
407    jmi_real_t remainder;
408    if (!jmi->atEvent || SURELY_LT_ZERO(t-offset) || jmi->atInitial) {
409        return JMI_FALSE;
410    }
411    remainder = jmi_dremainder(jmi, (t-offset),h);
412    if (jmi_abs(remainder) < jmi->time_events_epsilon)
413        return TRUE;
414    else
415        return FALSE;
416}
Note: See TracBrowser for help on using the repository browser.