src/language/expressions/operations.def

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