1 /* PSPP - a program for statistical analysis.
2 Copyright (C) 2008, 2010, 2011, 2015, 2016 Free Software Foundation, Inc.
4 This program is free software: you can redistribute it and/or modify
5 it under the terms of the GNU General Public License as published by
6 the Free Software Foundation, either version 3 of the License, or
7 (at your option) any later version.
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
14 You should have received a copy of the GNU General Public License
15 along with this program. If not, see <http://www.gnu.org/licenses/>. */
19 #include "language/expressions/helpers.h"
21 #include <gsl/gsl_roots.h>
22 #include <gsl/gsl_sf.h>
24 #include "language/expressions/private.h"
25 #include "libpspp/assertion.h"
26 #include "libpspp/pool.h"
28 #include "gl/minmax.h"
30 const struct substring empty_string = {NULL, 0};
33 expr_ymd_to_ofs (double year, double month, double day)
41 if (y != year || m != month || d != day)
43 msg (SE, _("One of the arguments to a DATE function is not an integer. "
44 "The result will be system-missing."));
48 ofs = calendar_gregorian_to_offset (y, m, d, &error);
51 msg (SE, "%s", error);
58 expr_ymd_to_date (double year, double month, double day)
60 double ofs = expr_ymd_to_ofs (year, month, day);
61 return ofs != SYSMIS ? ofs * DAY_S : SYSMIS;
65 expr_wkyr_to_date (double week, double year)
71 msg (SE, _("The week argument to DATE.WKYR is not an integer. "
72 "The result will be system-missing."));
75 else if (w < 1 || w > 53)
77 msg (SE, _("The week argument to DATE.WKYR is outside the acceptable "
79 "The result will be system-missing."));
84 double yr_1_1 = expr_ymd_to_ofs (year, 1, 1);
86 return DAY_S * (yr_1_1 + WEEK_DAY * (w - 1));
93 expr_yrday_to_date (double year, double yday)
99 msg (SE, _("The day argument to DATE.YRDAY is not an integer. "
100 "The result will be system-missing."));
103 else if (yd < 1 || yd > 366)
105 msg (SE, _("The day argument to DATE.YRDAY is outside the acceptable "
106 "range of 1 to 366. "
107 "The result will be system-missing."));
112 double yr_1_1 = expr_ymd_to_ofs (year, 1, 1);
113 if (yr_1_1 != SYSMIS)
114 return DAY_S * (yr_1_1 + yd - 1.);
121 expr_yrmoda (double year, double month, double day)
123 if (year >= 0 && year <= 99)
125 else if (year != (int) year && year > 47516)
127 msg (SE, _("The year argument to YRMODA is greater than 47516. "
128 "The result will be system-missing."));
132 return expr_ymd_to_ofs (year, month, day);
148 /* Stores in *UNIT the unit whose name is NAME.
150 static enum date_unit
151 recognize_unit (struct substring name, enum date_unit *unit)
156 const struct substring name;
158 static const struct unit_name unit_names[] =
160 { DATE_YEARS, SS_LITERAL_INITIALIZER ("years") },
161 { DATE_QUARTERS, SS_LITERAL_INITIALIZER ("quarters") },
162 { DATE_MONTHS, SS_LITERAL_INITIALIZER ("months") },
163 { DATE_WEEKS, SS_LITERAL_INITIALIZER ("weeks") },
164 { DATE_DAYS, SS_LITERAL_INITIALIZER ("days") },
165 { DATE_HOURS, SS_LITERAL_INITIALIZER ("hours") },
166 { DATE_MINUTES, SS_LITERAL_INITIALIZER ("minutes") },
167 { DATE_SECONDS, SS_LITERAL_INITIALIZER ("seconds") },
169 const int unit_name_cnt = sizeof unit_names / sizeof *unit_names;
171 const struct unit_name *un;
173 for (un = unit_names; un < &unit_names[unit_name_cnt]; un++)
174 if (ss_equals_case (un->name, name))
180 msg (SE, _("Unrecognized date unit `%.*s'. "
181 "Valid date units are `%s', `%s', `%s', "
182 "`%s', `%s', `%s', `%s', and `%s'."),
183 (int) ss_length (name), ss_data (name),
184 "years", "quarters", "months",
185 "weeks", "days", "hours", "minutes", "seconds");
190 /* Returns the number of whole years from DATE1 to DATE2,
191 where a year is defined as the same or later month, day, and
194 year_diff (double date1, double date2)
200 assert (date2 >= date1);
201 calendar_offset_to_gregorian (date1 / DAY_S, &y1, &m1, &d1, &yd1);
202 calendar_offset_to_gregorian (date2 / DAY_S, &y2, &m2, &d2, &yd2);
207 int yd1 = 32 * m1 + d1;
208 int yd2 = 32 * m2 + d2;
210 || (yd2 == yd1 && fmod (date2, DAY_S) < fmod (date1, DAY_S)))
216 /* Returns the number of whole months from DATE1 to DATE2,
217 where a month is defined as the same or later day and time of
220 month_diff (double date1, double date2)
226 assert (date2 >= date1);
227 calendar_offset_to_gregorian (date1 / DAY_S, &y1, &m1, &d1, &yd1);
228 calendar_offset_to_gregorian (date2 / DAY_S, &y2, &m2, &d2, &yd2);
230 diff = ((y2 * 12) + m2) - ((y1 * 12) + m1);
233 || (d2 == d1 && fmod (date2, DAY_S) < fmod (date1, DAY_S))))
238 /* Returns the number of whole quarter from DATE1 to DATE2,
239 where a quarter is defined as three months. */
241 quarter_diff (double date1, double date2)
243 return month_diff (date1, date2) / 3;
246 /* Returns the number of seconds in the given UNIT. */
248 date_unit_duration (enum date_unit unit)
272 /* Returns the span from DATE1 to DATE2 in terms of UNIT_NAME. */
274 expr_date_difference (double date1, double date2, struct substring unit_name)
278 if (!recognize_unit (unit_name, &unit))
284 return (date2 >= date1
285 ? year_diff (date1, date2)
286 : -year_diff (date2, date1));
289 return (date2 >= date1
290 ? quarter_diff (date1, date2)
291 : -quarter_diff (date2, date1));
294 return (date2 >= date1
295 ? month_diff (date1, date2)
296 : -month_diff (date2, date1));
303 return trunc ((date2 - date1) / date_unit_duration (unit));
309 /* How to deal with days out of range for a given month. */
312 SUM_ROLLOVER, /* Roll them over to the next month. */
313 SUM_CLOSEST /* Use the last day of the month. */
316 /* Stores in *METHOD the method whose name is NAME.
319 recognize_method (struct substring method_name, enum date_sum_method *method)
321 if (ss_equals_case (method_name, ss_cstr ("closest")))
323 *method = SUM_CLOSEST;
326 else if (ss_equals_case (method_name, ss_cstr ("rollover")))
328 *method = SUM_ROLLOVER;
333 msg (SE, _("Invalid DATESUM method. "
334 "Valid choices are `%s' and `%s'."), "closest", "rollover");
339 /* Returns DATE advanced by the given number of MONTHS, with
340 day-of-month overflow resolved using METHOD. */
342 add_months (double date, int months, enum date_sum_method method)
348 calendar_offset_to_gregorian (date / DAY_S, &y, &m, &d, &yd);
361 assert (m >= 1 && m <= 12);
363 if (method == SUM_CLOSEST && d > calendar_days_in_month (y, m))
364 d = calendar_days_in_month (y, m);
366 output = calendar_gregorian_to_offset (y, m, d, &error);
367 if (output != SYSMIS)
368 output = (output * DAY_S) + fmod (date, DAY_S);
371 msg (SE, "%s", error);
377 /* Returns DATE advanced by the given QUANTITY of units given in
378 UNIT_NAME, with day-of-month overflow resolved using
381 expr_date_sum (double date, double quantity, struct substring unit_name,
382 struct substring method_name)
385 enum date_sum_method method;
387 if (!recognize_unit (unit_name, &unit)
388 || !recognize_method (method_name, &method))
394 return add_months (date, trunc (quantity) * 12, method);
397 return add_months (date, trunc (quantity) * 3, method);
400 return add_months (date, trunc (quantity), method);
407 return date + quantity * date_unit_duration (unit);
414 compare_string_3way (const struct substring *a, const struct substring *b)
418 for (i = 0; i < a->length && i < b->length; i++)
419 if (a->string[i] != b->string[i])
420 return a->string[i] < b->string[i] ? -1 : 1;
421 for (; i < a->length; i++)
422 if (a->string[i] != ' ')
424 for (; i < b->length; i++)
425 if (b->string[i] != ' ')
431 count_valid (double *d, size_t d_cnt)
437 for (i = 0; i < d_cnt; i++)
438 valid_cnt += is_valid (d[i]);
443 alloc_string (struct expression *e, size_t length)
447 s.string = pool_alloc (e->eval_pool, length);
452 copy_string (struct expression *e, const char *old, size_t length)
454 struct substring s = alloc_string (e, length);
455 memcpy (s.string, old, length);
459 /* Returns the noncentral beta cumulative distribution function
460 value for the given arguments.
462 FIXME: The accuracy of this function is not entirely
463 satisfactory. We only match the example values given in AS
464 310 to the first 5 significant digits. */
466 ncdf_beta (double x, double a, double b, double lambda)
470 if (x <= 0. || x >= 1. || a <= 0. || b <= 0. || lambda <= 0.)
476 /* Algorithm AS 226. */
477 double x0, a0, beta, temp, gx, q, ax, sumq, sum;
478 double err_max = 2 * DBL_EPSILON;
483 x0 = floor (c - 5.0 * sqrt (c));
487 beta = (gsl_sf_lngamma (a0)
489 - gsl_sf_lngamma (a0 + b));
490 temp = gsl_sf_beta_inc (a0, b, x);
491 gx = exp (a0 * log (x) + b * log (1. - x) - beta - log (a0));
493 q = exp (-c + x0 * log (c)) - gsl_sf_lngamma (x0 + 1.);
505 gx = x * (a + b + iter - 1.) * gx / (a + iter);
511 err_bound = (temp - gx) * sumq;
513 while (iter < iter_max && err_bound > err_max);
519 /* Algorithm AS 310. */
521 int iter, iter_lower, iter_upper, iter1, iter2, j;
522 double t, q, r, psum, beta, s1, gx, fx, temp, ftemp, t0, s0, sum, s;
524 double err_max = 2 * DBL_EPSILON;
530 iter_lower = m - 5. * m_sqrt;
531 iter_upper = m + 5. * m_sqrt;
533 t = -c + m * log (c) - gsl_sf_lngamma (m + 1.);
537 beta = (gsl_sf_lngamma (a + m)
539 - gsl_sf_lngamma (a + m + b));
540 s1 = (a + m) * log (x) + b * log (1. - x) - log (a + m) - beta;
542 ftemp = temp = gsl_sf_beta_inc (a + m, b, x);
547 while (iter1 >= iter_lower && q >= err_max)
551 gx = (a + iter1) / (x * (a + b + iter1 - 1.)) * gx;
558 t0 = (gsl_sf_lngamma (a + b)
559 - gsl_sf_lngamma (a + 1.)
560 - gsl_sf_lngamma (b));
561 s0 = a * log (x) + b * log (1. - x);
564 for (j = 0; j < iter1; j++)
567 s += exp (t0 + s0 + j * log (x));
568 t1 = log (a + b + j) - log (a + 1. + j) + t0;
572 err_bound = (1. - gsl_sf_gamma_inc_P (iter1, c)) * (temp + s);
579 double ebd = err_bound + (1. - psum) * temp;
580 if (ebd < err_max || iter >= iter_upper)
588 gx = x * (a + b + iter2 - 1.) / (a + iter2) * gx;
597 cdf_bvnor (double x0, double x1, double r)
599 double z = pow2 (x0) - 2. * r * x0 * x1 + pow2 (x1);
600 return exp (-z / (2. * (1 - r * r))) * (2. * M_PI * sqrt (1 - r * r));
604 idf_fdist (double P, double df1, double df2)
606 double temp = gsl_cdf_beta_Pinv (P, df1 / 2, df2 / 2);
607 return temp * df2 / ((1. - temp) * df1);
611 * Mathlib : A C Library of Special Functions
612 * Copyright (C) 1998 Ross Ihaka
613 * Copyright (C) 2000 The R Development Core Team
615 * This program is free software; you can redistribute it and/or
617 * it under the terms of the GNU General Public License as
619 * the Free Software Foundation; either version 2 of the
621 * (at your option) any later version.
623 * This program is distributed in the hope that it will be
625 * but WITHOUT ANY WARRANTY; without even the implied warranty
627 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
628 * GNU General Public License for more details.
630 * You should have received a copy of the GNU General Public
632 * along with this program; if not, write to the Free Software
633 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
637 /* Returns the density of the noncentral beta distribution with
638 noncentrality parameter LAMBDA. */
640 npdf_beta (double x, double a, double b, double lambda)
642 if (lambda < 0. || a <= 0. || b <= 0.)
644 else if (lambda == 0.)
645 return gsl_ran_beta_pdf (x, a, b);
648 double max_error = 2 * DBL_EPSILON;
650 double term = gsl_ran_beta_pdf (x, a, b);
651 double lambda2 = 0.5 * lambda;
652 double weight = exp (-lambda2);
653 double sum = weight * term;
654 double psum = weight;
656 for (k = 1; k <= max_iter && 1 - psum < max_error; k++) {
657 weight *= lambda2 / k;
658 term *= x * (a + b) / a;
659 sum += weight * term;
668 round__ (double x, double mult, double fuzzbits, double adjustment)
671 fuzzbits = settings_get_fuzzbits ();
672 adjustment += exp2 (fuzzbits - DBL_MANT_DIG);
675 x = x >= 0. ? floor (x + adjustment) : -floor (-x + adjustment);
680 round_nearest (double x, double mult, double fuzzbits)
682 return round__ (x, mult, fuzzbits, .5);
686 round_zero (double x, double mult, double fuzzbits)
688 return round__ (x, mult, fuzzbits, 0);
692 replace_string (struct expression *e,
693 struct substring haystack,
694 struct substring needle,
695 struct substring replacement,
699 || haystack.length < needle.length
704 struct substring result = alloc_string (e, MAX_STRING);
708 while (i <= haystack.length - needle.length)
709 if (!memcmp (&haystack.string[i], needle.string, needle.length))
711 size_t copy_len = MIN (replacement.length, MAX_STRING - result.length);
712 memcpy (&result.string[result.length], replacement.string, copy_len);
713 result.length += copy_len;
721 if (result.length < MAX_STRING)
722 result.string[result.length++] = haystack.string[i];
725 while (i < haystack.length && result.length < MAX_STRING)
726 result.string[result.length++] = haystack.string[i++];
732 compare_doubles (const void *a_, const void *b_)
734 const double *ap = a_;
735 const double *bp = b_;
739 /* Sort SYSMIS to the end. */
747 median (double *a, size_t n)
749 /* Sort the array in-place, sorting SYSMIS to the end. */
750 qsort (a, n, sizeof *a, compare_doubles);
753 n = count_valid (a, n);
757 : (a[n / 2 - 1] + a[n / 2]) / 2.0);