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 "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)
108 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);
337 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;
412 for (i = 1; i <= n_args; i++)
413 if (!st_compare_pad (&sp[0].c[1], sp[0].c[0],
414 &sp[i].c[1], sp[i].c[0]))
430 for (i = 0; i < n_args; i++)
431 if (sp[i].f != SYSMIS)
435 sum[1] += sp[i].f * sp[i].f;
440 sp->f = calc_cfvar (sum, nv);
447 double max = -DBL_MAX;
450 for (i = 0; i < n_args; i++)
451 if (sp[i].f != SYSMIS)
471 for (i = 0; i < n_args; i++)
472 if (sp[i].f != SYSMIS)
480 sp->f = calc_mean (sum, nv);
487 double min = DBL_MAX;
490 for (i = 0; i < n_args; i++)
491 if (sp[i].f != SYSMIS)
509 for (i = 0; i < n_args; i++)
510 if (sp[i].f == SYSMIS)
521 for (i = 0; i < n_args; i++)
522 if (sp[i].f != SYSMIS)
535 for (i = 1; i <= n_args; i += 2)
536 if (sp[i].f == SYSMIS || sp[i + 1].f == SYSMIS)
538 else if (sp[0].f >= sp[i].f && sp[0].f <= sp[i + 1].f)
545 sp->f = sysmis ? SYSMIS : 0.0;
548 case OP_RANGE_STRING:
553 for (i = 1; i <= n_args; i += 2)
554 if (st_compare_pad (&sp[0].c[1], sp[0].c[0],
555 &sp[i].c[1], sp[i].c[0]) >= 0
556 && st_compare_pad (&sp[0].c[1], sp[0].c[0],
557 &sp[i + 1].c[1], sp[i + 1].c[0]) <= 0)
571 sum[0] = sum[1] = 0.0;
574 for (i = 0; i < n_args; i++)
575 if (sp[i].f != SYSMIS)
579 sum[1] += sp[i].f * sp[i].f;
584 sp->f = calc_stddev (calc_variance (sum, nv));
594 for (i = 0; i < n_args; i++)
595 if (sp[i].f != SYSMIS)
612 sum[0] = sum[1] = 0.0;
615 for (i = 0; i < n_args; i++)
616 if (sp[i].f != SYSMIS)
620 sum[1] += sp[i].f * sp[i].f;
625 sp->f = calc_variance (sum, nv);
629 /* Time construction function. */
632 if (sp[0].f == SYSMIS || sp[1].f == SYSMIS || sp[2].f == SYSMIS)
635 sp->f = 60. * (60. * sp[0].f + sp[1].f) + sp[2].f;
638 /* Date construction functions. */
641 sp->f = yrmoda (sp[2].f, sp[1].f, sp[0].f);
643 sp->f *= 60. * 60. * 24.;
647 sp->f = yrmoda (sp[2].f, sp[0].f, sp[1].f);
649 sp->f *= 60. * 60. * 24.;
653 sp->f = yrmoda (sp[1].f, sp[0].f, 1);
655 sp->f *= 60. * 60. * 24.;
659 if (sp[0].f == SYSMIS)
663 sp->f = yrmoda (sp[1].f, sp[0].f * 3 - 2, 1);
665 sp->f *= 60. * 60. * 24.;
670 if (sp[0].f == SYSMIS)
674 sp[1].f = yrmoda (sp[1].f, 1, 1);
676 sp[1].f = 60. * 60. * 24. * (sp[1].f + 7. * (floor (sp[0].f) - 1.));
682 if (sp[1].f == SYSMIS)
686 sp->f = yrmoda (sp[0].f, 1, 1);
688 sp->f = 60. * 60. * 24. * (sp->f + floor (sp[1].f) - 1);
693 sp->f = yrmoda (sp[0].f, sp[1].f, sp[2].f);
696 /* Date extraction functions. */
699 sp->f = floor (sp->f / 60. / 60. / 24.) * 60. * 60. * 24.;
703 sp->f = fmod (floor (sp->f / 60. / 60.), 24.);
707 sp->f = 86400. * julian_to_jday (sp->f / 86400.);
713 julian_to_calendar (sp->f / 86400., NULL, NULL, &day);
717 case OP_XDATE_MINUTE:
719 sp->f = fmod (floor (sp->f / 60.), 60.);
725 julian_to_calendar (sp->f / 86400., NULL, &month, NULL);
729 case OP_XDATE_QUARTER:
733 julian_to_calendar (sp->f / 86400., NULL, &month, NULL);
734 sp->f = (month - 1) / 3 + 1;
737 case OP_XDATE_SECOND:
739 sp->f = fmod (sp->f, 60.);
743 sp->f = floor (sp->f / 60. / 60. / 24.);
747 sp->f -= floor (sp->f / 60. / 60. / 24.) * 60. * 60. * 24.;
751 sp->f = (julian_to_jday (sp->f / 86400.) - 1) / 7 + 1;
755 sp->f = julian_to_wday (sp->f / 86400.);
761 julian_to_calendar (sp->f / 86400., &year, NULL, NULL);
766 /* String functions. */
772 dest = pool_alloc (e->pool, 256);
776 for (i = 0; i < n_args; i++)
779 if (sp[i].c[0] + dest[0] < 255)
781 memcpy (&dest[dest[0] + 1], &sp[i].c[1], sp[i].c[0]);
782 dest[0] += sp[i].c[0];
786 memcpy (&dest[dest[0] + 1], &sp[i].c[1], 255 - dest[0]);
800 int last = sp[0].c[0] - sp[1].c[0];
801 for (i = 0; i <= last; i++)
802 if (!memcmp (&sp[0].c[i + 1], &sp[0].c[1], sp[0].c[0]))
812 /* Length of each search string. */
813 int part_len = sp[2].f;
816 if (sp[1].c[0] == 0 || part_len <= 0 || sp[2].f == SYSMIS
817 || sp[1].c[0] % part_len != 0)
821 /* Last possible index. */
822 int last = sp[0].c[0] - part_len;
824 for (i = 0; i <= last; i++)
825 for (j = 0; j < sp[1].c[0]; j += part_len)
826 if (!memcmp (&sp[0].c[i], &sp[1].c[j], part_len))
841 for (i = sp[0].c[0] - sp[1].c[0]; i >= 0; i--)
842 if (!memcmp (&sp[0].c[i + 1], &sp[0].c[1], sp[0].c[0]))
852 /* Length of each search string. */
853 int part_len = sp[2].f;
856 if (sp[1].c[0] == 0 || part_len <= 0 || sp[2].f == SYSMIS
857 || sp[1].c[0] % part_len != 0)
861 for (i = sp[0].c[0] - part_len; i >= 0; i--)
862 for (j = 0; j < sp[1].c[0]; j += part_len)
863 if (!memcmp (&sp[0].c[i], &sp[1].c[j], part_len))
876 for (i = sp[0].c[0]; i >= 1; i--)
877 sp[0].c[i] = tolower ((unsigned char) (sp[0].c[i]));
880 for (i = sp[0].c[0]; i >= 1; i--)
881 sp[0].c[i] = toupper ((unsigned char) (sp[0].c[i]));
888 if (sp[1].f == SYSMIS || len < 0 || len > 255)
890 else if (len > sp[0].c[0])
894 dest = pool_alloc (e->pool, len + 1);
896 memset (&dest[1], ' ', len - sp->c[0]);
897 memcpy (&dest[len - sp->c[0] + 1], &sp->c[1], sp->c[0]);
907 if (sp[1].f == SYSMIS || len < 0 || len > 255 || sp[2].c[0] != 1)
909 else if (len > sp[0].c[0])
913 dest = pool_alloc (e->pool, len + 1);
915 memset (&dest[1], sp[2].c[1], len - sp->c[0]);
916 memcpy (&dest[len - sp->c[0] + 1], &sp->c[1], sp->c[0]);
926 if (sp[1].f == SYSMIS || len < 0 || len > 255)
928 else if (len > sp[0].c[0])
932 dest = pool_alloc (e->pool, len + 1);
934 memcpy (&dest[1], &sp->c[1], sp->c[0]);
935 memset (&dest[sp->c[0] + 1], ' ', len - sp->c[0]);
945 if (len < 0 || len > 255 || sp[2].c[0] != 1)
947 else if (len > sp[0].c[0])
951 dest = pool_alloc (e->pool, len + 1);
953 memcpy (&dest[1], &sp->c[1], sp->c[0]);
954 memset (&dest[sp->c[0] + 1], sp[2].c[1], len - sp->c[0]);
961 int len = sp[0].c[0];
964 while (i <= len && sp[0].c[i] == ' ')
968 sp[0].c[i] = sp[0].c[0] - i;
980 int len = sp[0].c[0];
981 int cmp = sp[1].c[1];
984 while (i <= len && sp[0].c[i] == cmp)
988 sp[0].c[i] = sp[0].c[0] - i;
996 while (sp[0].c[sp[0].c[0]] == ' ')
1001 if (sp[1].c[0] != 1)
1005 /* Note that NULs are not allowed in strings. This code
1006 needs to change if this decision is changed. */
1007 int cmp = sp[1].c[1];
1008 while (sp[0].c[sp[0].c[0]] == cmp)
1017 di.e = &sp->c[1] + sp->c[0];
1019 di.flags = DI_IGNORE_ERROR;
1021 di.format.type = FMT_F;
1022 di.format.w = sp->c[0];
1031 di.e = &sp->c[1] + sp->c[0];
1033 di.flags = DI_IGNORE_ERROR;
1035 di.format.type = *op++;
1036 di.format.w = *op++;
1037 di.format.d = *op++;
1044 unsigned char *dest;
1050 dest = pool_alloc (e->pool, f.w + 1);
1053 assert ((formats[f.type].cat & FCAT_STRING) == 0);
1054 data_out (&dest[1], &f, sp);
1064 if (index < 1 || index > sp[0].c[0])
1069 sp->c[index] = sp->c[0] - index;
1082 if (sp[1].f == SYSMIS || sp[2].f == SYSMIS || index < 1
1083 || index > sp[0].c[0] || n < 1)
1090 sp->c[index] = sp->c[0] - index;
1101 if (sp->f != SYSMIS)
1105 if (sp->f != SYSMIS)
1108 case OP_NUM_TO_BOOL:
1111 else if (sp->f == 1.0)
1113 else if (sp->f != SYSMIS)
1115 msg (SE, _("A number being treated as a Boolean in an "
1116 "expression was found to have a value other than "
1117 "0 (false), 1 (true), or the system-missing value. "
1118 "The result was forced to 0."));
1126 if (sp[0].f != SYSMIS)
1128 if (sp[1].f == SYSMIS)
1134 sp->f = fmod (sp[0].f, sp[1].f);
1138 if (sp->f != SYSMIS)
1139 sp->f *= rng_get_double_normal (pspp_rng ());
1142 if (sp->f != SYSMIS)
1143 sp->f *= rng_get_double (pspp_rng ());
1146 if (sp[0].f == SYSMIS || !finite (sp[0].f))
1151 case OP_VEC_ELEM_NUM:
1153 int rindx = sp[0].f + EPSILON;
1154 const struct vector *v = dict_get_vector (default_dict, *op++);
1156 if (sp[0].f == SYSMIS || rindx < 1 || rindx > v->cnt)
1158 if (sp[0].f == SYSMIS)
1159 msg (SE, _("SYSMIS is not a valid index value for vector "
1160 "%s. The result will be set to SYSMIS."),
1163 msg (SE, _("%g is not a valid index value for vector %s. "
1164 "The result will be set to SYSMIS."),
1169 sp->f = c->data[v->var[rindx - 1]->fv].f;
1172 case OP_VEC_ELEM_STR:
1174 int rindx = sp[0].f + EPSILON;
1175 const struct vector *vect = dict_get_vector (default_dict, *op++);
1178 if (sp[0].f == SYSMIS || rindx < 1 || rindx > vect->cnt)
1180 if (sp[0].f == SYSMIS)
1181 msg (SE, _("SYSMIS is not a valid index value for vector "
1182 "%s. The result will be set to the empty "
1186 msg (SE, _("%g is not a valid index value for vector %s. "
1187 "The result will be set to the empty string."),
1188 sp[0].f, vect->name);
1189 sp->c = pool_alloc (e->pool, 1);
1194 v = vect->var[rindx - 1];
1195 sp->c = pool_alloc (e->pool, v->width + 1);
1196 sp->c[0] = v->width;
1197 memcpy (&sp->c[1], c->data[v->fv].s, v->width);
1208 sp->c = pool_alloc (e->pool, *str + 1);
1209 memcpy (sp->c, str, *str + 1);
1214 sp->f = c->data[(*vars)->fv].f;
1215 if (is_num_user_missing (sp->f, *vars))
1221 int width = (*vars)->width;
1224 sp->c = pool_alloc (e->pool, width + 1);
1226 memcpy (&sp->c[1], &c->data[(*vars)->fv], width);
1232 struct ccase *c = lagged_case (*op++);
1239 sp->f = c->data[(*vars)->fv].f;
1240 if (is_num_user_missing (sp->f, *vars))
1248 struct ccase *c = lagged_case (*op++);
1249 int width = (*vars)->width;
1252 sp->c = pool_alloc (e->pool, width + 1);
1256 memset (sp->c, ' ', width);
1258 memcpy (&sp->c[1], &c->data[(*vars)->fv], width);
1265 sp->f = c->data[*op++].f == SYSMIS;
1269 sp->f = is_str_user_missing (c->data[(*vars)->fv].s, *vars);
1274 sp->f = c->data[*op++].f;
1278 sp->f = vfm_sink_info.ncases + 1;
1291 if (e->type != EX_STRING)
1293 double value = sp->f;
1294 if (!finite (value))