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                     break;
 741
 742                 default:
 743                     set_block_error(-2);
 744                     return;
 745             }
 746         }
 747 #if _MSC_VER
 748         if (!_finite(stack[bottom]) || _isnan(stack[bottom]))
 749         {
 750 #else
 751         if (isinf(stack[bottom]) || isnan(stack[bottom]))
 752         {
 753 #endif
 754             if (flag == 6)
 755             {
 756                 return;
 757             }
 758             set_block_error(-2);
 759             return;
 760         }
 761         else
 762         {
 763             block->outptr[0][0] = stack[bottom];
 764         }
 765     }
 766 }
 767 /*--------------------------------------------------------------------------*/