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
43 #include "julcal/julcal.h"
54 expr_evaluate (struct expression *e, struct ccase *c, union value *v)
56 unsigned char *op = e->op;
58 unsigned char *str = e->str;
59 struct variable **vars = e->var;
63 union value *sp = e->stack;
74 for (i = 1; i < *op; i++)
76 if (sp[i].f == SYSMIS)
89 for (i = 1; i < *op; i++)
91 if (sp[i].f == SYSMIS)
103 if (sp[0].f == SYSMIS)
105 if (approx_eq (sp[1].f, 0.0))
108 else if (sp[1].f == SYSMIS)
116 else if (approx_eq (sp[0].f, 0.0) && approx_eq (sp[1].f, 0.0))
119 sp->f = pow (sp[0].f, sp[1].f);
123 /* Note that the equality operator (==) may be used here
124 (instead of approx_eq) because booleans are always
125 *exactly* 0, 1, or SYSMIS.
127 Truth table (in order of detection):
139 1 and SYSMIS = SYSMIS
140 SYSMIS and SYSMIS = SYSMIS
144 SYSMIS and 1 = SYSMIS
148 if (sp[0].f == 0.0); /* 1 */
149 else if (sp[1].f == 0.0)
151 else if (sp[1].f == SYSMIS)
152 sp->f = SYSMIS; /* 3 */
155 /* Truth table (in order of detection):
168 SYSMIS or SYSMIS = SYSMIS
176 if (sp[0].f == 1.0); /* 1 */
177 else if (sp[1].f == 1.0)
179 else if (sp[1].f == SYSMIS)
180 sp->f = SYSMIS; /* 3 */
185 else if (sp[0].f == 1.0)
191 if (sp[0].f != SYSMIS)
193 if (sp[1].f == SYSMIS)
196 sp->f = approx_eq (sp[0].f, sp[1].f);
201 if (sp[0].f != SYSMIS)
203 if (sp[1].f == SYSMIS)
206 sp->f = approx_ge (sp[0].f, sp[1].f);
211 if (sp[0].f != SYSMIS)
213 if (sp[1].f == SYSMIS)
216 sp->f = approx_gt (sp[0].f, sp[1].f);
221 if (sp[0].f != SYSMIS)
223 if (sp[1].f == SYSMIS)
226 sp->f = approx_le (sp[0].f, sp[1].f);
231 if (sp[0].f != SYSMIS)
233 if (sp[1].f == SYSMIS)
236 sp->f = approx_lt (sp[0].f, sp[1].f);
241 if (sp[0].f != SYSMIS)
243 if (sp[1].f == SYSMIS)
246 sp->f = approx_ne (sp[0].f, sp[1].f);
250 /* String operators. */
253 sp[0].f = st_compare_pad (&sp[0].c[1], sp[0].c[0],
254 &sp[1].c[1], sp[1].c[0]) == 0;
258 sp[0].f = st_compare_pad (&sp[0].c[1], sp[0].c[0],
259 &sp[1].c[1], sp[1].c[0]) >= 0;
263 sp[0].f = st_compare_pad (&sp[0].c[1], sp[0].c[0],
264 &sp[1].c[1], sp[1].c[0]) > 0;
268 sp[0].f = st_compare_pad (&sp[0].c[1], sp[0].c[0],
269 &sp[1].c[1], sp[1].c[0]) <= 0;
273 sp[0].f = st_compare_pad (&sp[0].c[1], sp[0].c[0],
274 &sp[1].c[1], sp[1].c[0]) < 0;
278 sp[0].f = st_compare_pad (&sp[0].c[1], sp[0].c[0],
279 &sp[1].c[1], sp[1].c[0]) != 0;
282 /* Unary functions. */
289 sp->f = fabs (sp->f);
295 sp->f = acos (sp->f);
304 sp->f = asin (sp->f);
311 sp->f = atan (sp->f);
330 sp->f = log10 (sp->f);
339 sp->f = log10 (sp->f);
346 sp->f = fmod (sp->f, 10);
352 sp->f = floor (sp->f + 0.5);
354 sp->f = -floor (-sp->f + 0.5);
365 sp->f = sqrt (sp->f);
383 sp->f = floor (sp->f);
385 sp->f = -floor (-sp->f);
389 /* N-ary numeric functions. */
398 for (i = 1; i <= n_args; i++)
399 if (approx_eq (sp[0].f, sp[i].f))
404 else if (sp[i].f != SYSMIS)
406 sp->f = sysmis ? SYSMIS : 0.0;
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]))
432 for (i = 0; i < n_args; i++)
433 if (sp[i].f != SYSMIS)
437 sum[1] += sp[i].f * sp[i].f;
442 sp->f = calc_cfvar (sum, nv);
449 double max = -DBL_MAX;
452 for (i = 0; i < n_args; i++)
453 if (sp[i].f != SYSMIS)
473 for (i = 0; i < n_args; i++)
474 if (sp[i].f != SYSMIS)
482 sp->f = calc_mean (sum, nv);
489 double min = DBL_MAX;
492 for (i = 0; i < n_args; i++)
493 if (sp[i].f != SYSMIS)
511 for (i = 0; i < n_args; i++)
512 if (sp[i].f == SYSMIS)
523 for (i = 0; i < n_args; i++)
524 if (sp[i].f != SYSMIS)
537 for (i = 1; i <= n_args; i += 2)
538 if (sp[i].f == SYSMIS || sp[i + 1].f == SYSMIS)
540 else if (approx_ge (sp[0].f, sp[i].f)
541 && approx_le (sp[0].f, sp[i + 1].f))
548 sp->f = sysmis ? SYSMIS : 0.0;
551 case OP_RANGE_STRING:
556 for (i = 1; i <= n_args; i += 2)
557 if (st_compare_pad (&sp[0].c[1], sp[0].c[0],
558 &sp[i].c[1], sp[i].c[0]) >= 0
559 && st_compare_pad (&sp[0].c[1], sp[0].c[0],
560 &sp[i + 1].c[1], sp[i + 1].c[0]) <= 0)
574 sum[0] = sum[1] = 0.0;
577 for (i = 0; i < n_args; i++)
578 if (sp[i].f != SYSMIS)
582 sum[1] += sp[i].f * sp[i].f;
587 sp->f = calc_stddev (calc_variance (sum, nv));
597 for (i = 0; i < n_args; i++)
598 if (sp[i].f != SYSMIS)
615 sum[0] = sum[1] = 0.0;
618 for (i = 0; i < n_args; i++)
619 if (sp[i].f != SYSMIS)
623 sum[1] += sp[i].f * sp[i].f;
628 sp->f = calc_variance (sum, nv);
632 /* Time construction function. */
635 if (sp[0].f == SYSMIS || sp[1].f == SYSMIS || sp[2].f == SYSMIS)
638 sp->f = 60. * (60. * sp[0].f + sp[1].f) + sp[2].f;
641 /* Date construction functions. */
644 sp->f = yrmoda (sp[2].f, sp[1].f, sp[0].f);
646 sp->f *= 60. * 60. * 24.;
650 sp->f = yrmoda (sp[2].f, sp[0].f, sp[1].f);
652 sp->f *= 60. * 60. * 24.;
656 sp->f = yrmoda (sp[1].f, sp[0].f, 1);
658 sp->f *= 60. * 60. * 24.;
662 if (sp[0].f == SYSMIS)
666 sp->f = yrmoda (sp[1].f, sp[0].f * 3 - 2, 1);
668 sp->f *= 60. * 60. * 24.;
673 if (sp[0].f == SYSMIS)
677 sp[1].f = yrmoda (sp[1].f, 1, 1);
679 sp[1].f = 60. * 60. * 24. * (sp[1].f + 7. * (floor (sp[0].f) - 1.));
685 if (sp[1].f == SYSMIS)
689 sp->f = yrmoda (sp[0].f, 1, 1);
691 sp->f = 60. * 60. * 24. * (sp->f + floor (sp[1].f) - 1);
696 sp->f = yrmoda (sp[0].f, sp[1].f, sp[2].f);
699 /* Date extraction functions. */
702 sp->f = floor (sp->f / 60. / 60. / 24.) * 60. * 60. * 24.;
706 sp->f = fmod (floor (sp->f / 60. / 60.), 24.);
710 sp->f = 86400. * julian_to_jday (sp->f / 86400.);
716 julian_to_calendar (sp->f / 86400., NULL, NULL, &day);
720 case OP_XDATE_MINUTE:
722 sp->f = fmod (floor (sp->f / 60.), 60.);
728 julian_to_calendar (sp->f / 86400., NULL, &month, NULL);
732 case OP_XDATE_QUARTER:
736 julian_to_calendar (sp->f / 86400., NULL, &month, NULL);
737 sp->f = (month - 1) / 3 + 1;
740 case OP_XDATE_SECOND:
742 sp->f = fmod (sp->f, 60.);
746 sp->f = floor (sp->f / 60. / 60. / 24.);
750 sp->f -= floor (sp->f / 60. / 60. / 24.) * 60. * 60. * 24.;
754 sp->f = (julian_to_jday (sp->f / 86400.) - 1) / 7 + 1;
758 sp->f = julian_to_wday (sp->f / 86400.);
764 julian_to_calendar (sp->f / 86400., &year, NULL, NULL);
769 /* String functions. */
775 dest = pool_alloc (e->pool, 256);
779 for (i = 0; i < n_args; i++)
782 if (sp[i].c[0] + dest[0] < 255)
784 memcpy (&dest[dest[0] + 1], &sp[i].c[1], sp[i].c[0]);
785 dest[0] += sp[i].c[0];
789 memcpy (&dest[dest[0] + 1], &sp[i].c[1], 255 - dest[0]);
803 int last = sp[0].c[0] - sp[1].c[0];
804 for (i = 0; i <= last; i++)
805 if (!memcmp (&sp[0].c[i + 1], &sp[0].c[1], sp[0].c[0]))
815 /* Length of each search string. */
816 int part_len = sp[2].f;
819 if (sp[1].c[0] == 0 || part_len <= 0 || sp[2].f == SYSMIS
820 || sp[1].c[0] % part_len != 0)
824 /* Last possible index. */
825 int last = sp[0].c[0] - part_len;
827 for (i = 0; i <= last; i++)
828 for (j = 0; j < sp[1].c[0]; j += part_len)
829 if (!memcmp (&sp[0].c[i], &sp[1].c[j], part_len))
844 for (i = sp[0].c[0] - sp[1].c[0]; i >= 0; i--)
845 if (!memcmp (&sp[0].c[i + 1], &sp[0].c[1], sp[0].c[0]))
855 /* Length of each search string. */
856 int part_len = sp[2].f;
859 if (sp[1].c[0] == 0 || part_len <= 0 || sp[2].f == SYSMIS
860 || sp[1].c[0] % part_len != 0)
864 for (i = sp[0].c[0] - part_len; i >= 0; i--)
865 for (j = 0; j < sp[1].c[0]; j += part_len)
866 if (!memcmp (&sp[0].c[i], &sp[1].c[j], part_len))
879 for (i = sp[0].c[0]; i >= 1; i--)
880 sp[0].c[i] = tolower ((unsigned char) (sp[0].c[i]));
883 for (i = sp[0].c[0]; i >= 1; i--)
884 sp[0].c[i] = toupper ((unsigned char) (sp[0].c[i]));
891 if (sp[1].f == SYSMIS || len < 0 || len > 255)
893 else if (len > sp[0].c[0])
897 dest = pool_alloc (e->pool, len + 1);
899 memset (&dest[1], ' ', len - sp->c[0]);
900 memcpy (&dest[len - sp->c[0] + 1], &sp->c[1], sp->c[0]);
910 if (sp[1].f == SYSMIS || len < 0 || len > 255 || sp[2].c[0] != 1)
912 else if (len > sp[0].c[0])
916 dest = pool_alloc (e->pool, len + 1);
918 memset (&dest[1], sp[2].c[1], len - sp->c[0]);
919 memcpy (&dest[len - sp->c[0] + 1], &sp->c[1], sp->c[0]);
929 if (sp[1].f == SYSMIS || len < 0 || len > 255)
931 else if (len > sp[0].c[0])
935 dest = pool_alloc (e->pool, len + 1);
937 memcpy (&dest[1], &sp->c[1], sp->c[0]);
938 memset (&dest[sp->c[0] + 1], ' ', len - sp->c[0]);
948 if (len < 0 || len > 255 || sp[2].c[0] != 1)
950 else if (len > sp[0].c[0])
954 dest = pool_alloc (e->pool, len + 1);
956 memcpy (&dest[1], &sp->c[1], sp->c[0]);
957 memset (&dest[sp->c[0] + 1], sp[2].c[1], len - sp->c[0]);
964 int len = sp[0].c[0];
967 while (i <= len && sp[0].c[i] == ' ')
971 sp[0].c[i] = sp[0].c[0] - i;
983 int len = sp[0].c[0];
984 int cmp = sp[1].c[1];
987 while (i <= len && sp[0].c[i] == cmp)
991 sp[0].c[i] = sp[0].c[0] - i;
999 while (sp[0].c[sp[0].c[0]] == ' ')
1004 if (sp[1].c[0] != 1)
1008 /* Note that NULs are not allowed in strings. This code
1009 needs to change if this decision is changed. */
1010 int cmp = sp[1].c[1];
1011 while (sp[0].c[sp[0].c[0]] == cmp)
1020 di.e = &sp->c[1] + sp->c[0];
1022 di.flags = DI_IGNORE_ERROR;
1024 di.format.type = FMT_F;
1025 di.format.w = sp->c[0];
1034 di.e = &sp->c[1] + sp->c[0];
1036 di.flags = DI_IGNORE_ERROR;
1038 di.format.type = *op++;
1039 di.format.w = *op++;
1040 di.format.d = *op++;
1047 unsigned char *dest;
1053 dest = pool_alloc (e->pool, f.w + 1);
1056 assert ((formats[f.type].cat & FCAT_STRING) == 0);
1057 data_out (&dest[1], &f, sp);
1067 if (index < 1 || index > sp[0].c[0])
1072 sp->c[index] = sp->c[0] - index;
1085 if (sp[1].f == SYSMIS || sp[2].f == SYSMIS || index < 1
1086 || index > sp[0].c[0] || n < 1)
1093 sp->c[index] = sp->c[0] - index;
1104 if (sp->f != SYSMIS)
1108 if (sp->f != SYSMIS)
1111 case OP_NUM_TO_BOOL:
1112 if (approx_eq (sp->f, 0.0))
1114 else if (approx_eq (sp->f, 1.0))
1116 else if (sp->f != SYSMIS)
1118 msg (SE, _("A number being treated as a Boolean in an "
1119 "expression was found to have a value other than "
1120 "0 (false), 1 (true), or the system-missing value. "
1121 "The result was forced to 0."));
1129 if (sp[0].f != SYSMIS)
1131 if (sp[1].f == SYSMIS)
1133 if (approx_ne (sp[0].f, 0.0))
1137 sp->f = fmod (sp[0].f, sp[1].f);
1141 if (sp->f != SYSMIS)
1142 sp->f *= rng_get_double_normal (pspp_rng ());
1145 if (sp->f != SYSMIS)
1146 sp->f *= rng_get_double (pspp_rng ());
1149 if (sp[0].f == SYSMIS || !finite (sp[0].f))
1154 case OP_VEC_ELEM_NUM:
1156 int rindx = sp[0].f + EPSILON;
1157 const struct vector *v = dict_get_vector (default_dict, *op++);
1159 if (sp[0].f == SYSMIS || rindx < 1 || rindx > v->cnt)
1161 if (sp[0].f == SYSMIS)
1162 msg (SE, _("SYSMIS is not a valid index value for vector "
1163 "%s. The result will be set to SYSMIS."),
1166 msg (SE, _("%g is not a valid index value for vector %s. "
1167 "The result will be set to SYSMIS."),
1172 sp->f = c->data[v->var[rindx - 1]->fv].f;
1175 case OP_VEC_ELEM_STR:
1177 int rindx = sp[0].f + EPSILON;
1178 const struct vector *vect = dict_get_vector (default_dict, *op++);
1181 if (sp[0].f == SYSMIS || rindx < 1 || rindx > vect->cnt)
1183 if (sp[0].f == SYSMIS)
1184 msg (SE, _("SYSMIS is not a valid index value for vector "
1185 "%s. The result will be set to the empty "
1189 msg (SE, _("%g is not a valid index value for vector %s. "
1190 "The result will be set to the empty string."),
1191 sp[0].f, vect->name);
1192 sp->c = pool_alloc (e->pool, 1);
1197 v = vect->var[rindx - 1];
1198 sp->c = pool_alloc (e->pool, v->width + 1);
1199 sp->c[0] = v->width;
1200 memcpy (&sp->c[1], c->data[v->fv].s, v->width);
1211 sp->c = pool_alloc (e->pool, *str + 1);
1212 memcpy (sp->c, str, *str + 1);
1217 sp->f = c->data[(*vars)->fv].f;
1218 if (is_num_user_missing (sp->f, *vars))
1224 int width = (*vars)->width;
1227 sp->c = pool_alloc (e->pool, width + 1);
1229 memcpy (&sp->c[1], &c->data[(*vars)->fv], width);
1235 struct ccase *c = lagged_case (*op++);
1242 sp->f = c->data[(*vars)->fv].f;
1243 if (is_num_user_missing (sp->f, *vars))
1251 struct ccase *c = lagged_case (*op++);
1252 int width = (*vars)->width;
1255 sp->c = pool_alloc (e->pool, width + 1);
1259 memset (sp->c, ' ', width);
1261 memcpy (&sp->c[1], &c->data[(*vars)->fv], width);
1268 sp->f = c->data[*op++].f == SYSMIS;
1272 sp->f = is_str_user_missing (c->data[(*vars)->fv].s, *vars);
1277 sp->f = c->data[*op++].f;
1281 sp->f = vfm_sink_info.ncases + 1;
1294 if (e->type != EX_STRING)
1296 double value = sp->f;
1297 if (!finite (value))