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