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"
55 expr_evaluate (const struct expression *e, const struct ccase *c, int case_idx,
58 unsigned char *op = e->op;
60 unsigned char *str = e->str;
61 struct variable **vars = e->var;
65 union value *sp = e->stack;
75 if (sp[1].f == SYSMIS)
77 else if (sp[0].f != SYSMIS)
82 if (sp[1].f == SYSMIS)
84 else if (sp[0].f != SYSMIS)
89 if (sp[1].f == SYSMIS)
91 else if (sp[0].f != SYSMIS)
96 if (sp[1].f == SYSMIS || sp[1].f == 0.)
98 else if (sp[0].f != SYSMIS)
103 if (sp[0].f == SYSMIS)
108 else if (sp[1].f == SYSMIS)
115 else if (sp[0].f == 0.0 && sp[1].f <= 0.0)
118 sp->f = pow (sp[0].f, sp[1].f);
122 /* Note that booleans are always one of 0, 1, or SYSMIS.
124 Truth table (in order of detection):
136 1 and SYSMIS = SYSMIS
137 SYSMIS and SYSMIS = SYSMIS
141 SYSMIS and 1 = SYSMIS
145 if (sp[0].f == 0.0); /* 1 */
146 else if (sp[1].f == 0.0)
148 else if (sp[1].f == SYSMIS)
149 sp->f = SYSMIS; /* 3 */
152 /* Truth table (in order of detection):
165 SYSMIS or SYSMIS = SYSMIS
173 if (sp[0].f == 1.0); /* 1 */
174 else if (sp[1].f == 1.0)
176 else if (sp[1].f == SYSMIS)
177 sp->f = SYSMIS; /* 3 */
182 else if (sp[0].f == 1.0)
188 if (sp[0].f != SYSMIS)
190 if (sp[1].f == SYSMIS)
193 sp->f = sp[0].f == sp[1].f;
198 if (sp[0].f != SYSMIS)
200 if (sp[1].f == SYSMIS)
203 sp->f = sp[0].f >= sp[1].f;
208 if (sp[0].f != SYSMIS)
210 if (sp[1].f == SYSMIS)
213 sp->f = sp[0].f > sp[1].f;
218 if (sp[0].f != SYSMIS)
220 if (sp[1].f == SYSMIS)
223 sp->f = sp[0].f <= sp[1].f;
228 if (sp[0].f != SYSMIS)
230 if (sp[1].f == SYSMIS)
233 sp->f = sp[0].f < sp[1].f;
238 if (sp[0].f != SYSMIS)
240 if (sp[1].f == SYSMIS)
243 sp->f = sp[0].f != sp[1].f;
247 /* String operators. */
250 sp[0].f = st_compare_pad (&sp[0].c[1], sp[0].c[0],
251 &sp[1].c[1], sp[1].c[0]) == 0;
255 sp[0].f = st_compare_pad (&sp[0].c[1], sp[0].c[0],
256 &sp[1].c[1], sp[1].c[0]) >= 0;
260 sp[0].f = st_compare_pad (&sp[0].c[1], sp[0].c[0],
261 &sp[1].c[1], sp[1].c[0]) > 0;
265 sp[0].f = st_compare_pad (&sp[0].c[1], sp[0].c[0],
266 &sp[1].c[1], sp[1].c[0]) <= 0;
270 sp[0].f = st_compare_pad (&sp[0].c[1], sp[0].c[0],
271 &sp[1].c[1], sp[1].c[0]) < 0;
275 sp[0].f = st_compare_pad (&sp[0].c[1], sp[0].c[0],
276 &sp[1].c[1], sp[1].c[0]) != 0;
279 /* Unary functions. */
286 sp->f = fabs (sp->f);
292 sp->f = acos (sp->f);
301 sp->f = asin (sp->f);
308 sp->f = atan (sp->f);
327 sp->f = log10 (sp->f);
343 sp->f = fmod (sp->f, 10);
349 sp->f = floor (sp->f + 0.5);
351 sp->f = -floor (-sp->f + 0.5);
362 sp->f = sqrt (sp->f);
380 sp->f = floor (sp->f);
382 sp->f = -floor (-sp->f);
386 /* N-ary numeric functions. */
395 for (i = 1; i < n_args; i++)
396 if (sp[0].f == sp[i].f)
401 else if (sp[i].f != SYSMIS)
403 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]))
425 double weight, mean, variance;
429 moments_of_values (sp, n_args,
430 &weight, &mean, &variance, NULL, NULL);
432 if (weight < *op++ || mean == SYSMIS
433 || mean == 0 || variance == SYSMIS)
436 sp->f = sqrt (variance) / mean;
443 double max = -DBL_MAX;
446 for (i = 0; i < n_args; i++)
447 if (sp[i].f != SYSMIS)
465 for (i = 1; i < n_args; i++)
466 if (st_compare_pad (&sp[i].c[1], sp[i].c[0],
467 &sp[max_idx].c[1], sp[max_idx].c[0]) > 0)
469 sp[0].c = sp[max_idx].c;
479 moments_of_values (sp, n_args,
480 &weight, &mean, NULL, NULL, NULL);
481 sp->f = weight < *op++ ? SYSMIS : mean;
488 double min = DBL_MAX;
491 for (i = 0; i < n_args; i++)
492 if (sp[i].f != SYSMIS)
510 for (i = 1; i < n_args; i++)
511 if (st_compare_pad (&sp[i].c[1], sp[i].c[0],
512 &sp[min_idx].c[1], sp[min_idx].c[0]) < 0)
514 sp[0].c = sp[min_idx].c;
523 for (i = 0; i < n_args; i++)
524 if (sp[i].f == SYSMIS)
535 for (i = 0; i < n_args; i++)
536 if (sp[i].f != SYSMIS)
549 for (i = 1; i < n_args; i += 2)
550 if (sp[i].f == SYSMIS || sp[i + 1].f == SYSMIS)
552 else if (sp[0].f >= sp[i].f && sp[0].f <= sp[i + 1].f)
559 sp->f = sysmis ? SYSMIS : 0.0;
562 case OP_RANGE_STRING:
567 for (i = 1; i < n_args; i += 2)
568 if (st_compare_pad (&sp[0].c[1], sp[0].c[0],
569 &sp[i].c[1], sp[i].c[0]) >= 0
570 && st_compare_pad (&sp[0].c[1], sp[0].c[0],
571 &sp[i + 1].c[1], sp[i + 1].c[0]) <= 0)
582 double weight, variance;
585 moments_of_values (sp, n_args,
586 &weight, NULL, &variance, NULL, NULL);
587 sp->f = weight < *op++ ? SYSMIS : sqrt (variance);
597 for (i = 0; i < n_args; i++)
598 if (sp[i].f != SYSMIS)
612 double weight, variance;
615 moments_of_values (sp, n_args,
616 &weight, NULL, &variance, NULL, NULL);
617 sp->f = weight < *op++ ? SYSMIS : variance;
621 /* Time construction function. */
624 if (sp[0].f == SYSMIS || sp[1].f == SYSMIS || sp[2].f == SYSMIS)
629 min = min (sp[0].f, min (sp[1].f, sp[2].f));
630 max = max (sp[0].f, max (sp[1].f, sp[2].f));
631 if (min < 0. && max > 0.)
633 msg (SW, _("TIME.HMS cannot mix positive and negative "
634 "in its arguments."));
638 sp->f = 60. * (60. * sp[0].f + sp[1].f) + sp[2].f;
643 sp->f /= 60. * 60. * 24.;
649 case OP_CTIME_MINUTES:
655 sp->f *= 60. * 60. * 24.;
657 case OP_CTIME_SECONDS:
661 /* Date construction functions. */
664 sp->f = yrmoda (sp[2].f, sp[1].f, sp[0].f);
666 sp->f *= 60. * 60. * 24.;
670 sp->f = yrmoda (sp[2].f, sp[0].f, sp[1].f);
672 sp->f *= 60. * 60. * 24.;
676 sp->f = yrmoda (sp[1].f, sp[0].f, 1);
678 sp->f *= 60. * 60. * 24.;
682 if (sp[0].f == SYSMIS)
686 sp->f = yrmoda (sp[1].f, sp[0].f * 3 - 2, 1);
688 sp->f *= 60. * 60. * 24.;
693 if (sp[0].f == SYSMIS || sp[1].f == SYSMIS)
695 else if (sp[0].f < 0. || sp[0].f > 53.)
697 msg (SW, _("Week argument to WKYR must be in range 0 to 53."));
702 double result = yrmoda (sp[1].f, 1, 1);
703 if (result != SYSMIS)
704 result = 60. * 60. * 24. * (result + 7. * (sp->f - 1.));
710 if (sp[1].f == SYSMIS)
714 sp->f = yrmoda (sp[0].f, 1, 1);
716 sp->f = 60. * 60. * 24. * (sp->f + floor (sp[1].f) - 1);
721 sp->f = yrmoda (sp[0].f, sp[1].f, sp[2].f);
724 /* Date extraction functions. */
727 sp->f = floor (sp->f / 60. / 60. / 24.) * 60. * 60. * 24.;
731 sp->f = fmod (floor (sp->f / 60. / 60.), 24.);
735 sp->f = 86400. * julian_to_jday (sp->f / 86400.);
741 julian_to_calendar (sp->f / 86400., NULL, NULL, &day);
745 case OP_XDATE_MINUTE:
747 sp->f = fmod (floor (sp->f / 60.), 60.);
753 julian_to_calendar (sp->f / 86400., NULL, &month, NULL);
757 case OP_XDATE_QUARTER:
761 julian_to_calendar (sp->f / 86400., NULL, &month, NULL);
762 sp->f = (month - 1) / 3 + 1;
765 case OP_XDATE_SECOND:
767 sp->f = fmod (sp->f, 60.);
771 sp->f = floor (sp->f / 60. / 60. / 24.);
775 sp->f -= floor (sp->f / 60. / 60. / 24.) * 60. * 60. * 24.;
779 sp->f = (julian_to_jday (sp->f / 86400.) - 1) / 7 + 1;
783 sp->f = julian_to_wday (sp->f / 86400.);
789 julian_to_calendar (sp->f / 86400., &year, NULL, NULL);
794 /* String functions. */
800 dest = pool_alloc (e->pool, 256);
804 for (i = 0; i < n_args; i++)
807 if (sp[i].c[0] + dest[0] < 255)
809 memcpy (&dest[dest[0] + 1], &sp[i].c[1], sp[i].c[0]);
810 dest[0] += sp[i].c[0];
814 memcpy (&dest[dest[0] + 1], &sp[i].c[1], 255 - dest[0]);
828 int last = sp[0].c[0] - sp[1].c[0];
830 for (i = 0; i <= last; i++)
831 if (!memcmp (&sp[0].c[i + 1], &sp[1].c[1], sp[1].c[0]))
846 else if (sp[2].f == SYSMIS)
848 msg (SW, _("Argument 3 of RINDEX may not be system-missing."));
853 int part_len = sp[2].f;
855 if (part_len < 0 || part_len > sp[1].c[0]
856 || sp[1].c[0] % part_len != 0)
858 msg (SW, _("Argument 3 of RINDEX must be between 1 and "
859 "the length of argument 2, and it must "
860 "evenly divide the length of argument 2."));
866 int last = sp[0].c[0] - part_len;
867 for (i = 0; i <= last; i++)
868 for (j = 0; j < sp[1].c[0]; j += part_len)
869 if (!memcmp (&sp[0].c[i + 1], &sp[1].c[j + 1], part_len))
886 for (i = sp[0].c[0] - sp[1].c[0]; i >= 0; i--)
887 if (!memcmp (&sp[0].c[i + 1], &sp[1].c[1], sp[1].c[0]))
902 else if (sp[2].f == SYSMIS)
904 msg (SW, _("Argument 3 of RINDEX may not be system-missing."));
909 int part_len = sp[2].f;
911 if (part_len < 0 || part_len > sp[1].c[0]
912 || sp[1].c[0] % part_len != 0)
914 msg (SW, _("Argument 3 of RINDEX must be between 1 and "
915 "the length of argument 2, and it must "
916 "evenly divide the length of argument 2."));
921 for (i = sp[0].c[0] - part_len; i >= 0; i--)
922 for (j = 0; j < sp[1].c[0]; j += part_len)
923 if (!memcmp (&sp[0].c[i + 1], &sp[1].c[j + 1], part_len))
937 for (i = sp[0].c[0]; i >= 1; i--)
938 sp[0].c[i] = tolower ((unsigned char) (sp[0].c[i]));
941 for (i = sp[0].c[0]; i >= 1; i--)
942 sp[0].c[i] = toupper ((unsigned char) (sp[0].c[i]));
949 if (sp[1].f == SYSMIS || len < 0 || len > 255 || sp[2].c[0] != 1)
951 else if (len > sp[0].c[0])
955 dest = pool_alloc (e->pool, len + 1);
957 memset (&dest[1], sp[2].c[1], len - sp->c[0]);
958 memcpy (&dest[len - sp->c[0] + 1], &sp->c[1], sp->c[0]);
968 if (len < 0 || len > 255 || sp[2].c[0] != 1)
970 else if (len > sp[0].c[0])
974 dest = pool_alloc (e->pool, len + 1);
976 memcpy (&dest[1], &sp->c[1], sp->c[0]);
977 memset (&dest[sp->c[0] + 1], sp[2].c[1], len - sp->c[0]);
989 int len = sp[0].c[0];
990 int cmp = sp[1].c[1];
993 while (i <= len && sp[0].c[i] == cmp)
997 sp[0].c[i] = sp[0].c[0] - i;
1005 if (sp[1].c[0] != 1)
1009 int cmp = sp[1].c[1];
1010 while (sp[0].c[0] > 0 && sp[0].c[sp[0].c[0]] == cmp)
1022 di.format.type = *op++;
1023 di.format.w = *op++;
1024 di.format.d = *op++;
1025 di.e = &sp->c[1] + min (sp->c[0], di.format.w);
1033 unsigned char *dest;
1039 dest = pool_alloc (e->pool, f.w + 1);
1042 assert ((formats[f.type].cat & FCAT_STRING) == 0);
1043 data_out (&dest[1], &f, sp);
1053 if (index < 1 || index > sp[0].c[0])
1058 sp->c[index] = sp->c[0] - index;
1071 if (sp[1].f == SYSMIS || sp[2].f == SYSMIS || index < 1
1072 || index > sp[0].c[0] || n < 1)
1079 sp->c[index] = sp->c[0] - index;
1090 if (sp->f != SYSMIS)
1093 case OP_NUM_TO_BOOL:
1096 else if (sp->f == 1.0)
1098 else if (sp->f != SYSMIS)
1100 msg (SE, _("A number being treated as a Boolean in an "
1101 "expression was found to have a value other than "
1102 "0 (false), 1 (true), or the system-missing value. "
1103 "The result was forced to 0."));
1111 if (sp[0].f != SYSMIS)
1113 if (sp[1].f == SYSMIS)
1119 sp->f = fmod (sp[0].f, sp[1].f);
1123 if (sp->f != SYSMIS)
1124 sp->f *= rng_get_double_normal (pspp_rng ());
1127 if (sp->f != SYSMIS)
1128 sp->f *= rng_get_double (pspp_rng ());
1131 sp->f = sp->f == SYSMIS || !finite (sp->f);
1133 case OP_VEC_ELEM_NUM:
1135 int rindx = sp[0].f + EPSILON;
1136 const struct vector *v = dict_get_vector (default_dict, *op++);
1138 if (sp[0].f == SYSMIS || rindx < 1 || rindx > v->cnt)
1140 if (sp[0].f == SYSMIS)
1141 msg (SE, _("SYSMIS is not a valid index value for vector "
1142 "%s. The result will be set to SYSMIS."),
1145 msg (SE, _("%g is not a valid index value for vector %s. "
1146 "The result will be set to SYSMIS."),
1152 sp->f = case_num (c, v->var[rindx - 1]->fv);
1155 case OP_VEC_ELEM_STR:
1157 int rindx = sp[0].f + EPSILON;
1158 const struct vector *vect = dict_get_vector (default_dict, *op++);
1161 if (sp[0].f == SYSMIS || rindx < 1 || rindx > vect->cnt)
1163 if (sp[0].f == SYSMIS)
1164 msg (SE, _("SYSMIS is not a valid index value for vector "
1165 "%s. The result will be set to the empty "
1169 msg (SE, _("%g is not a valid index value for vector %s. "
1170 "The result will be set to the empty string."),
1171 sp[0].f, vect->name);
1172 sp->c = pool_alloc (e->pool, 1);
1177 v = vect->var[rindx - 1];
1178 sp->c = pool_alloc (e->pool, v->width + 1);
1179 sp->c[0] = v->width;
1181 memcpy (&sp->c[1], case_str (c, v->fv), v->width);
1192 sp->c = pool_alloc (e->pool, *str + 1);
1193 memcpy (sp->c, str, *str + 1);
1199 sp->f = case_num (c, (*vars)->fv);
1200 if (is_num_user_missing (sp->f, *vars))
1206 int width = (*vars)->width;
1209 sp->c = pool_alloc (e->pool, width + 1);
1212 memcpy (&sp->c[1], case_str (c, (*vars)->fv), width);
1218 struct ccase *c = lagged_case (*op++);
1225 sp->f = case_num (c, (*vars)->fv);
1226 if (is_num_user_missing (sp->f, *vars))
1234 struct ccase *c = lagged_case (*op++);
1235 int width = (*vars)->width;
1238 sp->c = pool_alloc (e->pool, width + 1);
1242 memset (sp->c, ' ', width);
1244 memcpy (&sp->c[1], case_str (c, (*vars)->fv), width);
1252 sp->f = case_num (c, *op++) == SYSMIS;
1257 sp->f = case_num (c, *op++);
1274 if (e->type != EXPR_STRING)
1276 double value = sp->f;
1277 if (!finite (value))