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