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