pintos-os.org Git - pspp/blob - src/language/expressions/operations.def

   1 // -*- c -*-
   2 //
   3 // PSPP - a program for statistical analysis.
   4 // Copyright (C) 2005, 2006, 2009, 2010, 2011, 2012, 2015, 2016 Free Software Foundation, Inc.
   5 //
   6 // This program is free software: you can redistribute it and/or modify
   7 // it under the terms of the GNU General Public License as published by
   8 // the Free Software Foundation, either version 3 of the License, or
   9 // (at your option) any later version.
  10 //
  11 // This program is distributed in the hope that it will be useful,
  12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14 // GNU General Public License for more details.
  15 //
  16 // You should have received a copy of the GNU General Public License
  17 // along with this program.  If not, see <http://www.gnu.org/licenses/>.
  18
  19 operator NEG (x) = -x;
  20
  21 operator ADD (a, b) = a + b;
  22 operator SUB (a, b) = a - b;
  23
  24 absorb_miss operator MUL (a, b)
  25 = (a == 0. || b == 0. ? 0.
  26    : a == SYSMIS || b == SYSMIS ? SYSMIS
  27    : a * b);
  28
  29 absorb_miss operator DIV (a, b)
  30 = (a == 0. ? 0.
  31    : a == SYSMIS || b == SYSMIS ? SYSMIS
  32    : a / b);
  33
  34 absorb_miss operator POW (a, b)
  35 = (a == SYSMIS ? (b == 0. ? 1. : a)
  36    : b == SYSMIS ? (a == 0. ? 0. : SYSMIS)
  37    : a == 0. && b <= 0. ? SYSMIS
  38    : pow (a, b));
  39
  40 absorb_miss boolean operator AND (boolean a, boolean b)
  41 = (a == 0. ? 0.
  42    : b == 0. ? 0.
  43    : b == SYSMIS ? SYSMIS
  44    : a);
  45
  46 absorb_miss boolean operator OR (boolean a, boolean b)
  47 = (a == 1. ? 1.
  48    : b == 1. ? 1.
  49    : b == SYSMIS ? SYSMIS
  50    : a);
  51
  52 boolean operator NOT (boolean a)
  53 = (a == 0. ? 1.
  54    : a == 1. ? 0.
  55    : SYSMIS);
  56
  57 // Numeric relational operators.
  58 boolean operator EQ (a, b) = a == b;
  59 boolean operator GE (a, b) = a >= b;
  60 boolean operator GT (a, b) = a > b;
  61 boolean operator LE (a, b) = a <= b;
  62 boolean operator LT (a, b) = a < b;
  63 boolean operator NE (a, b) = a != b;
  64
  65 // String relational operators.
  66 boolean operator EQ_STRING (string a, string b) = compare_string_3way (&a, &b) == 0;
  67 boolean operator GE_STRING (string a, string b) = compare_string_3way (&a, &b) >= 0;
  68 boolean operator GT_STRING (string a, string b) = compare_string_3way (&a, &b) > 0;
  69 boolean operator LE_STRING (string a, string b) = compare_string_3way (&a, &b) <= 0;
  70 boolean operator LT_STRING (string a, string b) = compare_string_3way (&a, &b) < 0;
  71 boolean operator NE_STRING (string a, string b) = compare_string_3way (&a, &b) != 0;
  72
  73 // Unary functions.
  74 function ABS (x) = fabs (x);
  75 extension function ACOS (x >= -1 && x <= 1) = acos (x);
  76 function ASIN (x >= -1 && x <= 1) = asin (x);
  77 function ATAN (x) = atan (x);
  78 extension function ARCOS (x >= -1 && x <= 1) = acos (x);
  79 function ARSIN (x >= -1 && x <= 1) = asin (x);
  80 function ARTAN (x) = atan (x);
  81 function COS (x) = cos (x);
  82 function EXP (x) = check_errno (exp (x));
  83 function LG10(x) = check_errno (log10 (x));
  84 function LN (x) = check_errno (log (x));
  85 function LNGAMMA (x >= 0) = gsl_sf_lngamma (x);
  86 function MOD10 (x) = fmod (x, 10);
  87 function RND (x) = round_nearest (x, 1, 0);
  88 function RND (x, mult != 0) = round_nearest (x, mult, 0);
  89 function RND (x, mult != 0, fuzzbits >= 0) = round_nearest (x, mult, fuzzbits);
  90 function SIN (x) = sin (x);
  91 function SQRT (x >= 0) = sqrt (x);
  92 function TAN (x) = check_errno (tan (x));
  93 function TRUNC (x) = round_zero (x, 1, 0);
  94 function TRUNC (x, mult != 0) = round_zero (x, mult, 0);
  95 function TRUNC (x, mult != 0, fuzzbits >= 0) = round_zero (x, mult, fuzzbits);
  96
  97 absorb_miss function MOD (n, d)
  98 {
  99   if (d != SYSMIS)
 100     return n != SYSMIS ? fmod (n, d) : SYSMIS;
 101   else
 102     return n != 0. ? SYSMIS : 0.;
 103 }
 104
 105 // N-ary numeric functions.
 106 absorb_miss boolean function ANY (x != SYSMIS, a[n])
 107 {
 108   int sysmis = 0;
 109   size_t i;
 110
 111   for (i = 0; i < n; i++)
 112     if (a[i] == x)
 113       return 1.;
 114     else if (a[i] == SYSMIS)
 115       sysmis = 1;
 116
 117   return sysmis ? SYSMIS : 0.;
 118 }
 119
 120 boolean function ANY (string x, string a[n])
 121 {
 122   size_t i;
 123
 124   for (i = 0; i < n; i++)
 125     if (!compare_string_3way (&x, &a[i]))
 126       return 1.;
 127   return 0.;
 128 }
 129
 130 function CFVAR.2 (a[n])
 131 {
 132   double mean, variance;
 133
 134   moments_of_doubles (a, n, NULL, &mean, &variance, NULL, NULL);
 135
 136   if (mean == SYSMIS || mean == 0 || variance == SYSMIS)
 137     return SYSMIS;
 138   else
 139     return sqrt (variance) / mean;
 140 }
 141
 142 function MAX.1 (a[n])
 143 {
 144   double max;
 145   size_t i;
 146
 147   max = -DBL_MAX;
 148   for (i = 0; i < n; i++)
 149     if (a[i] != SYSMIS && a[i] > max)
 150       max = a[i];
 151   return max;
 152 }
 153
 154 string function MAX (string a[n])
 155 {
 156   struct substring *max;
 157   size_t i;
 158
 159   max = &a[0];
 160   for (i = 1; i < n; i++)
 161     if (compare_string_3way (&a[i], max) > 0)
 162       max = &a[i];
 163   return *max;
 164 }
 165
 166 function MEAN.1 (a[n])
 167 {
 168   double mean;
 169   moments_of_doubles (a, n, NULL, &mean, NULL, NULL, NULL);
 170   return mean;
 171 }
 172
 173 function MEDIAN.1 (a[n])
 174 {
 175   return median (a, n);
 176 }
 177
 178 function MIN.1 (a[n])
 179 {
 180   double min;
 181   size_t i;
 182
 183   min = DBL_MAX;
 184   for (i = 0; i < n; i++)
 185     if (a[i] != SYSMIS && a[i] < min)
 186       min = a[i];
 187   return min;
 188 }
 189
 190 string function MIN (string a[n])
 191 {
 192   struct substring *min;
 193   size_t i;
 194
 195   min = &a[0];
 196   for (i = 1; i < n; i++)
 197     if (compare_string_3way (&a[i], min) < 0)
 198       min = &a[i];
 199   return *min;
 200 }
 201
 202 absorb_miss function NMISS (a[n])
 203 {
 204   size_t i;
 205   size_t missing_cnt = 0;
 206
 207   for (i = 0; i < n; i++)
 208     missing_cnt += a[i] == SYSMIS;
 209   return missing_cnt;
 210 }
 211
 212 absorb_miss function NVALID (a[n])
 213 {
 214   size_t i;
 215   size_t valid_cnt = 0;
 216
 217   for (i = 0; i < n; i++)
 218     valid_cnt += a[i] != SYSMIS;
 219   return valid_cnt;
 220 }
 221
 222 absorb_miss boolean function RANGE (x != SYSMIS, a[n*2])
 223 {
 224   size_t i;
 225   int sysmis = 0;
 226
 227   for (i = 0; i < n; i++)
 228     {
 229       double w = a[2 * i];
 230       double y = a[2 * i + 1];
 231       if (w != SYSMIS && y != SYSMIS)
 232         {
 233           if (w <= x && x <= y)
 234             return 1.0;
 235         }
 236       else
 237         sysmis = 1;
 238     }
 239   return sysmis ? SYSMIS : 0.;
 240 }
 241
 242 boolean function RANGE (string x, string a[n*2])
 243 {
 244   int i;
 245
 246   for (i = 0; i < n; i++)
 247     {
 248       struct substring *w = &a[2 * i];
 249       struct substring *y = &a[2 * i + 1];
 250       if (compare_string_3way (w, &x) <= 0 && compare_string_3way (&x, y) <= 0)
 251         return 1.;
 252     }
 253   return 0.;
 254 }
 255
 256 function SD.2 (a[n])
 257 {
 258   double variance;
 259   moments_of_doubles (a, n, NULL, NULL, &variance, NULL, NULL);
 260   return sqrt (variance);
 261 }
 262
 263 function SUM.1 (a[n])
 264 {
 265   double sum;
 266   size_t i;
 267
 268   sum = 0.;
 269   for (i = 0; i < n; i++)
 270     if (a[i] != SYSMIS)
 271       sum += a[i];
 272   return sum;
 273 }
 274
 275 function VARIANCE.2 (a[n])
 276 {
 277   double variance;
 278   moments_of_doubles (a, n, NULL, NULL, &variance, NULL, NULL);
 279   return variance;
 280 }
 281
 282 // Time construction & extraction functions.
 283 function TIME.HMS (h, m, s)
 284 {
 285   if ((h > 0. || m > 0. || s > 0.) && (h < 0. || m < 0. || s < 0.))
 286     {
 287       msg (SW, _("TIME.HMS cannot mix positive and negative arguments."));
 288       return SYSMIS;
 289     }
 290   else
 291     return H_S * h + MIN_S * m + s;
 292 }
 293 function TIME.DAYS (days) = days * DAY_S;
 294 function CTIME.DAYS (time) = time / DAY_S;
 295 function CTIME.HOURS (time) = time / H_S;
 296 function CTIME.MINUTES (time) = time / MIN_S;
 297 function CTIME.SECONDS (time) = time;
 298
 299 // Date construction functions.
 300 function DATE.DMY (d, m, y) = expr_ymd_to_date (y, m, d);
 301 function DATE.MDY (m, d, y) = expr_ymd_to_date (y, m, d);
 302 function DATE.MOYR (m, y) = expr_ymd_to_date (y, m, 1);
 303 function DATE.QYR (q, y)
 304 {
 305   if (q < 1.0 || q > 4.0 || q != (int) q)
 306     {
 307       msg (SW, _("The first argument to DATE.QYR must be 1, 2, 3, or 4."));
 308       return SYSMIS;
 309     }
 310    return expr_ymd_to_date (y, q * 3 - 2, 1);
 311 }
 312 function DATE.WKYR (w, y) = expr_wkyr_to_date (w, y);
 313 function DATE.YRDAY (y, yday) = expr_yrday_to_date (y, yday);
 314 function YRMODA (y, m, d) = expr_yrmoda (y, m, d);
 315
 316 // Date extraction functions.
 317 function XDATE.TDAY (date) = floor (date / DAY_S);
 318 function XDATE.HOUR (date) = fmod (floor (date / H_S), DAY_H);
 319 function XDATE.MINUTE (date) = fmod (floor (date / H_MIN), H_MIN);
 320 function XDATE.SECOND (date) = fmod (date, MIN_S);
 321 function XDATE.DATE (date) = floor (date / DAY_S) * DAY_S;
 322 function XDATE.TIME (date) = fmod (date, DAY_S);
 323
 324 function XDATE.JDAY (date >= DAY_S) = calendar_offset_to_yday (date / DAY_S);
 325 function XDATE.MDAY (date >= DAY_S) = calendar_offset_to_mday (date / DAY_S);
 326 function XDATE.MONTH (date >= DAY_S)
 327      = calendar_offset_to_month (date / DAY_S);
 328 function XDATE.QUARTER (date >= DAY_S)
 329     = (calendar_offset_to_month (date / DAY_S) - 1) / 3 + 1;
 330 function XDATE.WEEK (date >= DAY_S)
 331     = (calendar_offset_to_yday (date / DAY_S) - 1) / 7 + 1;
 332 function XDATE.WKDAY (date >= DAY_S) = calendar_offset_to_wday (date / DAY_S);
 333 function XDATE.YEAR (date >= DAY_S) = calendar_offset_to_year (date / DAY_S);
 334
 335 // Date arithmetic functions.
 336 no_abbrev function DATEDIFF (date2 >= DAY_S, date1 >= DAY_S, string unit)
 337      = expr_date_difference (date1, date2, unit);
 338 no_abbrev function DATESUM (date, quantity, string unit)
 339      = expr_date_sum (date, quantity, unit, ss_cstr ("closest"));
 340 no_abbrev function DATESUM (date, quantity, string unit, string method)
 341      = expr_date_sum (date, quantity, unit, method);
 342
 343
 344 // String functions.
 345 string function CONCAT (string a[n])
 346      expression e;
 347 {
 348   struct substring dst;
 349   size_t i;
 350
 351   dst = alloc_string (e, MAX_STRING);
 352   dst.length = 0;
 353   for (i = 0; i < n; i++)
 354     {
 355       struct substring *src = &a[i];
 356       size_t copy_len;
 357
 358       copy_len = src->length;
 359       if (dst.length + copy_len > MAX_STRING)
 360         copy_len = MAX_STRING - dst.length;
 361       memcpy (&dst.string[dst.length], src->string, copy_len);
 362       dst.length += copy_len;
 363     }
 364
 365   return dst;
 366 }
 367
 368 function INDEX (string haystack, string needle)
 369 {
 370   if (needle.length == 0)
 371     return SYSMIS;
 372   else
 373     {
 374       int limit = haystack.length - needle.length + 1;
 375       int i;
 376       for (i = 1; i <= limit; i++)
 377         if (!memcmp (&haystack.string[i - 1], needle.string, needle.length))
 378           return i;
 379       return 0;
 380     }
 381 }
 382
 383 function INDEX (string haystack, string needles, needle_len_d)
 384 {
 385   if (needle_len_d <= INT_MIN || needle_len_d >= INT_MAX
 386       || (int) needle_len_d != needle_len_d
 387       || needles.length == 0)
 388     return SYSMIS;
 389   else
 390     {
 391       int needle_len = needle_len_d;
 392       if (needle_len < 0 || needle_len > needles.length
 393           || needles.length % needle_len != 0)
 394         return SYSMIS;
 395       else
 396         {
 397           int limit = haystack.length - needle_len + 1;
 398           int i, j;
 399           for (i = 1; i <= limit; i++)
 400             for (j = 0; j < needles.length; j += needle_len)
 401               if (!memcmp (&haystack.string[i - 1], &needles.string[j],
 402                            needle_len))
 403                 return i;
 404           return 0;
 405         }
 406     }
 407 }
 408
 409 function RINDEX (string haystack, string needle)
 410 {
 411   if (needle.length == 0)
 412     return SYSMIS;
 413   else
 414     {
 415       int limit = haystack.length - needle.length + 1;
 416       int i;
 417       for (i = limit; i >= 1; i--)
 418         if (!memcmp (&haystack.string[i - 1], needle.string, needle.length))
 419           return i;
 420       return 0;
 421     }
 422 }
 423
 424 function RINDEX (string haystack, string needles, needle_len_d)
 425 {
 426   if (needle_len_d <= 0 || needle_len_d >= INT_MAX
 427       || (int) needle_len_d != needle_len_d
 428       || needles.length == 0)
 429     return SYSMIS;
 430   else
 431     {
 432       int needle_len = needle_len_d;
 433       if (needle_len < 0 || needle_len > needles.length
 434           || needles.length % needle_len != 0)
 435         return SYSMIS;
 436       else
 437         {
 438           int limit = haystack.length - needle_len + 1;
 439           int i, j;
 440           for (i = limit; i >= 1; i--)
 441             for (j = 0; j < needles.length; j += needle_len)
 442               if (!memcmp (&haystack.string[i - 1],
 443                            &needles.string[j], needle_len))
 444                 return i;
 445           return 0;
 446         }
 447     }
 448 }
 449
 450 function LENGTH (string s)
 451 {
 452   return s.length;
 453 }
 454
 455 string function LOWER (string s)
 456 {
 457   int i;
 458
 459   for (i = 0; i < s.length; i++)
 460     s.string[i] = tolower ((unsigned char) s.string[i]);
 461   return s;
 462 }
 463
 464 function MBLEN.BYTE (string s, idx)
 465 {
 466   if (idx < 0 || idx >= s.length || (int) idx != idx)
 467     return SYSMIS;
 468   else
 469     return 1;
 470 }
 471
 472 string function UPCASE (string s)
 473 {
 474   int i;
 475
 476   for (i = 0; i < s.length; i++)
 477     s.string[i] = toupper ((unsigned char) s.string[i]);
 478   return s;
 479 }
 480
 481 absorb_miss string function LPAD (string s, n)
 482      expression e;
 483 {
 484   if (n < 0 || n > MAX_STRING || (int) n != n)
 485     return empty_string;
 486   else if (s.length >= n)
 487     return s;
 488   else
 489     {
 490       struct substring t = alloc_string (e, n);
 491       memset (t.string, ' ', n - s.length);
 492       memcpy (&t.string[(int) n - s.length], s.string, s.length);
 493       return t;
 494     }
 495 }
 496
 497 absorb_miss string function LPAD (string s, n, string c)
 498      expression e;
 499 {
 500   if (n < 0 || n > MAX_STRING || (int) n != n || c.length != 1)
 501     return empty_string;
 502   else if (s.length >= n)
 503     return s;
 504   else
 505     {
 506       struct substring t = alloc_string (e, n);
 507       memset (t.string, c.string[0], n - s.length);
 508       memcpy (&t.string[(int) n - s.length], s.string, s.length);
 509       return t;
 510     }
 511 }
 512
 513 string function REPLACE (string haystack, string needle, string replacement)
 514     expression e;
 515   = replace_string (e, haystack, needle, replacement, DBL_MAX);
 516
 517 absorb_miss string function REPLACE (string haystack, string needle,
 518                                      string replacement, n)
 519     expression e;
 520   = replace_string (e, haystack, needle, replacement, n);
 521
 522 absorb_miss string function RPAD (string s, n)
 523      expression e;
 524 {
 525   if (n < 0 || n > MAX_STRING || (int) n != n)
 526     return empty_string;
 527   else if (s.length >= n)
 528     return s;
 529   else
 530     {
 531       struct substring t = alloc_string (e, n);
 532       memcpy (t.string, s.string, s.length);
 533       memset (&t.string[s.length], ' ', n - s.length);
 534       return t;
 535     }
 536 }
 537
 538 absorb_miss string function RPAD (string s, n, string c)
 539      expression e;
 540 {
 541   if (n < 0 || n > MAX_STRING || (int) n != n || c.length != 1)
 542     return empty_string;
 543   else if (s.length >= n)
 544     return s;
 545   else
 546     {
 547       struct substring t = alloc_string (e, n);
 548       memcpy (t.string, s.string, s.length);
 549       memset (&t.string[s.length], c.string[0], n - s.length);
 550       return t;
 551     }
 552 }
 553
 554 string function LTRIM (string s)
 555 {
 556   while (s.length > 0 && s.string[0] == ' ')
 557     {
 558       s.length--;
 559       s.string++;
 560     }
 561   return s;
 562 }
 563
 564 string function LTRIM (string s, string c)
 565 {
 566   if (c.length == 1)
 567     {
 568       while (s.length > 0 && s.string[0] == c.string[0])
 569         {
 570           s.length--;
 571           s.string++;
 572         }
 573       return s;
 574     }
 575   else
 576     return empty_string;
 577 }
 578
 579 string function RTRIM (string s)
 580 {
 581   while (s.length > 0 && s.string[s.length - 1] == ' ')
 582     s.length--;
 583   return s;
 584 }
 585
 586 string function RTRIM (string s, string c)
 587 {
 588   if (c.length == 1)
 589     {
 590       while (s.length > 0 && s.string[s.length - 1] == c.string[0])
 591         s.length--;
 592       return s;
 593     }
 594   else
 595     return empty_string;
 596 }
 597
 598 function NUMBER (string s, ni_format f)
 599 {
 600   union value out;
 601   char *error;
 602
 603   if (s.length > f->w)
 604     s.length = f->w;
 605   error = data_in (s, C_ENCODING, f->type, settings_get_fmt_settings (),
 606                    &out, 0, NULL);
 607   if (error == NULL)
 608     data_in_imply_decimals (s, C_ENCODING, f->type, f->d,
 609                             settings_get_fmt_settings (), &out);
 610   else
 611     {
 612       msg (SE, "Cannot parse `%.*s' as format %s: %s",
 613            (int) s.length, s.string, fmt_name (f->type), error);
 614       free (error);
 615     }
 616   return out.f;
 617 }
 618
 619 absorb_miss string function STRING (x, no_format f)
 620      expression e;
 621 {
 622   union value v;
 623   struct substring dst;
 624   char *s;
 625
 626   v.f = x;
 627
 628   assert (!fmt_is_string (f->type));
 629   s = data_out (&v, C_ENCODING, f, settings_get_fmt_settings ());
 630   dst = alloc_string (e, strlen (s));
 631   strcpy (dst.string, s);
 632   free (s);
 633   return dst;
 634 }
 635
 636 absorb_miss string function STRUNC (string s, n)
 637 {
 638   if (n < 1 || n == SYSMIS)
 639     return empty_string;
 640
 641   if (n < s.length)
 642     s.length = n;
 643   while (s.length > 0 && s.string[s.length - 1] == ' ')
 644     s.length--;
 645   return s;
 646 }
 647
 648 absorb_miss string function SUBSTR (string s, ofs)
 649      expression e;
 650 {
 651   if (ofs >= 1 && ofs <= s.length && (int) ofs == ofs)
 652     return copy_string (e, &s.string[(int) ofs - 1], s.length - ofs + 1);
 653   else
 654     return empty_string;
 655 }
 656
 657 absorb_miss string function SUBSTR (string s, ofs, cnt)
 658      expression e;
 659 {
 660   if (ofs >= 1 && ofs <= s.length && (int) ofs == ofs
 661       && cnt >= 1 && cnt <= INT_MAX && (int) cnt == cnt)
 662     {
 663       int cnt_max = s.length - (int) ofs + 1;
 664       return copy_string (e, &s.string[(int) ofs - 1],
 665                           cnt <= cnt_max ? cnt : cnt_max);
 666     }
 667   else
 668     return empty_string;
 669 }
 670
 671 absorb_miss no_opt no_abbrev string function VALUELABEL (var v)
 672      expression e;
 673      case c;
 674 {
 675   const char *label = var_lookup_value_label (v, case_data (c, v));
 676   if (label != NULL)
 677     return copy_string (e, label, strlen (label));
 678   else
 679     return empty_string;
 680 }
 681
 682 // Artificial.
 683 operator SQUARE (x) = x * x;
 684 boolean operator NUM_TO_BOOLEAN (x, string op_name)
 685 {
 686   if (x == 0. || x == 1. || x == SYSMIS)
 687     return x;
 688
 689   if (!ss_is_empty (op_name))
 690     msg (SE, _("An operand of the %.*s operator was found to have a value "
 691                "other than 0 (false), 1 (true), or the system-missing "
 692                "value.  The result was forced to 0."),
 693          (int) op_name.length, op_name.string);
 694   else
 695     msg (SE, _("A logical expression was found to have a value other than 0 "
 696                "(false), 1 (true), or the system-missing value.  The result "
 697                "was forced to 0."));
 698   return 0.;
 699 }
 700
 701 operator BOOLEAN_TO_NUM (boolean x) = x;
 702
 703 // Beta distribution.
 704 function PDF.BETA (x >= 0 && x <= 1, a > 0, b > 0)
 705      = gsl_ran_beta_pdf (x, a, b);
 706 function CDF.BETA (x >= 0 && x <= 1, a > 0, b > 0) = gsl_cdf_beta_P (x, a, b);
 707 function IDF.BETA (P >= 0 && P <= 1, a > 0, b > 0)
 708      = gsl_cdf_beta_Pinv (P, a, b);
 709 no_opt function RV.BETA (a > 0, b > 0) = gsl_ran_beta (get_rng (), a, b);
 710 function NCDF.BETA (x >= 0, a > 0, b > 0, lambda > 0)
 711      = ncdf_beta (x, a, b, lambda);
 712 function NPDF.BETA (x >= 0, a > 0, b > 0, lambda > 0)
 713      = npdf_beta (x, a, b, lambda);
 714
 715 // Bivariate normal distribution.
 716 function CDF.BVNOR (x0, x1, r >= -1 && r <= 1) = cdf_bvnor (x0, x1, r);
 717 function PDF.BVNOR (x0, x1, r >= -1 && r <= 1)
 718      = gsl_ran_bivariate_gaussian_pdf (x0, x1, 1, 1, r);
 719
 720 // Cauchy distribution.
 721 function CDF.CAUCHY (x, a, b > 0) = gsl_cdf_cauchy_P ((x - a) / b, 1);
 722 function IDF.CAUCHY (P > 0 && P < 1, a, b > 0)
 723      = a + b * gsl_cdf_cauchy_Pinv (P, 1);
 724 function PDF.CAUCHY (x, a, b > 0) = gsl_ran_cauchy_pdf ((x - a) / b, 1) / b;
 725 no_opt function RV.CAUCHY (a, b > 0) = a + b * gsl_ran_cauchy (get_rng (), 1);
 726
 727 // Chi-square distribution.
 728 function CDF.CHISQ (x >= 0, df > 0) = gsl_cdf_chisq_P (x, df);
 729 function IDF.CHISQ (P >= 0 && P < 1, df > 0) = gsl_cdf_chisq_Pinv (P, df);
 730 function PDF.CHISQ (x >= 0, df > 0) = gsl_ran_chisq_pdf (x, df);
 731 no_opt function RV.CHISQ (df > 0) = gsl_ran_chisq (get_rng (), df);
 732 function NCDF.CHISQ (x >= 0, df > 0, c) = unimplemented;
 733 function NPDF.CHISQ (x >= 0, df > 0, c) = unimplemented;
 734 function SIG.CHISQ (x >= 0, df > 0) = gsl_cdf_chisq_Q (x, df);
 735
 736 // Exponential distribution.
 737 function CDF.EXP (x >= 0, a > 0) = gsl_cdf_exponential_P (x, 1. / a);
 738 function IDF.EXP (P >= 0 && P < 1, a > 0)
 739      = gsl_cdf_exponential_Pinv (P, 1. / a);
 740 function PDF.EXP (x >= 0, a > 0) = gsl_ran_exponential_pdf (x, 1. / a);
 741 no_opt function RV.EXP (a > 0) = gsl_ran_exponential (get_rng (), 1. / a);
 742
 743 // Exponential power distribution.
 744 extension function PDF.XPOWER (x, a > 0, b >= 0)
 745      = gsl_ran_exppow_pdf (x, a, b);
 746 no_opt extension function RV.XPOWER (a > 0, b >= 0)
 747      = gsl_ran_exppow (get_rng (), a, b);
 748
 749 // F distribution.
 750 function CDF.F (x >= 0, df1 > 0, df2 > 0) = gsl_cdf_fdist_P (x, df1, df2);
 751 function IDF.F (P >= 0 && P < 1, df1 > 0, df2 > 0) = idf_fdist (P, df1, df2);
 752 function PDF.F (x >= 0, df1 > 0, df2 > 0) = gsl_ran_fdist_pdf (x, df1, df2);
 753 no_opt function RV.F (df1 > 0, df2 > 0) = gsl_ran_fdist (get_rng (), df1, df2);
 754 function NCDF.F (x >= 0, df1 > 0, df2 > 0, lambda >= 0) = unimplemented;
 755 function NPDF.F (x >= 0, df1 > 0, df2 > 0, lmabda >= 0) = unimplemented;
 756 function SIG.F (x >= 0, df1 > 0, df2 > 0) = gsl_cdf_fdist_Q (x, df1, df2);
 757
 758 // Gamma distribution.
 759 function CDF.GAMMA (x >= 0, a > 0, b > 0) = gsl_cdf_gamma_P (x, a, 1. / b);
 760 function IDF.GAMMA (P >= 0 && P <= 1, a > 0, b > 0)
 761      = gsl_cdf_gamma_Pinv (P, a, 1. / b);
 762 function PDF.GAMMA (x >= 0, a > 0, b > 0) = gsl_ran_gamma_pdf (x, a, 1. / b);
 763 no_opt function RV.GAMMA (a > 0, b > 0)
 764      = gsl_ran_gamma (get_rng (), a, 1. / b);
 765
 766 // Half-normal distribution.
 767 function CDF.HALFNRM (x, a, b > 0) = unimplemented;
 768 function IDF.HALFNRM (P > 0 && P < 1, a, b > 0) = unimplemented;
 769 function PDF.HALFNRM (x, a, b > 0) = unimplemented;
 770 no_opt function RV.HALFNRM (a, b > 0) = unimplemented;
 771
 772 // Inverse Gaussian distribution.
 773 function CDF.IGAUSS (x > 0, a > 0, b > 0) = unimplemented;
 774 function IDF.IGAUSS (P >= 0 && P < 1, a > 0, b > 0) = unimplemented;
 775 function PDF.IGAUSS (x > 0, a > 0, b > 0) = unimplemented;
 776 no_opt function RV.IGAUSS (a > 0, b > 0) = unimplemented;
 777
 778 // Landau distribution.
 779 extension function PDF.LANDAU (x) = gsl_ran_landau_pdf (x);
 780 no_opt extension function RV.LANDAU () = gsl_ran_landau (get_rng ());
 781
 782 // Laplace distribution.
 783 function CDF.LAPLACE (x, a, b > 0) = gsl_cdf_laplace_P ((x - a) / b, 1);
 784 function IDF.LAPLACE (P > 0 && P < 1, a, b > 0)
 785      = a + b * gsl_cdf_laplace_Pinv (P, 1);
 786 function PDF.LAPLACE (x, a, b > 0) = gsl_ran_laplace_pdf ((x - a) / b, 1) / b;
 787 no_opt function RV.LAPLACE (a, b > 0)
 788      = a + b * gsl_ran_laplace (get_rng (), 1);
 789
 790 // Levy alpha-stable distribution.
 791 no_opt extension function RV.LEVY (c, alpha > 0 && alpha <= 2)
 792      = gsl_ran_levy (get_rng (), c, alpha);
 793
 794 // Levy skew alpha-stable distribution.
 795 no_opt extension function RV.LVSKEW (c, alpha > 0 && alpha <= 2,
 796                                      beta >= -1 && beta <= 1)
 797      = gsl_ran_levy_skew (get_rng (), c, alpha, beta);
 798
 799 // Logistic distribution.
 800 function CDF.LOGISTIC (x, a, b > 0) = gsl_cdf_logistic_P ((x - a) / b, 1);
 801 function IDF.LOGISTIC (P > 0 && P < 1, a, b > 0)
 802      = a + b * gsl_cdf_logistic_Pinv (P, 1);
 803 function PDF.LOGISTIC (x, a, b > 0)
 804      = gsl_ran_logistic_pdf ((x - a) / b, 1) / b;
 805 no_opt function RV.LOGISTIC (a, b > 0)
 806      = a + b * gsl_ran_logistic (get_rng (), 1);
 807
 808 // Lognormal distribution.
 809 function CDF.LNORMAL (x >= 0, m > 0, s > 0)
 810      = gsl_cdf_lognormal_P (x, log (m), s);
 811 function IDF.LNORMAL (P >= 0 && P < 1, m > 0, s > 0)
 812      = gsl_cdf_lognormal_Pinv (P, log (m), s);
 813 function PDF.LNORMAL (x >= 0, m > 0, s > 0)
 814      = gsl_ran_lognormal_pdf (x, log (m), s);
 815 no_opt function RV.LNORMAL (m > 0, s > 0)
 816      = gsl_ran_lognormal (get_rng (), log (m), s);
 817
 818 // Normal distribution.
 819 function CDF.NORMAL (x, u, s > 0) = gsl_cdf_gaussian_P (x - u, s);
 820 function IDF.NORMAL (P > 0 && P < 1, u, s > 0)
 821      = u + gsl_cdf_gaussian_Pinv (P, s);
 822 function PDF.NORMAL (x, u, s > 0) = gsl_ran_gaussian_pdf ((x - u) / s, 1) / s;
 823 no_opt function RV.NORMAL (u, s > 0) = u + gsl_ran_gaussian (get_rng (), s);
 824 function CDFNORM (x) = gsl_cdf_ugaussian_P (x);
 825 function PROBIT (P > 0 && P < 1) = gsl_cdf_ugaussian_Pinv (P);
 826 no_opt function NORMAL (s > 0) = gsl_ran_gaussian (get_rng (), s);
 827
 828 // Normal tail distribution.
 829 function PDF.NTAIL (x, a > 0, sigma > 0)
 830      = gsl_ran_gaussian_tail_pdf (x, a, sigma);
 831 no_opt function RV.NTAIL (a > 0, sigma > 0)
 832      = gsl_ran_gaussian_tail (get_rng (), a, sigma);
 833
 834 // Pareto distribution.
 835 function CDF.PARETO (x >= a, a > 0, b > 0) = gsl_cdf_pareto_P (x, b, a);
 836 function IDF.PARETO (P >= 0 && P < 1, a > 0, b > 0)
 837      = gsl_cdf_pareto_Pinv (P, b, a);
 838 function PDF.PARETO (x >= a, a > 0, b > 0) = gsl_ran_pareto_pdf (x, b, a);
 839 no_opt function RV.PARETO (a > 0, b > 0) = gsl_ran_pareto (get_rng (), b, a);
 840
 841 // Rayleigh distribution.
 842 extension function CDF.RAYLEIGH (x, sigma > 0) = gsl_cdf_rayleigh_P (x, sigma);
 843 extension function IDF.RAYLEIGH (P >= 0 && P <= 1, sigma > 0)
 844      = gsl_cdf_rayleigh_Pinv (P, sigma);
 845 extension function PDF.RAYLEIGH (x, sigma > 0)
 846      = gsl_ran_rayleigh_pdf (x, sigma);
 847 no_opt extension function RV.RAYLEIGH (sigma > 0)
 848      = gsl_ran_rayleigh (get_rng (), sigma);
 849
 850 // Rayleigh tail distribution.
 851 extension function PDF.RTAIL (x, a, sigma)
 852      = gsl_ran_rayleigh_tail_pdf (x, a, sigma);
 853 no_opt extension function RV.RTAIL (a, sigma)
 854      = gsl_ran_rayleigh_tail (get_rng (), a, sigma);
 855
 856 // Studentized maximum modulus distribution.
 857 function CDF.SMOD (x > 0, a >= 1, b >= 1) = unimplemented;
 858 function IDF.SMOD (P >= 0 && P < 1, a >= 1, b >= 1) = unimplemented;
 859
 860 // Studentized range distribution.
 861 function CDF.SRANGE (x > 0, a >= 1, b >= 1) = unimplemented;
 862 function IDF.SRANGE (P >= 0 && P < 1, a >= 1, b >= 1) = unimplemented;
 863
 864 // Student t distribution.
 865 function CDF.T (x, df > 0) = gsl_cdf_tdist_P (x, df);
 866 function IDF.T (P > 0 && P < 1, df > 0) = gsl_cdf_tdist_Pinv (P, df);
 867 function PDF.T (x, df > 0) = gsl_ran_tdist_pdf (x, df);
 868 no_opt function RV.T (df > 0) = gsl_ran_tdist (get_rng (), df);
 869 function NCDF.T (x, df > 0, nc) = unimplemented;
 870 function NPDF.T (x, df > 0, nc) = unimplemented;
 871
 872 // Type-1 Gumbel distribution.
 873 extension function CDF.T1G (x, a, b) = gsl_cdf_gumbel1_P (x, a, b);
 874 extension function IDF.T1G (P >= 0 && P <= 1, a, b)
 875      = gsl_cdf_gumbel1_Pinv (P, a, b);
 876 extension function PDF.T1G (x, a, b) = gsl_ran_gumbel1_pdf (x, a, b);
 877 no_opt extension function RV.T1G (a, b) = gsl_ran_gumbel1 (get_rng (), a, b);
 878
 879 // Type-2 Gumbel distribution.
 880 extension function CDF.T2G (x, a, b) = gsl_cdf_gumbel2_P (x, a, b);
 881 extension function IDF.T2G (P >= 0 && P <= 1, a, b)
 882      = gsl_cdf_gumbel2_Pinv (P, a, b);
 883 extension function PDF.T2G (x, a, b) = gsl_ran_gumbel2_pdf (x, a, b);
 884 no_opt extension function RV.T2G (a, b) = gsl_ran_gumbel2 (get_rng (), a, b);
 885
 886 // Uniform distribution.
 887 function CDF.UNIFORM (x <= b, a <= x, b) = gsl_cdf_flat_P (x, a, b);
 888 function IDF.UNIFORM (P >= 0 && P <= 1, a <= b, b)
 889      = gsl_cdf_flat_Pinv (P, a, b);
 890 function PDF.UNIFORM (x <= b, a <= x, b) = gsl_ran_flat_pdf (x, a, b);
 891 no_opt function RV.UNIFORM (a <= b, b) = gsl_ran_flat (get_rng (), a, b);
 892 no_opt function UNIFORM (b >= 0) = gsl_ran_flat (get_rng (), 0, b);
 893
 894 // Weibull distribution.
 895 function CDF.WEIBULL (x >= 0, a > 0, b > 0) = gsl_cdf_weibull_P (x, a, b);
 896 function IDF.WEIBULL (P >= 0 && P < 1, a > 0, b > 0)
 897      = gsl_cdf_weibull_Pinv (P, a, b);
 898 function PDF.WEIBULL (x >= 0, a > 0, b > 0) = gsl_ran_weibull_pdf (x, a, b);
 899 no_opt function RV.WEIBULL (a > 0, b > 0) = gsl_ran_weibull (get_rng (), a, b);
 900
 901 // Bernoulli distribution.
 902 function CDF.BERNOULLI (k == 0 || k == 1, p >= 0 && p <= 1)
 903      = k ? 1 : 1 - p;
 904 function PDF.BERNOULLI (k == 0 || k == 1, p >= 0 && p <= 1)
 905      = gsl_ran_bernoulli_pdf (k, p);
 906 no_opt function RV.BERNOULLI (p >= 0 && p <= 1)
 907      = gsl_ran_bernoulli (get_rng (), p);
 908
 909 // Binomial distribution.
 910 function CDF.BINOM (k, n > 0 && n == floor (n), p >= 0 && p <= 1)
 911      = gsl_cdf_binomial_P (k, p, n);
 912 function PDF.BINOM (k >= 0 && k == floor (k) && k <= n,
 913                     n > 0 && n == floor (n),
 914                     p >= 0 && p <= 1)
 915      = gsl_ran_binomial_pdf (k, p, n);
 916 no_opt function RV.BINOM (p > 0 && p == floor (p), n >= 0 && n <= 1)
 917      = gsl_ran_binomial (get_rng (), p, n);
 918
 919 // Geometric distribution.
 920 function CDF.GEOM (k >= 1 && k == floor (k), p >= 0 && p <= 1)
 921      = gsl_cdf_geometric_P (k, p);
 922 function PDF.GEOM (k >= 1 && k == floor (k),
 923                    p >= 0 && p <= 1)
 924      = gsl_ran_geometric_pdf (k, p);
 925 no_opt function RV.GEOM (p >= 0 && p <= 1) = gsl_ran_geometric (get_rng (), p);
 926
 927 // Hypergeometric distribution.
 928 function CDF.HYPER (k >= 0 && k == floor (k) && k <= c,
 929                     a > 0 && a == floor (a),
 930                     b > 0 && b == floor (b) && b <= a,
 931                     c > 0 && c == floor (c) && c <= a)
 932      = gsl_cdf_hypergeometric_P (k, c, a - c, b);
 933 function PDF.HYPER (k >= 0 && k == floor (k) && k <= c,
 934                     a > 0 && a == floor (a),
 935                     b > 0 && b == floor (b) && b <= a,
 936                     c > 0 && c == floor (c) && c <= a)
 937      = gsl_ran_hypergeometric_pdf (k, c, a - c, b);
 938 no_opt function RV.HYPER (a > 0 && a == floor (a),
 939                           b > 0 && b == floor (b) && b <= a,
 940                           c > 0 && c == floor (c) && c <= a)
 941      = gsl_ran_hypergeometric (get_rng (), c, a - c, b);
 942
 943 // Logarithmic distribution.
 944 extension function PDF.LOG (k >= 1, p > 0 && p <= 1)
 945      = gsl_ran_logarithmic_pdf (k, p);
 946 no_opt extension function RV.LOG (p > 0 && p <= 1)
 947      = gsl_ran_logarithmic (get_rng (), p);
 948
 949 // Negative binomial distribution.
 950 function CDF.NEGBIN (k >= 1, n == floor (n), p > 0 && p <= 1)
 951      = gsl_cdf_negative_binomial_P (k, p, n);
 952 function PDF.NEGBIN (k >= 1, n == floor (n), p > 0 && p <= 1)
 953      = gsl_ran_negative_binomial_pdf (k, p, n);
 954 no_opt function RV.NEGBIN (n == floor (n), p > 0 && p <= 1)
 955      = gsl_ran_negative_binomial (get_rng (), p, n);
 956
 957 // Poisson distribution.
 958 function CDF.POISSON (k >= 0 && k == floor (k), mu > 0)
 959      = gsl_cdf_poisson_P (k, mu);
 960 function PDF.POISSON (k >= 0 && k == floor (k), mu > 0)
 961      = gsl_ran_poisson_pdf (k, mu);
 962 no_opt function RV.POISSON (mu > 0) = gsl_ran_poisson (get_rng (), mu);
 963
 964 // Weirdness.
 965 absorb_miss boolean function MISSING (x) = x == SYSMIS || !isfinite (x);
 966 absorb_miss boolean function SYSMIS (x) = x == SYSMIS || !isfinite (x);
 967 no_opt boolean function SYSMIS (num_var v)
 968      case c;
 969 {
 970   return case_num (c, v) == SYSMIS;
 971 }
 972 no_opt boolean function VALUE (num_var v)
 973      case c;
 974 {
 975   return case_num (c, v);
 976 }
 977
 978 no_opt operator VEC_ELEM_NUM (idx)
 979      vector v;
 980      case c;
 981 {
 982   if (idx >= 1 && idx <= vector_get_var_cnt (v))
 983     {
 984       const struct variable *var = vector_get_var (v, (size_t) idx - 1);
 985       double value = case_num (c, var);
 986       return !var_is_num_missing (var, value, MV_USER) ? value : SYSMIS;
 987     }
 988   else
 989     {
 990       if (idx == SYSMIS)
 991         msg (SE, _("SYSMIS is not a valid index value for vector "
 992                    "%s.  The result will be set to SYSMIS."),
 993              vector_get_name (v));
 994       else
 995         msg (SE, _("%g is not a valid index value for vector %s.  "
 996                    "The result will be set to SYSMIS."),
 997              idx, vector_get_name (v));
 998       return SYSMIS;
 999     }
1000 }
1001
1002 absorb_miss no_opt string operator VEC_ELEM_STR (idx)
1003      expression e;
1004      vector v;
1005      case c;
1006 {
1007   if (idx >= 1 && idx <= vector_get_var_cnt (v))
1008     {
1009       struct variable *var = vector_get_var (v, (size_t) idx - 1);
1010       return copy_string (e, CHAR_CAST_BUG (char *, case_str (c, var)),
1011                           var_get_width (var));
1012     }
1013   else
1014     {
1015       if (idx == SYSMIS)
1016         msg (SE, _("SYSMIS is not a valid index value for vector "
1017                    "%s.  The result will be set to the empty string."),
1018              vector_get_name (v));
1019       else
1020         msg (SE, _("%g is not a valid index value for vector %s.  "
1021                    "The result will be set to the empty string."),
1022              idx, vector_get_name (v));
1023       return empty_string;
1024     }
1025 }
1026
1027 // Terminals.
1028
1029 no_opt operator NUM_VAR ()
1030      case c;
1031      num_var v;
1032 {
1033   double d = case_num (c, v);
1034   return !var_is_num_missing (v, d, MV_USER) ? d : SYSMIS;
1035 }
1036
1037 no_opt string operator STR_VAR ()
1038      case c;
1039      expression e;
1040      str_var v;
1041 {
1042   struct substring s = alloc_string (e, var_get_width (v));
1043   memcpy (s.string, case_str (c, v), var_get_width (v));
1044   return s;
1045 }
1046
1047 no_opt perm_only function LAG (num_var v, pos_int n_before)
1048     dataset ds;
1049 {
1050   const struct ccase *c = lagged_case (ds, n_before);
1051   if (c != NULL)
1052     {
1053       double x = case_num (c, v);
1054       return !var_is_num_missing (v, x, MV_USER) ? x : SYSMIS;
1055     }
1056   else
1057     return SYSMIS;
1058 }
1059
1060 no_opt perm_only function LAG (num_var v)
1061     dataset ds;
1062 {
1063   const struct ccase *c = lagged_case (ds, 1);
1064   if (c != NULL)
1065     {
1066       double x = case_num (c, v);
1067       return !var_is_num_missing (v, x, MV_USER) ? x : SYSMIS;
1068     }
1069   else
1070     return SYSMIS;
1071 }
1072
1073 no_opt perm_only string function LAG (str_var v, pos_int n_before)
1074      expression e;
1075      dataset ds;
1076 {
1077   const struct ccase *c = lagged_case (ds, n_before);
1078   if (c != NULL)
1079     return copy_string (e, CHAR_CAST_BUG (char *, case_str (c, v)),
1080                         var_get_width (v));
1081   else
1082     return empty_string;
1083 }
1084
1085 no_opt perm_only string function LAG (str_var v)
1086      expression e;
1087      dataset ds;
1088 {
1089   const struct ccase *c = lagged_case (ds, 1);
1090   if (c != NULL)
1091     return copy_string (e, CHAR_CAST_BUG (char *, case_str (c, v)),
1092                         var_get_width (v));
1093   else
1094     return empty_string;
1095 }
1096
1097 no_opt operator NUM_SYS ()
1098      case c;
1099      num_var v;
1100 {
1101   return case_num (c, v) == SYSMIS;
1102 }
1103
1104 no_opt operator NUM_VAL ()
1105      case c;
1106      num_var v;
1107 {
1108   return case_num (c, v);
1109 }
1110
1111 no_opt operator CASENUM ()
1112      case_idx idx;
1113 {
1114   return idx;
1115 }