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