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