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