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