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