src/expr-evl.c

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