Fix assertion for proper Huffman merge pattern: 0 == 1 modulo 1.
[pspp] / 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 }