src/expr-evl.c

   1 /* PSPP - computes sample statistics.
   2    Copyright (C) 1997-9, 2000 Free Software Foundation, Inc.
   3    Written by Ben Pfaff <blp@gnu.org>.
   4
   5    This program is free software; you can redistribute it and/or
   6    modify it under the terms of the GNU General Public License as
   7    published by the Free Software Foundation; either version 2 of the
   8    License, or (at your option) any later version.
   9
  10    This program is distributed in the hope that it will be useful, but
  11    WITHOUT ANY WARRANTY; without even the implied warranty of
  12    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13    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., 59 Temple Place - Suite 330, Boston, MA
  18    02111-1307, USA. */
  19
  20 #include <config.h>
  21
  22 #if TIME_WITH_SYS_TIME
  23 #include <sys/time.h>
  24 #include <time.h>
  25 #else
  26 #if HAVE_SYS_TIME_H
  27 #include <sys/time.h>
  28 #else
  29 #include <time.h>
  30 #endif
  31 #endif
  32
  33 #include <ctype.h>
  34 #include "expr.h"
  35 #include "exprP.h"
  36 #include <assert.h>
  37 #include <math.h>
  38 #include <errno.h>
  39 #include <stdio.h>
  40 #include "approx.h"
  41 #include "data-in.h"
  42 #include "error.h"
  43 #include "julcal/julcal.h"
  44 #include "magic.h"
  45 #include "pool.h"
  46 #include "random.h"
  47 #include "stats.h"
  48 #include "str.h"
  49 #include "var.h"
  50 #include "vfm.h"
  51 #include "vfmP.h"
  52
  53 double
  54 expr_evaluate (struct expression *e, struct ccase *c, union value *v)
  55 {
  56   unsigned char *op = e->op;
  57   double *dbl = e->num;
  58   unsigned char *str = e->str;
  59   struct variable **vars = e->var;
  60   int i, j;
  61
  62   /* Stack pointer. */
  63   union value *sp = e->stack;
  64
  65   pool_clear (e->pool);
  66
  67   for (;;)
  68     {
  69       switch (*op++)
  70         {
  71         case OP_PLUS:
  72           sp -= *op - 1;
  73           if (sp->f != SYSMIS)
  74             for (i = 1; i < *op; i++)
  75               {
  76                 if (sp[i].f == SYSMIS)
  77                   {
  78                     sp->f = SYSMIS;
  79                     break;
  80                   }
  81                 else
  82                   sp->f += sp[i].f;
  83               }
  84           op++;
  85           break;
  86         case OP_MUL:
  87           sp -= *op - 1;
  88           if (sp->f != SYSMIS)
  89             for (i = 1; i < *op; i++)
  90               {
  91                 if (sp[i].f == SYSMIS)
  92                   {
  93                     sp->f = SYSMIS;
  94                     break;
  95                   }
  96                 else
  97                   sp->f *= sp[i].f;
  98               }
  99           op++;
 100           break;
 101         case OP_POW:
 102           sp--;
 103           if (sp[0].f == SYSMIS)
 104             {
 105               if (approx_eq (sp[1].f, 0.0))
 106                 sp->f = 1.0;
 107             }
 108           else if (sp[1].f == SYSMIS)
 109             {
 110               if (sp[0].f == 0.0)
 111                 /* SYSMIS**0 */
 112                 sp->f = 0.0;
 113               else
 114                 sp->f = SYSMIS;
 115             }
 116           else if (approx_eq (sp[0].f, 0.0) && approx_eq (sp[1].f, 0.0))
 117             sp->f = SYSMIS;
 118           else
 119             sp->f = pow (sp[0].f, sp[1].f);
 120           break;
 121
 122         case OP_AND:
 123           /* Note that the equality operator (==) may be used here
 124              (instead of approx_eq) because booleans are always
 125              *exactly* 0, 1, or SYSMIS.
 126
 127              Truth table (in order of detection):
 128
 129              1:
 130              0 and 0 = 0
 131              0 and 1 = 0
 132              0 and SYSMIS = 0
 133
 134              2:
 135              1 and 0 = 0
 136              SYSMIS and 0 = 0
 137
 138              3:
 139              1 and SYSMIS = SYSMIS
 140              SYSMIS and SYSMIS = SYSMIS
 141
 142              4:
 143              1 and 1 = 1
 144              SYSMIS and 1 = SYSMIS
 145
 146            */
 147           sp--;
 148           if (sp[0].f == 0.0);  /* 1 */
 149           else if (sp[1].f == 0.0)
 150             sp->f = 0.0;        /* 2 */
 151           else if (sp[1].f == SYSMIS)
 152             sp->f = SYSMIS;     /* 3 */
 153           break;
 154         case OP_OR:
 155           /* Truth table (in order of detection):
 156
 157              1:
 158              1 or 1 = 1
 159              1 or 0 = 1
 160              1 or SYSMIS = 1
 161
 162              2:
 163              0 or 1 = 1
 164              SYSMIS or 1 = 1
 165
 166              3:
 167              0 or SYSMIS = SYSMIS
 168              SYSMIS or SYSMIS = SYSMIS
 169
 170              4:
 171              0 or 0 = 0
 172              SYSMIS or 0 = SYSMIS
 173
 174            */
 175           sp--;
 176           if (sp[0].f == 1.0);  /* 1 */
 177           else if (sp[1].f == 1.0)
 178             sp->f = 1.0;        /* 2 */
 179           else if (sp[1].f == SYSMIS)
 180             sp->f = SYSMIS;     /* 3 */
 181           break;
 182         case OP_NOT:
 183           if (sp[0].f == 0.0)
 184             sp->f = 1.0;
 185           else if (sp[0].f == 1.0)
 186             sp->f = 0.0;
 187           break;
 188
 189         case OP_EQ:
 190           sp--;
 191           if (sp[0].f != SYSMIS)
 192             {
 193               if (sp[1].f == SYSMIS)
 194                 sp->f = SYSMIS;
 195               else
 196                 sp->f = approx_eq (sp[0].f, sp[1].f);
 197             }
 198           break;
 199         case OP_GE:
 200           sp--;
 201           if (sp[0].f != SYSMIS)
 202             {
 203               if (sp[1].f == SYSMIS)
 204                 sp->f = SYSMIS;
 205               else
 206                 sp->f = approx_ge (sp[0].f, sp[1].f);
 207             }
 208           break;
 209         case OP_GT:
 210           sp--;
 211           if (sp[0].f != SYSMIS)
 212             {
 213               if (sp[1].f == SYSMIS)
 214                 sp->f = SYSMIS;
 215               else
 216                 sp->f = approx_gt (sp[0].f, sp[1].f);
 217             }
 218           break;
 219         case OP_LE:
 220           sp--;
 221           if (sp[0].f != SYSMIS)
 222             {
 223               if (sp[1].f == SYSMIS)
 224                 sp->f = SYSMIS;
 225               else
 226                 sp->f = approx_le (sp[0].f, sp[1].f);
 227             }
 228           break;
 229         case OP_LT:
 230           sp--;
 231           if (sp[0].f != SYSMIS)
 232             {
 233               if (sp[1].f == SYSMIS)
 234                 sp->f = SYSMIS;
 235               else
 236                 sp->f = approx_lt (sp[0].f, sp[1].f);
 237             }
 238           break;
 239         case OP_NE:
 240           sp--;
 241           if (sp[0].f != SYSMIS)
 242             {
 243               if (sp[1].f == SYSMIS)
 244                 sp->f = SYSMIS;
 245               else
 246                 sp->f = approx_ne (sp[0].f, sp[1].f);
 247             }
 248           break;
 249
 250           /* String operators. */
 251         case OP_STRING_EQ:
 252           sp--;
 253           sp[0].f = st_compare_pad (&sp[0].c[1], sp[0].c[0],
 254                                     &sp[1].c[1], sp[1].c[0]) == 0;
 255           break;
 256         case OP_STRING_GE:
 257           sp--;
 258           sp[0].f = st_compare_pad (&sp[0].c[1], sp[0].c[0],
 259                                     &sp[1].c[1], sp[1].c[0]) >= 0;
 260           break;
 261         case OP_STRING_GT:
 262           sp--;
 263           sp[0].f = st_compare_pad (&sp[0].c[1], sp[0].c[0],
 264                                     &sp[1].c[1], sp[1].c[0]) > 0;
 265           break;
 266         case OP_STRING_LE:
 267           sp--;
 268           sp[0].f = st_compare_pad (&sp[0].c[1], sp[0].c[0],
 269                                     &sp[1].c[1], sp[1].c[0]) <= 0;
 270           break;
 271         case OP_STRING_LT:
 272           sp--;
 273           sp[0].f = st_compare_pad (&sp[0].c[1], sp[0].c[0],
 274                                     &sp[1].c[1], sp[1].c[0]) < 0;
 275           break;
 276         case OP_STRING_NE:
 277           sp--;
 278           sp[0].f = st_compare_pad (&sp[0].c[1], sp[0].c[0],
 279                                     &sp[1].c[1], sp[1].c[0]) != 0;
 280           break;
 281
 282           /* Unary functions. */
 283         case OP_NEG:
 284           if (sp->f != SYSMIS)
 285             sp->f = -sp->f;
 286           break;
 287         case OP_ABS:
 288           if (sp->f != SYSMIS)
 289             sp->f = fabs (sp->f);
 290           break;
 291         case OP_ARCOS:
 292           if (sp->f != SYSMIS)
 293             {
 294               errno = 0;
 295               sp->f = acos (sp->f);
 296               if (errno)
 297                 sp->f = SYSMIS;
 298             }
 299           break;
 300         case OP_ARSIN:
 301           if (sp->f != SYSMIS)
 302             {
 303               errno = 0;
 304               sp->f = asin (sp->f);
 305               if (errno)
 306                 sp->f = SYSMIS;
 307             }
 308           break;
 309         case OP_ARTAN:
 310           if (sp->f != SYSMIS)
 311             sp->f = atan (sp->f);
 312           break;
 313         case OP_COS:
 314           if (sp->f != SYSMIS)
 315             sp->f = cos (sp->f);
 316           break;
 317         case OP_EXP:
 318           if (sp->f != SYSMIS)
 319             {
 320               errno = 0;
 321               sp->f = exp (sp->f);
 322               if (errno)
 323                 sp->f = SYSMIS;
 324             }
 325           break;
 326         case OP_LG10:
 327           if (sp->f != SYSMIS)
 328             {
 329               errno = 0;
 330               sp->f = log10 (sp->f);
 331               if (errno)
 332                 sp->f = SYSMIS;
 333             }
 334           break;
 335         case OP_LN:
 336           if (sp->f != SYSMIS)
 337             {
 338               errno = 0;
 339               sp->f = log10 (sp->f);
 340               if (errno)
 341                 sp->f = SYSMIS;
 342             }
 343           break;
 344         case OP_MOD10:
 345           if (sp->f != SYSMIS)
 346             sp->f = fmod (sp->f, 10);
 347           break;
 348         case OP_RND:
 349           if (sp->f != SYSMIS)
 350             {
 351               if (sp->f >= 0.0)
 352                 sp->f = floor (sp->f + 0.5);
 353               else
 354                 sp->f = -floor (-sp->f + 0.5);
 355             }
 356           break;
 357         case OP_SIN:
 358           if (sp->f != SYSMIS)
 359             sp->f = sin (sp->f);
 360           break;
 361         case OP_SQRT:
 362           if (sp->f != SYSMIS)
 363             {
 364               errno = 0;
 365               sp->f = sqrt (sp->f);
 366               if (errno)
 367                 sp->f = SYSMIS;
 368             }
 369           break;
 370         case OP_TAN:
 371           if (sp->f != SYSMIS)
 372             {
 373               errno = 0;
 374               sp->f = tan (sp->f);
 375               if (errno)
 376                 sp->f = SYSMIS;
 377             }
 378           break;
 379         case OP_TRUNC:
 380           if (sp->f != SYSMIS)
 381             {
 382               if (sp->f >= 0.0)
 383                 sp->f = floor (sp->f);
 384               else
 385                 sp->f = -floor (-sp->f);
 386             }
 387           break;
 388
 389           /* N-ary numeric functions. */
 390         case OP_ANY:
 391           {
 392             int n_args = *op++;
 393             int sysmis = 1;
 394
 395             sp -= n_args - 1;
 396             if (sp->f == SYSMIS)
 397               break;
 398             for (i = 1; i <= n_args; i++)
 399               if (approx_eq (sp[0].f, sp[i].f))
 400                 {
 401                   sp->f = 1.0;
 402                   goto main_loop;
 403                 }
 404               else if (sp[i].f != SYSMIS)
 405                 sysmis = 0;
 406             sp->f = sysmis ? SYSMIS : 0.0;
 407           }
 408           break;
 409         case OP_ANY_STRING:
 410           {
 411             int n_args = *op++;
 412
 413             sp -= n_args - 1;
 414             for (i = 1; i <= n_args; i++)
 415               if (!st_compare_pad (&sp[0].c[1], sp[0].c[0],
 416                                    &sp[i].c[1], sp[i].c[0]))
 417                 {
 418                   sp->f = 1.0;
 419                   goto main_loop;
 420                 }
 421             sp->f = 0.0;
 422           }
 423           break;
 424         case OP_CFVAR:
 425           {
 426             int n_args = *op++;
 427             int nv = 0;
 428             double sum[2] =
 429             {0.0, 0.0};
 430
 431             sp -= n_args - 1;
 432             for (i = 0; i < n_args; i++)
 433               if (sp[i].f != SYSMIS)
 434                 {
 435                   nv++;
 436                   sum[0] += sp[i].f;
 437                   sum[1] += sp[i].f * sp[i].f;
 438                 }
 439             if (nv < *op++)
 440               sp->f = SYSMIS;
 441             else
 442               sp->f = calc_cfvar (sum, nv);
 443           }
 444           break;
 445         case OP_MAX:
 446           {
 447             int n_args = *op++;
 448             int nv = 0;
 449             double max = -DBL_MAX;
 450
 451             sp -= n_args - 1;
 452             for (i = 0; i < n_args; i++)
 453               if (sp[i].f != SYSMIS)
 454                 {
 455                   nv++;
 456                   if (sp[i].f > max)
 457                     max = sp[i].f;
 458                 }
 459             if (nv < *op++)
 460               sp->f = SYSMIS;
 461             else
 462               sp->f = max;
 463           }
 464           break;
 465         case OP_MEAN:
 466           {
 467             int n_args = *op++;
 468             int nv = 0;
 469             double sum[1] =
 470             {0.0};
 471
 472             sp -= n_args - 1;
 473             for (i = 0; i < n_args; i++)
 474               if (sp[i].f != SYSMIS)
 475                 {
 476                   nv++;
 477                   sum[0] += sp[i].f;
 478                 }
 479             if (nv < *op++)
 480               sp->f = SYSMIS;
 481             else
 482               sp->f = calc_mean (sum, nv);
 483           }
 484           break;
 485         case OP_MIN:
 486           {
 487             int n_args = *op++;
 488             int nv = 0;
 489             double min = DBL_MAX;
 490
 491             sp -= n_args - 1;
 492             for (i = 0; i < n_args; i++)
 493               if (sp[i].f != SYSMIS)
 494                 {
 495                   nv++;
 496                   if (sp[i].f < min)
 497                     min = sp[i].f;
 498                 }
 499             if (nv < *op++)
 500               sp->f = SYSMIS;
 501             else
 502               sp->f = min;
 503           }
 504           break;
 505         case OP_NMISS:
 506           {
 507             int n_args = *op++;
 508             int n_missing = 0;
 509
 510             sp -= n_args - 1;
 511             for (i = 0; i < n_args; i++)
 512               if (sp[i].f == SYSMIS)
 513                 n_missing++;
 514             sp->f = n_missing;
 515           }
 516           break;
 517         case OP_NVALID:
 518           {
 519             int n_args = *op++;
 520             int n_valid = 0;
 521
 522             sp -= n_args - 1;
 523             for (i = 0; i < n_args; i++)
 524               if (sp[i].f != SYSMIS)
 525                 n_valid++;
 526             sp->f = n_valid;
 527           }
 528           break;
 529         case OP_RANGE:
 530           {
 531             int n_args = *op++;
 532             int sysmis = 1;
 533
 534             sp -= n_args - 1;
 535             if (sp->f == SYSMIS)
 536               break;
 537             for (i = 1; i <= n_args; i += 2)
 538               if (sp[i].f == SYSMIS || sp[i + 1].f == SYSMIS)
 539                 continue;
 540               else if (approx_ge (sp[0].f, sp[i].f)
 541                        && approx_le (sp[0].f, sp[i + 1].f))
 542                 {
 543                   sp->f = 1.0;
 544                   goto main_loop;
 545                 }
 546               else
 547                 sysmis = 0;
 548             sp->f = sysmis ? SYSMIS : 0.0;
 549           }
 550           break;
 551         case OP_RANGE_STRING:
 552           {
 553             int n_args = *op++;
 554
 555             sp -= n_args - 1;
 556             for (i = 1; i <= n_args; i += 2)
 557               if (st_compare_pad (&sp[0].c[1], sp[0].c[0],
 558                                   &sp[i].c[1], sp[i].c[0]) >= 0
 559                   && st_compare_pad (&sp[0].c[1], sp[0].c[0],
 560                                      &sp[i + 1].c[1], sp[i + 1].c[0]) <= 0)
 561                 {
 562                   sp->f = 1.0;
 563                   goto main_loop;
 564                 }
 565             sp->f = 0.0;
 566           }
 567           break;
 568         case OP_SD:
 569           {
 570             int n_args = *op++;
 571             int nv = 0;
 572             double sum[2];
 573
 574             sum[0] = sum[1] = 0.0;
 575
 576             sp -= n_args - 1;
 577             for (i = 0; i < n_args; i++)
 578               if (sp[i].f != SYSMIS)
 579                 {
 580                   nv++;
 581                   sum[0] += sp[i].f;
 582                   sum[1] += sp[i].f * sp[i].f;
 583                 }
 584             if (nv < *op++)
 585               sp->f = SYSMIS;
 586             else
 587               sp->f = calc_stddev (calc_variance (sum, nv));
 588           }
 589           break;
 590         case OP_SUM:
 591           {
 592             int n_args = *op++;
 593             int nv = 0;
 594             double sum = 0.0;
 595
 596             sp -= n_args - 1;
 597             for (i = 0; i < n_args; i++)
 598               if (sp[i].f != SYSMIS)
 599                 {
 600                   nv++;
 601                   sum += sp[i].f;
 602                 }
 603             if (nv < *op++)
 604               sp->f = SYSMIS;
 605             else
 606               sp->f = sum;
 607           }
 608           break;
 609         case OP_VARIANCE:
 610           {
 611             int n_args = *op++;
 612             int nv = 0;
 613             double sum[2];
 614
 615             sum[0] = sum[1] = 0.0;
 616
 617             sp -= n_args - 1;
 618             for (i = 0; i < n_args; i++)
 619               if (sp[i].f != SYSMIS)
 620                 {
 621                   nv++;
 622                   sum[0] += sp[i].f;
 623                   sum[1] += sp[i].f * sp[i].f;
 624                 }
 625             if (nv < *op++)
 626               sp->f = SYSMIS;
 627             else
 628               sp->f = calc_variance (sum, nv);
 629           }
 630           break;
 631
 632           /* Time construction function. */
 633         case OP_TIME_HMS:
 634           sp -= 2;
 635           if (sp[0].f == SYSMIS || sp[1].f == SYSMIS || sp[2].f == SYSMIS)
 636             sp->f = SYSMIS;
 637           else
 638             sp->f = 60. * (60. * sp[0].f + sp[1].f) + sp[2].f;
 639           break;
 640
 641           /* Date construction functions. */
 642         case OP_DATE_DMY:
 643           sp -= 2;
 644           sp->f = yrmoda (sp[2].f, sp[1].f, sp[0].f);
 645           if (sp->f != SYSMIS)
 646             sp->f *= 60. * 60. * 24.;
 647           break;
 648         case OP_DATE_MDY:
 649           sp -= 2;
 650           sp->f = yrmoda (sp[2].f, sp[0].f, sp[1].f);
 651           if (sp->f != SYSMIS)
 652             sp->f *= 60. * 60. * 24.;
 653           break;
 654         case OP_DATE_MOYR:
 655           sp--;
 656           sp->f = yrmoda (sp[1].f, sp[0].f, 1);
 657           if (sp->f != SYSMIS)
 658             sp->f *= 60. * 60. * 24.;
 659           break;
 660         case OP_DATE_QYR:
 661           sp--;
 662           if (sp[0].f == SYSMIS)
 663             sp->f = SYSMIS;
 664           else
 665             {
 666               sp->f = yrmoda (sp[1].f, sp[0].f * 3 - 2, 1);
 667               if (sp->f != SYSMIS)
 668                 sp->f *= 60. * 60. * 24.;
 669             }
 670           break;
 671         case OP_DATE_WKYR:
 672           sp--;
 673           if (sp[0].f == SYSMIS)
 674             sp->f = SYSMIS;
 675           else
 676             {
 677               sp[1].f = yrmoda (sp[1].f, 1, 1);
 678               if (sp->f != SYSMIS)
 679                 sp[1].f = 60. * 60. * 24. * (sp[1].f + 7. * (floor (sp[0].f) - 1.));
 680               sp->f = sp[1].f;
 681             }
 682           break;
 683         case OP_DATE_YRDAY:
 684           sp--;
 685           if (sp[1].f == SYSMIS)
 686             sp->f = SYSMIS;
 687           else
 688             {
 689               sp->f = yrmoda (sp[0].f, 1, 1);
 690               if (sp->f != SYSMIS)
 691                 sp->f = 60. * 60. * 24. * (sp->f + floor (sp[1].f) - 1);
 692             }
 693           break;
 694         case OP_YRMODA:
 695           sp -= 2;
 696           sp->f = yrmoda (sp[0].f, sp[1].f, sp[2].f);
 697           break;
 698
 699           /* Date extraction functions. */
 700         case OP_XDATE_DATE:
 701           if (sp->f != SYSMIS)
 702             sp->f = floor (sp->f / 60. / 60. / 24.) * 60. * 60. * 24.;
 703           break;
 704         case OP_XDATE_HOUR:
 705           if (sp->f != SYSMIS)
 706             sp->f = fmod (floor (sp->f / 60. / 60.), 24.);
 707           break;
 708         case OP_XDATE_JDAY:
 709           if (sp->f != SYSMIS)
 710             sp->f = 86400. * julian_to_jday (sp->f / 86400.);
 711           break;
 712         case OP_XDATE_MDAY:
 713           if (sp->f != SYSMIS)
 714             {
 715               int day;
 716               julian_to_calendar (sp->f / 86400., NULL, NULL, &day);
 717               sp->f = day;
 718             }
 719           break;
 720         case OP_XDATE_MINUTE:
 721           if (sp->f != SYSMIS)
 722             sp->f = fmod (floor (sp->f / 60.), 60.);
 723           break;
 724         case OP_XDATE_MONTH:
 725           if (sp->f != SYSMIS)
 726             {
 727               int month;
 728               julian_to_calendar (sp->f / 86400., NULL, &month, NULL);
 729               sp->f = month;
 730             }
 731           break;
 732         case OP_XDATE_QUARTER:
 733           if (sp->f != SYSMIS)
 734             {
 735               int month;
 736               julian_to_calendar (sp->f / 86400., NULL, &month, NULL);
 737               sp->f = (month - 1) / 3 + 1;
 738             }
 739           break;
 740         case OP_XDATE_SECOND:
 741           if (sp->f != SYSMIS)
 742             sp->f = fmod (sp->f, 60.);
 743           break;
 744         case OP_XDATE_TDAY:
 745           if (sp->f != SYSMIS)
 746             sp->f = floor (sp->f / 60. / 60. / 24.);
 747           break;
 748         case OP_XDATE_TIME:
 749           if (sp->f != SYSMIS)
 750             sp->f -= floor (sp->f / 60. / 60. / 24.) * 60. * 60. * 24.;
 751           break;
 752         case OP_XDATE_WEEK:
 753           if (sp->f != SYSMIS)
 754             sp->f = (julian_to_jday (sp->f / 86400.) - 1) / 7 + 1;
 755           break;
 756         case OP_XDATE_WKDAY:
 757           if (sp->f != SYSMIS)
 758             sp->f = julian_to_wday (sp->f / 86400.);
 759           break;
 760         case OP_XDATE_YEAR:
 761           if (sp->f != SYSMIS)
 762             {
 763               int year;
 764               julian_to_calendar (sp->f / 86400., &year, NULL, NULL);
 765               sp->f = year;
 766             }
 767           break;
 768
 769           /* String functions. */
 770         case OP_CONCAT:
 771           {
 772             int n_args = *op++;
 773             unsigned char *dest;
 774
 775             dest = pool_alloc (e->pool, 256);
 776             dest[0] = 0;
 777
 778             sp -= n_args - 1;
 779             for (i = 0; i < n_args; i++)
 780               if (sp[i].c[0] != 0)
 781                 {
 782                   if (sp[i].c[0] + dest[0] < 255)
 783                     {
 784                       memcpy (&dest[dest[0] + 1], &sp[i].c[1], sp[i].c[0]);
 785                       dest[0] += sp[i].c[0];
 786                     }
 787                   else
 788                     {
 789                       memcpy (&dest[dest[0] + 1], &sp[i].c[1], 255 - dest[0]);
 790                       dest[0] = 255;
 791                       break;
 792                     }
 793                 }
 794             sp[0].c = dest;
 795           }
 796           break;
 797         case OP_INDEX:
 798           sp--;
 799           if (sp[1].c[0] == 0)
 800             sp->f = SYSMIS;
 801           else
 802             {
 803               int last = sp[0].c[0] - sp[1].c[0];
 804               for (i = 0; i <= last; i++)
 805                 if (!memcmp (&sp[0].c[i + 1], &sp[0].c[1], sp[0].c[0]))
 806                   {
 807                     sp->f = i + 1;
 808                     goto main_loop;
 809                   }
 810               sp->f = 0.0;
 811             }
 812           break;
 813         case OP_INDEX_OPT:
 814           {
 815             /* Length of each search string. */
 816             int part_len = sp[2].f;
 817
 818             sp -= 2;
 819             if (sp[1].c[0] == 0 || part_len <= 0 || sp[2].f == SYSMIS
 820                 || sp[1].c[0] % part_len != 0)
 821               sp->f = SYSMIS;
 822             else
 823               {
 824                 /* Last possible index. */
 825                 int last = sp[0].c[0] - part_len;
 826
 827                 for (i = 0; i <= last; i++)
 828                   for (j = 0; j < sp[1].c[0]; j += part_len)
 829                     if (!memcmp (&sp[0].c[i], &sp[1].c[j], part_len))
 830                       {
 831                         sp->f = i + 1;
 832                         goto main_loop;
 833                       }
 834                 sp->f = 0.0;
 835               }
 836           }
 837           break;
 838         case OP_RINDEX:
 839           sp--;
 840           if (sp[1].c[0] == 0)
 841             sp->f = SYSMIS;
 842           else
 843             {
 844               for (i = sp[0].c[0] - sp[1].c[0]; i >= 0; i--)
 845                 if (!memcmp (&sp[0].c[i + 1], &sp[0].c[1], sp[0].c[0]))
 846                   {
 847                     sp->f = i + 1;
 848                     goto main_loop;
 849                   }
 850               sp->f = 0.0;
 851             }
 852           break;
 853         case OP_RINDEX_OPT:
 854           {
 855             /* Length of each search string. */
 856             int part_len = sp[2].f;
 857
 858             sp -= 2;
 859             if (sp[1].c[0] == 0 || part_len <= 0 || sp[2].f == SYSMIS
 860                 || sp[1].c[0] % part_len != 0)
 861               sp->f = SYSMIS;
 862             else
 863               {
 864                 for (i = sp[0].c[0] - part_len; i >= 0; i--)
 865                   for (j = 0; j < sp[1].c[0]; j += part_len)
 866                     if (!memcmp (&sp[0].c[i], &sp[1].c[j], part_len))
 867                       {
 868                         sp->f = i + 1;
 869                         goto main_loop;
 870                       }
 871                 sp->f = 0.0;
 872               }
 873           }
 874           break;
 875         case OP_LENGTH:
 876           sp->f = sp[0].c[0];
 877           break;
 878         case OP_LOWER:
 879           for (i = sp[0].c[0]; i >= 1; i--)
 880             sp[0].c[i] = tolower ((unsigned char) (sp[0].c[i]));
 881           break;
 882         case OP_UPPER:
 883           for (i = sp[0].c[0]; i >= 1; i--)
 884             sp[0].c[i] = toupper ((unsigned char) (sp[0].c[i]));
 885           break;
 886         case OP_LPAD:
 887           {
 888             int len;
 889             sp--;
 890             len = sp[1].f;
 891             if (sp[1].f == SYSMIS || len < 0 || len > 255)
 892               sp->c[0] = 0;
 893             else if (len > sp[0].c[0])
 894               {
 895                 unsigned char *dest;
 896
 897                 dest = pool_alloc (e->pool, len + 1);
 898                 dest[0] = len;
 899                 memset (&dest[1], ' ', len - sp->c[0]);
 900                 memcpy (&dest[len - sp->c[0] + 1], &sp->c[1], sp->c[0]);
 901                 sp->c = dest;
 902               }
 903           }
 904           break;
 905         case OP_LPAD_OPT:
 906           {
 907             int len;
 908             sp -= 2;
 909             len = sp[1].f;
 910             if (sp[1].f == SYSMIS || len < 0 || len > 255 || sp[2].c[0] != 1)
 911               sp->c[0] = 0;
 912             else if (len > sp[0].c[0])
 913               {
 914                 unsigned char *dest;
 915
 916                 dest = pool_alloc (e->pool, len + 1);
 917                 dest[0] = len;
 918                 memset (&dest[1], sp[2].c[1], len - sp->c[0]);
 919                 memcpy (&dest[len - sp->c[0] + 1], &sp->c[1], sp->c[0]);
 920                 sp->c = dest;
 921               }
 922           }
 923           break;
 924         case OP_RPAD:
 925           {
 926             int len;
 927             sp--;
 928             len = sp[1].f;
 929             if (sp[1].f == SYSMIS || len < 0 || len > 255)
 930               sp->c[0] = 0;
 931             else if (len > sp[0].c[0])
 932               {
 933                 unsigned char *dest;
 934
 935                 dest = pool_alloc (e->pool, len + 1);
 936                 dest[0] = len;
 937                 memcpy (&dest[1], &sp->c[1], sp->c[0]);
 938                 memset (&dest[sp->c[0] + 1], ' ', len - sp->c[0]);
 939                 sp->c = dest;
 940               }
 941           }
 942           break;
 943         case OP_RPAD_OPT:
 944           {
 945             int len;
 946             sp -= 2;
 947             len = sp[1].f;
 948             if (len < 0 || len > 255 || sp[2].c[0] != 1)
 949               sp->c[0] = 0;
 950             else if (len > sp[0].c[0])
 951               {
 952                 unsigned char *dest;
 953
 954                 dest = pool_alloc (e->pool, len + 1);
 955                 dest[0] = len;
 956                 memcpy (&dest[1], &sp->c[1], sp->c[0]);
 957                 memset (&dest[sp->c[0] + 1], sp[2].c[1], len - sp->c[0]);
 958                 sp->c = dest;
 959               }
 960           }
 961           break;
 962         case OP_LTRIM:
 963           {
 964             int len = sp[0].c[0];
 965
 966             i = 1;
 967             while (i <= len && sp[0].c[i] == ' ')
 968               i++;
 969             if (--i)
 970               {
 971                 sp[0].c[i] = sp[0].c[0] - i;
 972                 sp->c = &sp[0].c[i];
 973               }
 974           }
 975           break;
 976         case OP_LTRIM_OPT:
 977           {
 978             sp--;
 979             if (sp[1].c[0] != 1)
 980               sp[0].c[0] = 0;
 981             else
 982               {
 983                 int len = sp[0].c[0];
 984                 int cmp = sp[1].c[1];
 985
 986                 i = 1;
 987                 while (i <= len && sp[0].c[i] == cmp)
 988                   i++;
 989                 if (--i)
 990                   {
 991                     sp[0].c[i] = sp[0].c[0] - i;
 992                     sp->c = &sp[0].c[i];
 993                   }
 994               }
 995           }
 996           break;
 997         case OP_RTRIM:
 998           assert (' ' != 0);
 999           while (sp[0].c[sp[0].c[0]] == ' ')
1000             sp[0].c[0]--;
1001           break;
1002         case OP_RTRIM_OPT:
1003           sp--;
1004           if (sp[1].c[0] != 1)
1005             sp[0].c[0] = 0;
1006           else
1007             {
1008               /* Note that NULs are not allowed in strings.  This code
1009                  needs to change if this decision is changed. */
1010               int cmp = sp[1].c[1];
1011               while (sp[0].c[sp[0].c[0]] == cmp)
1012                 sp[0].c[0]--;
1013             }
1014           break;
1015         case OP_NUMBER:
1016           {
1017             struct data_in di;
1018
1019             di.s = &sp->c[1];
1020             di.e = &sp->c[1] + sp->c[0];
1021             di.v = sp;
1022             di.flags = DI_IGNORE_ERROR;
1023             di.f1 = 1;
1024             di.format.type = FMT_F;
1025             di.format.w = sp->c[0];
1026             di.format.d = 0;
1027             data_in (&di);
1028           }
1029           break;
1030         case OP_NUMBER_OPT:
1031           {
1032             struct data_in di;
1033             di.s = &sp->c[1];
1034             di.e = &sp->c[1] + sp->c[0];
1035             di.v = sp;
1036             di.flags = DI_IGNORE_ERROR;
1037             di.f1 = 1;
1038             di.format.type = *op++;
1039             di.format.w = *op++;
1040             di.format.d = *op++;
1041             data_in (&di);
1042           }
1043           break;
1044         case OP_STRING:
1045           {
1046             struct fmt_spec f;
1047             unsigned char *dest;
1048
1049             f.type = *op++;
1050             f.w = *op++;
1051             f.d = *op++;
1052
1053             dest = pool_alloc (e->pool, f.w + 1);
1054             dest[0] = f.w;
1055
1056             assert ((formats[f.type].cat & FCAT_STRING) == 0);
1057             data_out (&dest[1], &f, sp);
1058             sp->c = dest;
1059           }
1060           break;
1061         case OP_SUBSTR:
1062           {
1063             int index;
1064
1065             sp--;
1066             index = sp[1].f;
1067             if (index < 1 || index > sp[0].c[0])
1068               sp->c[0] = 0;
1069             else if (index > 1)
1070               {
1071                 index--;
1072                 sp->c[index] = sp->c[0] - index;
1073                 sp->c += index;
1074               }
1075           }
1076           break;
1077         case OP_SUBSTR_OPT:
1078           {
1079             int index;
1080             int n;
1081
1082             sp -= 2;
1083             index = sp[1].f;
1084             n = sp[2].f;
1085             if (sp[1].f == SYSMIS || sp[2].f == SYSMIS || index < 1
1086                 || index > sp[0].c[0] || n < 1)
1087               sp->c[0] = 0;
1088             else
1089               {
1090                 if (index > 1)
1091                   {
1092                     index--;
1093                     sp->c[index] = sp->c[0] - index;
1094                     sp->c += index;
1095                   }
1096                 if (sp->c[0] > n)
1097                   sp->c[0] = n;
1098               }
1099           }
1100           break;
1101
1102           /* Artificial. */
1103         case OP_INV:
1104           if (sp->f != SYSMIS)
1105             sp->f = 1. / sp->f;
1106           break;
1107         case OP_SQUARE:
1108           if (sp->f != SYSMIS)
1109             sp->f *= sp->f;
1110           break;
1111         case OP_NUM_TO_BOOL:
1112           if (approx_eq (sp->f, 0.0))
1113             sp->f = 0.0;
1114           else if (approx_eq (sp->f, 1.0))
1115             sp->f = 1.0;
1116           else if (sp->f != SYSMIS)
1117             {
1118               msg (SE, _("A number being treated as a Boolean in an "
1119                          "expression was found to have a value other than "
1120                          "0 (false), 1 (true), or the system-missing value.  "
1121                          "The result was forced to 0."));
1122               sp->f = 0.0;
1123             }
1124           break;
1125
1126           /* Weirdness. */
1127         case OP_MOD:
1128           sp--;
1129           if (sp[0].f != SYSMIS)
1130             {
1131               if (sp[1].f == SYSMIS)
1132                 {
1133                   if (approx_ne (sp[0].f, 0.0))
1134                     sp->f = SYSMIS;
1135                 }
1136               else
1137                 sp->f = fmod (sp[0].f, sp[1].f);
1138             }
1139           break;
1140         case OP_NORMAL:
1141           if (sp->f != SYSMIS)
1142             sp->f *= rng_get_double_normal (pspp_rng ());
1143           break;
1144         case OP_UNIFORM:
1145           if (sp->f != SYSMIS)
1146             sp->f *= rng_get_double (pspp_rng ());
1147           break;
1148         case OP_SYSMIS:
1149           if (sp[0].f == SYSMIS || !finite (sp[0].f))
1150             sp->f = 1.0;
1151           else
1152             sp->f = 0.0;
1153           break;
1154         case OP_VEC_ELEM_NUM:
1155           {
1156             int rindx = sp[0].f + EPSILON;
1157             const struct vector *v = dict_get_vector (default_dict, *op++);
1158
1159             if (sp[0].f == SYSMIS || rindx < 1 || rindx > v->cnt)
1160               {
1161                 if (sp[0].f == SYSMIS)
1162                   msg (SE, _("SYSMIS is not a valid index value for vector "
1163                              "%s.  The result will be set to SYSMIS."),
1164                        v->name);
1165                 else
1166                   msg (SE, _("%g is not a valid index value for vector %s.  "
1167                              "The result will be set to SYSMIS."),
1168                        sp[0].f, v->name);
1169                 sp->f = SYSMIS;
1170                 break;
1171               }
1172             sp->f = c->data[v->var[rindx - 1]->fv].f;
1173           }
1174           break;
1175         case OP_VEC_ELEM_STR:
1176           {
1177             int rindx = sp[0].f + EPSILON;
1178             const struct vector *vect = dict_get_vector (default_dict, *op++);
1179             struct variable *v;
1180
1181             if (sp[0].f == SYSMIS || rindx < 1 || rindx > vect->cnt)
1182               {
1183                 if (sp[0].f == SYSMIS)
1184                   msg (SE, _("SYSMIS is not a valid index value for vector "
1185                              "%s.  The result will be set to the empty "
1186                              "string."),
1187                        vect->name);
1188                 else
1189                   msg (SE, _("%g is not a valid index value for vector %s.  "
1190                              "The result will be set to the empty string."),
1191                        sp[0].f, vect->name);
1192                 sp->c = pool_alloc (e->pool, 1);
1193                 sp->c[0] = 0;
1194                 break;
1195               }
1196
1197             v = vect->var[rindx - 1];
1198             sp->c = pool_alloc (e->pool, v->width + 1);
1199             sp->c[0] = v->width;
1200             memcpy (&sp->c[1], c->data[v->fv].s, v->width);
1201           }
1202           break;
1203
1204           /* Terminals. */
1205         case OP_NUM_CON:
1206           sp++;
1207           sp->f = *dbl++;
1208           break;
1209         case OP_STR_CON:
1210           sp++;
1211           sp->c = pool_alloc (e->pool, *str + 1);
1212           memcpy (sp->c, str, *str + 1);
1213           str += *str + 1;
1214           break;
1215         case OP_NUM_VAR:
1216           sp++;
1217           sp->f = c->data[(*vars)->fv].f;
1218           if (is_num_user_missing (sp->f, *vars))
1219             sp->f = SYSMIS;
1220           vars++;
1221           break;
1222         case OP_STR_VAR:
1223           {
1224             int width = (*vars)->width;
1225
1226             sp++;
1227             sp->c = pool_alloc (e->pool, width + 1);
1228             sp->c[0] = width;
1229             memcpy (&sp->c[1], &c->data[(*vars)->fv], width);
1230             vars++;
1231           }
1232           break;
1233         case OP_NUM_LAG:
1234           {
1235             struct ccase *c = lagged_case (*op++);
1236
1237             sp++;
1238             if (c == NULL)
1239               sp->f = SYSMIS;
1240             else
1241               {
1242                 sp->f = c->data[(*vars)->fv].f;
1243                 if (is_num_user_missing (sp->f, *vars))
1244                   sp->f = SYSMIS;
1245               }
1246             vars++;
1247             break;
1248           }
1249         case OP_STR_LAG:
1250           {
1251             struct ccase *c = lagged_case (*op++);
1252             int width = (*vars)->width;
1253
1254             sp++;
1255             sp->c = pool_alloc (e->pool, width + 1);
1256             sp->c[0] = width;
1257
1258             if (c == NULL)
1259               memset (sp->c, ' ', width);
1260             else
1261               memcpy (&sp->c[1], &c->data[(*vars)->fv], width);
1262
1263             vars++;
1264           }
1265           break;
1266         case OP_NUM_SYS:
1267           sp++;
1268           sp->f = c->data[*op++].f == SYSMIS;
1269           break;
1270         case OP_STR_MIS:
1271           sp++;
1272           sp->f = is_str_user_missing (c->data[(*vars)->fv].s, *vars);
1273           vars++;
1274           break;
1275         case OP_NUM_VAL:
1276           sp++;
1277           sp->f = c->data[*op++].f;
1278           break;
1279         case OP_CASENUM:
1280           sp++;
1281           sp->f = vfm_sink_info.ncases + 1;
1282           break;
1283
1284         case OP_SENTINEL:
1285           goto finished;
1286
1287         default:
1288           assert (0);
1289         }
1290
1291     main_loop: ;
1292     }
1293 finished:
1294   if (e->type != EX_STRING)
1295     {
1296       double value = sp->f;
1297       if (!finite (value))
1298         value = SYSMIS;
1299       if (v)
1300         v->f = value;
1301       return value;
1302     }
1303   else
1304     {
1305       assert (v);
1306
1307       v->c = sp->c;
1308
1309       return 0.0;
1310     }
1311 }