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 (const struct expression *e, const struct ccase *c, int case_num,
57 unsigned char *op = e->op;
59 unsigned char *str = e->str;
60 struct variable **vars = e->var;
64 union value *sp = e->stack;
74 if (sp[1].f == SYSMIS)
76 else if (sp[0].f != SYSMIS)
81 if (sp[1].f == SYSMIS)
83 else if (sp[0].f != SYSMIS)
88 if (sp[1].f == SYSMIS)
90 else if (sp[0].f != SYSMIS)
95 if (sp[1].f == SYSMIS || sp[1].f == 0.)
97 else if (sp[0].f != SYSMIS)
102 if (sp[0].f == SYSMIS)
107 else if (sp[1].f == SYSMIS)
114 else if (sp[0].f == 0.0 && sp[1].f <= 0.0)
117 sp->f = pow (sp[0].f, sp[1].f);
121 /* Note that booleans are always one of 0, 1, or SYSMIS.
123 Truth table (in order of detection):
135 1 and SYSMIS = SYSMIS
136 SYSMIS and SYSMIS = SYSMIS
140 SYSMIS and 1 = SYSMIS
144 if (sp[0].f == 0.0); /* 1 */
145 else if (sp[1].f == 0.0)
147 else if (sp[1].f == SYSMIS)
148 sp->f = SYSMIS; /* 3 */
151 /* Truth table (in order of detection):
164 SYSMIS or SYSMIS = SYSMIS
172 if (sp[0].f == 1.0); /* 1 */
173 else if (sp[1].f == 1.0)
175 else if (sp[1].f == SYSMIS)
176 sp->f = SYSMIS; /* 3 */
181 else if (sp[0].f == 1.0)
187 if (sp[0].f != SYSMIS)
189 if (sp[1].f == SYSMIS)
192 sp->f = sp[0].f == sp[1].f;
197 if (sp[0].f != SYSMIS)
199 if (sp[1].f == SYSMIS)
202 sp->f = sp[0].f >= sp[1].f;
207 if (sp[0].f != SYSMIS)
209 if (sp[1].f == SYSMIS)
212 sp->f = sp[0].f > sp[1].f;
217 if (sp[0].f != SYSMIS)
219 if (sp[1].f == SYSMIS)
222 sp->f = sp[0].f <= sp[1].f;
227 if (sp[0].f != SYSMIS)
229 if (sp[1].f == SYSMIS)
232 sp->f = sp[0].f < sp[1].f;
237 if (sp[0].f != SYSMIS)
239 if (sp[1].f == SYSMIS)
242 sp->f = sp[0].f != sp[1].f;
246 /* String operators. */
249 sp[0].f = st_compare_pad (&sp[0].c[1], sp[0].c[0],
250 &sp[1].c[1], sp[1].c[0]) == 0;
254 sp[0].f = st_compare_pad (&sp[0].c[1], sp[0].c[0],
255 &sp[1].c[1], sp[1].c[0]) >= 0;
259 sp[0].f = st_compare_pad (&sp[0].c[1], sp[0].c[0],
260 &sp[1].c[1], sp[1].c[0]) > 0;
264 sp[0].f = st_compare_pad (&sp[0].c[1], sp[0].c[0],
265 &sp[1].c[1], sp[1].c[0]) <= 0;
269 sp[0].f = st_compare_pad (&sp[0].c[1], sp[0].c[0],
270 &sp[1].c[1], sp[1].c[0]) < 0;
274 sp[0].f = st_compare_pad (&sp[0].c[1], sp[0].c[0],
275 &sp[1].c[1], sp[1].c[0]) != 0;
278 /* Unary functions. */
285 sp->f = fabs (sp->f);
291 sp->f = acos (sp->f);
300 sp->f = asin (sp->f);
307 sp->f = atan (sp->f);
326 sp->f = log10 (sp->f);
342 sp->f = fmod (sp->f, 10);
348 sp->f = floor (sp->f + 0.5);
350 sp->f = -floor (-sp->f + 0.5);
361 sp->f = sqrt (sp->f);
379 sp->f = floor (sp->f);
381 sp->f = -floor (-sp->f);
385 /* N-ary numeric functions. */
394 for (i = 1; i < n_args; i++)
395 if (sp[0].f == sp[i].f)
400 else if (sp[i].f != SYSMIS)
402 sp->f = sysmis ? SYSMIS : 0.0;
411 for (i = 1; i < n_args; i++)
412 if (!st_compare_pad (&sp[0].c[1], sp[0].c[0],
413 &sp[i].c[1], sp[i].c[0]))
429 for (i = 0; i < n_args; i++)
430 if (sp[i].f != SYSMIS)
434 sum[1] += sp[i].f * sp[i].f;
439 sp->f = calc_cfvar (sum, nv);
446 double max = -DBL_MAX;
449 for (i = 0; i < n_args; i++)
450 if (sp[i].f != SYSMIS)
468 for (i = 1; i < n_args; i++)
469 if (st_compare_pad (&sp[i].c[1], sp[i].c[0],
470 &sp[max_idx].c[1], sp[max_idx].c[0]) > 0)
472 sp[0].c = sp[max_idx].c;
483 for (i = 0; i < n_args; i++)
484 if (sp[i].f != SYSMIS)
492 sp->f = calc_mean (sum, nv);
499 double min = DBL_MAX;
502 for (i = 0; i < n_args; i++)
503 if (sp[i].f != SYSMIS)
521 for (i = 1; i < n_args; i++)
522 if (st_compare_pad (&sp[i].c[1], sp[i].c[0],
523 &sp[min_idx].c[1], sp[min_idx].c[0]) < 0)
525 sp[0].c = sp[min_idx].c;
534 for (i = 0; i < n_args; i++)
535 if (sp[i].f == SYSMIS)
546 for (i = 0; i < n_args; i++)
547 if (sp[i].f != SYSMIS)
560 for (i = 1; i < n_args; i += 2)
561 if (sp[i].f == SYSMIS || sp[i + 1].f == SYSMIS)
563 else if (sp[0].f >= sp[i].f && sp[0].f <= sp[i + 1].f)
570 sp->f = sysmis ? SYSMIS : 0.0;
573 case OP_RANGE_STRING:
578 for (i = 1; i < n_args; i += 2)
579 if (st_compare_pad (&sp[0].c[1], sp[0].c[0],
580 &sp[i].c[1], sp[i].c[0]) >= 0
581 && st_compare_pad (&sp[0].c[1], sp[0].c[0],
582 &sp[i + 1].c[1], sp[i + 1].c[0]) <= 0)
596 sum[0] = sum[1] = 0.0;
599 for (i = 0; i < n_args; i++)
600 if (sp[i].f != SYSMIS)
604 sum[1] += sp[i].f * sp[i].f;
609 sp->f = calc_stddev (calc_variance (sum, nv));
619 for (i = 0; i < n_args; i++)
620 if (sp[i].f != SYSMIS)
637 sum[0] = sum[1] = 0.0;
640 for (i = 0; i < n_args; i++)
641 if (sp[i].f != SYSMIS)
645 sum[1] += sp[i].f * sp[i].f;
650 sp->f = calc_variance (sum, nv);
654 /* Time construction function. */
657 if (sp[0].f == SYSMIS || sp[1].f == SYSMIS || sp[2].f == SYSMIS)
662 min = min (sp[0].f, min (sp[1].f, sp[2].f));
663 max = max (sp[0].f, max (sp[1].f, sp[2].f));
664 if (min < 0. && max > 0.)
666 msg (SW, _("TIME.HMS cannot mix positive and negative "
667 "in its arguments."));
671 sp->f = 60. * (60. * sp[0].f + sp[1].f) + sp[2].f;
676 sp->f /= 60. * 60. * 24.;
682 case OP_CTIME_MINUTES:
688 sp->f *= 60. * 60. * 24.;
690 case OP_CTIME_SECONDS:
694 /* Date construction functions. */
697 sp->f = yrmoda (sp[2].f, sp[1].f, sp[0].f);
699 sp->f *= 60. * 60. * 24.;
703 sp->f = yrmoda (sp[2].f, sp[0].f, sp[1].f);
705 sp->f *= 60. * 60. * 24.;
709 sp->f = yrmoda (sp[1].f, sp[0].f, 1);
711 sp->f *= 60. * 60. * 24.;
715 if (sp[0].f == SYSMIS)
719 sp->f = yrmoda (sp[1].f, sp[0].f * 3 - 2, 1);
721 sp->f *= 60. * 60. * 24.;
726 if (sp[0].f == SYSMIS || sp[1].f == SYSMIS)
728 else if (sp[0].f < 0. || sp[0].f > 53.)
730 msg (SW, _("Week argument to WKYR must be in range 0 to 53."));
735 double result = yrmoda (sp[1].f, 1, 1);
736 if (result != SYSMIS)
737 result = 60. * 60. * 24. * (result + 7. * (sp->f - 1.));
743 if (sp[1].f == SYSMIS)
747 sp->f = yrmoda (sp[0].f, 1, 1);
749 sp->f = 60. * 60. * 24. * (sp->f + floor (sp[1].f) - 1);
754 sp->f = yrmoda (sp[0].f, sp[1].f, sp[2].f);
757 /* Date extraction functions. */
760 sp->f = floor (sp->f / 60. / 60. / 24.) * 60. * 60. * 24.;
764 sp->f = fmod (floor (sp->f / 60. / 60.), 24.);
768 sp->f = 86400. * julian_to_jday (sp->f / 86400.);
774 julian_to_calendar (sp->f / 86400., NULL, NULL, &day);
778 case OP_XDATE_MINUTE:
780 sp->f = fmod (floor (sp->f / 60.), 60.);
786 julian_to_calendar (sp->f / 86400., NULL, &month, NULL);
790 case OP_XDATE_QUARTER:
794 julian_to_calendar (sp->f / 86400., NULL, &month, NULL);
795 sp->f = (month - 1) / 3 + 1;
798 case OP_XDATE_SECOND:
800 sp->f = fmod (sp->f, 60.);
804 sp->f = floor (sp->f / 60. / 60. / 24.);
808 sp->f -= floor (sp->f / 60. / 60. / 24.) * 60. * 60. * 24.;
812 sp->f = (julian_to_jday (sp->f / 86400.) - 1) / 7 + 1;
816 sp->f = julian_to_wday (sp->f / 86400.);
822 julian_to_calendar (sp->f / 86400., &year, NULL, NULL);
827 /* String functions. */
833 dest = pool_alloc (e->pool, 256);
837 for (i = 0; i < n_args; i++)
840 if (sp[i].c[0] + dest[0] < 255)
842 memcpy (&dest[dest[0] + 1], &sp[i].c[1], sp[i].c[0]);
843 dest[0] += sp[i].c[0];
847 memcpy (&dest[dest[0] + 1], &sp[i].c[1], 255 - dest[0]);
861 int last = sp[0].c[0] - sp[1].c[0];
863 for (i = 0; i <= last; i++)
864 if (!memcmp (&sp[0].c[i + 1], &sp[1].c[1], sp[1].c[0]))
879 else if (sp[2].f == SYSMIS)
881 msg (SW, _("Argument 3 of RINDEX may not be system-missing."));
886 int part_len = sp[2].f;
888 if (part_len < 0 || part_len > sp[1].c[0]
889 || sp[1].c[0] % part_len != 0)
891 msg (SW, _("Argument 3 of RINDEX must be between 1 and "
892 "the length of argument 2, and it must "
893 "evenly divide the length of argument 2."));
899 int last = sp[0].c[0] - part_len;
900 for (i = 0; i <= last; i++)
901 for (j = 0; j < sp[1].c[0]; j += part_len)
902 if (!memcmp (&sp[0].c[i + 1], &sp[1].c[j + 1], part_len))
919 for (i = sp[0].c[0] - sp[1].c[0]; i >= 0; i--)
920 if (!memcmp (&sp[0].c[i + 1], &sp[1].c[1], sp[1].c[0]))
935 else if (sp[2].f == SYSMIS)
937 msg (SW, _("Argument 3 of RINDEX may not be system-missing."));
942 int part_len = sp[2].f;
944 if (part_len < 0 || part_len > sp[1].c[0]
945 || sp[1].c[0] % part_len != 0)
947 msg (SW, _("Argument 3 of RINDEX must be between 1 and "
948 "the length of argument 2, and it must "
949 "evenly divide the length of argument 2."));
954 for (i = sp[0].c[0] - part_len; i >= 0; i--)
955 for (j = 0; j < sp[1].c[0]; j += part_len)
956 if (!memcmp (&sp[0].c[i + 1], &sp[1].c[j + 1], part_len))
970 for (i = sp[0].c[0]; i >= 1; i--)
971 sp[0].c[i] = tolower ((unsigned char) (sp[0].c[i]));
974 for (i = sp[0].c[0]; i >= 1; i--)
975 sp[0].c[i] = toupper ((unsigned char) (sp[0].c[i]));
982 if (sp[1].f == SYSMIS || len < 0 || len > 255 || sp[2].c[0] != 1)
984 else if (len > sp[0].c[0])
988 dest = pool_alloc (e->pool, len + 1);
990 memset (&dest[1], sp[2].c[1], len - sp->c[0]);
991 memcpy (&dest[len - sp->c[0] + 1], &sp->c[1], sp->c[0]);
1001 if (len < 0 || len > 255 || sp[2].c[0] != 1)
1003 else if (len > sp[0].c[0])
1005 unsigned char *dest;
1007 dest = pool_alloc (e->pool, len + 1);
1009 memcpy (&dest[1], &sp->c[1], sp->c[0]);
1010 memset (&dest[sp->c[0] + 1], sp[2].c[1], len - sp->c[0]);
1018 if (sp[1].c[0] != 1)
1022 int len = sp[0].c[0];
1023 int cmp = sp[1].c[1];
1026 while (i <= len && sp[0].c[i] == cmp)
1030 sp[0].c[i] = sp[0].c[0] - i;
1031 sp->c = &sp[0].c[i];
1038 if (sp[1].c[0] != 1)
1042 int cmp = sp[1].c[1];
1043 while (sp[0].c[0] > 0 && sp[0].c[sp[0].c[0]] == cmp)
1055 di.format.type = *op++;
1056 di.format.w = *op++;
1057 di.format.d = *op++;
1058 di.e = &sp->c[1] + min (sp->c[0], di.format.w);
1066 unsigned char *dest;
1072 dest = pool_alloc (e->pool, f.w + 1);
1075 assert ((formats[f.type].cat & FCAT_STRING) == 0);
1076 data_out (&dest[1], &f, sp);
1086 if (index < 1 || index > sp[0].c[0])
1091 sp->c[index] = sp->c[0] - index;
1104 if (sp[1].f == SYSMIS || sp[2].f == SYSMIS || index < 1
1105 || index > sp[0].c[0] || n < 1)
1112 sp->c[index] = sp->c[0] - index;
1123 if (sp->f != SYSMIS)
1126 case OP_NUM_TO_BOOL:
1129 else if (sp->f == 1.0)
1131 else if (sp->f != SYSMIS)
1133 msg (SE, _("A number being treated as a Boolean in an "
1134 "expression was found to have a value other than "
1135 "0 (false), 1 (true), or the system-missing value. "
1136 "The result was forced to 0."));
1144 if (sp[0].f != SYSMIS)
1146 if (sp[1].f == SYSMIS)
1152 sp->f = fmod (sp[0].f, sp[1].f);
1156 if (sp->f != SYSMIS)
1157 sp->f *= rng_get_double_normal (pspp_rng ());
1160 if (sp->f != SYSMIS)
1161 sp->f *= rng_get_double (pspp_rng ());
1164 sp->f = sp->f == SYSMIS || !finite (sp->f);
1166 case OP_VEC_ELEM_NUM:
1168 int rindx = sp[0].f + EPSILON;
1169 const struct vector *v = dict_get_vector (default_dict, *op++);
1171 if (sp[0].f == SYSMIS || rindx < 1 || rindx > v->cnt)
1173 if (sp[0].f == SYSMIS)
1174 msg (SE, _("SYSMIS is not a valid index value for vector "
1175 "%s. The result will be set to SYSMIS."),
1178 msg (SE, _("%g is not a valid index value for vector %s. "
1179 "The result will be set to SYSMIS."),
1185 sp->f = c->data[v->var[rindx - 1]->fv].f;
1188 case OP_VEC_ELEM_STR:
1190 int rindx = sp[0].f + EPSILON;
1191 const struct vector *vect = dict_get_vector (default_dict, *op++);
1194 if (sp[0].f == SYSMIS || rindx < 1 || rindx > vect->cnt)
1196 if (sp[0].f == SYSMIS)
1197 msg (SE, _("SYSMIS is not a valid index value for vector "
1198 "%s. The result will be set to the empty "
1202 msg (SE, _("%g is not a valid index value for vector %s. "
1203 "The result will be set to the empty string."),
1204 sp[0].f, vect->name);
1205 sp->c = pool_alloc (e->pool, 1);
1210 v = vect->var[rindx - 1];
1211 sp->c = pool_alloc (e->pool, v->width + 1);
1212 sp->c[0] = v->width;
1214 memcpy (&sp->c[1], c->data[v->fv].s, v->width);
1225 sp->c = pool_alloc (e->pool, *str + 1);
1226 memcpy (sp->c, str, *str + 1);
1232 sp->f = c->data[(*vars)->fv].f;
1233 if (is_num_user_missing (sp->f, *vars))
1239 int width = (*vars)->width;
1242 sp->c = pool_alloc (e->pool, width + 1);
1245 memcpy (&sp->c[1], &c->data[(*vars)->fv], width);
1251 struct ccase *c = lagged_case (*op++);
1258 sp->f = c->data[(*vars)->fv].f;
1259 if (is_num_user_missing (sp->f, *vars))
1267 struct ccase *c = lagged_case (*op++);
1268 int width = (*vars)->width;
1271 sp->c = pool_alloc (e->pool, width + 1);
1275 memset (sp->c, ' ', width);
1277 memcpy (&sp->c[1], &c->data[(*vars)->fv], width);
1285 sp->f = c->data[*op++].f == SYSMIS;
1290 sp->f = c->data[*op++].f;
1307 if (e->type != EXPR_STRING)
1309 double value = sp->f;
1310 if (!finite (value))