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