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