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)->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 = rand_normal (sp->f);
1205           break;
1206         case OP_UNIFORM:
1207           if (sp->f != SYSMIS)
1208             sp->f = rand_uniform (sp->f);
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             struct vector *v = &vec[*op++];
1220
1221             if (sp[0].f == SYSMIS || rindx < 1 || rindx > v->nv)
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->v[rindx - 1]->fv].f;
1235           }
1236           break;
1237         case OP_VEC_ELEM_STR:
1238           {
1239             int rindx = sp[0].f + EPSILON;
1240             struct vector *vect = &vec[*op++];
1241             struct variable *v;
1242
1243             if (sp[0].f == SYSMIS || rindx < 1 || rindx > vect->nv)
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->v[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 #if __CHECKER__
1355           /* This case prevents Checker from choking. */
1356         case 42000:
1357           assert (0);
1358 #endif
1359
1360         default:
1361 #if GLOBAL_DEBUGGING
1362           printf (_("evaluate_expression(): not implemented: %s\n"),
1363                   ops[op[-1]].name);
1364 #else
1365           printf (_("evaluate_expression(): not implemented: %d\n"), op[-1]);
1366 #endif
1367           assert (0);
1368         }
1369
1370     main_loop: ;
1371     }
1372 finished:
1373   if (e->type != EX_STRING)
1374     {
1375       double value = sp->f;
1376       if (!finite (value))
1377         value = SYSMIS;
1378       if (v)
1379         v->f = value;
1380       return value;
1381     }
1382   else
1383     {
1384       assert (v);
1385
1386 #if PAGED_STACK
1387       memcpy (e->str_stack, sp->c, sp->c[0] + 1);
1388       v->c = e->str_stack;
1389 #else
1390       v->c = sp->c;
1391 #endif
1392
1393       return 0.0;
1394     }
1395 }