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