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