1 /* PSPP - computes sample statistics.
2 Copyright (C) 1997-9, 2000 Free Software Foundation, Inc.
3 Written by Ben Pfaff <blp@gnu.org>.
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.
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.
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
22 #if TIME_WITH_SYS_TIME
42 #include "dictionary.h"
44 #include "julcal/julcal.h"
56 expr_evaluate (const struct expression *e, const struct ccase *c, int case_idx,
59 unsigned char *op = e->op;
61 unsigned char *str = e->str;
62 struct variable **vars = e->var;
66 union value *sp = e->stack;
76 if (sp[1].f == SYSMIS)
78 else if (sp[0].f != SYSMIS)
83 if (sp[1].f == SYSMIS)
85 else if (sp[0].f != SYSMIS)
90 if (sp[1].f == SYSMIS)
92 else if (sp[0].f != SYSMIS)
97 if (sp[1].f == SYSMIS || sp[1].f == 0.)
99 else if (sp[0].f != SYSMIS)
104 if (sp[0].f == SYSMIS)
109 else if (sp[1].f == SYSMIS)
116 else if (sp[0].f == 0.0 && sp[1].f <= 0.0)
119 sp->f = pow (sp[0].f, sp[1].f);
123 /* Note that booleans are always one of 0, 1, or SYSMIS.
125 Truth table (in order of detection):
137 1 and SYSMIS = SYSMIS
138 SYSMIS and SYSMIS = SYSMIS
142 SYSMIS and 1 = SYSMIS
146 if (sp[0].f == 0.0); /* 1 */
147 else if (sp[1].f == 0.0)
149 else if (sp[1].f == SYSMIS)
150 sp->f = SYSMIS; /* 3 */
153 /* Truth table (in order of detection):
166 SYSMIS or SYSMIS = SYSMIS
174 if (sp[0].f == 1.0); /* 1 */
175 else if (sp[1].f == 1.0)
177 else if (sp[1].f == SYSMIS)
178 sp->f = SYSMIS; /* 3 */
183 else if (sp[0].f == 1.0)
189 if (sp[0].f != SYSMIS)
191 if (sp[1].f == SYSMIS)
194 sp->f = sp[0].f == sp[1].f;
199 if (sp[0].f != SYSMIS)
201 if (sp[1].f == SYSMIS)
204 sp->f = sp[0].f >= sp[1].f;
209 if (sp[0].f != SYSMIS)
211 if (sp[1].f == SYSMIS)
214 sp->f = sp[0].f > sp[1].f;
219 if (sp[0].f != SYSMIS)
221 if (sp[1].f == SYSMIS)
224 sp->f = sp[0].f <= sp[1].f;
229 if (sp[0].f != SYSMIS)
231 if (sp[1].f == SYSMIS)
234 sp->f = sp[0].f < sp[1].f;
239 if (sp[0].f != SYSMIS)
241 if (sp[1].f == SYSMIS)
244 sp->f = sp[0].f != sp[1].f;
248 /* String operators. */
251 sp[0].f = st_compare_pad (&sp[0].c[1], sp[0].c[0],
252 &sp[1].c[1], sp[1].c[0]) == 0;
256 sp[0].f = st_compare_pad (&sp[0].c[1], sp[0].c[0],
257 &sp[1].c[1], sp[1].c[0]) >= 0;
261 sp[0].f = st_compare_pad (&sp[0].c[1], sp[0].c[0],
262 &sp[1].c[1], sp[1].c[0]) > 0;
266 sp[0].f = st_compare_pad (&sp[0].c[1], sp[0].c[0],
267 &sp[1].c[1], sp[1].c[0]) <= 0;
271 sp[0].f = st_compare_pad (&sp[0].c[1], sp[0].c[0],
272 &sp[1].c[1], sp[1].c[0]) < 0;
276 sp[0].f = st_compare_pad (&sp[0].c[1], sp[0].c[0],
277 &sp[1].c[1], sp[1].c[0]) != 0;
280 /* Unary functions. */
287 sp->f = fabs (sp->f);
293 sp->f = acos (sp->f);
302 sp->f = asin (sp->f);
309 sp->f = atan (sp->f);
328 sp->f = log10 (sp->f);
344 sp->f = fmod (sp->f, 10);
350 sp->f = floor (sp->f + 0.5);
352 sp->f = -floor (-sp->f + 0.5);
363 sp->f = sqrt (sp->f);
381 sp->f = floor (sp->f);
383 sp->f = -floor (-sp->f);
387 /* N-ary numeric functions. */
396 for (i = 1; i < n_args; i++)
397 if (sp[0].f == sp[i].f)
402 else if (sp[i].f != SYSMIS)
404 sp->f = sysmis ? SYSMIS : 0.0;
413 for (i = 1; i < n_args; i++)
414 if (!st_compare_pad (&sp[0].c[1], sp[0].c[0],
415 &sp[i].c[1], sp[i].c[0]))
426 double weight, mean, variance;
430 moments_of_values (sp, n_args,
431 &weight, &mean, &variance, NULL, NULL);
433 if (weight < *op++ || mean == SYSMIS
434 || mean == 0 || variance == SYSMIS)
437 sp->f = sqrt (variance) / mean;
444 double max = -DBL_MAX;
447 for (i = 0; i < n_args; i++)
448 if (sp[i].f != SYSMIS)
466 for (i = 1; i < n_args; i++)
467 if (st_compare_pad (&sp[i].c[1], sp[i].c[0],
468 &sp[max_idx].c[1], sp[max_idx].c[0]) > 0)
470 sp[0].c = sp[max_idx].c;
480 moments_of_values (sp, n_args,
481 &weight, &mean, NULL, NULL, NULL);
482 sp->f = weight < *op++ ? SYSMIS : mean;
489 double min = DBL_MAX;
492 for (i = 0; i < n_args; i++)
493 if (sp[i].f != SYSMIS)
511 for (i = 1; i < n_args; i++)
512 if (st_compare_pad (&sp[i].c[1], sp[i].c[0],
513 &sp[min_idx].c[1], sp[min_idx].c[0]) < 0)
515 sp[0].c = sp[min_idx].c;
524 for (i = 0; i < n_args; i++)
525 if (sp[i].f == SYSMIS)
536 for (i = 0; i < n_args; i++)
537 if (sp[i].f != SYSMIS)
550 for (i = 1; i < n_args; i += 2)
551 if (sp[i].f == SYSMIS || sp[i + 1].f == SYSMIS)
553 else if (sp[0].f >= sp[i].f && sp[0].f <= sp[i + 1].f)
560 sp->f = sysmis ? SYSMIS : 0.0;
563 case OP_RANGE_STRING:
568 for (i = 1; i < n_args; i += 2)
569 if (st_compare_pad (&sp[0].c[1], sp[0].c[0],
570 &sp[i].c[1], sp[i].c[0]) >= 0
571 && st_compare_pad (&sp[0].c[1], sp[0].c[0],
572 &sp[i + 1].c[1], sp[i + 1].c[0]) <= 0)
583 double weight, variance;
586 moments_of_values (sp, n_args,
587 &weight, NULL, &variance, NULL, NULL);
588 sp->f = weight < *op++ ? SYSMIS : sqrt (variance);
598 for (i = 0; i < n_args; i++)
599 if (sp[i].f != SYSMIS)
613 double weight, variance;
616 moments_of_values (sp, n_args,
617 &weight, NULL, &variance, NULL, NULL);
618 sp->f = weight < *op++ ? SYSMIS : variance;
622 /* Time construction function. */
625 if (sp[0].f == SYSMIS || sp[1].f == SYSMIS || sp[2].f == SYSMIS)
630 min = min (sp[0].f, min (sp[1].f, sp[2].f));
631 max = max (sp[0].f, max (sp[1].f, sp[2].f));
632 if (min < 0. && max > 0.)
634 msg (SW, _("TIME.HMS cannot mix positive and negative "
635 "in its arguments."));
639 sp->f = 60. * (60. * sp[0].f + sp[1].f) + sp[2].f;
644 sp->f /= 60. * 60. * 24.;
650 case OP_CTIME_MINUTES:
656 sp->f *= 60. * 60. * 24.;
658 case OP_CTIME_SECONDS:
662 /* Date construction functions. */
665 sp->f = yrmoda (sp[2].f, sp[1].f, sp[0].f);
667 sp->f *= 60. * 60. * 24.;
671 sp->f = yrmoda (sp[2].f, sp[0].f, sp[1].f);
673 sp->f *= 60. * 60. * 24.;
677 sp->f = yrmoda (sp[1].f, sp[0].f, 1);
679 sp->f *= 60. * 60. * 24.;
683 if (sp[0].f == SYSMIS)
687 sp->f = yrmoda (sp[1].f, sp[0].f * 3 - 2, 1);
689 sp->f *= 60. * 60. * 24.;
694 if (sp[0].f == SYSMIS || sp[1].f == SYSMIS)
696 else if (sp[0].f < 0. || sp[0].f > 53.)
698 msg (SW, _("Week argument to WKYR must be in range 0 to 53."));
703 double result = yrmoda (sp[1].f, 1, 1);
704 if (result != SYSMIS)
705 result = 60. * 60. * 24. * (result + 7. * (sp->f - 1.));
711 if (sp[1].f == SYSMIS)
715 sp->f = yrmoda (sp[0].f, 1, 1);
717 sp->f = 60. * 60. * 24. * (sp->f + floor (sp[1].f) - 1);
722 sp->f = yrmoda (sp[0].f, sp[1].f, sp[2].f);
725 /* Date extraction functions. */
728 sp->f = floor (sp->f / 60. / 60. / 24.) * 60. * 60. * 24.;
732 sp->f = fmod (floor (sp->f / 60. / 60.), 24.);
736 sp->f = 86400. * julian_to_jday (sp->f / 86400.);
742 julian_to_calendar (sp->f / 86400., NULL, NULL, &day);
746 case OP_XDATE_MINUTE:
748 sp->f = fmod (floor (sp->f / 60.), 60.);
754 julian_to_calendar (sp->f / 86400., NULL, &month, NULL);
758 case OP_XDATE_QUARTER:
762 julian_to_calendar (sp->f / 86400., NULL, &month, NULL);
763 sp->f = (month - 1) / 3 + 1;
766 case OP_XDATE_SECOND:
768 sp->f = fmod (sp->f, 60.);
772 sp->f = floor (sp->f / 60. / 60. / 24.);
776 sp->f -= floor (sp->f / 60. / 60. / 24.) * 60. * 60. * 24.;
780 sp->f = (julian_to_jday (sp->f / 86400.) - 1) / 7 + 1;
784 sp->f = julian_to_wday (sp->f / 86400.);
790 julian_to_calendar (sp->f / 86400., &year, NULL, NULL);
795 /* String functions. */
801 dest = pool_alloc (e->pool, 256);
805 for (i = 0; i < n_args; i++)
808 if (sp[i].c[0] + dest[0] < 255)
810 memcpy (&dest[dest[0] + 1], &sp[i].c[1], sp[i].c[0]);
811 dest[0] += sp[i].c[0];
815 memcpy (&dest[dest[0] + 1], &sp[i].c[1], 255 - dest[0]);
829 int last = sp[0].c[0] - sp[1].c[0];
831 for (i = 0; i <= last; i++)
832 if (!memcmp (&sp[0].c[i + 1], &sp[1].c[1], sp[1].c[0]))
847 else if (sp[2].f == SYSMIS)
849 msg (SW, _("Argument 3 of RINDEX may not be system-missing."));
854 int part_len = sp[2].f;
856 if (part_len < 0 || part_len > sp[1].c[0]
857 || sp[1].c[0] % part_len != 0)
859 msg (SW, _("Argument 3 of RINDEX must be between 1 and "
860 "the length of argument 2, and it must "
861 "evenly divide the length of argument 2."));
867 int last = sp[0].c[0] - part_len;
868 for (i = 0; i <= last; i++)
869 for (j = 0; j < sp[1].c[0]; j += part_len)
870 if (!memcmp (&sp[0].c[i + 1], &sp[1].c[j + 1], part_len))
887 for (i = sp[0].c[0] - sp[1].c[0]; i >= 0; i--)
888 if (!memcmp (&sp[0].c[i + 1], &sp[1].c[1], sp[1].c[0]))
903 else if (sp[2].f == SYSMIS)
905 msg (SW, _("Argument 3 of RINDEX may not be system-missing."));
910 int part_len = sp[2].f;
912 if (part_len < 0 || part_len > sp[1].c[0]
913 || sp[1].c[0] % part_len != 0)
915 msg (SW, _("Argument 3 of RINDEX must be between 1 and "
916 "the length of argument 2, and it must "
917 "evenly divide the length of argument 2."));
922 for (i = sp[0].c[0] - part_len; i >= 0; i--)
923 for (j = 0; j < sp[1].c[0]; j += part_len)
924 if (!memcmp (&sp[0].c[i + 1], &sp[1].c[j + 1], part_len))
938 for (i = sp[0].c[0]; i >= 1; i--)
939 sp[0].c[i] = tolower ((unsigned char) (sp[0].c[i]));
942 for (i = sp[0].c[0]; i >= 1; i--)
943 sp[0].c[i] = toupper ((unsigned char) (sp[0].c[i]));
950 if (sp[1].f == SYSMIS || len < 0 || len > 255 || sp[2].c[0] != 1)
952 else if (len > sp[0].c[0])
956 dest = pool_alloc (e->pool, len + 1);
958 memset (&dest[1], sp[2].c[1], len - sp->c[0]);
959 memcpy (&dest[len - sp->c[0] + 1], &sp->c[1], sp->c[0]);
969 if (len < 0 || len > 255 || sp[2].c[0] != 1)
971 else if (len > sp[0].c[0])
975 dest = pool_alloc (e->pool, len + 1);
977 memcpy (&dest[1], &sp->c[1], sp->c[0]);
978 memset (&dest[sp->c[0] + 1], sp[2].c[1], len - sp->c[0]);
990 int len = sp[0].c[0];
991 int cmp = sp[1].c[1];
994 while (i <= len && sp[0].c[i] == cmp)
998 sp[0].c[i] = sp[0].c[0] - i;
1006 if (sp[1].c[0] != 1)
1010 int cmp = sp[1].c[1];
1011 while (sp[0].c[0] > 0 && sp[0].c[sp[0].c[0]] == cmp)
1023 di.format.type = *op++;
1024 di.format.w = *op++;
1025 di.format.d = *op++;
1026 di.e = &sp->c[1] + min (sp->c[0], di.format.w);
1034 unsigned char *dest;
1040 dest = pool_alloc (e->pool, f.w + 1);
1043 assert ((formats[f.type].cat & FCAT_STRING) == 0);
1044 data_out (&dest[1], &f, sp);
1054 if (index < 1 || index > sp[0].c[0])
1059 sp->c[index] = sp->c[0] - index;
1072 if (sp[1].f == SYSMIS || sp[2].f == SYSMIS || index < 1
1073 || index > sp[0].c[0] || n < 1)
1080 sp->c[index] = sp->c[0] - index;
1091 if (sp->f != SYSMIS)
1094 case OP_NUM_TO_BOOL:
1097 else if (sp->f == 1.0)
1099 else if (sp->f != SYSMIS)
1101 msg (SE, _("A number being treated as a Boolean in an "
1102 "expression was found to have a value other than "
1103 "0 (false), 1 (true), or the system-missing value. "
1104 "The result was forced to 0."));
1112 if (sp[0].f != SYSMIS)
1114 if (sp[1].f == SYSMIS)
1120 sp->f = fmod (sp[0].f, sp[1].f);
1124 if (sp->f != SYSMIS)
1125 sp->f *= rng_get_double_normal (pspp_rng ());
1128 if (sp->f != SYSMIS)
1129 sp->f *= rng_get_double (pspp_rng ());
1132 sp->f = sp->f == SYSMIS || !finite (sp->f);
1134 case OP_VEC_ELEM_NUM:
1136 int rindx = sp[0].f + EPSILON;
1137 const struct vector *v = dict_get_vector (default_dict, *op++);
1139 if (sp[0].f == SYSMIS || rindx < 1 || rindx > v->cnt)
1141 if (sp[0].f == SYSMIS)
1142 msg (SE, _("SYSMIS is not a valid index value for vector "
1143 "%s. The result will be set to SYSMIS."),
1146 msg (SE, _("%g is not a valid index value for vector %s. "
1147 "The result will be set to SYSMIS."),
1153 sp->f = case_num (c, v->var[rindx - 1]->fv);
1156 case OP_VEC_ELEM_STR:
1158 int rindx = sp[0].f + EPSILON;
1159 const struct vector *vect = dict_get_vector (default_dict, *op++);
1162 if (sp[0].f == SYSMIS || rindx < 1 || rindx > vect->cnt)
1164 if (sp[0].f == SYSMIS)
1165 msg (SE, _("SYSMIS is not a valid index value for vector "
1166 "%s. The result will be set to the empty "
1170 msg (SE, _("%g is not a valid index value for vector %s. "
1171 "The result will be set to the empty string."),
1172 sp[0].f, vect->name);
1173 sp->c = pool_alloc (e->pool, 1);
1178 v = vect->var[rindx - 1];
1179 sp->c = pool_alloc (e->pool, v->width + 1);
1180 sp->c[0] = v->width;
1182 memcpy (&sp->c[1], case_str (c, v->fv), v->width);
1193 sp->c = pool_alloc (e->pool, *str + 1);
1194 memcpy (sp->c, str, *str + 1);
1200 sp->f = case_num (c, (*vars)->fv);
1201 if (is_num_user_missing (sp->f, *vars))
1207 int width = (*vars)->width;
1210 sp->c = pool_alloc (e->pool, width + 1);
1213 memcpy (&sp->c[1], case_str (c, (*vars)->fv), width);
1219 struct ccase *c = lagged_case (*op++);
1226 sp->f = case_num (c, (*vars)->fv);
1227 if (is_num_user_missing (sp->f, *vars))
1235 struct ccase *c = lagged_case (*op++);
1236 int width = (*vars)->width;
1239 sp->c = pool_alloc (e->pool, width + 1);
1243 memset (sp->c, ' ', width);
1245 memcpy (&sp->c[1], case_str (c, (*vars)->fv), width);
1253 sp->f = case_num (c, *op++) == SYSMIS;
1258 sp->f = case_num (c, *op++);
1275 if (e->type != EXPR_STRING)
1277 double value = sp->f;
1278 if (!finite (value))