* Bug #13513 fixed - EXPRESSION block with "u1" as expression failed with a
[scilab.git] / scilab / modules / scicos_blocks / src / c / evaluate_expr.c
1 /*  Scicos
2 *
3 *  Copyright (C) INRIA - METALAU Project <scicos@inria.fr>
4 *
5 * This program is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation; either version 2 of the License, or
8 * (at your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13 * GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program; if not, write to the Free Software
17 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
18 *
19 * See the file ./license.txt
20 */
21 /*--------------------------------------------------------------------------*/
22 #include <math.h>
23 #include <stdlib.h>
24 #if _MSC_VER
25 #include <float.h>
26 #endif
27
28 #ifdef solaris
29 #include <ieeefp.h>
30 #endif
31
32 #include "core_math.h"
33 #include "scicos_block.h"
34 #include "machine.h" /* isinf */
35 #include "dynlib_scicos_blocks.h"
36 /*--------------------------------------------------------------------------*/
37 #if _MSC_VER
38 /*
39 arcsinh z = log (z+sqrt(1+z2))
40 */
41 double asinh(double x)
42 {
43     return log(x + sqrt(x * x + 1));
44 }
45 /*--------------------------------------------------------------------------*/
46 double acosh(double x)
47 {
48     return log(x + sqrt(x * x - 1));
49 }
50 /*--------------------------------------------------------------------------*/
51 /*
52 Inverse hyperbolic tangent (Atanh(x)) Log((1 + x) / (1 – x)) / 2
53 */
54 double atanh(double x)
55 {
56     return (double)(log ((1. + x) / (1. - x)) / 2);
57 }
58 #endif
59 /*--------------------------------------------------------------------------*/
60 SCICOS_BLOCKS_IMPEXP void evaluate_expr(scicos_block *block, int flag)
61 {
62     static double stack [1000];
63     static int count = 0, bottom = 0, nzcr = 0, i = 0, phase = 0;
64     int j = 0;
65     if (flag == 1 || flag == 6 || flag == 9)
66     {
67         phase = get_phase_simulation();
68         bottom = -1;
69         count = -1;
70         nzcr = -1;
71         while (count < block->nipar - 1)
72         {
73             count = count + 1;
74             switch (block->ipar[count])
75             {
76                 case 2:
77                     count = count + 1;
78                     bottom = bottom + 1;
79                     if (bottom > 999)
80                     {
81                         set_block_error(-16);
82                         return;
83                     }
84                     if (block->nin > 1)
85                     {
86                         stack[bottom] = block->inptr[block->ipar[count] - 1][0];
87                     }
88                     else
89                     {
90                         j = block->ipar[count] - 1;
91                         if (j < block->insz[0])
92                         {
93                             stack[bottom] = block->inptr[0][block->ipar[count] - 1];
94                         }
95                         else
96                         {
97                             stack[bottom] = 0.;
98                         }
99                     }
100                     break;
101                 case 6:
102                     count = count + 1;
103                     bottom = bottom + 1;
104                     if (bottom > 999)
105                     {
106                         set_block_error(-16);
107                         return;
108                     }
109                     stack[bottom] = block->rpar[block->ipar[count] - 1];
110                     break;
111                 case 5:
112                     count = count + 1;
113                     /* invalid script : call a function without lhs */
114                     if (bottom < 0)
115                     {
116                         set_block_error(-2);
117                         return;
118                     }
119
120                     switch (block->ipar[count])
121                     {
122                         case 1:
123                             stack[bottom - 1] = stack[bottom - 1] + stack[bottom];
124                             bottom = bottom - 1;
125                             break;
126                         case 2:
127                             stack[bottom - 1] = stack[bottom - 1] - stack[bottom];
128                             bottom = bottom - 1;
129                             break;
130                         case 3:
131                             stack[bottom - 1] = stack[bottom - 1] * stack[bottom];
132                             bottom = bottom - 1;
133                             break;
134                         case 7:
135                             stack[bottom - 1] = stack[bottom - 1] / stack[bottom];
136                             bottom = bottom - 1;
137                             break;
138                         case 15:
139                             stack[bottom - 1] = pow(stack[bottom - 1], stack[bottom]);
140                             bottom = bottom - 1;
141                             break;
142                         case 16: /* case == */
143                             if (block->ng > 0)
144                             {
145                                 nzcr = nzcr + 1;
146                             }
147                             if (flag == 9)
148                             {
149                                 block->g[nzcr] = stack[bottom - 1] - stack[bottom];
150                                 if (phase == 1)
151                                 {
152                                     block->mode[nzcr] = (stack[bottom - 1] == stack[bottom]);
153                                 }
154                             }
155                             if (phase == 1 || block->ng == 0)
156                             {
157                                 i = (stack[bottom - 1] == stack[bottom]);
158                             }
159                             else
160                             {
161                                 i = block->mode[nzcr];
162                             }
163                             stack[bottom - 1] = (double)i;
164                             bottom = bottom - 1;
165                             break;
166
167                         case 17:
168                             if (block->ng > 0)
169                             {
170                                 nzcr = nzcr + 1;
171                             }
172                             if (flag == 9)
173                             {
174                                 block->g[nzcr] = stack[bottom - 1] - stack[bottom];
175                                 if (phase == 1)
176                                 {
177                                     block->mode[nzcr] = (stack[bottom - 1] < stack[bottom]);
178                                 }
179                             }
180                             if (phase == 1 || block->ng == 0)
181                             {
182                                 i = (stack[bottom - 1] < stack[bottom]);
183                             }
184                             else
185                             {
186                                 i = block->mode[nzcr];
187                             }
188                             stack[bottom - 1] = (double)i;
189                             bottom = bottom - 1;
190                             break;
191                         case 18:
192                             if (block->ng > 0)
193                             {
194                                 nzcr = nzcr + 1;
195                             }
196                             if (flag == 9)
197                             {
198                                 block->g[nzcr] = stack[bottom - 1] - stack[bottom];
199                                 if (phase == 1)
200                                 {
201                                     block->mode[nzcr] = (stack[bottom - 1] > stack[bottom]);
202                                 }
203                             }
204                             if (phase == 1 || block->ng == 0)
205                             {
206                                 i = (stack[bottom - 1] > stack[bottom]);
207                             }
208                             else
209                             {
210                                 i = block->mode[nzcr];
211                             }
212                             stack[bottom - 1] = (double)i;
213                             bottom = bottom - 1;
214                             break;
215                         case 19:
216                             if (block->ng > 0)
217                             {
218                                 nzcr = nzcr + 1;
219                             }
220                             if (flag == 9)
221                             {
222                                 block->g[nzcr] = stack[bottom - 1] - stack[bottom];
223                                 if (phase == 1)
224                                 {
225                                     block->mode[nzcr] = (stack[bottom - 1] <= stack[bottom]);
226                                 }
227                             }
228                             if (phase == 1 || block->ng == 0)
229                             {
230                                 i = (stack[bottom - 1] <= stack[bottom]);
231                             }
232                             else
233                             {
234                                 i = block->mode[nzcr];
235                             }
236                             stack[bottom - 1] = (double)i;
237                             bottom = bottom - 1;
238                             break;
239                         case 20:
240                             if (block->ng > 0)
241                             {
242                                 nzcr = nzcr + 1;
243                             }
244                             if (flag == 9)
245                             {
246                                 block->g[nzcr] = stack[bottom - 1] - stack[bottom];
247                                 if (phase == 1)
248                                 {
249                                     block->mode[nzcr] = (stack[bottom - 1] >= stack[bottom]);
250                                 }
251                             }
252                             if (phase == 1 || block->ng == 0)
253                             {
254                                 i = (stack[bottom - 1] >= stack[bottom]);
255                             }
256                             else
257                             {
258                                 i = block->mode[nzcr];
259                             }
260                             stack[bottom - 1] = (double)i;
261                             bottom = bottom - 1;
262                             break;
263                         case 21:
264                             if (block->ng > 0)
265                             {
266                                 nzcr = nzcr + 1;
267                             }
268                             if (flag == 9)
269                             {
270                                 block->g[nzcr] = stack[bottom - 1] - stack[bottom];
271                                 if (phase == 1)
272                                 {
273                                     block->mode[nzcr] = (stack[bottom - 1] != stack[bottom]);
274                                 }
275                             }
276                             if (phase == 1 || block->ng == 0)
277                             {
278                                 i = (stack[bottom - 1] != stack[bottom]);
279                             }
280                             else
281                             {
282                                 i = block->mode[nzcr];
283                             }
284                             stack[bottom - 1] = (double)i;
285                             bottom = bottom - 1;
286                             break;
287                         case 28:
288                             if (block->ng > 0)
289                             {
290                                 nzcr = nzcr + 1;
291                             }
292                             if (flag == 9)
293                             {
294                                 block->g[nzcr] = stack[bottom - 1] - stack[bottom];
295                                 if (phase == 1)
296                                 {
297                                     block->mode[nzcr] = ((int)stack[bottom - 1] || (int)stack[bottom]);
298                                 }
299                             }
300                             if (phase == 1 || block->ng == 0)
301                             {
302                                 i = ((int)stack[bottom - 1] || (int)stack[bottom]);
303                             }
304                             else
305                             {
306                                 i = block->mode[nzcr];
307                             }
308                             stack[bottom - 1] = (double)i;
309                             bottom = bottom - 1;
310                             break;
311                         case 29:
312                             if (block->ng > 0)
313                             {
314                                 nzcr = nzcr + 1;
315                             }
316                             if (flag == 9)
317                             {
318                                 block->g[nzcr] = stack[bottom - 1] - stack[bottom];
319                                 if (phase == 1)
320                                 {
321                                     block->mode[nzcr] = ((int)stack[bottom - 1] && (int)stack[bottom]);
322                                 }
323                             }
324                             if (phase == 1 || block->ng == 0)
325                             {
326                                 i = ((int)stack[bottom - 1] && (int)stack[bottom]);
327                             }
328                             else
329                             {
330                                 i = block->mode[nzcr];
331                             }
332                             stack[bottom - 1] = (double)i;
333                             bottom = bottom - 1;
334                             break;
335
336                         case 30:
337                             if (block->ng > 0)
338                             {
339                                 nzcr = nzcr + 1;
340                             }
341                             if (flag == 9)
342                             {
343                                 block->g[nzcr] = stack[bottom];
344                                 if (phase == 1)
345                                 {
346                                     block->mode[nzcr] = (0.0 == stack[bottom]);
347                                 }
348                             }
349                             if (phase == 1 || block->ng == 0)
350                             {
351                                 i = (stack[bottom] == 0.0);
352                             }
353                             else
354                             {
355                                 i = block->mode[nzcr];
356                             }
357                             if (i)
358                             {
359                                 stack[bottom] = 1.0;
360                             }
361                             else
362                             {
363                                 stack[bottom] = 0.0;
364                             }
365                             break;
366                         case 99:
367                             stack[bottom] = -stack[bottom];
368                             break;
369                         case 101:
370                             stack[bottom] = sin(stack[bottom]);
371                             break;
372                         case 102:
373                             stack[bottom] = cos(stack[bottom]);
374                             break;
375                         case 103:
376                             stack[bottom] = tan(stack[bottom]);
377                             break;
378                         case 104:
379                             stack[bottom] = exp(stack[bottom]);
380                             break;
381                         case 105:
382                             stack[bottom] = log(stack[bottom]);
383                             break;
384                         case 106:
385                             stack[bottom] = sinh(stack[bottom]);
386                             break;
387                         case 107:
388                             stack[bottom] = cosh(stack[bottom]);
389                             break;
390                         case 108:
391                             stack[bottom] = tanh(stack[bottom]);
392                             break;
393
394                         case 109:
395                             if (block->ng > 0)
396                             {
397                                 nzcr = nzcr + 1;
398                             }
399                             if (flag == 9)
400                             {
401                                 if (stack[bottom] > 0)
402                                 {
403                                     i = (int)floor(stack[bottom]);
404                                 }
405                                 else
406                                 {
407                                     i = (int)ceil(stack[bottom]);
408                                 }
409                                 if (i == 0)
410                                 {
411                                     block->g[nzcr] = (stack[bottom] - 1) * (stack[bottom] + 1);
412                                 }
413                                 else if (i > 0)
414                                 {
415                                     block->g[nzcr] = (stack[bottom] - i - 1.) * (stack[bottom] - i);
416                                 }
417                                 else
418                                 {
419                                     block->g[nzcr] = (stack[bottom] - i) * (stack[bottom] - i + 1);
420                                 }
421                                 if (i % 2)
422                                 {
423                                     block->g[nzcr] = -block->g[nzcr];
424                                 }
425                                 if (phase == 1)
426                                 {
427                                     block->mode[nzcr] = i;
428                                 }
429                             }
430                             if (phase == 1 || block->ng == 0)
431                             {
432                                 if (stack[bottom] > 0)
433                                 {
434                                     stack[bottom] = floor(stack[bottom]);
435                                 }
436                                 else
437                                 {
438                                     stack[bottom] = ceil(stack[bottom]);
439                                 }
440                             }
441                             else
442                             {
443                                 stack[bottom] = (double) block->mode[nzcr];
444                             }
445                             break;
446                             /*
447                             if (stack[bottom]>0) {
448                               stack[bottom]=floor(stack[bottom]);
449                             }else{
450                               stack[bottom]=ceil(stack[bottom]);
451                               }*/
452                             break;
453                         case 110:
454                             if (block->ng > 0)
455                             {
456                                 nzcr = nzcr + 1;
457                             }
458                             if (flag == 9)
459                             {
460                                 if (stack[bottom] > 0)
461                                 {
462                                     i = (int)floor(stack[bottom] + .5);
463                                 }
464                                 else
465                                 {
466                                     i = (int)ceil(stack[bottom] - .5);
467                                 }
468                                 block->g[nzcr] = (stack[bottom] - i - .5) * (stack[bottom] - i + .5);
469                                 if (i % 2)
470                                 {
471                                     block->g[nzcr] = -block->g[nzcr];
472                                 }
473                                 if (phase == 1)
474                                 {
475                                     block->mode[nzcr] = i;
476                                 }
477                             }
478                             if (phase == 1 || block->ng == 0)
479                             {
480                                 if (stack[bottom] > 0)
481                                 {
482                                     stack[bottom] = floor(stack[bottom] + .5);
483                                 }
484                                 else
485                                 {
486                                     stack[bottom] = ceil(stack[bottom] - .5);
487                                 }
488                             }
489                             else
490                             {
491                                 stack[bottom] = (double) block->mode[nzcr];
492                             }
493                             break;
494                         /*  if (stack[bottom]>0) {
495                           stack[bottom]=floor(stack[bottom]+.5);
496                         }else{
497                           stack[bottom]=ceil(stack[bottom]-.5);
498                         }*/
499                         case 111:
500                             if (block->ng > 0)
501                             {
502                                 nzcr = nzcr + 1;
503                             }
504                             if (flag == 9)
505                             {
506                                 i = (int)ceil(stack[bottom]);
507                                 block->g[nzcr] = (stack[bottom] - i) * (stack[bottom] - i + 1);
508                                 if (i % 2)
509                                 {
510                                     block->g[nzcr] = -block->g[nzcr];
511                                 }
512                                 if (phase == 1)
513                                 {
514                                     block->mode[nzcr] = i;
515                                 }
516                             }
517                             if (phase == 1 || block->ng == 0)
518                             {
519                                 stack[bottom] = ceil(stack[bottom]);
520                             }
521                             else
522                             {
523                                 stack[bottom] = (double) block->mode[nzcr];
524                             }
525                             break;
526                         case 112:
527                             if (block->ng > 0)
528                             {
529                                 nzcr = nzcr + 1;
530                             }
531                             if (flag == 9)
532                             {
533                                 i = (int)floor(stack[bottom]);
534                                 block->g[nzcr] = (stack[bottom] - i - 1) * (stack[bottom] - i);
535                                 if (i % 2)
536                                 {
537                                     block->g[nzcr] = -block->g[nzcr];
538                                 }
539                                 if (phase == 1)
540                                 {
541                                     block->mode[nzcr] = i;
542                                 }
543                             }
544                             if (phase == 1 || block->ng == 0)
545                             {
546                                 stack[bottom] = floor(stack[bottom]);
547                             }
548                             else
549                             {
550                                 stack[bottom] = (double) block->mode[nzcr];
551                             }
552                             break;
553                         case 113:
554                             if (block->ng > 0)
555                             {
556                                 nzcr = nzcr + 1;
557                             }
558                             if (flag == 9)
559                             {
560                                 if (stack[bottom] > 0)
561                                 {
562                                     i = 1;
563                                 }
564                                 else if (stack[bottom] < 0)
565                                 {
566                                     i = -1;
567                                 }
568                                 else
569                                 {
570                                     i = 0;
571                                 }
572                                 block->g[nzcr] = stack[bottom];
573                                 if (phase == 1)
574                                 {
575                                     block->mode[nzcr] = i;
576                                 }
577                             }
578                             if (phase == 1 || block->ng == 0)
579                             {
580                                 if (stack[bottom] > 0)
581                                 {
582                                     stack[bottom] = 1.0;
583                                 }
584                                 else if (stack[bottom] < 0)
585                                 {
586                                     stack[bottom] = -1.0;
587                                 }
588                                 else
589                                 {
590                                     stack[bottom] = 0.0;
591                                 }
592                             }
593                             else
594                             {
595                                 stack[bottom] = (double) block->mode[nzcr];
596                             }
597                             break;
598                         /* if (stack[bottom]>0) {
599                           stack[bottom]=1.0;
600                         }else if(stack[bottom]<0){
601                           stack[bottom]=-1.0;
602                         }else{
603                           stack[bottom]=0.0;
604                           }*/
605                         case 114:  /* abs */
606                             if (block->ng > 0)
607                             {
608                                 nzcr = nzcr + 1;
609                             }
610                             if (flag == 9)
611                             {
612                                 if (stack[bottom] > 0)
613                                 {
614                                     i = 1;
615                                 }
616                                 else if (stack[bottom] < 0)
617                                 {
618                                     i = -1;
619                                 }
620                                 else
621                                 {
622                                     i = 0;
623                                 }
624                                 block->g[nzcr] = stack[bottom];
625                                 if (phase == 1)
626                                 {
627                                     block->mode[nzcr] = i;
628                                 }
629                             }
630                             if (phase == 1 || block->ng == 0)
631                             {
632                                 if (stack[bottom] > 0)
633                                 {
634                                     /* stack[bottom] = stack[bottom]; */
635                                 }
636                                 else
637                                 {
638                                     stack[bottom] = -stack[bottom];
639                                 }
640                             }
641                             else
642                             {
643                                 stack[bottom] = stack[bottom] * (block->mode[nzcr]);
644                             }
645                             break;
646                         /* if (stack[bottom]>0) {
647                           stack[bottom]=stack[bottom];
648                         }else {
649                           stack[bottom]=-stack[bottom];
650                           }*/
651                         case 115:
652                             if (block->ng > 0)
653                             {
654                                 nzcr = nzcr + 1;
655                             }
656                             if (flag == 9)
657                             {
658                                 if (stack[bottom] > stack[bottom - 1])
659                                 {
660                                     i = 0;
661                                 }
662                                 else
663                                 {
664                                     i = 1;
665                                 }
666                                 block->g[nzcr] = stack[bottom] - stack[bottom - 1];
667                                 if (phase == 1)
668                                 {
669                                     block->mode[nzcr] = i;
670                                 }
671                             }
672                             if (phase == 1 || block->ng == 0)
673                             {
674                                 stack[bottom - 1] = Max(stack[bottom - 1], stack[bottom]);
675                             }
676                             else
677                             {
678                                 stack[bottom - 1] = stack[bottom - block->mode[nzcr]];
679                             }
680                             bottom = bottom - 1;
681                             break;
682                         case 116:
683                             if (block->ng > 0)
684                             {
685                                 nzcr = nzcr + 1;
686                             }
687                             if (flag == 9)
688                             {
689                                 if (stack[bottom] < stack[bottom - 1])
690                                 {
691                                     i = 0;
692                                 }
693                                 else
694                                 {
695                                     i = 1;
696                                 }
697                                 block->g[nzcr] = stack[bottom] - stack[bottom - 1];
698                                 if (phase == 1)
699                                 {
700                                     block->mode[nzcr] = i;
701                                 }
702                             }
703                             if (phase == 1 || block->ng == 0)
704                             {
705                                 stack[bottom - 1] = Min(stack[bottom - 1], stack[bottom]);
706                             }
707                             else
708                             {
709                                 stack[bottom - 1] = stack[bottom - block->mode[nzcr]];
710                             }
711                             bottom = bottom - 1;
712                             break;
713                         case 117:
714                             stack[bottom] = asin(stack[bottom]);
715                             break;
716                         case 118:
717                             stack[bottom] = acos(stack[bottom]);
718                             break;
719                         case 119:
720                             stack[bottom] = atan(stack[bottom]);
721                             break;
722                         case 120:
723                             stack[bottom] = asinh(stack[bottom]);
724                             break;
725                         case 121:
726                             stack[bottom] = acosh(stack[bottom]);
727                             break;
728                         case 122:
729                             stack[bottom] = atanh(stack[bottom]);
730                             break;
731                         case 123:
732                             stack[bottom - 1] = atan2(stack[bottom - 1], stack[bottom]);
733                             bottom = bottom - 1;
734                             break;
735
736                         case 124:
737                             stack[bottom] = log10(stack[bottom]);
738                             break;
739                     }
740             }
741         }
742 #if _MSC_VER
743         if (!_finite(stack[bottom]) || _isnan(stack[bottom]))
744         {
745 #else
746         if (isinf(stack[bottom]) || isnan(stack[bottom]))
747         {
748 #endif
749             if (flag == 6)
750             {
751                 return;
752             }
753             set_block_error(-2);
754             return;
755         }
756         else
757         {
758             block->outptr[0][0] = stack[bottom];
759         }
760     }
761 }
762 /*--------------------------------------------------------------------------*/