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 /* Returns the noncentral beta cumulative distribution function
 462    value for the given arguments.
 463
 464    FIXME: The accuracy of this function is not entirely
 465    satisfactory.  We only match the example values given in AS
 466    310 to the first 5 significant digits. */
 467 double
 468 ncdf_beta (double x, double a, double b, double lambda)
 469 {
 470   double c;
 471
 472   if (x <= 0. || x >= 1. || a <= 0. || b <= 0. || lambda <= 0.)
 473     return SYSMIS;
 474
 475   c = lambda / 2.;
 476   if (lambda < 54.)
 477     {
 478       /* Algorithm AS 226. */
 479       double x0, a0, beta, temp, gx, q, ax, sumq, sum;
 480       double err_max = 2 * DBL_EPSILON;
 481       double err_bound;
 482       int iter_max = 100;
 483       int iter;
 484
 485       x0 = floor (c - 5.0 * sqrt (c));
 486       if (x0 < 0.)
 487         x0 = 0.;
 488       a0 = a + x0;
 489       beta = (gsl_sf_lngamma (a0)
 490               + gsl_sf_lngamma (b)
 491               - gsl_sf_lngamma (a0 + b));
 492       temp = gsl_sf_beta_inc (a0, b, x);
 493       gx = exp (a0 * log (x) + b * log (1. - x) - beta - log (a0));
 494       if (a0 >= a)
 495         q = exp (-c + x0 * log (c)) - gsl_sf_lngamma (x0 + 1.);
 496       else
 497         q = exp (-c);
 498       ax = q * temp;
 499       sumq = 1. - q;
 500       sum = ax;
 501
 502       iter = 0;
 503       do
 504         {
 505           iter++;
 506           temp -= gx;
 507           gx = x * (a + b + iter - 1.) * gx / (a + iter);
 508           q *= c / iter;
 509           sumq -= q;
 510           ax = temp * q;
 511           sum += ax;
 512
 513           err_bound = (temp - gx) * sumq;
 514         }
 515       while (iter < iter_max && err_bound > err_max);
 516
 517       return sum;
 518     }
 519   else
 520     {
 521       /* Algorithm AS 310. */
 522       double m, m_sqrt;
 523       int iter, iter_lower, iter_upper, iter1, iter2, j;
 524       double t, q, r, psum, beta, s1, gx, fx, temp, ftemp, t0, s0, sum, s;
 525       double err_bound;
 526       double err_max = 2 * DBL_EPSILON;
 527
 528       iter = 0;
 529
 530       m = floor (c + .5);
 531       m_sqrt = sqrt (m);
 532       iter_lower = m - 5. * m_sqrt;
 533       iter_upper = m + 5. * m_sqrt;
 534
 535       t = -c + m * log (c) - gsl_sf_lngamma (m + 1.);
 536       q = exp (t);
 537       r = q;
 538       psum = q;
 539       beta = (gsl_sf_lngamma (a + m)
 540               + gsl_sf_lngamma (b)
 541               - gsl_sf_lngamma (a + m + b));
 542       s1 = (a + m) * log (x) + b * log (1. - x) - log (a + m) - beta;
 543       fx = gx = exp (s1);
 544       ftemp = temp = gsl_sf_beta_inc (a + m, b, x);
 545       iter++;
 546       sum = q * temp;
 547       iter1 = m;
 548
 549       while (iter1 >= iter_lower && q >= err_max)
 550         {
 551           q = q * iter1 / c;
 552           iter++;
 553           gx = (a + iter1) / (x * (a + b + iter1 - 1.)) * gx;
 554           iter1--;
 555           temp += gx;
 556           psum += q;
 557           sum += q * temp;
 558         }
 559
 560       t0 = (gsl_sf_lngamma (a + b)
 561             - gsl_sf_lngamma (a + 1.)
 562             - gsl_sf_lngamma (b));
 563       s0 = a * log (x) + b * log (1. - x);
 564
 565       s = 0.;
 566       for (j = 0; j < iter1; j++)
 567         {
 568           double t1;
 569           s += exp (t0 + s0 + j * log (x));
 570           t1 = log (a + b + j) - log (a + 1. + j) + t0;
 571           t0 = t1;
 572         }
 573
 574       err_bound = (1. - gsl_sf_gamma_inc_P (iter1, c)) * (temp + s);
 575       q = r;
 576       temp = ftemp;
 577       gx = fx;
 578       iter2 = m;
 579       for (;;)
 580         {
 581           double ebd = err_bound + (1. - psum) * temp;
 582           if (ebd < err_max || iter >= iter_upper)
 583             break;
 584
 585           iter2++;
 586           iter++;
 587           q = q * c / iter2;
 588           psum += q;
 589           temp -= gx;
 590           gx = x * (a + b + iter2 - 1.) / (a + iter2) * gx;
 591           sum += q * temp;
 592         }
 593
 594       return sum;
 595     }
 596 }
 597
 598 double
 599 cdf_bvnor (double x0, double x1, double r)
 600 {
 601   double z = pow2 (x0) - 2. * r * x0 * x1 + pow2 (x1);
 602   return exp (-z / (2. * (1 - r * r))) * (2. * M_PI * sqrt (1 - r * r));
 603 }
 604
 605 double
 606 idf_fdist (double P, double df1, double df2)
 607 {
 608   double temp = gsl_cdf_beta_Pinv (P, df1 / 2, df2 / 2);
 609   return temp * df2 / ((1. - temp) * df1);
 610 }
 611
 612 /*
 613  *  Mathlib : A C Library of Special Functions
 614  *  Copyright (C) 1998 Ross Ihaka
 615  *  Copyright (C) 2000 The R Development Core Team
 616  *
 617  *  This program is free software; you can redistribute it and/or
 618  *  modify
 619  *  it under the terms of the GNU General Public License as
 620  *  published by
 621  *  the Free Software Foundation; either version 2 of the
 622  *  License, or
 623  *  (at your option) any later version.
 624  *
 625  *  This program is distributed in the hope that it will be
 626  *  useful,
 627  *  but WITHOUT ANY WARRANTY; without even the implied warranty
 628  *  of
 629  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 630  *  GNU General Public License for more details.
 631  *
 632  *  You should have received a copy of the GNU General Public
 633  *  License
 634  *  along with this program; if not, write to the Free Software
 635  *  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
 636  *  02110-1301 USA.
 637  */
 638
 639 /* Returns the density of the noncentral beta distribution with
 640    noncentrality parameter LAMBDA. */
 641 double
 642 npdf_beta (double x, double a, double b, double lambda)
 643 {
 644   if (lambda < 0. || a <= 0. || b <= 0.)
 645     return SYSMIS;
 646   else if (lambda == 0.)
 647     return gsl_ran_beta_pdf (x, a, b);
 648   else
 649     {
 650       double max_error = 2 * DBL_EPSILON;
 651       int max_iter = 200;
 652       double term = gsl_ran_beta_pdf (x, a, b);
 653       double lambda2 = 0.5 * lambda;
 654       double weight = exp (-lambda2);
 655       double sum = weight * term;
 656       double psum = weight;
 657       int k;
 658       for (k = 1; k <= max_iter && 1 - psum < max_error; k++) {
 659         weight *= lambda2 / k;
 660         term *= x * (a + b) / a;
 661         sum += weight * term;
 662         psum += weight;
 663         a += 1;
 664       }
 665       return sum;
 666     }
 667 }
 668
 669 static double
 670 round__ (double x, double mult, double fuzzbits, double adjustment)
 671 {
 672   if (fuzzbits <= 0)
 673     fuzzbits = settings_get_fuzzbits ();
 674   adjustment += exp2 (fuzzbits - DBL_MANT_DIG);
 675
 676   x /= mult;
 677   x = x >= 0. ? floor (x + adjustment) : -floor (-x + adjustment);
 678   return x * mult;
 679 }
 680
 681 double
 682 round_nearest (double x, double mult, double fuzzbits)
 683 {
 684   return round__ (x, mult, fuzzbits, .5);
 685 }
 686
 687 double
 688 round_zero (double x, double mult, double fuzzbits)
 689 {
 690   return round__ (x, mult, fuzzbits, 0);
 691 }
 692
 693 struct substring
 694 replace_string (struct expression *e,
 695                 struct substring haystack,
 696                 struct substring needle,
 697                 struct substring replacement,
 698                 double n)
 699 {
 700   if (!needle.length
 701       || haystack.length < needle.length
 702       || n <= 0
 703       || n == SYSMIS)
 704     return haystack;
 705
 706   struct substring result = alloc_string (e, MAX_STRING);
 707   result.length = 0;
 708
 709   size_t i = 0;
 710   while (i <= haystack.length - needle.length)
 711     if (!memcmp (&haystack.string[i], needle.string, needle.length))
 712       {
 713         size_t copy_len = MIN (replacement.length, MAX_STRING - result.length);
 714         memcpy (&result.string[result.length], replacement.string, copy_len);
 715         result.length += copy_len;
 716         i += needle.length;
 717
 718         if (--n < 1)
 719           break;
 720       }
 721     else
 722       {
 723         if (result.length < MAX_STRING)
 724           result.string[result.length++] = haystack.string[i];
 725         i++;
 726       }
 727   while (i < haystack.length && result.length < MAX_STRING)
 728     result.string[result.length++] = haystack.string[i++];
 729
 730   return result;
 731 }
 732
 733 static int
 734 compare_doubles (const void *a_, const void *b_)
 735 {
 736   const double *ap = a_;
 737   const double *bp = b_;
 738   double a = *ap;
 739   double b = *bp;
 740
 741   /* Sort SYSMIS to the end. */
 742   return (a == b ? 0
 743           : a == SYSMIS ? 1
 744           : b == SYSMIS ? -1
 745           : a > b ? 1 : -1);
 746 }
 747
 748 double
 749 median (double *a, size_t n)
 750 {
 751   /* Sort the array in-place, sorting SYSMIS to the end. */
 752   qsort (a, n, sizeof *a, compare_doubles);
 753
 754   /* Drop SYSMIS. */
 755   n = count_valid (a, n);
 756
 757   return (!n ? SYSMIS
 758           : n % 2 ? a[n / 2]
 759           : (a[n / 2 - 1] + a[n / 2]) / 2.0);
 760 }