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