Fix missing @clicksequence problem with older Texinfo versions.
[pspp-builds.git] / expressions / helpers.c
1 #include <config.h>
2 #include "helpers.h"
3 #include <gsl/gsl_roots.h>
4 #include <gsl/gsl_sf.h>
5 #include <libpspp/assertion.h>
6 #include <libpspp/pool.h>
7 #include "private.h"
8
9 const struct substring empty_string = {NULL, 0};
10
11 static void
12 expr_error (void *aux UNUSED, const char *format, ...)
13 {
14   struct msg m;
15   va_list args;
16
17   m.category = MSG_SYNTAX;
18   m.severity = MSG_ERROR;
19   va_start (args, format);
20   m.text = xvasprintf (format, args);
21   va_end (args);
22
23   msg_emit (&m);
24 }
25
26 double
27 expr_ymd_to_ofs (double year, double month, double day)
28 {
29   int y = year;
30   int m = month;
31   int d = day;
32
33   if (y != year || m != month || d != day)
34     {
35       msg (SE, _("One of the arguments to a DATE function is not an integer.  "
36                  "The result will be system-missing."));
37       return SYSMIS;
38     }
39
40   return calendar_gregorian_to_offset (y, m, d, expr_error, NULL);
41 }
42
43 double
44 expr_ymd_to_date (double year, double month, double day)
45 {
46   double ofs = expr_ymd_to_ofs (year, month, day);
47   return ofs != SYSMIS ? ofs * DAY_S : SYSMIS;
48 }
49
50 double
51 expr_wkyr_to_date (double week, double year)
52 {
53   int w = week;
54
55   if (w != week)
56     {
57       msg (SE, _("The week argument to DATE.WKYR is not an integer.  "
58                  "The result will be system-missing."));
59       return SYSMIS;
60     }
61   else if (w < 1 || w > 53)
62     {
63       msg (SE, _("The week argument to DATE.WKYR is outside the acceptable "
64                  "range of 1 to 53.  "
65                  "The result will be system-missing."));
66       return SYSMIS;
67     }
68   else
69     {
70       double yr_1_1 = expr_ymd_to_ofs (year, 1, 1);
71       if (yr_1_1 != SYSMIS)
72         return DAY_S * (yr_1_1 + WEEK_DAY * (w - 1));
73       else
74         return SYSMIS;
75     }
76 }
77
78 double
79 expr_yrday_to_date (double year, double yday)
80 {
81   int yd = yday;
82
83   if (yd != yday)
84     {
85       msg (SE, _("The day argument to DATE.YRDAY is not an integer.  "
86                  "The result will be system-missing."));
87       return SYSMIS;
88     }
89   else if (yd < 1 || yd > 366)
90     {
91       msg (SE, _("The day argument to DATE.YRDAY is outside the acceptable "
92                  "range of 1 to 366.  "
93                  "The result will be system-missing."));
94       return SYSMIS;
95     }
96   else
97     {
98       double yr_1_1 = expr_ymd_to_ofs (year, 1, 1);
99       if (yr_1_1 != SYSMIS)
100         return DAY_S * (yr_1_1 + yd - 1.);
101       else
102         return SYSMIS;
103     }
104 }
105
106 double
107 expr_yrmoda (double year, double month, double day)
108 {
109   if (year >= 0 && year <= 99)
110     year += 1900;
111   else if (year != (int) year && year > 47516)
112     {
113       msg (SE, _("The year argument to YRMODA is greater than 47516.  "
114                  "The result will be system-missing."));
115       return SYSMIS;
116     }
117
118   return expr_ymd_to_ofs (year, month, day);
119 }
120 \f
121 /* A date unit. */
122 enum date_unit
123   {
124     DATE_YEARS,
125     DATE_QUARTERS,
126     DATE_MONTHS,
127     DATE_WEEKS,
128     DATE_DAYS,
129     DATE_HOURS,
130     DATE_MINUTES,
131     DATE_SECONDS
132   };
133
134 /* Stores in *UNIT the unit whose name is NAME.
135    Return success. */
136 static enum date_unit
137 recognize_unit (struct substring name, enum date_unit *unit)
138 {
139   struct unit_name
140     {
141       enum date_unit unit;
142       const struct substring name;
143     };
144   static const struct unit_name unit_names[] =
145     {
146       { DATE_YEARS, SS_LITERAL_INITIALIZER ("years") },
147       { DATE_QUARTERS, SS_LITERAL_INITIALIZER ("quarters") },
148       { DATE_MONTHS, SS_LITERAL_INITIALIZER ("months") },
149       { DATE_WEEKS, SS_LITERAL_INITIALIZER ("weeks") },
150       { DATE_DAYS, SS_LITERAL_INITIALIZER ("days") },
151       { DATE_HOURS, SS_LITERAL_INITIALIZER ("hours") },
152       { DATE_MINUTES, SS_LITERAL_INITIALIZER ("minutes") },
153       { DATE_SECONDS, SS_LITERAL_INITIALIZER ("seconds") },
154     };
155   const int unit_name_cnt = sizeof unit_names / sizeof *unit_names;
156
157   const struct unit_name *un;
158
159   for (un = unit_names; un < &unit_names[unit_name_cnt]; un++)
160     if (ss_equals_case (un->name, name))
161       {
162         *unit = un->unit;
163         return true;
164       }
165
166   msg (SE, _("Unrecognized date unit \"%.*s\".  "
167              "Valid date units are \"years\", \"quarters\", \"months\", "
168              "\"weeks\", \"days\", \"hours\", \"minutes\", and \"seconds\"."),
169        (int) ss_length (name), ss_data (name));
170   return false;
171 }
172
173 /* Returns the number of whole years from DATE1 to DATE2,
174    where a year is defined as the same or later month, day, and
175    time of day. */
176 static int
177 year_diff (double date1, double date2)
178 {
179   int y1, m1, d1, yd1;
180   int y2, m2, d2, yd2;
181   int diff;
182
183   assert (date2 >= date1);
184   calendar_offset_to_gregorian (date1 / DAY_S, &y1, &m1, &d1, &yd1);
185   calendar_offset_to_gregorian (date2 / DAY_S, &y2, &m2, &d2, &yd2);
186
187   diff = y2 - y1;
188   if (diff > 0)
189     {
190       int yd1 = 32 * m1 + d1;
191       int yd2 = 32 * m2 + d2;
192       if (yd2 < yd1
193           || (yd2 == yd1 && fmod (date2, DAY_S) < fmod (date1, DAY_S)))
194         diff--;
195     }
196   return diff;
197 }
198
199 /* Returns the number of whole months from DATE1 to DATE2,
200    where a month is defined as the same or later day and time of
201    day. */
202 static int
203 month_diff (double date1, double date2)
204 {
205   int y1, m1, d1, yd1;
206   int y2, m2, d2, yd2;
207   int diff;
208
209   assert (date2 >= date1);
210   calendar_offset_to_gregorian (date1 / DAY_S, &y1, &m1, &d1, &yd1);
211   calendar_offset_to_gregorian (date2 / DAY_S, &y2, &m2, &d2, &yd2);
212
213   diff = ((y2 * 12) + m2) - ((y1 * 12) + m1);
214   if (diff > 0
215       && (d2 < d1
216           || (d2 == d1 && fmod (date2, DAY_S) < fmod (date1, DAY_S))))
217     diff--;
218   return diff;
219 }
220
221 /* Returns the number of whole quarter from DATE1 to DATE2,
222    where a quarter is defined as three months. */
223 static int
224 quarter_diff (double date1, double date2)
225 {
226   return month_diff (date1, date2) / 3;
227 }
228
229 /* Returns the number of seconds in the given UNIT. */
230 static int
231 date_unit_duration (enum date_unit unit)
232 {
233   switch (unit)
234     {
235     case DATE_WEEKS:
236       return WEEK_S;
237
238     case DATE_DAYS:
239       return DAY_S;
240
241     case DATE_HOURS:
242       return H_S;
243
244     case DATE_MINUTES:
245       return MIN_S;
246
247     case DATE_SECONDS:
248       return 1;
249
250     default:
251       NOT_REACHED ();
252     }
253 }
254
255 /* Returns the span from DATE1 to DATE2 in terms of UNIT_NAME. */
256 double
257 expr_date_difference (double date1, double date2, struct substring unit_name)
258 {
259   enum date_unit unit;
260
261   if (!recognize_unit (unit_name, &unit))
262     return SYSMIS;
263
264   switch (unit)
265     {
266     case DATE_YEARS:
267       return (date2 >= date1
268               ? year_diff (date1, date2)
269               : -year_diff (date2, date1));
270
271     case DATE_QUARTERS:
272       return (date2 >= date1
273               ? quarter_diff (date1, date2)
274               : -quarter_diff (date2, date1));
275
276     case DATE_MONTHS:
277       return (date2 >= date1
278               ? month_diff (date1, date2)
279               : -month_diff (date2, date1));
280
281     case DATE_WEEKS:
282     case DATE_DAYS:
283     case DATE_HOURS:
284     case DATE_MINUTES:
285     case DATE_SECONDS:
286       return trunc ((date2 - date1) / date_unit_duration (unit));
287     }
288
289   NOT_REACHED ();
290 }
291
292 /* How to deal with days out of range for a given month. */
293 enum date_sum_method
294   {
295     SUM_ROLLOVER,       /* Roll them over to the next month. */
296     SUM_CLOSEST         /* Use the last day of the month. */
297   };
298
299 /* Stores in *METHOD the method whose name is NAME.
300    Return success. */
301 static bool
302 recognize_method (struct substring method_name, enum date_sum_method *method)
303 {
304   if (ss_equals_case (method_name, ss_cstr ("closest")))
305     {
306       *method = SUM_CLOSEST;
307       return true;
308     }
309   else if (ss_equals_case (method_name, ss_cstr ("rollover")))
310     {
311       *method = SUM_ROLLOVER;
312       return true;
313     }
314   else
315     {
316       msg (SE, _("Invalid DATESUM method.  "
317                  "Valid choices are \"closest\" and \"rollover\"."));
318       return false;
319     }
320 }
321
322 /* Returns DATE advanced by the given number of MONTHS, with
323    day-of-month overflow resolved using METHOD. */
324 static double
325 add_months (double date, int months, enum date_sum_method method)
326 {
327   int y, m, d, yd;
328   double output;
329
330   calendar_offset_to_gregorian (date / DAY_S, &y, &m, &d, &yd);
331   y += months / 12;
332   m += months % 12;
333   if (m < 1)
334     {
335       m += 12;
336       y--;
337     }
338   else if (m > 12)
339     {
340       m -= 12;
341       y++;
342     }
343   assert (m >= 1 && m <= 12);
344
345   if (method == SUM_CLOSEST && d > calendar_days_in_month (y, m))
346     d = calendar_days_in_month (y, m);
347
348   output = calendar_gregorian_to_offset (y, m, d, expr_error, NULL);
349   if (output != SYSMIS)
350     output = (output * DAY_S) + fmod (date, DAY_S);
351   return output;
352 }
353
354 /* Returns DATE advanced by the given QUANTITY of units given in
355    UNIT_NAME, with day-of-month overflow resolved using
356    METHOD_NAME. */
357 double
358 expr_date_sum (double date, double quantity, struct substring unit_name,
359                struct substring method_name)
360 {
361   enum date_unit unit;
362   enum date_sum_method method;
363
364   if (!recognize_unit (unit_name, &unit)
365       || !recognize_method (method_name, &method))
366     return SYSMIS;
367
368   switch (unit)
369     {
370     case DATE_YEARS:
371       return add_months (date, trunc (quantity) * 12, method);
372
373     case DATE_QUARTERS:
374       return add_months (date, trunc (quantity) * 3, method);
375
376     case DATE_MONTHS:
377       return add_months (date, trunc (quantity), method);
378
379     case DATE_WEEKS:
380     case DATE_DAYS:
381     case DATE_HOURS:
382     case DATE_MINUTES:
383     case DATE_SECONDS:
384       return date + quantity * date_unit_duration (unit);
385     }
386
387   NOT_REACHED ();
388 }
389
390 int
391 compare_string (const struct substring *a, const struct substring *b)
392 {
393   size_t i;
394
395   for (i = 0; i < a->length && i < b->length; i++)
396     if (a->string[i] != b->string[i])
397       return a->string[i] < b->string[i] ? -1 : 1;
398   for (; i < a->length; i++)
399     if (a->string[i] != ' ')
400       return 1;
401   for (; i < b->length; i++)
402     if (b->string[i] != ' ')
403       return -1;
404   return 0;
405 }
406
407 size_t
408 count_valid (double *d, size_t d_cnt)
409 {
410   size_t valid_cnt;
411   size_t i;
412
413   valid_cnt = 0;
414   for (i = 0; i < d_cnt; i++)
415     valid_cnt += is_valid (d[i]);
416   return valid_cnt;
417 }
418
419 struct substring
420 alloc_string (struct expression *e, size_t length)
421 {
422   struct substring s;
423   s.length = length;
424   s.string = pool_alloc (e->eval_pool, length);
425   return s;
426 }
427
428 struct substring
429 copy_string (struct expression *e, const char *old, size_t length)
430 {
431   struct substring s = alloc_string (e, length);
432   memcpy (s.string, old, length);
433   return s;
434 }
435
436 /* Returns the noncentral beta cumulative distribution function
437    value for the given arguments.
438
439    FIXME: The accuracy of this function is not entirely
440    satisfactory.  We only match the example values given in AS
441    310 to the first 5 significant digits. */
442 double
443 ncdf_beta (double x, double a, double b, double lambda)
444 {
445   double c;
446
447   if (x <= 0. || x >= 1. || a <= 0. || b <= 0. || lambda <= 0.)
448     return SYSMIS;
449
450   c = lambda / 2.;
451   if (lambda < 54.)
452     {
453       /* Algorithm AS 226. */
454       double x0, a0, beta, temp, gx, q, ax, sumq, sum;
455       double err_max = 2 * DBL_EPSILON;
456       double err_bound;
457       int iter_max = 100;
458       int iter;
459
460       x0 = floor (c - 5.0 * sqrt (c));
461       if (x0 < 0.)
462         x0 = 0.;
463       a0 = a + x0;
464       beta = (gsl_sf_lngamma (a0)
465               + gsl_sf_lngamma (b)
466               - gsl_sf_lngamma (a0 + b));
467       temp = gsl_sf_beta_inc (a0, b, x);
468       gx = exp (a0 * log (x) + b * log (1. - x) - beta - log (a0));
469       if (a0 >= a)
470         q = exp (-c + x0 * log (c)) - gsl_sf_lngamma (x0 + 1.);
471       else
472         q = exp (-c);
473       ax = q * temp;
474       sumq = 1. - q;
475       sum = ax;
476
477       iter = 0;
478       do
479         {
480           iter++;
481           temp -= gx;
482           gx = x * (a + b + iter - 1.) * gx / (a + iter);
483           q *= c / iter;
484           sumq -= q;
485           ax = temp * q;
486           sum += ax;
487
488           err_bound = (temp - gx) * sumq;
489         }
490       while (iter < iter_max && err_bound > err_max);
491
492       return sum;
493     }
494   else
495     {
496       /* Algorithm AS 310. */
497       double m, m_sqrt;
498       int iter, iter_lower, iter_upper, iter1, iter2, j;
499       double t, q, r, psum, beta, s1, gx, fx, temp, ftemp, t0, s0, sum, s;
500       double err_bound;
501       double err_max = 2 * DBL_EPSILON;
502
503       iter = 0;
504
505       m = floor (c + .5);
506       m_sqrt = sqrt (m);
507       iter_lower = m - 5. * m_sqrt;
508       iter_upper = m + 5. * m_sqrt;
509
510       t = -c + m * log (c) - gsl_sf_lngamma (m + 1.);
511       q = exp (t);
512       r = q;
513       psum = q;
514       beta = (gsl_sf_lngamma (a + m)
515               + gsl_sf_lngamma (b)
516               - gsl_sf_lngamma (a + m + b));
517       s1 = (a + m) * log (x) + b * log (1. - x) - log (a + m) - beta;
518       fx = gx = exp (s1);
519       ftemp = temp = gsl_sf_beta_inc (a + m, b, x);
520       iter++;
521       sum = q * temp;
522       iter1 = m;
523
524       while (iter1 >= iter_lower && q >= err_max)
525         {
526           q = q * iter1 / c;
527           iter++;
528           gx = (a + iter1) / (x * (a + b + iter1 - 1.)) * gx;
529           iter1--;
530           temp += gx;
531           psum += q;
532           sum += q * temp;
533         }
534
535       t0 = (gsl_sf_lngamma (a + b)
536             - gsl_sf_lngamma (a + 1.)
537             - gsl_sf_lngamma (b));
538       s0 = a * log (x) + b * log (1. - x);
539
540       s = 0.;
541       for (j = 0; j < iter1; j++)
542         {
543           double t1;
544           s += exp (t0 + s0 + j * log (x));
545           t1 = log (a + b + j) - log (a + 1. + j) + t0;
546           t0 = t1;
547         }
548
549       err_bound = (1. - gsl_sf_gamma_inc_P (iter1, c)) * (temp + s);
550       q = r;
551       temp = ftemp;
552       gx = fx;
553       iter2 = m;
554       for (;;)
555         {
556           double ebd = err_bound + (1. - psum) * temp;
557           if (ebd < err_max || iter >= iter_upper)
558             break;
559
560           iter2++;
561           iter++;
562           q = q * c / iter2;
563           psum += q;
564           temp -= gx;
565           gx = x * (a + b + iter2 - 1.) / (a + iter2) * gx;
566           sum += q * temp;
567         }
568
569       return sum;
570     }
571 }
572
573 double
574 cdf_bvnor (double x0, double x1, double r)
575 {
576   double z = x0 * x0 - 2. * r * x0 * x1 + x1 * x1;
577   return exp (-z / (2. * (1 - r * r))) * (2. * M_PI * sqrt (1 - r * r));
578 }
579
580 double
581 idf_fdist (double P, double df1, double df2)
582 {
583   double temp = gslextras_cdf_beta_Pinv (P, df1 / 2, df2 / 2);
584   return temp * df2 / ((1. - temp) * df1);
585 }
586
587 /*
588  *  Mathlib : A C Library of Special Functions
589  *  Copyright (C) 1998 Ross Ihaka
590  *  Copyright (C) 2000 The R Development Core Team
591  *
592  *  This program is free software; you can redistribute it and/or
593  *  modify
594  *  it under the terms of the GNU General Public License as
595  *  published by
596  *  the Free Software Foundation; either version 2 of the
597  *  License, or
598  *  (at your option) any later version.
599  *
600  *  This program is distributed in the hope that it will be
601  *  useful,
602  *  but WITHOUT ANY WARRANTY; without even the implied warranty
603  *  of
604  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
605  *  GNU General Public License for more details.
606  *
607  *  You should have received a copy of the GNU General Public
608  *  License
609  *  along with this program; if not, write to the Free Software
610  *  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
611  *  02110-1301 USA.
612  */
613
614 /* Returns the density of the noncentral beta distribution with
615    noncentrality parameter LAMBDA. */
616 double
617 npdf_beta (double x, double a, double b, double lambda)
618 {
619   if (lambda < 0. || a <= 0. || b <= 0.)
620     return SYSMIS;
621   else if (lambda == 0.)
622     return gsl_ran_beta_pdf (x, a, b);
623   else
624     {
625       double max_error = 2 * DBL_EPSILON;
626       int max_iter = 200;
627       double term = gsl_ran_beta_pdf (x, a, b);
628       double lambda2 = 0.5 * lambda;
629       double weight = exp (-lambda2);
630       double sum = weight * term;
631       double psum = weight;
632       int k;
633       for (k = 1; k <= max_iter && 1 - psum < max_error; k++) {
634         weight *= lambda2 / k;
635         term *= x * (a + b) / a;
636         sum += weight * term;
637         psum += weight;
638         a += 1;
639       }
640       return sum;
641     }
642 }