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, 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;
75 for (i = 1; i < *op; i++)
77 if (sp[i].f == SYSMIS)
90 for (i = 1; i < *op; i++)
92 if (sp[i].f == SYSMIS)
104 if (sp[0].f == SYSMIS)
109 else if (sp[1].f == SYSMIS)
117 else if (sp[0].f == 0.0 && sp[1].f == 0.0)
120 sp->f = pow (sp[0].f, sp[1].f);
124 /* Note that booleans are always one of 0, 1, or SYSMIS.
126 Truth table (in order of detection):
138 1 and SYSMIS = SYSMIS
139 SYSMIS and SYSMIS = SYSMIS
143 SYSMIS and 1 = SYSMIS
147 if (sp[0].f == 0.0); /* 1 */
148 else if (sp[1].f == 0.0)
150 else if (sp[1].f == SYSMIS)
151 sp->f = SYSMIS; /* 3 */
154 /* Truth table (in order of detection):
167 SYSMIS or SYSMIS = SYSMIS
175 if (sp[0].f == 1.0); /* 1 */
176 else if (sp[1].f == 1.0)
178 else if (sp[1].f == SYSMIS)
179 sp->f = SYSMIS; /* 3 */
184 else if (sp[0].f == 1.0)
190 if (sp[0].f != SYSMIS)
192 if (sp[1].f == SYSMIS)
195 sp->f = sp[0].f == sp[1].f;
200 if (sp[0].f != SYSMIS)
202 if (sp[1].f == SYSMIS)
205 sp->f = sp[0].f >= sp[1].f;
210 if (sp[0].f != SYSMIS)
212 if (sp[1].f == SYSMIS)
215 sp->f = sp[0].f > sp[1].f;
220 if (sp[0].f != SYSMIS)
222 if (sp[1].f == SYSMIS)
225 sp->f = sp[0].f <= sp[1].f;
230 if (sp[0].f != SYSMIS)
232 if (sp[1].f == SYSMIS)
235 sp->f = sp[0].f < sp[1].f;
240 if (sp[0].f != SYSMIS)
242 if (sp[1].f == SYSMIS)
245 sp->f = sp[0].f != sp[1].f;
249 /* String operators. */
252 sp[0].f = st_compare_pad (&sp[0].c[1], sp[0].c[0],
253 &sp[1].c[1], sp[1].c[0]) == 0;
257 sp[0].f = st_compare_pad (&sp[0].c[1], sp[0].c[0],
258 &sp[1].c[1], sp[1].c[0]) >= 0;
262 sp[0].f = st_compare_pad (&sp[0].c[1], sp[0].c[0],
263 &sp[1].c[1], sp[1].c[0]) > 0;
267 sp[0].f = st_compare_pad (&sp[0].c[1], sp[0].c[0],
268 &sp[1].c[1], sp[1].c[0]) <= 0;
272 sp[0].f = st_compare_pad (&sp[0].c[1], sp[0].c[0],
273 &sp[1].c[1], sp[1].c[0]) < 0;
277 sp[0].f = st_compare_pad (&sp[0].c[1], sp[0].c[0],
278 &sp[1].c[1], sp[1].c[0]) != 0;
281 /* Unary functions. */
288 sp->f = fabs (sp->f);
294 sp->f = acos (sp->f);
303 sp->f = asin (sp->f);
310 sp->f = atan (sp->f);
329 sp->f = log10 (sp->f);
338 sp->f = log10 (sp->f);
345 sp->f = fmod (sp->f, 10);
351 sp->f = floor (sp->f + 0.5);
353 sp->f = -floor (-sp->f + 0.5);
364 sp->f = sqrt (sp->f);
382 sp->f = floor (sp->f);
384 sp->f = -floor (-sp->f);
388 /* N-ary numeric functions. */
397 for (i = 1; i <= n_args; i++)
398 if (sp[0].f == sp[i].f)
403 else if (sp[i].f != SYSMIS)
405 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]))
431 for (i = 0; i < n_args; i++)
432 if (sp[i].f != SYSMIS)
436 sum[1] += sp[i].f * sp[i].f;
441 sp->f = calc_cfvar (sum, nv);
448 double max = -DBL_MAX;
451 for (i = 0; i < n_args; i++)
452 if (sp[i].f != SYSMIS)
472 for (i = 0; i < n_args; i++)
473 if (sp[i].f != SYSMIS)
481 sp->f = calc_mean (sum, nv);
488 double min = DBL_MAX;
491 for (i = 0; i < n_args; i++)
492 if (sp[i].f != SYSMIS)
510 for (i = 0; i < n_args; i++)
511 if (sp[i].f == SYSMIS)
522 for (i = 0; i < n_args; i++)
523 if (sp[i].f != SYSMIS)
536 for (i = 1; i <= n_args; i += 2)
537 if (sp[i].f == SYSMIS || sp[i + 1].f == SYSMIS)
539 else if (sp[0].f >= sp[i].f && sp[0].f <= sp[i + 1].f)
546 sp->f = sysmis ? SYSMIS : 0.0;
549 case OP_RANGE_STRING:
554 for (i = 1; i <= n_args; i += 2)
555 if (st_compare_pad (&sp[0].c[1], sp[0].c[0],
556 &sp[i].c[1], sp[i].c[0]) >= 0
557 && st_compare_pad (&sp[0].c[1], sp[0].c[0],
558 &sp[i + 1].c[1], sp[i + 1].c[0]) <= 0)
572 sum[0] = sum[1] = 0.0;
575 for (i = 0; i < n_args; i++)
576 if (sp[i].f != SYSMIS)
580 sum[1] += sp[i].f * sp[i].f;
585 sp->f = calc_stddev (calc_variance (sum, nv));
595 for (i = 0; i < n_args; i++)
596 if (sp[i].f != SYSMIS)
613 sum[0] = sum[1] = 0.0;
616 for (i = 0; i < n_args; i++)
617 if (sp[i].f != SYSMIS)
621 sum[1] += sp[i].f * sp[i].f;
626 sp->f = calc_variance (sum, nv);
630 /* Time construction function. */
633 if (sp[0].f == SYSMIS || sp[1].f == SYSMIS || sp[2].f == SYSMIS)
636 sp->f = 60. * (60. * sp[0].f + sp[1].f) + sp[2].f;
639 /* Date construction functions. */
642 sp->f = yrmoda (sp[2].f, sp[1].f, sp[0].f);
644 sp->f *= 60. * 60. * 24.;
648 sp->f = yrmoda (sp[2].f, sp[0].f, sp[1].f);
650 sp->f *= 60. * 60. * 24.;
654 sp->f = yrmoda (sp[1].f, sp[0].f, 1);
656 sp->f *= 60. * 60. * 24.;
660 if (sp[0].f == SYSMIS)
664 sp->f = yrmoda (sp[1].f, sp[0].f * 3 - 2, 1);
666 sp->f *= 60. * 60. * 24.;
671 if (sp[0].f == SYSMIS)
675 sp[1].f = yrmoda (sp[1].f, 1, 1);
677 sp[1].f = 60. * 60. * 24. * (sp[1].f + 7. * (floor (sp[0].f) - 1.));
683 if (sp[1].f == SYSMIS)
687 sp->f = yrmoda (sp[0].f, 1, 1);
689 sp->f = 60. * 60. * 24. * (sp->f + floor (sp[1].f) - 1);
694 sp->f = yrmoda (sp[0].f, sp[1].f, sp[2].f);
697 /* Date extraction functions. */
700 sp->f = floor (sp->f / 60. / 60. / 24.) * 60. * 60. * 24.;
704 sp->f = fmod (floor (sp->f / 60. / 60.), 24.);
708 sp->f = 86400. * julian_to_jday (sp->f / 86400.);
714 julian_to_calendar (sp->f / 86400., NULL, NULL, &day);
718 case OP_XDATE_MINUTE:
720 sp->f = fmod (floor (sp->f / 60.), 60.);
726 julian_to_calendar (sp->f / 86400., NULL, &month, NULL);
730 case OP_XDATE_QUARTER:
734 julian_to_calendar (sp->f / 86400., NULL, &month, NULL);
735 sp->f = (month - 1) / 3 + 1;
738 case OP_XDATE_SECOND:
740 sp->f = fmod (sp->f, 60.);
744 sp->f = floor (sp->f / 60. / 60. / 24.);
748 sp->f -= floor (sp->f / 60. / 60. / 24.) * 60. * 60. * 24.;
752 sp->f = (julian_to_jday (sp->f / 86400.) - 1) / 7 + 1;
756 sp->f = julian_to_wday (sp->f / 86400.);
762 julian_to_calendar (sp->f / 86400., &year, NULL, NULL);
767 /* String functions. */
773 dest = pool_alloc (e->pool, 256);
777 for (i = 0; i < n_args; i++)
780 if (sp[i].c[0] + dest[0] < 255)
782 memcpy (&dest[dest[0] + 1], &sp[i].c[1], sp[i].c[0]);
783 dest[0] += sp[i].c[0];
787 memcpy (&dest[dest[0] + 1], &sp[i].c[1], 255 - dest[0]);
801 int last = sp[0].c[0] - sp[1].c[0];
802 for (i = 0; i <= last; i++)
803 if (!memcmp (&sp[0].c[i + 1], &sp[0].c[1], sp[0].c[0]))
813 /* Length of each search string. */
814 int part_len = sp[2].f;
817 if (sp[1].c[0] == 0 || part_len <= 0 || sp[2].f == SYSMIS
818 || sp[1].c[0] % part_len != 0)
822 /* Last possible index. */
823 int last = sp[0].c[0] - part_len;
825 for (i = 0; i <= last; i++)
826 for (j = 0; j < sp[1].c[0]; j += part_len)
827 if (!memcmp (&sp[0].c[i], &sp[1].c[j], part_len))
842 for (i = sp[0].c[0] - sp[1].c[0]; i >= 0; i--)
843 if (!memcmp (&sp[0].c[i + 1], &sp[0].c[1], sp[0].c[0]))
853 /* Length of each search string. */
854 int part_len = sp[2].f;
857 if (sp[1].c[0] == 0 || part_len <= 0 || sp[2].f == SYSMIS
858 || sp[1].c[0] % part_len != 0)
862 for (i = sp[0].c[0] - part_len; i >= 0; i--)
863 for (j = 0; j < sp[1].c[0]; j += part_len)
864 if (!memcmp (&sp[0].c[i], &sp[1].c[j], part_len))
877 for (i = sp[0].c[0]; i >= 1; i--)
878 sp[0].c[i] = tolower ((unsigned char) (sp[0].c[i]));
881 for (i = sp[0].c[0]; i >= 1; i--)
882 sp[0].c[i] = toupper ((unsigned char) (sp[0].c[i]));
889 if (sp[1].f == SYSMIS || len < 0 || len > 255)
891 else if (len > sp[0].c[0])
895 dest = pool_alloc (e->pool, len + 1);
897 memset (&dest[1], ' ', len - sp->c[0]);
898 memcpy (&dest[len - sp->c[0] + 1], &sp->c[1], sp->c[0]);
908 if (sp[1].f == SYSMIS || len < 0 || len > 255 || sp[2].c[0] != 1)
910 else if (len > sp[0].c[0])
914 dest = pool_alloc (e->pool, len + 1);
916 memset (&dest[1], sp[2].c[1], len - sp->c[0]);
917 memcpy (&dest[len - sp->c[0] + 1], &sp->c[1], sp->c[0]);
927 if (sp[1].f == SYSMIS || len < 0 || len > 255)
929 else if (len > sp[0].c[0])
933 dest = pool_alloc (e->pool, len + 1);
935 memcpy (&dest[1], &sp->c[1], sp->c[0]);
936 memset (&dest[sp->c[0] + 1], ' ', len - sp->c[0]);
946 if (len < 0 || len > 255 || sp[2].c[0] != 1)
948 else if (len > sp[0].c[0])
952 dest = pool_alloc (e->pool, len + 1);
954 memcpy (&dest[1], &sp->c[1], sp->c[0]);
955 memset (&dest[sp->c[0] + 1], sp[2].c[1], len - sp->c[0]);
962 int len = sp[0].c[0];
965 while (i <= len && sp[0].c[i] == ' ')
969 sp[0].c[i] = sp[0].c[0] - i;
981 int len = sp[0].c[0];
982 int cmp = sp[1].c[1];
985 while (i <= len && sp[0].c[i] == cmp)
989 sp[0].c[i] = sp[0].c[0] - i;
997 while (sp[0].c[sp[0].c[0]] == ' ')
1002 if (sp[1].c[0] != 1)
1006 /* Note that NULs are not allowed in strings. This code
1007 needs to change if this decision is changed. */
1008 int cmp = sp[1].c[1];
1009 while (sp[0].c[sp[0].c[0]] == cmp)
1018 di.e = &sp->c[1] + sp->c[0];
1020 di.flags = DI_IGNORE_ERROR;
1022 di.format.type = FMT_F;
1023 di.format.w = sp->c[0];
1032 di.e = &sp->c[1] + sp->c[0];
1034 di.flags = DI_IGNORE_ERROR;
1036 di.format.type = *op++;
1037 di.format.w = *op++;
1038 di.format.d = *op++;
1045 unsigned char *dest;
1051 dest = pool_alloc (e->pool, f.w + 1);
1054 assert ((formats[f.type].cat & FCAT_STRING) == 0);
1055 data_out (&dest[1], &f, sp);
1065 if (index < 1 || index > sp[0].c[0])
1070 sp->c[index] = sp->c[0] - index;
1083 if (sp[1].f == SYSMIS || sp[2].f == SYSMIS || index < 1
1084 || index > sp[0].c[0] || n < 1)
1091 sp->c[index] = sp->c[0] - index;
1102 if (sp->f != SYSMIS)
1106 if (sp->f != SYSMIS)
1109 case OP_NUM_TO_BOOL:
1112 else if (sp->f == 1.0)
1114 else if (sp->f != SYSMIS)
1116 msg (SE, _("A number being treated as a Boolean in an "
1117 "expression was found to have a value other than "
1118 "0 (false), 1 (true), or the system-missing value. "
1119 "The result was forced to 0."));
1127 if (sp[0].f != SYSMIS)
1129 if (sp[1].f == SYSMIS)
1135 sp->f = fmod (sp[0].f, sp[1].f);
1139 if (sp->f != SYSMIS)
1140 sp->f *= rng_get_double_normal (pspp_rng ());
1143 if (sp->f != SYSMIS)
1144 sp->f *= rng_get_double (pspp_rng ());
1147 if (sp[0].f == SYSMIS || !finite (sp[0].f))
1152 case OP_VEC_ELEM_NUM:
1154 int rindx = sp[0].f + EPSILON;
1155 const struct vector *v = dict_get_vector (default_dict, *op++);
1157 if (sp[0].f == SYSMIS || rindx < 1 || rindx > v->cnt)
1159 if (sp[0].f == SYSMIS)
1160 msg (SE, _("SYSMIS is not a valid index value for vector "
1161 "%s. The result will be set to SYSMIS."),
1164 msg (SE, _("%g is not a valid index value for vector %s. "
1165 "The result will be set to SYSMIS."),
1171 sp->f = c->data[v->var[rindx - 1]->fv].f;
1174 case OP_VEC_ELEM_STR:
1176 int rindx = sp[0].f + EPSILON;
1177 const struct vector *vect = dict_get_vector (default_dict, *op++);
1180 if (sp[0].f == SYSMIS || rindx < 1 || rindx > vect->cnt)
1182 if (sp[0].f == SYSMIS)
1183 msg (SE, _("SYSMIS is not a valid index value for vector "
1184 "%s. The result will be set to the empty "
1188 msg (SE, _("%g is not a valid index value for vector %s. "
1189 "The result will be set to the empty string."),
1190 sp[0].f, vect->name);
1191 sp->c = pool_alloc (e->pool, 1);
1196 v = vect->var[rindx - 1];
1197 sp->c = pool_alloc (e->pool, v->width + 1);
1198 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);
1218 sp->f = c->data[(*vars)->fv].f;
1219 if (is_num_user_missing (sp->f, *vars))
1225 int width = (*vars)->width;
1228 sp->c = pool_alloc (e->pool, width + 1);
1231 memcpy (&sp->c[1], &c->data[(*vars)->fv], width);
1237 struct ccase *c = lagged_case (*op++);
1244 sp->f = c->data[(*vars)->fv].f;
1245 if (is_num_user_missing (sp->f, *vars))
1253 struct ccase *c = lagged_case (*op++);
1254 int width = (*vars)->width;
1257 sp->c = pool_alloc (e->pool, width + 1);
1261 memset (sp->c, ' ', width);
1263 memcpy (&sp->c[1], &c->data[(*vars)->fv], width);
1271 sp->f = c->data[*op++].f == SYSMIS;
1276 sp->f = is_str_user_missing (c->data[(*vars)->fv].s, *vars);
1282 sp->f = c->data[*op++].f;
1299 if (e->type != EX_STRING)
1301 double value = sp->f;
1302 if (!finite (value))