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