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