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