src/expr-opt.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 #include <assert.h>
  22 #include <math.h>
  23 #include <ctype.h>
  24 #include <errno.h>
  25 #include <stdlib.h>
  26 #include "alloc.h"
  27 #include "approx.h"
  28 #include "data-in.h"
  29 #include "error.h"
  30 #include "expr.h"
  31 #include "exprP.h"
  32 #include "julcal/julcal.h"
  33 #include "misc.h"
  34 #include "stats.h"
  35 #include "str.h"
  36 #include "var.h"
  37
  38 /*
  39    Expression "optimizer"
  40
  41    Operates on the tree representation of expressions.
  42    optimize_expression() performs the optimizations listed below:
  43
  44    1. Constant folding
  45    Any operation with constant operands is replaced by its value.
  46    (Exception: random-number-generator functions.)
  47
  48    2. Strength reduction (x is any expression; a is a numeric constant)
  49    x/0 => SYSMIS
  50    x*0 => 0
  51    x**0 => 1
  52    x**1, x+0, x-0, x*1 => x
  53    x**2 => sqr(x)
  54    x/a => x*(1/a)   (where 1/a is evaluated at optimization time)
  55
  56    I thought about adding additional optimizations but decided that what
  57    is here could already be considered overkill.
  58  */
  59
  60 static struct nonterm_node *evaluate_tree (struct nonterm_node * n);
  61 static struct nonterm_node *optimize_tree (struct nonterm_node * n);
  62
  63 struct nonterm_node *
  64 optimize_expression (struct nonterm_node * n)
  65 {
  66   int i;
  67
  68   /* Set to 1 if a child is nonconstant. */
  69   int nonconst = 0;
  70
  71   /* Number of system-missing children. */
  72   int sysmis = 0;
  73
  74   /* We can't optimize a terminal node. */
  75   if (n->type > OP_TERMINAL)
  76     return n;
  77
  78   /* Start by optimizing all the children. */
  79   for (i = 0; i < n->n; i++)
  80     {
  81       n->arg[i] = ((union any_node *)
  82                    optimize_expression ((struct nonterm_node *) n->arg[i]));
  83       if (n->arg[i]->type == OP_NUM_CON)
  84         {
  85           if (n->arg[i]->num_con.value == SYSMIS)
  86             sysmis++;
  87         }
  88       else if (n->arg[i]->type != OP_STR_CON)
  89         nonconst = 1;
  90     }
  91
  92   if (sysmis && !(ops[n->type].flags & OP_ABSORB_MISS))
  93     /* Just about any operation produces SYSMIS when given any SYSMIS
  94        arguments. */
  95     {
  96       struct num_con_node *num = xmalloc (sizeof *num);
  97       free_node ((union any_node *) n);
  98       num->type = OP_NUM_CON;
  99       num->value = SYSMIS;
 100       n = (struct nonterm_node *) num;
 101     }
 102   else if (!nonconst)
 103     /* If all the children of this node are constants, then there are
 104        obvious optimizations. */
 105     n = evaluate_tree (n);
 106   else
 107     /* Otherwise, we may be able to make certain optimizations
 108        anyway. */
 109     n = optimize_tree (n);
 110   return n;
 111 }
 112
 113 static struct nonterm_node *repl_num_con (struct nonterm_node *, double);
 114 static struct nonterm_node *force_repl_num_con (struct nonterm_node *, double);
 115 static struct nonterm_node *repl_str_con (struct nonterm_node *, char *, int);
 116
 117 #define n0 n->arg[0]->num_con.value
 118 #define n1 n->arg[1]->num_con.value
 119 #define n2 n->arg[2]->num_con.value
 120
 121 #define s0 n->arg[0]->str_con.s
 122 #define s0l n->arg[0]->str_con.len
 123 #define s1 n->arg[1]->str_con.s
 124 #define s1l n->arg[1]->str_con.len
 125 #define s2 n->arg[2]->str_con.s
 126 #define s2l n->arg[2]->str_con.len
 127 #define s(X) n->arg[X]->str_con.s
 128 #define sl(X) n->arg[X]->str_con.len
 129
 130 static struct nonterm_node *
 131 optimize_tree (struct nonterm_node * n)
 132 {
 133   int i;
 134
 135   errno = 0;
 136   if (n->type == OP_PLUS || n->type == OP_MUL)
 137     {
 138       /* Default constant value. */
 139       double def = n->type == OP_MUL ? 1.0 : 0.0;
 140
 141       /* Total value of all the constants. */
 142       double cval = def;
 143
 144       /* Number of nonconst arguments. */
 145       int nvar = 0;
 146
 147       /* New node. */
 148       struct nonterm_node *m;
 149
 150       /* Argument copying counter. */
 151       int c;
 152
 153       /* 1=SYSMIS encountered */
 154       int sysmis = 0;
 155
 156       for (i = 0; i < n->n; i++)
 157         if (n->arg[i]->type == OP_NUM_CON)
 158           {
 159             if (n->arg[i]->num_con.value != SYSMIS)
 160               {
 161                 if (n->type == OP_MUL)
 162                   cval *= n->arg[i]->num_con.value;
 163                 else
 164                   cval += n->arg[i]->num_con.value;
 165               }
 166             else
 167               sysmis++;
 168           }
 169         else
 170           nvar++;
 171
 172       /* 0*SYSMIS=0, 0/SYSMIS=0; otherwise, SYSMIS and infinities
 173          produce SYSMIS. */
 174       if (approx_eq (cval, 0.0) && n->type == OP_MUL)
 175         nvar = 0;
 176       else if (sysmis || !finite (cval))
 177         {
 178           nvar = 0;
 179           cval = SYSMIS;
 180         }
 181
 182       /* If no nonconstant terms, replace with a constant node. */
 183       if (nvar == 0)
 184         return force_repl_num_con (n, cval);
 185
 186       if (nvar == 1 && cval == def)
 187         {
 188           /* If there is exactly one nonconstant term and no constant
 189              terms, replace with the nonconstant term. */
 190           for (i = 0; i < n->n; i++)
 191             if (n->arg[i]->type != OP_NUM_CON)
 192               m = (struct nonterm_node *) n->arg[i];
 193             else
 194               free_node (n->arg[i]);
 195         }
 196       else
 197         {
 198           /* Otherwise consolidate all the nonconstant terms. */
 199           m = xmalloc (sizeof (struct nonterm_node)
 200                        + ((nvar + approx_ne (cval, def) - 1)
 201                           * sizeof (union any_node *)));
 202           for (i = c = 0; i < n->n; i++)
 203             if (n->arg[i]->type != OP_NUM_CON)
 204               m->arg[c++] = n->arg[i];
 205             else
 206               free_node (n->arg[i]);
 207
 208           if (approx_ne (cval, def))
 209             {
 210               m->arg[c] = xmalloc (sizeof (struct num_con_node));
 211               m->arg[c]->num_con.type = OP_NUM_CON;
 212               m->arg[c]->num_con.value = cval;
 213               c++;
 214             }
 215
 216           m->type = n->type;
 217           m->n = c;
 218         }
 219       free (n);
 220       n = m;
 221     }
 222   else if (n->type == OP_POW)
 223     {
 224       if (n->arg[1]->type == OP_NUM_CON)
 225         {
 226           if (approx_eq (n1, 1.0))
 227             {
 228               struct nonterm_node *m = (struct nonterm_node *) n->arg[0];
 229
 230               free_node (n->arg[1]);
 231               free (n);
 232               return m;
 233             }
 234           else if (approx_eq (n1, 2.0))
 235             {
 236               n = xrealloc (n, sizeof (struct nonterm_node));
 237               n->type = OP_SQUARE;
 238               n->n = 1;
 239             }
 240         }
 241     }
 242   return n;
 243 }
 244
 245 #define rnc(D)                                  \
 246         (n = repl_num_con (n, D))
 247
 248 #define frnc(D)                                 \
 249         (n = force_repl_num_con (n, D))
 250
 251 /* Finds the first NEEDLE of length NEEDLE_LEN in a HAYSTACK of length
 252    HAYSTACK_LEN.  Returns a 1-based index, 0 on failure. */
 253 static inline int
 254 str_search (char *haystack, int haystack_len, char *needle, int needle_len)
 255 {
 256   char *p = memmem (haystack, haystack_len, needle, needle_len);
 257   return p ? p - haystack + 1 : 0;
 258 }
 259
 260 /* Finds the last NEEDLE of length NEEDLE_LEN in a HAYSTACK of length
 261    HAYSTACK_LEN.  Returns a 1-based index, 0 on failure. */
 262 static inline int
 263 str_rsearch (char *haystack, int haystack_len, char *needle, int needle_len)
 264 {
 265   char *p = mm_find_reverse (haystack, haystack_len, needle, needle_len);
 266   return p ? p - haystack + 1 : 0;
 267 }
 268
 269 static struct nonterm_node *
 270 evaluate_tree (struct nonterm_node * n)
 271 {
 272   static char *strbuf;
 273   int add;
 274   int len;
 275   int i;
 276
 277   if (!strbuf)
 278     strbuf = xmalloc (256);
 279   errno = 0;
 280
 281   switch (n->type)
 282     {
 283     case OP_PLUS:
 284     case OP_MUL:
 285       return optimize_tree (n);
 286
 287     case OP_POW:
 288       if (approx_eq (n0, 0.0) && approx_eq (n1, 0.0))
 289         frnc (SYSMIS);
 290       else if (n0 == SYSMIS && n1 == 0.0)
 291         frnc (1.0);
 292       else if (n0 == 0.0 && n1 == SYSMIS)
 293         frnc (0.0);
 294       else
 295         rnc (pow (n0, n1));
 296       break;
 297
 298     case OP_AND:
 299       if (n0 == 0.0 || n1 == 0.0)
 300         frnc (0.0);
 301       else if (n0 == SYSMIS || n1 == SYSMIS)
 302         frnc (SYSMIS);
 303       else
 304         frnc (1.0);
 305       break;
 306     case OP_OR:
 307       if (n0 == 1.0 || n1 == 1.0)
 308         frnc (1.0);
 309       else if (n0 == SYSMIS || n1 == SYSMIS)
 310         frnc (SYSMIS);
 311       else
 312         frnc (0.0);
 313       break;
 314     case OP_NOT:
 315       rnc (n0 == 0.0 ? 1.0 : 0.0);
 316       break;
 317
 318     case OP_EQ:
 319       rnc (approx_eq (n0, n1));
 320       break;
 321     case OP_GE:
 322       rnc (approx_ge (n0, n1));
 323       break;
 324     case OP_GT:
 325       rnc (approx_gt (n0, n1));
 326       break;
 327     case OP_LE:
 328       rnc (approx_le (n0, n1));
 329       break;
 330     case OP_LT:
 331       rnc (approx_lt (n0, n1));
 332       break;
 333     case OP_NE:
 334       rnc (approx_ne (n0, n1));
 335       break;
 336
 337       /* String operators. */
 338     case OP_STRING_EQ:
 339       rnc (st_compare_pad (s0, s0l, s1, s1l) == 0);
 340       break;
 341     case OP_STRING_GE:
 342       rnc (st_compare_pad (s0, s0l, s1, s1l) >= 0);
 343       break;
 344     case OP_STRING_GT:
 345       rnc (st_compare_pad (s0, s0l, s1, s1l) > 0);
 346       break;
 347     case OP_STRING_LE:
 348       rnc (st_compare_pad (s0, s0l, s1, s1l) <= 0);
 349       break;
 350     case OP_STRING_LT:
 351       rnc (st_compare_pad (s0, s0l, s1, s1l) < 0);
 352       break;
 353     case OP_STRING_NE:
 354       rnc (st_compare_pad (s0, s0l, s1, s1l) != 0);
 355       break;
 356
 357       /* Unary functions. */
 358     case OP_NEG:
 359       rnc (-n0);
 360       break;
 361     case OP_ABS:
 362       rnc (fabs (n0));
 363       break;
 364     case OP_ARCOS:
 365       rnc (acos (n0));
 366       break;
 367     case OP_ARSIN:
 368       rnc (asin (n0));
 369       break;
 370     case OP_ARTAN:
 371       rnc (atan (n0));
 372       break;
 373     case OP_COS:
 374       rnc (cos (n0));
 375       break;
 376     case OP_EXP:
 377       rnc (exp (n0));
 378       break;
 379     case OP_LG10:
 380       rnc (log10 (n0));
 381       break;
 382     case OP_LN:
 383       rnc (log (n0));
 384       break;
 385     case OP_MOD10:
 386       rnc (fmod (n0, 10));
 387       break;
 388     case OP_RND:
 389       rnc (n0 >= 0.0 ? floor (n0 + 0.5) : -floor (-n0 + 0.5));
 390       break;
 391     case OP_SIN:
 392       rnc (sin (n0));
 393       break;
 394     case OP_SQRT:
 395       rnc (sqrt (n0));
 396       break;
 397     case OP_TAN:
 398       rnc (tan (n0));
 399       break;
 400     case OP_TRUNC:
 401       rnc (n0 >= 0.0 ? floor (n0) : -floor (-n0));
 402       break;
 403
 404       /* N-ary numeric functions. */
 405     case OP_ANY:
 406       if (n0 == SYSMIS)
 407         frnc (SYSMIS);
 408       else
 409         {
 410           int sysmis = 1;
 411           double ni;
 412
 413           for (i = 1; i < n->n; i++)
 414             {
 415               ni = n->arg[i]->num_con.value;
 416               if (approx_eq (n0, ni))
 417                 {
 418                   frnc (1.0);
 419                   goto any_done;
 420                 }
 421               if (ni != SYSMIS)
 422                 sysmis = 0;
 423             }
 424           frnc (sysmis ? SYSMIS : 0.0);
 425         }
 426     any_done:
 427       break;
 428     case OP_ANY_STRING:
 429       for (i = 1; i < n->n; i++)
 430         if (!st_compare_pad (n->arg[0]->str_con.s, n->arg[0]->str_con.len,
 431                              n->arg[i]->str_con.s, n->arg[i]->str_con.len))
 432           {
 433             frnc (1.0);
 434             goto any_string_done;
 435           }
 436       frnc (0.0);
 437     any_string_done:
 438       break;
 439
 440     case OP_CFVAR:
 441     case OP_MAX:
 442     case OP_MEAN:
 443     case OP_MIN:
 444     case OP_NMISS:
 445     case OP_NVALID:
 446     case OP_SD:
 447     case OP_SUM:
 448     case OP_VARIANCE:
 449       {
 450         double d[2] =
 451         {0.0, 0.0};             /* sum, sum of squares */
 452         double min = DBL_MAX;   /* minimum value */
 453         double max = -DBL_MAX;  /* maximum value */
 454         double ni;              /* value of i'th argument */
 455         int nv = 0;             /* number of valid arguments */
 456
 457         for (i = 0; i < n->n; i++)
 458           {
 459             ni = n->arg[i]->num_con.value;
 460             if (ni != SYSMIS)
 461               {
 462                 nv++;
 463                 d[0] += ni;
 464                 d[1] += ni * ni;
 465                 if (ni < min)
 466                   min = ni;
 467                 if (ni > max)
 468                   max = ni;
 469               }
 470           }
 471         if (n->type == OP_NMISS)
 472           frnc (i - nv);
 473         else if (n->type == OP_NVALID)
 474           frnc (nv);
 475         else if (nv >= (int) n->arg[i])
 476           {
 477             switch (n->type)
 478               {
 479               case OP_CFVAR:
 480                 frnc (calc_cfvar (d, nv));
 481                 break;
 482               case OP_MAX:
 483                 frnc (max);
 484                 break;
 485               case OP_MEAN:
 486                 frnc (calc_mean (d, nv));
 487                 break;
 488               case OP_MIN:
 489                 frnc (min);
 490                 break;
 491               case OP_SD:
 492                 frnc (calc_stddev (calc_variance (d, nv)));
 493                 break;
 494               case OP_SUM:
 495                 frnc (d[0]);
 496                 break;
 497               case OP_VARIANCE:
 498                 frnc (calc_variance (d, nv));
 499                 break;
 500               }
 501           }
 502         else
 503           frnc (SYSMIS);
 504       }
 505       break;
 506     case OP_RANGE:
 507       if (n0 == SYSMIS)
 508         frnc (SYSMIS);
 509       else
 510         {
 511           double min, max;
 512           int sysmis = 1;
 513
 514           for (i = 1; i < n->n; i += 2)
 515             {
 516               min = n->arg[i]->num_con.value;
 517               max = n->arg[i + 1]->num_con.value;
 518               if (min == SYSMIS || max == SYSMIS)
 519                 continue;
 520               sysmis = 0;
 521               if (approx_ge (n0, min) && approx_le (n0, max))
 522                 {
 523                   frnc (1.0);
 524                   goto range_done;
 525                 }
 526             }
 527           frnc (sysmis ? SYSMIS : 0.0);
 528         }
 529     range_done:
 530       break;
 531     case OP_RANGE_STRING:
 532       for (i = 1; i < n->n; i += 2)
 533         if (st_compare_pad (n->arg[0]->str_con.s, n->arg[0]->str_con.len,
 534                             n->arg[i]->str_con.s, n->arg[i]->str_con.len) >= 0
 535             && st_compare_pad (n->arg[0]->str_con.s, n->arg[0]->str_con.len,
 536                                n->arg[i + 1]->str_con.s,
 537                                n->arg[i + 1]->str_con.len) <= 0)
 538           {
 539             frnc (1.0);
 540             goto range_str_done;
 541           }
 542       frnc (0.0);
 543     range_str_done:
 544       break;
 545
 546       /* Time function. */
 547     case OP_TIME_HMS:
 548       rnc (60. * (60. * n0 + n1) + n2);
 549       break;
 550
 551       /* Date construction functions. */
 552     case OP_DATE_DMY:
 553       rnc (60. * 60. * 24. * yrmoda (n2, n1, n0));
 554       break;
 555     case OP_DATE_MDY:
 556       rnc (60. * 60. * 24. * yrmoda (n2, n0, n1));
 557       break;
 558     case OP_DATE_MOYR:
 559       rnc (60. * 60. * 24. * yrmoda (n1, n0, 1));
 560       break;
 561     case OP_DATE_QYR:
 562       rnc (60. * 60. * 24. * yrmoda (n1, 3 * (int) n0 - 2, 1));
 563       break;
 564     case OP_DATE_WKYR:
 565       {
 566         double t = yrmoda (n1, 1, 1);
 567         if (t != SYSMIS)
 568           t = 60. * 60. * 24. * (t + 7. * (n0 - 1));
 569         rnc (t);
 570       }
 571       break;
 572     case OP_DATE_YRDAY:
 573       {
 574         double t = yrmoda (n0, 1, 1);
 575         if (t != SYSMIS)
 576           t = 60. * 60. * 24. * (t + n0 - 1);
 577         rnc (t);
 578       }
 579       break;
 580     case OP_YRMODA:
 581       rnc (yrmoda (n0, n1, n2));
 582       break;
 583       /* Date extraction functions. */
 584     case OP_XDATE_DATE:
 585       rnc (floor (n0 / 60. / 60. / 24.) * 60. * 60. * 24.);
 586       break;
 587     case OP_XDATE_HOUR:
 588       rnc (fmod (floor (n0 / 60. / 60.), 24.));
 589       break;
 590     case OP_XDATE_JDAY:
 591       rnc (julian_to_jday (n0 / 86400.));
 592       break;
 593     case OP_XDATE_MDAY:
 594       {
 595         int day;
 596         julian_to_calendar (n0 / 86400., NULL, NULL, &day);
 597         rnc (day);
 598       }
 599       break;
 600     case OP_XDATE_MINUTE:
 601       rnc (fmod (floor (n0 / 60.), 60.));
 602       break;
 603     case OP_XDATE_MONTH:
 604       {
 605         int month;
 606         julian_to_calendar (n0 / 86400., NULL, &month, NULL);
 607         rnc (month);
 608       }
 609       break;
 610     case OP_XDATE_QUARTER:
 611       {
 612         int month;
 613         julian_to_calendar (n0 / 86400., NULL, &month, NULL);
 614         rnc ((month - 1) / 3 + 1);
 615       }
 616       break;
 617     case OP_XDATE_SECOND:
 618       rnc (fmod (n0, 60.));
 619       break;
 620     case OP_XDATE_TDAY:
 621       rnc (floor (n0 / 60. / 60. / 24.));
 622       break;
 623     case OP_XDATE_TIME:
 624       rnc (n0 - floor (n0 / 60. / 60. / 24.) * 60. * 60. * 24.);
 625       break;
 626     case OP_XDATE_WEEK:
 627       rnc ((julian_to_jday (n0) - 1) / 7 + 1);
 628       break;
 629     case OP_XDATE_WKDAY:
 630       rnc (julian_to_wday (n0));
 631       break;
 632     case OP_XDATE_YEAR:
 633       {
 634         int year;
 635         julian_to_calendar (n0 / 86400., &year, NULL, NULL);
 636         rnc (year);
 637       }
 638       break;
 639
 640       /* String functions. */
 641     case OP_CONCAT:
 642       {
 643         len = s0l;
 644         memcpy (strbuf, s0, len);
 645         for (i = 1; i < n->n; i++)
 646           {
 647             add = sl (i);
 648             if (add + len > 255)
 649               add = 255 - len;
 650             memcpy (&strbuf[len], s (i), add);
 651             len += add;
 652           }
 653         n = repl_str_con (n, strbuf, len);
 654       }
 655       break;
 656     case OP_INDEX:
 657       rnc (s1l ? str_search (s0, s0l, s1, s1l) : SYSMIS);
 658       break;
 659     case OP_INDEX_OPT:
 660       if (n2 == SYSMIS || (int) n2 <= 0 || s1l % (int) n2)
 661         {
 662           msg (SW, _("While optimizing a constant expression, there was "
 663                "a bad value for the third argument to INDEX."));
 664           frnc (SYSMIS);
 665         }
 666       else
 667         {
 668           int pos = 0;
 669           int c = s1l / (int) n2;
 670           int r;
 671
 672           for (i = 0; i < c; i++)
 673             {
 674               r = str_search (s0, s0l, s (i), sl (i));
 675               if (r < pos || pos == 0)
 676                 pos = r;
 677             }
 678           frnc (pos);
 679         }
 680       break;
 681     case OP_RINDEX:
 682       rnc (str_rsearch (s0, s0l, s1, s1l));
 683       break;
 684     case OP_RINDEX_OPT:
 685       if (n2 == SYSMIS || (int) n2 <= 0 || s1l % (int) n2)
 686         {
 687           msg (SE, _("While optimizing a constant expression, there was "
 688                "a bad value for the third argument to RINDEX."));
 689           frnc (SYSMIS);
 690         }
 691       else
 692         {
 693           int pos = 0;
 694           int c = s1l / n2;
 695           int r;
 696
 697           for (i = 0; i < c; i++)
 698             {
 699               r = str_rsearch (s0, s0l, s (i), sl (i));
 700               if (r > pos)
 701                 pos = r;
 702             }
 703           frnc (pos);
 704         }
 705       break;
 706     case OP_LENGTH:
 707       frnc (s0l);
 708       break;
 709     case OP_LOWER:
 710       {
 711         char *cp;
 712         for (cp = &s0[s0l]; cp >= s0; cp--)
 713           *cp = tolower ((unsigned char) (*cp));
 714         n = repl_str_con (n, s0, s0l);
 715       }
 716       break;
 717     case OP_UPPER:
 718       {
 719         char *cp;
 720         for (cp = &s0[s0[0] + 1]; cp > s0; cp--)
 721           *cp = toupper ((unsigned char) (*cp));
 722         n = repl_str_con (n, s0, s0l);
 723       }
 724       break;
 725     case OP_LPAD:
 726     case OP_LPAD_OPT:
 727     case OP_RPAD:
 728     case OP_RPAD_OPT:
 729       {
 730         int c;
 731
 732         if (n1 == SYSMIS)
 733           {
 734             n = repl_str_con (n, NULL, 0);
 735             break;
 736           }
 737         len = n1;
 738         len = range (len, 1, 255);
 739         add = max (n1 - s0l, 0);
 740
 741         if (n->type == OP_LPAD_OPT || n->type == OP_RPAD_OPT)
 742           {
 743             if (s2l < 1)
 744               {
 745                 c = n->type == OP_LPAD_OPT ? 'L' : 'R';
 746                 msg (SE, _("Third argument to %cPAD() must be at least one "
 747                      "character in length."), c);
 748                 c = ' ';
 749               }
 750             else
 751               c = s2[0];
 752           }
 753         else
 754           c = ' ';
 755
 756         if (n->type == OP_LPAD || n->type == OP_LPAD_OPT)
 757           memmove (&s0[add], s0, len);
 758         if (n->type == OP_LPAD || n->type == OP_LPAD_OPT)
 759           memset (s0, c, add);
 760         else
 761           memset (&s0[s0l], c, add);
 762
 763         n = repl_str_con (n, s0, len);
 764       }
 765       break;
 766     case OP_LTRIM:
 767     case OP_LTRIM_OPT:
 768     case OP_RTRIM:
 769     case OP_RTRIM_OPT:
 770       {
 771         int c;
 772         char *cp = s0;
 773
 774         if (n->type == OP_LTRIM_OPT || n->type == OP_RTRIM_OPT)
 775           {
 776             if (s1l < 1)
 777               {
 778                 c = n->type == OP_LTRIM_OPT ? 'L' : 'R';
 779                 msg (SE, _("Second argument to %cTRIM() must be at least one "
 780                      "character in length."), c);
 781                 c = ' ';
 782               }
 783             else
 784               c = s1[0];
 785           }
 786         len = s0l;
 787         if (n->type == OP_LTRIM || n->type == OP_LTRIM_OPT)
 788           {
 789             while (*cp == c && cp < &s0[len])
 790               cp++;
 791             len -= cp - s0;
 792           }
 793         else
 794           while (len > 0 && s0[len - 1] == c)
 795             len--;
 796         n = repl_str_con (n, cp, len);
 797       }
 798       break;
 799     case OP_NUMBER:
 800     case OP_NUMBER_OPT:
 801       {
 802         union value v;
 803         struct data_in di;
 804
 805         di.s = s0;
 806         di.e = s0 + s0l;
 807         di.v = &v;
 808         di.flags = DI_IGNORE_ERROR;
 809         di.f1 = 1;
 810
 811         if (n->type == OP_NUMBER_OPT)
 812           {
 813             di.format.type = (int) n->arg[1];
 814             di.format.w = (int) n->arg[2];
 815             di.format.d = (int) n->arg[3];
 816           }
 817         else
 818           {
 819             di.format.type = FMT_F;
 820             di.format.w = s0l;
 821             di.format.d = 0;
 822           }
 823
 824         data_in (&di);
 825         frnc (v.f);
 826       }
 827       break;
 828     case OP_STRING:
 829       {
 830         union value v;
 831         struct fmt_spec f;
 832         f.type = (int) n->arg[1];
 833         f.w = (int) n->arg[2];
 834         f.d = (int) n->arg[3];
 835         v.f = n0;
 836
 837         data_out (strbuf, &f, &v);
 838         n = repl_str_con (n, strbuf, f.w);
 839       }
 840       break;
 841     case OP_SUBSTR:
 842     case OP_SUBSTR_OPT:
 843       {
 844         int pos = (int) n1;
 845         if (pos > s0l || pos <= 0 || n1 == SYSMIS
 846             || (n->type == OP_SUBSTR_OPT && n2 == SYSMIS))
 847           n = repl_str_con (n, NULL, 0);
 848         else
 849           {
 850             if (n->type == OP_SUBSTR_OPT)
 851               {
 852                 len = (int) n2;
 853                 if (len + pos - 1 > s0l)
 854                   len = s0l - pos + 1;
 855               }
 856             else
 857               len = s0l - pos + 1;
 858             n = repl_str_con (n, &s0[pos - 1], len);
 859           }
 860       }
 861       break;
 862
 863       /* Weirdness. */
 864     case OP_INV:
 865       rnc (1.0 / n0);
 866       break;
 867     case OP_MOD:
 868       if (approx_eq (n0, 0.0) && n1 == SYSMIS)
 869         frnc (0.0);
 870       else
 871         rnc (fmod (n0, n1));
 872       break;
 873     case OP_NUM_TO_BOOL:
 874       if (approx_eq (n0, 0.0))
 875         n0 = 0.0;
 876       else if (approx_eq (n0, 1.0))
 877         n0 = 1.0;
 878       else if (n0 != SYSMIS)
 879         {
 880           msg (SE, _("When optimizing a constant expression, an integer "
 881                "that was being used as an Boolean value was found "
 882                "to have a constant value other than 0, 1, or SYSMIS."));
 883           n0 = 0.0;
 884         }
 885       rnc (n0);
 886       break;
 887
 888 #if __CHECKER__
 889       /* This case prevents Checker from choking. */
 890     case 42000:
 891       assert (0);
 892 #endif
 893     }
 894   return n;
 895 }
 896
 897 #undef n0
 898 #undef n1
 899 #undef n2
 900
 901 #undef s0
 902 #undef s0l
 903 #undef s1
 904 #undef s1l
 905 #undef s2
 906 #undef s2l
 907 #undef s
 908 #undef sl
 909
 910 #undef rnc
 911 #undef frnc
 912
 913 static struct nonterm_node *
 914 repl_num_con (struct nonterm_node * n, double d)
 915 {
 916   int i;
 917   if (!finite (d) || errno)
 918     d = SYSMIS;
 919   else
 920     for (i = 0; i < n->n; i++)
 921       if (n->arg[i]->type == OP_NUM_CON && n->arg[i]->num_con.value == SYSMIS)
 922         {
 923           d = SYSMIS;
 924           break;
 925         }
 926   return force_repl_num_con (n, d);
 927 }
 928
 929 static struct nonterm_node *
 930 force_repl_num_con (struct nonterm_node * n, double d)
 931 {
 932   struct num_con_node *num;
 933
 934   if (!finite (d) || errno)
 935     d = SYSMIS;
 936   free_node ((union any_node *) n);
 937   num = xmalloc (sizeof *num);
 938   num->type = OP_NUM_CON;
 939   num->value = d;
 940   return (struct nonterm_node *) num;
 941 }
 942
 943 static struct nonterm_node *
 944 repl_str_con (struct nonterm_node * n, char *s, int len)
 945 {
 946   struct str_con_node *str;
 947
 948   /* The ordering here is important since the source string may be
 949      part of a subnode of n. */
 950   str = xmalloc (sizeof *str + len - 1);
 951   str->type = OP_STR_CON;
 952   str->len = len;
 953   memcpy (str->s, s, len);
 954   free_node ((union any_node *) n);
 955   return (struct nonterm_node *) str;
 956 }
 957
 958 /* Returns the number of days since 10 Oct 1582 for the date
 959    YEAR/MONTH/DAY, where YEAR is in range 0..199 or 1582..19999, MONTH
 960    is in 1..12, and DAY is in 1..31. */
 961 double
 962 yrmoda (double year, double month, double day)
 963 {
 964   if (year == SYSMIS || month == SYSMIS || day == SYSMIS)
 965     return SYSMIS;
 966
 967   /* The addition of EPSILON avoids converting, for example,
 968      1991.9999997=>1991. */
 969   year = floor (year + EPSILON);
 970   month = floor (month + EPSILON);
 971   day = floor (day + EPSILON);
 972
 973   if (year >= 0. && year <= 199.)
 974     year += 1900.;
 975   if ((year < 1582. || year > 19999.)
 976       || (year == 1582. && (month < 10. || (month == 10. && day < 15.)))
 977       || (month < -1 || month > 13)
 978       || (day < -1 || day > 32))
 979     return SYSMIS;
 980   return calendar_to_julian (year, month, day);
 981 }
 982 \f
 983 /* Expression dumper. */
 984
 985 static struct expression *e;
 986 static int nop, mop;
 987 static int ndbl, mdbl;
 988 static int nstr, mstr;
 989 static int nvars, mvars;
 990
 991 static void dump_node (union any_node * n);
 992 static void emit (int);
 993 static void emit_num_con (double);
 994 static void emit_str_con (char *, int);
 995 static void emit_var (struct variable *);
 996
 997 void
 998 dump_expression (union any_node * n, struct expression * expr)
 999 {
1000   unsigned char *o;
1001
1002   int height = 0;
1003
1004   int max_height = 0;
1005
1006   e = expr;
1007   e->op = NULL;
1008   e->num = NULL;
1009   e->str = NULL;
1010   e->var = NULL;
1011   nop = mop = 0;
1012   ndbl = mdbl = 0;
1013   nstr = mstr = 0;
1014   nvars = mvars = 0;
1015   dump_node (n);
1016   emit (OP_SENTINEL);
1017
1018   /* Now compute the stack height needed to evaluate the expression. */
1019   for (o = e->op; *o != OP_SENTINEL; o++)
1020     {
1021       if (ops[*o].flags & OP_VAR_ARGS)
1022         height += 1 - o[1];
1023       else
1024         height += ops[*o].height;
1025       o += ops[*o].skip;
1026       if (height > max_height)
1027         max_height = height;
1028     }
1029
1030   /* ANSI says we have to waste space for one `value' since pointers
1031      are not guaranteed to be able to point to a spot *before* a
1032      block.  If only all the world were a VAX... */
1033   max_height++;
1034
1035   e->stack = xmalloc (max_height * sizeof *e->stack);
1036
1037 #if PAGED_STACK
1038   e->str_stack = e->type == EX_STRING ? xmalloc (256) : NULL;
1039 #else
1040   e->str_stack = xmalloc (256);
1041   e->str_size = 256;
1042 #endif
1043 }
1044
1045 static void
1046 dump_node (union any_node * n)
1047 {
1048   if (n->type == OP_AND || n->type == OP_OR)
1049     {
1050       int i;
1051
1052       dump_node (n->nonterm.arg[0]);
1053       for (i = 1; i < n->nonterm.n; i++)
1054         {
1055           dump_node (n->nonterm.arg[i]);
1056           emit (n->type);
1057         }
1058       return;
1059     }
1060   else if (n->type < OP_TERMINAL)
1061     {
1062       int i;
1063       for (i = 0; i < n->nonterm.n; i++)
1064         dump_node (n->nonterm.arg[i]);
1065       emit (n->type);
1066       if (ops[n->type].flags & OP_VAR_ARGS)
1067         emit (n->nonterm.n);
1068       if (ops[n->type].flags & OP_MIN_ARGS)
1069         emit ((int) n->nonterm.arg[n->nonterm.n]);
1070       if (ops[n->type].flags & OP_FMT_SPEC)
1071         {
1072           emit ((int) n->nonterm.arg[n->nonterm.n]);
1073           emit ((int) n->nonterm.arg[n->nonterm.n + 1]);
1074           emit ((int) n->nonterm.arg[n->nonterm.n + 2]);
1075         }
1076       return;
1077     }
1078
1079   emit (n->type);
1080   if (n->type == OP_NUM_CON)
1081     emit_num_con (n->num_con.value);
1082   else if (n->type == OP_STR_CON)
1083     emit_str_con (n->str_con.s, n->str_con.len);
1084   else if (n->type == OP_NUM_VAR || n->type == OP_STR_VAR
1085            || n->type == OP_STR_MIS)
1086     emit_var (n->var.v);
1087   else if (n->type == OP_NUM_LAG || n->type == OP_STR_LAG)
1088     {
1089       emit_var (n->lag.v);
1090       emit (n->lag.lag);
1091     }
1092   else if (n->type == OP_NUM_SYS || n->type == OP_NUM_VAL)
1093     emit (n->var.v->fv);
1094   else
1095     assert (n->type == OP_CASENUM);
1096 }
1097
1098 static void
1099 emit (int op)
1100 {
1101   if (nop >= mop)
1102     {
1103       mop += 16;
1104       e->op = xrealloc (e->op, mop * sizeof *e->op);
1105     }
1106   e->op[nop++] = op;
1107 }
1108
1109 static void
1110 emit_num_con (double dbl)
1111 {
1112   if (ndbl >= mdbl)
1113     {
1114       mdbl += 16;
1115       e->num = xrealloc (e->num, mdbl * sizeof *e->num);
1116     }
1117   e->num[ndbl++] = dbl;
1118 }
1119
1120 static void
1121 emit_str_con (char *str, int len)
1122 {
1123   if (nstr + len + 1 > mstr)
1124     {
1125       mstr += 256;
1126       e->str = xrealloc (e->str, mstr);
1127     }
1128   e->str[nstr++] = len;
1129   memcpy (&e->str[nstr], str, len);
1130   nstr += len;
1131 }
1132
1133 static void
1134 emit_var (struct variable * v)
1135 {
1136   if (nvars >= mvars)
1137     {
1138       mvars += 16;
1139       e->var = xrealloc (e->var, mvars * sizeof *e->var);
1140     }
1141   e->var[nvars++] = v;
1142 }