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