source: trunk/RuntimeLibrary/src/jmi/jmi_math.c @ 10859

Last change on this file since 10859 was 10859, checked in by amartensen, 20 months ago

#5576 merge to trunk

File size: 13.2 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 cathegory_name[], const char func_name[], const char msg[], const char val[]) {
25    if (func_name != NULL) {
26        char buf[64];
27        sprintf(buf, "%s%s", cathegory_name, "InFunc");
28        jmi_log_node(jmi->log, logWarning, buf, "<func: %s, exp: %s, val:%s>", func_name, msg, val);
29    } else {
30        jmi_log_node(jmi->log, logWarning, cathegory_name, "<exp:%s, val: %s>", msg, val);
31    }
32}
33
34int jmi_check_nan(jmi_t *jmi, jmi_real_t* val, size_t n_val, jmi_int_t* index_of_nan) {
35    size_t i = 0;
36    for (i = 0; i < n_val; i++) {
37        if ( val[i] - val[i] != 0) {
38            *index_of_nan = i;
39            return JMI_ERROR;
40        }
41    }
42    return JMI_OK;
43}
44
45void jmi_inf_log(jmi_t *jmi, const char func_name[], const char msg[], jmi_real_t res, jmi_real_t x) {
46    if (((res - res) != 0)) {
47        if (jmi == NULL) jmi = jmi_get_current();
48       
49        if (res > 0) {
50            /* res is +inf */
51            char val[64];
52            sprintf(val, "%.14E", x);
53            jmi_log_func_or_eq(jmi, "RangeError", func_name, msg, val);
54        } else if (res < 0){
55            /* res is -inf */
56            char val[64];
57            sprintf(val, "%.14E", x);
58            jmi_log_func_or_eq(jmi, "RangeError", func_name, msg, val);
59        }
60    }
61}
62
63/*Some of these functions return types are a temporary remnant of CppAD*/
64jmi_real_t jmi_divide(jmi_t *jmi, const char func_name[], jmi_real_t num, jmi_real_t den, const char msg[]) {
65    if (den == 0) {
66        char val[64];
67        sprintf(val, "%.14E, %.14E", num, den);
68       
69        if (jmi == NULL) jmi = jmi_get_current();
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        if (jmi == NULL) jmi = jmi_get_current();
155        jmi_log_func_or_eq(jmi, "IllegalAtan2Input", func_name, msg, val);
156    }
157   
158    return to_return;
159}
160
161jmi_real_t jmi_atan2_function(const char func_name[], jmi_real_t x, jmi_real_t y, const char msg[]) {
162    return jmi_atan2(NULL, func_name, x, y, msg);
163}
164
165jmi_real_t jmi_atan2_equation(jmi_t *jmi, jmi_real_t x, jmi_real_t y, const char msg[]) {
166    return jmi_atan2(jmi, NULL, x, y, msg);
167}
168
169jmi_real_t jmi_pow(jmi_t *jmi, const char func_name[], jmi_real_t x, jmi_real_t y, const char msg[]) {
170    jmi_real_t to_return = pow(x, y);
171
172    if ((to_return - to_return) != 0) {
173        /* The returned value is not a number */
174        if (jmi == NULL)
175            jmi = jmi_get_current();
176
177        /* Check that the inputs are in the domain of the function*/
178        if (x > 0 || (x == 0 && y > 0) || (x < 0 && (int) y == y)) {
179            /* Range problem, will return JMI_INF or -JMI_INF */
180            char val[64];
181            sprintf(val, "%.14E, %.14E", x, y);
182            jmi_log_func_or_eq(jmi, "RangeError", func_name, msg, val);
183        } else if (x == 0 && y < 0) {
184            /* Pole error */
185            char val[64];
186            sprintf(val, "%.14E, %.14E", x, y);
187            jmi_log_func_or_eq(jmi, "DivideByZero", func_name, msg, val);
188        }
189    }
190    /* jmi_inf_log(jmi, func_name, msg, to_return); */
191    return to_return;
192}
193
194jmi_real_t jmi_pow_function(const char func_name[], jmi_real_t x, jmi_real_t y, const char msg[]) {
195    return jmi_pow(NULL, func_name, x, y, msg);
196}
197
198jmi_real_t jmi_pow_equation(jmi_t *jmi, jmi_real_t x, jmi_real_t y, const char msg[]) {
199    return jmi_pow(jmi, NULL, x, y, msg);
200}
201
202jmi_real_t jmi_exp(jmi_t *jmi, const char func_name[], jmi_real_t x, const char msg[]) {
203
204    jmi_real_t to_return = exp(x);
205    jmi_inf_log(jmi, func_name, msg, to_return, x);
206    return to_return;
207}
208
209jmi_real_t jmi_exp_function(const char func_name[], jmi_real_t x, const char msg[]) {
210    return jmi_exp(NULL, func_name, x, msg);
211}
212
213jmi_real_t jmi_exp_equation(jmi_t *jmi, jmi_real_t x, const char msg[]) {
214    return jmi_exp(jmi, NULL, x, msg);
215}
216
217jmi_real_t jmi_log(jmi_t *jmi, const char func_name[], jmi_real_t x, const char msg[]) {
218
219    jmi_real_t to_return = log(x);
220   
221    if ((to_return - to_return) != 0) {
222        /* The returned value is not a number */
223        if (jmi == NULL) jmi = jmi_get_current();
224       
225        if (x == 0) {
226            /* Pole problem, will return -JMI_INF */
227            char val[64];
228            sprintf(val, "%.14E", x);
229            jmi_log_func_or_eq(jmi, "LogarithmOfZero", func_name, msg, val);
230        } else if (x > 0) {
231            /* Range problem, will return JMI_INF */
232            char val[64];
233            sprintf(val, "%.14E", x);
234            jmi_log_func_or_eq(jmi, "LogarithmOfInf", func_name, msg, val);
235        }
236    }
237    return to_return;
238}
239
240jmi_real_t jmi_log_function(const char func_name[], jmi_real_t x, const char msg[]) {
241    return jmi_log(NULL, func_name, x, msg);
242}
243
244jmi_real_t jmi_log_equation(jmi_t *jmi, jmi_real_t x, const char msg[]) {
245    return jmi_log(jmi, NULL, x, msg);
246}
247
248jmi_real_t jmi_log10(jmi_t *jmi, const char func_name[], jmi_real_t x, const char msg[]) {
249
250    jmi_real_t to_return = log10(x);
251   
252    if ((to_return - to_return) != 0) {
253        /* The returned value is not a number */
254        if (jmi == NULL) jmi = jmi_get_current();
255       
256        if (x == 0) {
257            /* Pole problem, will return -JMI_INF */
258            char val[64];
259            sprintf(val, "%.14E", x);
260            jmi_log_func_or_eq(jmi, "LogarithmOfZero", func_name, msg, val);
261        } else if (x > 0) {
262            /* Infinity problem, will return JMI_INF */
263            char val[64];
264            sprintf(val, "%.14E", x);
265            jmi_log_func_or_eq(jmi, "LogarithmOfInf", func_name, msg, val);
266        }
267    }
268    return to_return;
269}
270
271jmi_real_t jmi_log10_function(const char func_name[], jmi_real_t x, const char msg[]) {
272    return jmi_log10(NULL, func_name, x, msg);
273}
274
275jmi_real_t jmi_log10_equation(jmi_t *jmi, jmi_real_t x, const char msg[]) {
276    return jmi_log10(jmi, NULL, x, msg);
277}
278
279jmi_real_t jmi_sinh(jmi_t *jmi, const char func_name[], jmi_real_t x, const char msg[]) {
280
281    jmi_real_t to_return = sinh(x);
282    jmi_inf_log(jmi, func_name, msg, to_return, x);
283    return to_return;
284}
285
286jmi_real_t jmi_sinh_function(const char func_name[], jmi_real_t x, const char msg[]) {
287    return jmi_sinh(NULL, func_name, x, msg);
288}
289
290jmi_real_t jmi_sinh_equation(jmi_t *jmi, jmi_real_t x, const char msg[]) {
291    return jmi_sinh(jmi, NULL, x, msg);
292}
293
294jmi_real_t jmi_cosh(jmi_t *jmi, const char func_name[], jmi_real_t x, const char msg[]) {
295
296    jmi_real_t to_return = cosh(x);
297    jmi_inf_log(jmi, func_name, msg, to_return, x);
298    return to_return;
299}
300
301jmi_real_t jmi_cosh_function(const char func_name[], jmi_real_t x, const char msg[]) {
302    return jmi_cosh(NULL, func_name, x, msg);
303}
304
305jmi_real_t jmi_cosh_equation(jmi_t *jmi, jmi_real_t x, const char msg[]) {
306    return jmi_cosh(jmi, NULL, x, msg);
307}
308
309jmi_real_t jmi_tan(jmi_t *jmi, const char func_name[], jmi_real_t x, const char msg[]) {
310
311    jmi_real_t to_return = tan(x);
312    jmi_inf_log(jmi, func_name, msg, to_return, x);
313    return to_return;
314}
315
316jmi_real_t jmi_tan_function(const char func_name[], jmi_real_t x, const char msg[]) {
317    return jmi_tan(NULL, func_name, x, msg);
318}
319
320jmi_real_t jmi_tan_equation(jmi_t *jmi, jmi_real_t x, const char msg[]) {
321    return jmi_tan(jmi, NULL, x, msg);
322}
323
324jmi_real_t jmi_sin(jmi_t *jmi, const char func_name[], jmi_real_t x, const char msg[]) {
325
326    jmi_real_t to_return = sin(x);
327    jmi_inf_log(jmi, func_name, msg, to_return, x);
328    return to_return;
329}
330
331jmi_real_t jmi_sin_function(const char func_name[], jmi_real_t x, const char msg[]) {
332    return jmi_sin(NULL, func_name, x, msg);
333}
334
335jmi_real_t jmi_sin_equation(jmi_t *jmi, jmi_real_t x, const char msg[]) {
336    return jmi_sin(jmi, NULL, x, msg);
337}
338
339jmi_real_t jmi_cos(jmi_t *jmi, const char func_name[], jmi_real_t x, const char msg[]) {
340
341    jmi_real_t to_return = cos(x);
342    jmi_inf_log(jmi, func_name, msg, to_return, x);
343    return to_return;
344}
345
346jmi_real_t jmi_cos_function(const char func_name[], jmi_real_t x, const char msg[]) {
347    return jmi_cos(NULL, func_name, x, msg);
348}
349
350jmi_real_t jmi_cos_equation(jmi_t *jmi, jmi_real_t x, const char msg[]) {
351    return jmi_cos(jmi, NULL, x, msg);
352}
353
354jmi_real_t jmi_atan(jmi_t *jmi, const char func_name[], jmi_real_t x, const char msg[]) {
355
356    jmi_real_t to_return = atan(x);
357    jmi_inf_log(jmi, func_name, msg, to_return, x);
358    return to_return;
359}
360
361jmi_real_t jmi_atan_function(const char func_name[], jmi_real_t x, const char msg[]) {
362    return jmi_atan(NULL, func_name, x, msg);
363}
364
365jmi_real_t jmi_atan_equation(jmi_t *jmi, jmi_real_t x, const char msg[]) {
366    return jmi_atan(jmi, NULL, x, msg);
367}
368
369jmi_real_t jmi_tanh(jmi_t *jmi, const char func_name[], jmi_real_t x, const char msg[]) {
370
371    jmi_real_t to_return = tanh(x);
372    jmi_inf_log(jmi, func_name, msg, to_return, x);
373    return to_return;
374}
375
376jmi_real_t jmi_tanh_function(const char func_name[], jmi_real_t x, const char msg[]) {
377    return jmi_tanh(NULL, func_name, x, msg);
378}
379
380jmi_real_t jmi_tanh_equation(jmi_t *jmi, jmi_real_t x, const char msg[]) {
381    return jmi_tanh(jmi, NULL, x, msg);
382}
383
384jmi_real_t jmi_abs(jmi_real_t v) {
385    return COND_EXP_GE(v, 0.0, v, -v);
386}
387
388jmi_real_t jmi_sign(jmi_real_t v) {
389    return COND_EXP_GT(v, 0.0, 1.0, COND_EXP_LT(v, 0.0, -1.0, 0.0));
390}
391
392jmi_real_t jmi_min(jmi_real_t x, jmi_real_t y) {
393    return COND_EXP_LT(x, y, x ,y);
394}
395
396jmi_real_t jmi_max(jmi_real_t x, jmi_real_t y) {
397    return COND_EXP_GT(x, y, x ,y);
398}
399
400jmi_real_t jmi_dround(jmi_real_t x) {
401    return (x >= 0)? floor(x + 0.5) : floor(x - 0.5);
402}
403
404jmi_real_t jmi_dremainder(jmi_t* jmi, jmi_real_t x, jmi_real_t y) {
405    jmi_real_t res = fmod(x,y);
406    jmi_real_t scaling = jmi_max(1.0, jmi_max(x,y));
407    return ((jmi_abs(res-y)/scaling)<jmi->time_events_epsilon)? (res-y)/scaling : res/scaling;
408}
409
410jmi_real_t jmi_sample(jmi_t* jmi, jmi_real_t offset, jmi_real_t h) {
411    jmi_real_t t = jmi_get_t(jmi)[0];
412    jmi_real_t remainder;
413    if (!jmi->atEvent || SURELY_LT_ZERO(t-offset) || jmi->atInitial) {
414        return JMI_FALSE;
415    }
416    remainder = jmi_dremainder(jmi, (t-offset),h);
417    if (jmi_abs(remainder) < jmi->time_events_epsilon)
418        return TRUE;
419    else
420        return FALSE;
421}
Note: See TracBrowser for help on using the repository browser.