pintos-os.org Git - pspp/blob - src/language/expressions/helpers.c

   1 /* PSPP - a program for statistical analysis.
   2    Copyright (C) 2008, 2010, 2011, 2015, 2016 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
  19 #include "language/expressions/helpers.h"
  20
  21 #include <gsl/gsl_roots.h>
  22 #include <gsl/gsl_sf.h>
  23
  24 #include "language/expressions/private.h"
  25 #include "libpspp/assertion.h"
  26 #include "libpspp/pool.h"
  27
  28 #include "gl/minmax.h"
  29
  30 const struct substring empty_string = {NULL, 0};
  31
  32 double
  33 expr_ymd_to_ofs (int y, int m, int d,
  34                  const struct expression *e, const struct expr_node *node,
  35                  int ya, int ma, int da)
  36 {
  37   int *error = calendar_gregorian_adjust (&y, &m, &d,
  38                                           settings_get_fmt_settings ());
  39   if (!error)
  40     return calendar_raw_gregorian_to_offset (y, m, d);
  41   else
  42     {
  43       msg_at (SE, expr_location (e, node),
  44               _("Invalid arguments to %s function."),
  45               operations[node->type].name);
  46
  47       if (error == &y && ya > 0)
  48         msg_at (SN, expr_location (e, y < 1582 ? node->args[ya - 1] : node),
  49                 _("Date %04d-%d-%d is before the earliest supported date "
  50                   "1582-10-15."), y, m, d);
  51       else if (error == &m && ma > 0)
  52         msg_at (SN, expr_location (e, node->args[ma - 1]),
  53                 _("Month %d is not in the acceptable range of 0 to 13."), m);
  54       else if (error == &d && da > 0)
  55         msg_at (SN, expr_location (e, node->args[da - 1]),
  56                 _("Day %d is not in the acceptable range of 0 to 31."), d);
  57       return SYSMIS;
  58     }
  59 }
  60
  61 double
  62 expr_ymd_to_date (int y, int m, int d,
  63                   const struct expression *e, const struct expr_node *n,
  64                   int ya, int ma, int da)
  65 {
  66   double ofs = expr_ymd_to_ofs (y, m, d, e, n, ya, ma, da);
  67   return ofs != SYSMIS ? ofs * DAY_S : SYSMIS;
  68 }
  69 \f
  70 /* A date unit. */
  71 enum date_unit
  72   {
  73     DATE_YEARS,
  74     DATE_QUARTERS,
  75     DATE_MONTHS,
  76     DATE_WEEKS,
  77     DATE_DAYS,
  78     DATE_HOURS,
  79     DATE_MINUTES,
  80     DATE_SECONDS
  81   };
  82
  83 /* Stores in *UNIT the unit whose name is NAME.
  84    Return success. */
  85 static enum date_unit
  86 recognize_unit (struct substring name, const struct expression *e,
  87                 const struct expr_node *n, enum date_unit *unit)
  88 {
  89   struct unit_name
  90     {
  91       enum date_unit unit;
  92       const struct substring name;
  93     };
  94   static const struct unit_name unit_names[] =
  95     {
  96       { DATE_YEARS, SS_LITERAL_INITIALIZER ("years") },
  97       { DATE_QUARTERS, SS_LITERAL_INITIALIZER ("quarters") },
  98       { DATE_MONTHS, SS_LITERAL_INITIALIZER ("months") },
  99       { DATE_WEEKS, SS_LITERAL_INITIALIZER ("weeks") },
 100       { DATE_DAYS, SS_LITERAL_INITIALIZER ("days") },
 101       { DATE_HOURS, SS_LITERAL_INITIALIZER ("hours") },
 102       { DATE_MINUTES, SS_LITERAL_INITIALIZER ("minutes") },
 103       { DATE_SECONDS, SS_LITERAL_INITIALIZER ("seconds") },
 104     };
 105   const int n_unit_names = sizeof unit_names / sizeof *unit_names;
 106
 107   const struct unit_name *un;
 108
 109   for (un = unit_names; un < &unit_names[n_unit_names]; un++)
 110     if (ss_equals_case (un->name, name))
 111       {
 112         *unit = un->unit;
 113         return true;
 114       }
 115
 116   msg_at (SE, expr_location (e, n),
 117           _("Unrecognized date unit `%.*s'.  "
 118             "Valid date units are `%s', `%s', `%s', "
 119             "`%s', `%s', `%s', `%s', and `%s'."),
 120           (int) ss_length (name), ss_data (name),
 121           "years", "quarters", "months",
 122           "weeks", "days", "hours", "minutes", "seconds");
 123
 124   return false;
 125 }
 126
 127 /* Returns the number of whole years from DATE1 to DATE2,
 128    where a year is defined as the same or later month, day, and
 129    time of day. */
 130 static int
 131 year_diff (double date1, double date2)
 132 {
 133   int y1, m1, d1, yd1;
 134   int y2, m2, d2, yd2;
 135   int diff;
 136
 137   assert (date2 >= date1);
 138   calendar_offset_to_gregorian (date1 / DAY_S, &y1, &m1, &d1, &yd1);
 139   calendar_offset_to_gregorian (date2 / DAY_S, &y2, &m2, &d2, &yd2);
 140
 141   diff = y2 - y1;
 142   if (diff > 0)
 143     {
 144       int yd1 = 32 * m1 + d1;
 145       int yd2 = 32 * m2 + d2;
 146       if (yd2 < yd1
 147           || (yd2 == yd1 && fmod (date2, DAY_S) < fmod (date1, DAY_S)))
 148         diff--;
 149     }
 150   return diff;
 151 }
 152
 153 /* Returns the number of whole months from DATE1 to DATE2,
 154    where a month is defined as the same or later day and time of
 155    day. */
 156 static int
 157 month_diff (double date1, double date2)
 158 {
 159   int y1, m1, d1, yd1;
 160   int y2, m2, d2, yd2;
 161   int diff;
 162
 163   assert (date2 >= date1);
 164   calendar_offset_to_gregorian (date1 / DAY_S, &y1, &m1, &d1, &yd1);
 165   calendar_offset_to_gregorian (date2 / DAY_S, &y2, &m2, &d2, &yd2);
 166
 167   diff = ((y2 * 12) + m2) - ((y1 * 12) + m1);
 168   if (diff > 0
 169       && (d2 < d1
 170           || (d2 == d1 && fmod (date2, DAY_S) < fmod (date1, DAY_S))))
 171     diff--;
 172   return diff;
 173 }
 174
 175 /* Returns the number of whole quarter from DATE1 to DATE2,
 176    where a quarter is defined as three months. */
 177 static int
 178 quarter_diff (double date1, double date2)
 179 {
 180   return month_diff (date1, date2) / 3;
 181 }
 182
 183 /* Returns the number of seconds in the given UNIT. */
 184 static int
 185 date_unit_duration (enum date_unit unit)
 186 {
 187   switch (unit)
 188     {
 189     case DATE_WEEKS:
 190       return WEEK_S;
 191
 192     case DATE_DAYS:
 193       return DAY_S;
 194
 195     case DATE_HOURS:
 196       return H_S;
 197
 198     case DATE_MINUTES:
 199       return MIN_S;
 200
 201     case DATE_SECONDS:
 202       return 1;
 203
 204     default:
 205       NOT_REACHED ();
 206     }
 207 }
 208
 209 /* Returns the span from DATE1 to DATE2 in terms of UNIT_NAME. */
 210 double
 211 expr_date_difference (double date1, double date2, struct substring unit_name,
 212                       const struct expression *e, const struct expr_node *n)
 213 {
 214   enum date_unit unit;
 215   if (!recognize_unit (unit_name, e, n->args[2], &unit))
 216     return SYSMIS;
 217
 218   switch (unit)
 219     {
 220     case DATE_YEARS:
 221       return (date2 >= date1
 222               ? year_diff (date1, date2)
 223               : -year_diff (date2, date1));
 224
 225     case DATE_QUARTERS:
 226       return (date2 >= date1
 227               ? quarter_diff (date1, date2)
 228               : -quarter_diff (date2, date1));
 229
 230     case DATE_MONTHS:
 231       return (date2 >= date1
 232               ? month_diff (date1, date2)
 233               : -month_diff (date2, date1));
 234
 235     case DATE_WEEKS:
 236     case DATE_DAYS:
 237     case DATE_HOURS:
 238     case DATE_MINUTES:
 239     case DATE_SECONDS:
 240       return trunc ((date2 - date1) / date_unit_duration (unit));
 241     }
 242
 243   NOT_REACHED ();
 244 }
 245
 246 /* How to deal with days out of range for a given month. */
 247 enum date_sum_method
 248   {
 249     SUM_ROLLOVER,       /* Roll them over to the next month. */
 250     SUM_CLOSEST         /* Use the last day of the month. */
 251   };
 252
 253 /* Stores in *METHOD the method whose name is NAME.
 254    Return success. */
 255 static bool
 256 recognize_method (struct substring method_name,
 257                   const struct expression *e, const struct expr_node *n,
 258                   enum date_sum_method *method)
 259 {
 260   if (ss_equals_case (method_name, ss_cstr ("closest")))
 261     {
 262       *method = SUM_CLOSEST;
 263       return true;
 264     }
 265   else if (ss_equals_case (method_name, ss_cstr ("rollover")))
 266     {
 267       *method = SUM_ROLLOVER;
 268       return true;
 269     }
 270   else
 271     {
 272       msg_at (SE, expr_location (e, n),
 273               _("Invalid DATESUM method.  "
 274                 "Valid choices are `%s' and `%s'."), "closest", "rollover");
 275       return false;
 276     }
 277 }
 278
 279 /* Returns DATE advanced by the given number of MONTHS, with
 280    day-of-month overflow resolved using METHOD. */
 281 static double
 282 add_months (double date, int months, enum date_sum_method method,
 283             const struct expression *e, const struct expr_node *n)
 284 {
 285   int y, m, d, yd;
 286   double output;
 287   char *error;
 288
 289   calendar_offset_to_gregorian (date / DAY_S, &y, &m, &d, &yd);
 290   y += months / 12;
 291   m += months % 12;
 292   if (m < 1)
 293     {
 294       m += 12;
 295       y--;
 296     }
 297   else if (m > 12)
 298     {
 299       m -= 12;
 300       y++;
 301     }
 302   assert (m >= 1 && m <= 12);
 303
 304   if (method == SUM_CLOSEST && d > calendar_days_in_month (y, m))
 305     d = calendar_days_in_month (y, m);
 306
 307   output = calendar_gregorian_to_offset (y, m, d, settings_get_fmt_settings (),
 308                                          &error);
 309   if (output != SYSMIS)
 310     output = (output * DAY_S) + fmod (date, DAY_S);
 311   else
 312     {
 313       msg_at (SE, expr_location (e, n), "%s", error);
 314       free (error);
 315     }
 316   return output;
 317 }
 318
 319 /* Returns DATE advanced by the given QUANTITY of units given in
 320    UNIT_NAME, with day-of-month overflow resolved using
 321    METHOD_NAME. */
 322 static double
 323 expr_date_sum__ (double date, double quantity, struct substring unit_name,
 324                  enum date_sum_method method,
 325                  const struct expression *e, const struct expr_node *n)
 326 {
 327   enum date_unit unit;
 328   if (!recognize_unit (unit_name, e, n->args[2], &unit))
 329     return SYSMIS;
 330
 331   switch (unit)
 332     {
 333     case DATE_YEARS:
 334       return add_months (date, trunc (quantity) * 12, method, e, n);
 335
 336     case DATE_QUARTERS:
 337       return add_months (date, trunc (quantity) * 3, method, e, n);
 338
 339     case DATE_MONTHS:
 340       return add_months (date, trunc (quantity), method, e, n);
 341
 342     case DATE_WEEKS:
 343     case DATE_DAYS:
 344     case DATE_HOURS:
 345     case DATE_MINUTES:
 346     case DATE_SECONDS:
 347       return date + quantity * date_unit_duration (unit);
 348     }
 349
 350   NOT_REACHED ();
 351 }
 352
 353 /* Returns DATE advanced by the given QUANTITY of units given in
 354    UNIT_NAME, with day-of-month overflow resolved using
 355    METHOD_NAME. */
 356 double
 357 expr_date_sum (double date, double quantity, struct substring unit_name,
 358                struct substring method_name,
 359                const struct expression *e, const struct expr_node *n)
 360 {
 361   enum date_sum_method method;
 362   if (!recognize_method (method_name, e, n->args[3], &method))
 363     return SYSMIS;
 364
 365   return expr_date_sum__ (date, quantity, unit_name, method, e, n);
 366 }
 367
 368 /* Returns DATE advanced by the given QUANTITY of units given in
 369    UNIT_NAME, with day-of-month overflow resolved using
 370    METHOD_NAME. */
 371 double
 372 expr_date_sum_closest (double date, double quantity, struct substring unit_name,
 373                        const struct expression *e, const struct expr_node *n)
 374 {
 375   return expr_date_sum__ (date, quantity, unit_name, SUM_CLOSEST, e, n);
 376 }
 377
 378 int
 379 compare_string_3way (const struct substring *a, const struct substring *b)
 380 {
 381   size_t i;
 382
 383   for (i = 0; i < a->length && i < b->length; i++)
 384     if (a->string[i] != b->string[i])
 385       return a->string[i] < b->string[i] ? -1 : 1;
 386   for (; i < a->length; i++)
 387     if (a->string[i] != ' ')
 388       return 1;
 389   for (; i < b->length; i++)
 390     if (b->string[i] != ' ')
 391       return -1;
 392   return 0;
 393 }
 394
 395 size_t
 396 count_valid (double *d, size_t n)
 397 {
 398   size_t n_valid = 0;
 399   for (size_t i = 0; i < n; i++)
 400     n_valid += is_valid (d[i]);
 401   return n_valid;
 402 }
 403
 404 struct substring
 405 alloc_string (struct expression *e, size_t length)
 406 {
 407   struct substring s;
 408   s.length = length;
 409   s.string = pool_alloc (e->eval_pool, length);
 410   return s;
 411 }
 412
 413 struct substring
 414 copy_string (struct expression *e, const char *old, size_t length)
 415 {
 416   struct substring s = alloc_string (e, length);
 417   memcpy (s.string, old, length);
 418   return s;
 419 }
 420
 421 static double
 422 round__ (double x, double mult, double fuzzbits, double adjustment)
 423 {
 424   if (fuzzbits <= 0)
 425     fuzzbits = settings_get_fuzzbits ();
 426   adjustment += exp2 (fuzzbits - DBL_MANT_DIG);
 427
 428   x /= mult;
 429   x = x >= 0. ? floor (x + adjustment) : -floor (-x + adjustment);
 430   return x * mult;
 431 }
 432
 433 double
 434 round_nearest (double x, double mult, double fuzzbits)
 435 {
 436   return round__ (x, mult, fuzzbits, .5);
 437 }
 438
 439 double
 440 round_zero (double x, double mult, double fuzzbits)
 441 {
 442   return round__ (x, mult, fuzzbits, 0);
 443 }
 444
 445 struct substring
 446 replace_string (struct expression *e,
 447                 struct substring haystack,
 448                 struct substring needle,
 449                 struct substring replacement,
 450                 int n)
 451 {
 452   if (!needle.length || haystack.length < needle.length || n <= 0)
 453     return haystack;
 454
 455   struct substring result = alloc_string (e, MAX_STRING);
 456   result.length = 0;
 457
 458   size_t i = 0;
 459   while (i <= haystack.length - needle.length)
 460     if (!memcmp (&haystack.string[i], needle.string, needle.length))
 461       {
 462         size_t copy_len = MIN (replacement.length, MAX_STRING - result.length);
 463         memcpy (&result.string[result.length], replacement.string, copy_len);
 464         result.length += copy_len;
 465         i += needle.length;
 466
 467         if (--n < 1)
 468           break;
 469       }
 470     else
 471       {
 472         if (result.length < MAX_STRING)
 473           result.string[result.length++] = haystack.string[i];
 474         i++;
 475       }
 476   while (i < haystack.length && result.length < MAX_STRING)
 477     result.string[result.length++] = haystack.string[i++];
 478
 479   return result;
 480 }
 481
 482 static int
 483 compare_doubles (const void *a_, const void *b_)
 484 {
 485   const double *ap = a_;
 486   const double *bp = b_;
 487   double a = *ap;
 488   double b = *bp;
 489
 490   /* Sort SYSMIS to the end. */
 491   return (a == b ? 0
 492           : a == SYSMIS ? 1
 493           : b == SYSMIS ? -1
 494           : a > b ? 1 : -1);
 495 }
 496
 497 double
 498 median (double *a, size_t n)
 499 {
 500   /* Sort the array in-place, sorting SYSMIS to the end. */
 501   qsort (a, n, sizeof *a, compare_doubles);
 502
 503   /* Drop SYSMIS. */
 504   n = count_valid (a, n);
 505
 506   return (!n ? SYSMIS
 507           : n % 2 ? a[n / 2]
 508           : (a[n / 2 - 1] + a[n / 2]) / 2.0);
 509 }
 510
 511 const struct variable *
 512 expr_index_vector (const struct expression *e, const struct expr_node *n,
 513                    const struct vector *v, double idx)
 514 {
 515   if (idx >= 1 && idx <= vector_get_n_vars (v))
 516     return vector_get_var (v, idx - 1);
 517
 518   msg_at (SE, expr_location (e, n),
 519           _("Index outside valid range 1 to %zu, inclusive, for vector %s.  "
 520             "The value will be treated as system-missing."),
 521           vector_get_n_vars (v), vector_get_name (v));
 522   if (idx == SYSMIS)
 523     msg_at (SN, expr_location (e, n->args[0]),
 524             _("The index is system-missing."));
 525   else
 526     msg_at (SN, expr_location (e, n->args[0]),
 527             _("The index has value %g."), idx);
 528   return NULL;
 529 }