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