src/language/expressions/operations.def

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