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