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