// -*- c -*-
//
// PSPP - a program for statistical analysis.
// Copyright (C) 2005, 2006, 2009, 2010, 2011, 2012, 2015, 2016 Free Software Foundation, Inc.
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program.  If not, see <http://www.gnu.org/licenses/>.

operator NEG (x) = -x;

operator ADD (a, b) = a + b;
operator SUB (a, b) = a - b;

absorb_miss operator MUL (a, b)
= (a == 0. || b == 0. ? 0.
   : a == SYSMIS || b == SYSMIS ? SYSMIS
   : a * b);

absorb_miss operator DIV (a, b)
= (a == 0. ? 0.
   : a == SYSMIS || b == SYSMIS ? SYSMIS
   : a / b);

absorb_miss operator POW (a, b)
= (a == SYSMIS ? (b == 0. ? 1. : a)
   : b == SYSMIS ? (a == 0. ? 0. : SYSMIS)
   : a == 0. && b <= 0. ? SYSMIS
   : pow (a, b));

absorb_miss boolean operator AND (boolean a, boolean b)
= (a == 0. ? 0.
   : b == 0. ? 0.
   : b == SYSMIS ? SYSMIS
   : a);

absorb_miss boolean operator OR (boolean a, boolean b)
= (a == 1. ? 1.
   : b == 1. ? 1.
   : b == SYSMIS ? SYSMIS
   : a);

boolean operator NOT (boolean a)
= (a == 0. ? 1.
   : a == 1. ? 0.
   : SYSMIS);

// Numeric relational operators.
boolean operator EQ (a, b) = a == b;
boolean operator GE (a, b) = a >= b;
boolean operator GT (a, b) = a > b;
boolean operator LE (a, b) = a <= b;
boolean operator LT (a, b) = a < b;
boolean operator NE (a, b) = a != b;

// String relational operators.
boolean operator EQ_STRING (string a, string b) = compare_string_3way (&a, &b) == 0;
boolean operator GE_STRING (string a, string b) = compare_string_3way (&a, &b) >= 0;
boolean operator GT_STRING (string a, string b) = compare_string_3way (&a, &b) > 0;
boolean operator LE_STRING (string a, string b) = compare_string_3way (&a, &b) <= 0;
boolean operator LT_STRING (string a, string b) = compare_string_3way (&a, &b) < 0;
boolean operator NE_STRING (string a, string b) = compare_string_3way (&a, &b) != 0;

// Unary functions.
function ABS (x) = fabs (x);
extension function ACOS (x >= -1 && x <= 1) = acos (x);
function ASIN (x >= -1 && x <= 1) = asin (x);
function ATAN (x) = atan (x);
extension function ARCOS (x >= -1 && x <= 1) = acos (x);
function ARSIN (x >= -1 && x <= 1) = asin (x);
function ARTAN (x) = atan (x);
function COS (x) = cos (x);
function EXP (x) = check_errno (exp (x));
function LG10(x) = check_errno (log10 (x));
function LN (x) = check_errno (log (x));
function LNGAMMA (x >= 0) = gsl_sf_lngamma (x);
function MOD10 (x) = fmod (x, 10);
function RND (x) = round_nearest (x, 1, 0);
function RND (x, mult != 0) = round_nearest (x, mult, 0);
function RND (x, mult != 0, fuzzbits >= 0) = round_nearest (x, mult, fuzzbits);
function SIN (x) = sin (x);
function SQRT (x >= 0) = sqrt (x);
function TAN (x) = check_errno (tan (x));
function TRUNC (x) = round_zero (x, 1, 0);
function TRUNC (x, mult != 0) = round_zero (x, mult, 0);
function TRUNC (x, mult != 0, fuzzbits >= 0) = round_zero (x, mult, fuzzbits);

absorb_miss function MOD (n, d)
{
  if (d != SYSMIS)
    return n != SYSMIS ? fmod (n, d) : SYSMIS;
  else
    return n != 0. ? SYSMIS : 0.;
}

// N-ary numeric functions.
absorb_miss boolean function ANY (x, a[n])
{
  double retval = SYSMIS;
  if (x != SYSMIS)
    {
      for (size_t i = 0; i < n; i++)
        if (a[i] == x)
          return 1.;
        else if (a[i] != SYSMIS)
          retval = 0.;
    }
  return retval;
}

boolean function ANY (string x, string a[n])
{
  size_t i;

  for (i = 0; i < n; i++)
    if (!compare_string_3way (&x, &a[i]))
      return 1.;
  return 0.;
}

function CFVAR.2 (a[n])
{
  double mean, variance;

  moments_of_doubles (a, n, NULL, &mean, &variance, NULL, NULL);

  if (mean == SYSMIS || mean == 0 || variance == SYSMIS)
    return SYSMIS;
  else
    return sqrt (variance) / mean;
}

function MAX.1 (a[n])
{
  double max;
  size_t i;

  max = -DBL_MAX;
  for (i = 0; i < n; i++)
    if (a[i] != SYSMIS && a[i] > max)
      max = a[i];
  return max;
}

string function MAX (string a[n])
{
  struct substring *max;
  size_t i;

  max = &a[0];
  for (i = 1; i < n; i++)
    if (compare_string_3way (&a[i], max) > 0)
      max = &a[i];
  return *max;
}

function MEAN.1 (a[n])
{
  double mean;
  moments_of_doubles (a, n, NULL, &mean, NULL, NULL, NULL);
  return mean;
}

function MEDIAN.1 (a[n])
{
  return median (a, n);
}

function MIN.1 (a[n])
{
  double min;
  size_t i;

  min = DBL_MAX;
  for (i = 0; i < n; i++)
    if (a[i] != SYSMIS && a[i] < min)
      min = a[i];
  return min;
}

string function MIN (string a[n])
{
  struct substring *min;
  size_t i;

  min = &a[0];
  for (i = 1; i < n; i++)
    if (compare_string_3way (&a[i], min) < 0)
      min = &a[i];
  return *min;
}

absorb_miss function NMISS (a[n])
{
  size_t n_missings = 0;
  for (size_t i = 0; i < n; i++)
    n_missings += a[i] == SYSMIS;
  return n_missings;
}

absorb_miss function NVALID (a[n])
{
  size_t n_valids = 0;
  for (size_t i = 0; i < n; i++)
    n_valids += a[i] != SYSMIS;
  return n_valids;
}

absorb_miss boolean function RANGE (x != SYSMIS, a[n*2])
{
  bool found = false;
  bool valid = false;
  for (size_t i = 0; i < n; i++)
    {
      double w = a[2 * i];
      double y = a[2 * i + 1];
      if (w != SYSMIS && y != SYSMIS)
        {
          if (w <= x && x <= y)
            found = true;
          else if (w <= y)
            valid = true;
          else
            return SYSMIS;
        }
    }
  return found ? true : valid ? false : SYSMIS;
}

boolean function RANGE (string x, string a[n*2])
{
  bool found = false;
  for (size_t i = 0; i < n; i++)
    {
      struct substring *w = &a[2 * i];
      struct substring *y = &a[2 * i + 1];
      if (compare_string_3way (w, &x) <= 0 && compare_string_3way (&x, y) <= 0)
        found = true;
      else if (compare_string_3way (w, y) > 0)
        return SYSMIS;
    }
  return found;
}

function SD.2 (a[n])
{
  double variance;
  moments_of_doubles (a, n, NULL, NULL, &variance, NULL, NULL);
  return sqrt (variance);
}

function SUM.1 (a[n])
{
  double sum;
  size_t i;

  sum = 0.;
  for (i = 0; i < n; i++)
    if (a[i] != SYSMIS)
      sum += a[i];
  return sum;
}

function VARIANCE.2 (a[n])
{
  double variance;
  moments_of_doubles (a, n, NULL, NULL, &variance, NULL, NULL);
  return variance;
}

// Time construction & extraction functions.
function TIME.HMS (h, m, s)
  expression e;
  expr_node n;
{
  if ((h > 0. || m > 0. || s > 0.) && (h < 0. || m < 0. || s < 0.))
    {
      msg_at (SW, expr_location (e, n),
              _("TIME.HMS cannot accept a mix of positive and negative "
                "arguments."));
      double args[] = { h, m, s };
      for (size_t i = 0; i < 3; i++)
        if (args[i] > 0)
          msg_at (SN, expr_location (e, n->args[i]),
                  _("This argument has positive value %g."), args[i]);
        else if (args[i] < 0)
          msg_at (SN, expr_location (e, n->args[i]),
                  _("This argument has negative value %g."), args[i]);
      return SYSMIS;
    }
  else
    return H_S * h + MIN_S * m + s;
}
function TIME.DAYS (days) = days * DAY_S;
function CTIME.DAYS (time) = time / DAY_S;
function CTIME.HOURS (time) = time / H_S;
function CTIME.MINUTES (time) = time / MIN_S;
function CTIME.SECONDS (time) = time;

// Date construction functions.
function DATE.DMY (integer d, integer m, integer y)
  expression e;
  expr_node n;
= expr_ymd_to_date (y, m, d, e, n, 3, 2, 1);

function DATE.MDY (integer m, integer d, integer y)
  expression e;
  expr_node n;
= expr_ymd_to_date (y, m, d, e, n, 3, 1, 2);

function DATE.MOYR (integer m, integer y)
  expression e;
  expr_node n;
= expr_ymd_to_date (y, m, 1, e, n, 2, 1, 0);

function DATE.QYR (integer q, integer y)
  expression e;
  expr_node n;
{
  if (q < 1 || q > 4)
    {
      msg_at (SW, expr_location (e, n->args[0]),
              _("Argument 1 to DATE.QYR must be 1, 2, 3, or 4 (not %d)."), q);
      return SYSMIS;
    }
  return expr_ymd_to_date (y, q * 3 - 2, 1, e, n, 2, 0, 0);
}

function DATE.WKYR (integer w, integer y)
  expression e;
  expr_node n;
{
  if (w < 1 || w > 53)
    {
      msg_at (SE, expr_location (e, n->args[0]),
              _("The week argument to DATE.WKYR is outside the acceptable "
                "range of 1 to 53.  The result will be system-missing."));
      return SYSMIS;
    }
  else
    {
      double yr_1_1 = expr_ymd_to_ofs (y, 1, 1, e, n, 2, 0, 0);
      if (yr_1_1 != SYSMIS)
        return DAY_S * (yr_1_1 + WEEK_DAY * (w - 1));
      else
        return SYSMIS;
    }
}

function DATE.YRDAY (integer y, integer yd)
  expression e;
  expr_node n;
{
  if (yd < 1 || yd > 366)
    {
      msg_at (SE, expr_location (e, n->args[1]),
              _("The value %d as day argument to DATE.YRDAY is outside the "
                "acceptable range of 1 to 366.  "
                "The result will be system-missing."), yd);
      return SYSMIS;
    }
  else
    {
      double yr_1_1 = expr_ymd_to_ofs (y, 1, 1, e, n, 1, 0, 0);
      if (yr_1_1 != SYSMIS)
        return DAY_S * (yr_1_1 + yd - 1.);
      else
        return SYSMIS;
    }
}

function YRMODA (integer y, integer m, integer d)
  expression e;
  expr_node n;
{
  if (y >= 0 && y <= 99)
    y += 1900;
  else if (y > 47516)
    {
      msg_at (SE, expr_location (e, n->args[0]),
              _("The year argument to YRMODA is greater than 47516.  "
                "The result will be system-missing."));
      return SYSMIS;
    }

  return expr_ymd_to_ofs (y, m, d, e, n, 1, 2, 3);
}

// Date extraction functions.
function XDATE.TDAY (date) = floor (date / DAY_S);
function XDATE.HOUR (date) = fmod (floor (date / H_S), DAY_H);
function XDATE.MINUTE (date) = fmod (floor (date / H_MIN), H_MIN);
function XDATE.SECOND (date) = fmod (date, MIN_S);
function XDATE.DATE (date) = floor (date / DAY_S) * DAY_S;
function XDATE.TIME (date) = fmod (date, DAY_S);

function XDATE.JDAY (date >= DAY_S) = calendar_offset_to_yday (date / DAY_S);
function XDATE.MDAY (date >= DAY_S) = calendar_offset_to_mday (date / DAY_S);
function XDATE.MONTH (date >= DAY_S)
     = calendar_offset_to_month (date / DAY_S);
function XDATE.QUARTER (date >= DAY_S)
    = (calendar_offset_to_month (date / DAY_S) - 1) / 3 + 1;
function XDATE.WEEK (date >= DAY_S)
    = (calendar_offset_to_yday (date / DAY_S) - 1) / 7 + 1;
function XDATE.WKDAY (date >= DAY_S) = calendar_offset_to_wday (date / DAY_S);
function XDATE.YEAR (date >= DAY_S) = calendar_offset_to_year (date / DAY_S);

// Date arithmetic functions.
no_abbrev function DATEDIFF (date2 >= DAY_S, date1 >= DAY_S, string unit)
  expression e;
  expr_node n;
= expr_date_difference (date1, date2, unit, e, n);

no_abbrev function DATESUM (date, quantity, string unit)
  expression e;
  expr_node n;
= expr_date_sum_closest (date, quantity, unit, e, n);
no_abbrev function DATESUM (date, quantity, string unit, string method)
  expression e;
  expr_node n;
= expr_date_sum (date, quantity, unit, method, e, n);


// String functions.
string function CONCAT (string a[n])
     expression e;
{
  struct substring dst;
  size_t i;

  dst = alloc_string (e, MAX_STRING);
  dst.length = 0;
  for (i = 0; i < n; i++)
    {
      struct substring *src = &a[i];
      size_t copy_len;

      copy_len = src->length;
      if (dst.length + copy_len > MAX_STRING)
        copy_len = MAX_STRING - dst.length;
      memcpy (&dst.string[dst.length], src->string, copy_len);
      dst.length += copy_len;
    }

  return dst;
}

function INDEX (string haystack, string needle)
{
  if (haystack.length >= needle.length)
    {
      size_t limit = haystack.length - needle.length + 1;
      for (size_t i = 1; i <= limit; i++)
        if (!memcmp (&haystack.string[i - 1], needle.string, needle.length))
          return i;
    }
  return 0;
}

function INDEX (string haystack, string needles, integer needle_len)
  expression e;
  expr_node n;
{
  if (needle_len <= 0 || needles.length % needle_len != 0)
    {
      msg_at (SE, expr_location (e, n),
              _("INDEX needle length argument must evenly divide the "
                "length of the needles argument."));
      msg_at (SN, expr_location (e, n->args[1]),
              _("The needles argument has length %zu."), needles.length);
      msg_at (SN, expr_location (e, n->args[2]),
              _("The needle length argument has value %d."), needle_len);
      return SYSMIS;
    }

  if (haystack.length >= needle_len)
    {
      size_t limit = haystack.length - needle_len + 1;
      for (size_t i = 1; i <= limit; i++)
        for (size_t j = 0; j < needles.length; j += needle_len)
          if (!memcmp (&haystack.string[i - 1], &needles.string[j], needle_len))
            return i;
    }

  return 0;
}

function RINDEX (string haystack, string needle)
{
  if (haystack.length >= needle.length)
    {
      size_t limit = haystack.length - needle.length + 1;
      for (size_t i = limit; i >= 1; i--)
        if (!memcmp (&haystack.string[i - 1], needle.string, needle.length))
          return i;
    }

  return 0;
}

function RINDEX (string haystack, string needles, integer needle_len)
  expression e;
  expr_node n;
{
  if (needle_len <= 0 || needles.length % needle_len != 0)
    {
      msg_at (SE, expr_location (e, n),
              _("RINDEX needle length argument must evenly divide the "
                "length of the needles argument."));
      msg_at (SN, expr_location (e, n->args[1]),
              _("The needles argument has length %zu."), needles.length);
      msg_at (SN, expr_location (e, n->args[2]),
              _("The needle length argument has value %d."), needle_len);
      return SYSMIS;
    }

  if (haystack.length >= needle_len)
    {
      size_t limit = haystack.length - needle_len + 1;
      for (size_t i = limit; i >= 1; i--)
        for (size_t j = 0; j < needles.length; j += needle_len)
          if (!memcmp (&haystack.string[i - 1], &needles.string[j], needle_len))
            return i;
    }

  return 0;
}

function LENGTH (string s)
{
  return s.length;
}

string function LOWER (string s)
{
  int i;

  for (i = 0; i < s.length; i++)
    s.string[i] = tolower ((unsigned char) s.string[i]);
  return s;
}

function MBLEN.BYTE (string s, idx)
{
  if (idx < 0 || idx >= s.length || (int) idx != idx)
    return SYSMIS;
  else
    return 1;
}

string function UPCASE (string s)
{
  int i;

  for (i = 0; i < s.length; i++)
    s.string[i] = toupper ((unsigned char) s.string[i]);
  return s;
}

absorb_miss string function LPAD (string s, integer n)
     expression e;
     expr_node node;
{
  if (n < 0 || n > MAX_STRING)
    {
      if (n != INT_MIN)
        {
          msg_at (SE, expr_location (e, node),
                  _("The length argument to LPAD must be between 0 and %d."),
                  MAX_STRING);
          msg_at (SN, expr_location (e, node->args[1]),
                  _("The length argument is %d."), n);
        }

      return s;
    }
  else if (s.length >= n)
    return s;
  else
    {
      struct substring t = alloc_string (e, n);
      size_t pad = n - s.length;
      memset (t.string, ' ', pad);
      memcpy (&t.string[pad], s.string, s.length);
      return t;
    }
}

absorb_miss string function LPAD (string s, integer n, string c)
     expression e;
     expr_node node;
{
  if (n < 0 || n > MAX_STRING)
    {
      if (n != INT_MIN)
        {
          msg_at (SE, expr_location (e, node),
                  _("The length argument to LPAD must be between 0 and %d."),
                  MAX_STRING);
          msg_at (SN, expr_location (e, node->args[1]),
                  _("The length argument is %d."), n);
        }

      return s;
    }
  else if (s.length >= n)
    return s;
  else if (c.length == 0)
    {
      msg_at (SE, expr_location (e, node),
              _("The padding argument to LPAD must not be an empty string."));
      return s;
    }
  else
    {
      size_t n_pad = (n - s.length) / c.length;
      if (!n_pad)
        return s;

      struct substring t = alloc_string (e, n);
      t.length = 0;
      for (size_t i = 0; i < n_pad; i++)
        {
          memcpy (t.string + t.length, c.string, c.length);
          t.length += c.length;
        }
      memcpy (t.string + t.length, s.string, s.length);
      t.length += s.length;
      return t;
    }
}

string function REPLACE (string haystack, string needle, string replacement)
    expression e;
  = replace_string (e, haystack, needle, replacement, INT_MAX);

absorb_miss string function REPLACE (string haystack, string needle,
                         string replacement, integer n)
    expression e;
  = replace_string (e, haystack, needle, replacement, n);

absorb_miss string function RPAD (string s, integer n)
     expression e;
     expr_node node;
{
  if (n < 0 || n > MAX_STRING)
    {
      if (n != INT_MIN)
        {
          msg_at (SE, expr_location (e, node),
                  _("The length argument to RPAD must be between 0 and %d."),
                  MAX_STRING);
          msg_at (SN, expr_location (e, node->args[1]),
                  _("The length argument is %d."), n);
        }

      return s;
    }
  else if (s.length >= n)
    return s;
  else
    {
      struct substring t = alloc_string (e, n);
      size_t pad = n - s.length;
      memcpy (t.string, s.string, s.length);
      memset (t.string + s.length, ' ', pad);
      return t;
    }
}

absorb_miss string function RPAD (string s, integer n, string c)
     expression e;
     expr_node node;
{
  if (n < 0 || n > MAX_STRING)
    {
      if (n != INT_MIN)
        {
          msg_at (SE, expr_location (e, node),
                  _("The length argument to RPAD must be between 0 and %d."),
                  MAX_STRING);
          msg_at (SN, expr_location (e, node->args[1]),
                  _("The length argument is %d."), n);
        }

      return s;
    }
  else if (s.length >= n)
    return s;
  else if (c.length == 0)
    {
      msg_at (SE, expr_location (e, node),
              _("The padding argument to RPAD must not be an empty string."));
      return s;
    }
  else
    {
      size_t n_pad = (n - s.length) / c.length;
      if (!n_pad)
        return s;

      struct substring t = alloc_string (e, n);
      memcpy (t.string, s.string, s.length);
      t.length = s.length;
      for (size_t i = 0; i < n_pad; i++)
        {
          memcpy (t.string + t.length, c.string, c.length);
          t.length += c.length;
        }
      return t;
    }
}

string function LTRIM (string s)
{
  while (s.length > 0 && s.string[0] == ' ')
    {
      s.length--;
      s.string++;
    }
  return s;
}

string function LTRIM (string s, string c)
{
  if (c.length > 0)
    while (s.length >= c.length && !memcmp (s.string, c.string, c.length))
      {
        s.length -= c.length;
        s.string += c.length;
      }
  return s;
}

string function RTRIM (string s)
{
  while (s.length > 0 && s.string[s.length - 1] == ' ')
    s.length--;
  return s;
}

string function RTRIM (string s, string c)
{
  if (c.length > 0)
    while (s.length >= c.length
           && !memcmp (&s.string[s.length - c.length], c.string, c.length))
      s.length -= c.length;
  return s;
}

function NUMBER (string s, ni_format f)
  expression e;
  expr_node n;
{
  if (s.length > f->w)
    s.length = f->w;

  union value out;
  char *error = data_in (s, C_ENCODING, f->type, settings_get_fmt_settings (),
                         &out, 0, NULL);
  if (error == NULL)
    data_in_imply_decimals (s, C_ENCODING, f->type, f->d,
                            settings_get_fmt_settings (), &out);
  else
    {
      msg_at (SE, expr_location (e, n->args[0]),
              _("Cannot parse \"%.*s\" as format %s: %s"),
              (int) s.length, s.string, fmt_name (f->type), error);
      free (error);
    }
  return out.f;
}

absorb_miss string function STRING (x, no_format f)
     expression e;
{
  union value v;
  struct substring dst;
  char *s;

  v.f = x;

  assert (!fmt_is_string (f->type));
  s = data_out (&v, C_ENCODING, f, settings_get_fmt_settings ());
  dst = alloc_string (e, strlen (s));
  strcpy (dst.string, s);
  free (s);
  return dst;
}

absorb_miss string function STRUNC (string s, integer n)
{
  if (n < 1)
    return n == INT_MIN ? s : empty_string;

  if (n < s.length)
    s.length = n;
  while (s.length > 0 && s.string[s.length - 1] == ' ')
    s.length--;
  return s;
}

absorb_miss string function SUBSTR (string s, integer ofs)
{
  return (ofs >= 1 && ofs <= s.length
          ? ss_substr (s, ofs - 1, SIZE_MAX)
          : empty_string);
}

absorb_miss string function SUBSTR (string s, integer ofs, integer cnt)
{
  return (ofs >= 1 && cnt >= 1
          ? ss_substr (s, ofs - 1, cnt)
          : empty_string);
}

absorb_miss no_opt no_abbrev string function VALUELABEL (var v)
     expression e;
     case c;
{
  const char *label = var_lookup_value_label (v, case_data (c, v));
  if (label != NULL)
    return copy_string (e, label, strlen (label));
  else
    return empty_string;
}

// Artificial.
operator SQUARE (x) = x * x;

absorb_miss boolean operator OPERAND_TO_BOOLEAN (x, expr_node parent)
  expression e;
  expr_node n;
{
  if (x == 0. || x == 1. || x == SYSMIS)
    return x;

  switch (parent->n_args)
    {
    case 2:
      msg_at (SE, expr_location (e, parent),
              /* TRANSLATORS: There are exactly two operands. */
              _("The operands of %s must have value 0 or 1."),
              operations[parent->type].name);
      break;

    case 1:
      msg_at (SE, expr_location (e, parent),
              _("The operand of %s must have value 0 or 1."),
              operations[parent->type].name);
      break;

    default:
      NOT_REACHED ();
    }

  msg_at (SN, expr_location (e, n),
          _("This operand with unexpected value %g will be treated as 0."), x);
  return 0.;
}

absorb_miss boolean operator EXPR_TO_BOOLEAN (x)
  expression e;
  expr_node n;
{
  if (x == 0. || x == 1. || x == SYSMIS)
    return x;

  msg_at (SE, expr_location (e, n),
          _("This expression, which must be 0 or 1, evaluated to %g.  "
            "It will be treated as 0."), x);
  return 0.;
}

operator NUM_TO_INTEGER (x)
  expression e;
  expr_node n;
{
  if (x == floor (x) && x > INT_MIN && x <= INT_MAX)
    return x;

  msg_at (SE, expr_location (e, n),
          _("Treating unexpected non-integer value %g as missing."), x);
  return SYSMIS;
}

operator BOOLEAN_TO_NUM (boolean x) = x;

// Beta distribution.
function PDF.BETA (x >= 0 && x <= 1, a > 0, b > 0)
     = gsl_ran_beta_pdf (x, a, b);
function CDF.BETA (x >= 0 && x <= 1, a > 0, b > 0) = gsl_cdf_beta_P (x, a, b);
function IDF.BETA (P >= 0 && P <= 1, a > 0, b > 0)
     = gsl_cdf_beta_Pinv (P, a, b);
no_opt function RV.BETA (a > 0, b > 0) = gsl_ran_beta (get_rng (), a, b);
function NCDF.BETA (x >= 0, a > 0, b > 0, lambda > 0)
     = ncdf_beta (x, a, b, lambda);
function NPDF.BETA (x >= 0, a > 0, b > 0, lambda > 0)
     = npdf_beta (x, a, b, lambda);

// Bivariate normal distribution.
function CDF.BVNOR (x0, x1, r >= -1 && r <= 1) = cdf_bvnor (x0, x1, r);
function PDF.BVNOR (x0, x1, r >= -1 && r <= 1)
     = gsl_ran_bivariate_gaussian_pdf (x0, x1, 1, 1, r);

// Cauchy distribution.
function CDF.CAUCHY (x, a, b > 0) = gsl_cdf_cauchy_P ((x - a) / b, 1);
function IDF.CAUCHY (P > 0 && P < 1, a, b > 0)
     = a + b * gsl_cdf_cauchy_Pinv (P, 1);
function PDF.CAUCHY (x, a, b > 0) = gsl_ran_cauchy_pdf ((x - a) / b, 1) / b;
no_opt function RV.CAUCHY (a, b > 0) = a + b * gsl_ran_cauchy (get_rng (), 1);

// Chi-square distribution.
function CDF.CHISQ (x >= 0, df > 0) = gsl_cdf_chisq_P (x, df);
function IDF.CHISQ (P >= 0 && P < 1, df > 0) = gsl_cdf_chisq_Pinv (P, df);
function PDF.CHISQ (x >= 0, df > 0) = gsl_ran_chisq_pdf (x, df);
no_opt function RV.CHISQ (df > 0) = gsl_ran_chisq (get_rng (), df);
function NCDF.CHISQ (x >= 0, df > 0, c) = unimplemented;
function NPDF.CHISQ (x >= 0, df > 0, c) = unimplemented;
function SIG.CHISQ (x >= 0, df > 0) = gsl_cdf_chisq_Q (x, df);

// Exponential distribution.
function CDF.EXP (x >= 0, a > 0) = gsl_cdf_exponential_P (x, 1. / a);
function IDF.EXP (P >= 0 && P < 1, a > 0)
     = gsl_cdf_exponential_Pinv (P, 1. / a);
function PDF.EXP (x >= 0, a > 0) = gsl_ran_exponential_pdf (x, 1. / a);
no_opt function RV.EXP (a > 0) = gsl_ran_exponential (get_rng (), 1. / a);

// Exponential power distribution.
extension function PDF.XPOWER (x, a > 0, b >= 0)
     = gsl_ran_exppow_pdf (x, a, b);
no_opt extension function RV.XPOWER (a > 0, b >= 0)
     = gsl_ran_exppow (get_rng (), a, b);

// F distribution.
function CDF.F (x >= 0, df1 > 0, df2 > 0) = gsl_cdf_fdist_P (x, df1, df2);
function IDF.F (P >= 0 && P < 1, df1 > 0, df2 > 0) = idf_fdist (P, df1, df2);
function PDF.F (x >= 0, df1 > 0, df2 > 0) = gsl_ran_fdist_pdf (x, df1, df2);
no_opt function RV.F (df1 > 0, df2 > 0) = gsl_ran_fdist (get_rng (), df1, df2);
function NCDF.F (x >= 0, df1 > 0, df2 > 0, lambda >= 0) = unimplemented;
function NPDF.F (x >= 0, df1 > 0, df2 > 0, lmabda >= 0) = unimplemented;
function SIG.F (x >= 0, df1 > 0, df2 > 0) = gsl_cdf_fdist_Q (x, df1, df2);

// Gamma distribution.
function CDF.GAMMA (x >= 0, a > 0, b > 0) = gsl_cdf_gamma_P (x, a, 1. / b);
function IDF.GAMMA (P >= 0 && P <= 1, a > 0, b > 0)
     = gsl_cdf_gamma_Pinv (P, a, 1. / b);
function PDF.GAMMA (x >= 0, a > 0, b > 0) = gsl_ran_gamma_pdf (x, a, 1. / b);
no_opt function RV.GAMMA (a > 0, b > 0)
     = gsl_ran_gamma (get_rng (), a, 1. / b);

// Half-normal distribution.
function CDF.HALFNRM (x, a, b > 0) = unimplemented;
function IDF.HALFNRM (P > 0 && P < 1, a, b > 0) = unimplemented;
function PDF.HALFNRM (x, a, b > 0) = unimplemented;
no_opt function RV.HALFNRM (a, b > 0) = unimplemented;

// Inverse Gaussian distribution.
function CDF.IGAUSS (x > 0, a > 0, b > 0) = unimplemented;
function IDF.IGAUSS (P >= 0 && P < 1, a > 0, b > 0) = unimplemented;
function PDF.IGAUSS (x > 0, a > 0, b > 0) = unimplemented;
no_opt function RV.IGAUSS (a > 0, b > 0) = unimplemented;

// Landau distribution.
extension function PDF.LANDAU (x) = gsl_ran_landau_pdf (x);
no_opt extension function RV.LANDAU () = gsl_ran_landau (get_rng ());

// Laplace distribution.
function CDF.LAPLACE (x, a, b > 0) = gsl_cdf_laplace_P ((x - a) / b, 1);
function IDF.LAPLACE (P > 0 && P < 1, a, b > 0)
     = a + b * gsl_cdf_laplace_Pinv (P, 1);
function PDF.LAPLACE (x, a, b > 0) = gsl_ran_laplace_pdf ((x - a) / b, 1) / b;
no_opt function RV.LAPLACE (a, b > 0)
     = a + b * gsl_ran_laplace (get_rng (), 1);

// Levy alpha-stable distribution.
no_opt extension function RV.LEVY (c, alpha > 0 && alpha <= 2)
     = gsl_ran_levy (get_rng (), c, alpha);

// Levy skew alpha-stable distribution.
no_opt extension function RV.LVSKEW (c, alpha > 0 && alpha <= 2,
                                     beta >= -1 && beta <= 1)
     = gsl_ran_levy_skew (get_rng (), c, alpha, beta);

// Logistic distribution.
function CDF.LOGISTIC (x, a, b > 0) = gsl_cdf_logistic_P ((x - a) / b, 1);
function IDF.LOGISTIC (P > 0 && P < 1, a, b > 0)
     = a + b * gsl_cdf_logistic_Pinv (P, 1);
function PDF.LOGISTIC (x, a, b > 0)
     = gsl_ran_logistic_pdf ((x - a) / b, 1) / b;
no_opt function RV.LOGISTIC (a, b > 0)
     = a + b * gsl_ran_logistic (get_rng (), 1);

// Lognormal distribution.
function CDF.LNORMAL (x >= 0, m > 0, s > 0)
     = gsl_cdf_lognormal_P (x, log (m), s);
function IDF.LNORMAL (P >= 0 && P < 1, m > 0, s > 0)
     = gsl_cdf_lognormal_Pinv (P, log (m), s);
function PDF.LNORMAL (x >= 0, m > 0, s > 0)
     = gsl_ran_lognormal_pdf (x, log (m), s);
no_opt function RV.LNORMAL (m > 0, s > 0)
     = gsl_ran_lognormal (get_rng (), log (m), s);

// Normal distribution.
function CDF.NORMAL (x, u, s > 0) = gsl_cdf_gaussian_P (x - u, s);
function IDF.NORMAL (P > 0 && P < 1, u, s > 0)
     = u + gsl_cdf_gaussian_Pinv (P, s);
function PDF.NORMAL (x, u, s > 0) = gsl_ran_gaussian_pdf ((x - u) / s, 1) / s;
no_opt function RV.NORMAL (u, s > 0) = u + gsl_ran_gaussian (get_rng (), s);
function CDFNORM (x) = gsl_cdf_ugaussian_P (x);
function PROBIT (P > 0 && P < 1) = gsl_cdf_ugaussian_Pinv (P);
no_opt function NORMAL (s > 0) = gsl_ran_gaussian (get_rng (), s);

// Normal tail distribution.
function PDF.NTAIL (x, a > 0, sigma > 0)
     = gsl_ran_gaussian_tail_pdf (x, a, sigma);
no_opt function RV.NTAIL (a > 0, sigma > 0)
     = gsl_ran_gaussian_tail (get_rng (), a, sigma);

// Pareto distribution.
function CDF.PARETO (x >= a, a > 0, b > 0) = gsl_cdf_pareto_P (x, b, a);
function IDF.PARETO (P >= 0 && P < 1, a > 0, b > 0)
     = gsl_cdf_pareto_Pinv (P, b, a);
function PDF.PARETO (x >= a, a > 0, b > 0) = gsl_ran_pareto_pdf (x, b, a);
no_opt function RV.PARETO (a > 0, b > 0) = gsl_ran_pareto (get_rng (), b, a);

// Rayleigh distribution.
extension function CDF.RAYLEIGH (x, sigma > 0) = gsl_cdf_rayleigh_P (x, sigma);
extension function IDF.RAYLEIGH (P >= 0 && P <= 1, sigma > 0)
     = gsl_cdf_rayleigh_Pinv (P, sigma);
extension function PDF.RAYLEIGH (x, sigma > 0)
     = gsl_ran_rayleigh_pdf (x, sigma);
no_opt extension function RV.RAYLEIGH (sigma > 0)
     = gsl_ran_rayleigh (get_rng (), sigma);

// Rayleigh tail distribution.
extension function PDF.RTAIL (x, a, sigma)
     = gsl_ran_rayleigh_tail_pdf (x, a, sigma);
no_opt extension function RV.RTAIL (a, sigma)
     = gsl_ran_rayleigh_tail (get_rng (), a, sigma);

// Studentized maximum modulus distribution.
function CDF.SMOD (x > 0, a >= 1, b >= 1) = unimplemented;
function IDF.SMOD (P >= 0 && P < 1, a >= 1, b >= 1) = unimplemented;

// Studentized range distribution.
function CDF.SRANGE (x > 0, a >= 1, b >= 1) = unimplemented;
function IDF.SRANGE (P >= 0 && P < 1, a >= 1, b >= 1) = unimplemented;

// Student t distribution.
function CDF.T (x, df > 0) = gsl_cdf_tdist_P (x, df);
function IDF.T (P > 0 && P < 1, df > 0) = gsl_cdf_tdist_Pinv (P, df);
function PDF.T (x, df > 0) = gsl_ran_tdist_pdf (x, df);
no_opt function RV.T (df > 0) = gsl_ran_tdist (get_rng (), df);
function NCDF.T (x, df > 0, nc) = unimplemented;
function NPDF.T (x, df > 0, nc) = unimplemented;

// Type-1 Gumbel distribution.
extension function CDF.T1G (x, a, b) = gsl_cdf_gumbel1_P (x, a, b);
extension function IDF.T1G (P >= 0 && P <= 1, a, b)
     = gsl_cdf_gumbel1_Pinv (P, a, b);
extension function PDF.T1G (x, a, b) = gsl_ran_gumbel1_pdf (x, a, b);
no_opt extension function RV.T1G (a, b) = gsl_ran_gumbel1 (get_rng (), a, b);

// Type-2 Gumbel distribution.
extension function CDF.T2G (x, a, b) = gsl_cdf_gumbel2_P (x, a, b);
extension function IDF.T2G (P >= 0 && P <= 1, a, b)
     = gsl_cdf_gumbel2_Pinv (P, a, b);
extension function PDF.T2G (x, a, b) = gsl_ran_gumbel2_pdf (x, a, b);
no_opt extension function RV.T2G (a, b) = gsl_ran_gumbel2 (get_rng (), a, b);

// Uniform distribution.
function CDF.UNIFORM (x <= b, a <= x, b) = gsl_cdf_flat_P (x, a, b);
function IDF.UNIFORM (P >= 0 && P <= 1, a <= b, b)
     = gsl_cdf_flat_Pinv (P, a, b);
function PDF.UNIFORM (x <= b, a <= x, b) = gsl_ran_flat_pdf (x, a, b);
no_opt function RV.UNIFORM (a <= b, b) = gsl_ran_flat (get_rng (), a, b);
no_opt function UNIFORM (b >= 0) = gsl_ran_flat (get_rng (), 0, b);

// Weibull distribution.
function CDF.WEIBULL (x >= 0, a > 0, b > 0) = gsl_cdf_weibull_P (x, a, b);
function IDF.WEIBULL (P >= 0 && P < 1, a > 0, b > 0)
     = gsl_cdf_weibull_Pinv (P, a, b);
function PDF.WEIBULL (x >= 0, a > 0, b > 0) = gsl_ran_weibull_pdf (x, a, b);
no_opt function RV.WEIBULL (a > 0, b > 0) = gsl_ran_weibull (get_rng (), a, b);

// Bernoulli distribution.
function CDF.BERNOULLI (k == 0 || k == 1, p >= 0 && p <= 1)
     = k ? 1 : 1 - p;
function PDF.BERNOULLI (k == 0 || k == 1, p >= 0 && p <= 1)
     = gsl_ran_bernoulli_pdf (k, p);
no_opt function RV.BERNOULLI (p >= 0 && p <= 1)
     = gsl_ran_bernoulli (get_rng (), p);

// Binomial distribution.
function CDF.BINOM (k, n > 0 && n == floor (n), p >= 0 && p <= 1)
     = gsl_cdf_binomial_P (k, p, n);
function PDF.BINOM (k >= 0 && k == floor (k) && k <= n,
                    n > 0 && n == floor (n),
                    p >= 0 && p <= 1)
     = gsl_ran_binomial_pdf (k, p, n);
no_opt function RV.BINOM (p > 0 && p == floor (p), n >= 0 && n <= 1)
     = gsl_ran_binomial (get_rng (), p, n);

// Geometric distribution.
function CDF.GEOM (k >= 1 && k == floor (k), p >= 0 && p <= 1)
     = gsl_cdf_geometric_P (k, p);
function PDF.GEOM (k >= 1 && k == floor (k),
                   p >= 0 && p <= 1)
     = gsl_ran_geometric_pdf (k, p);
no_opt function RV.GEOM (p >= 0 && p <= 1) = gsl_ran_geometric (get_rng (), p);

// Hypergeometric distribution.
function CDF.HYPER (k >= 0 && k == floor (k) && k <= c,
                    a > 0 && a == floor (a),
                    b > 0 && b == floor (b) && b <= a,
                    c > 0 && c == floor (c) && c <= a)
     = gsl_cdf_hypergeometric_P (k, c, a - c, b);
function PDF.HYPER (k >= 0 && k == floor (k) && k <= c,
                    a > 0 && a == floor (a),
                    b > 0 && b == floor (b) && b <= a,
                    c > 0 && c == floor (c) && c <= a)
     = gsl_ran_hypergeometric_pdf (k, c, a - c, b);
no_opt function RV.HYPER (a > 0 && a == floor (a),
                          b > 0 && b == floor (b) && b <= a,
                          c > 0 && c == floor (c) && c <= a)
     = gsl_ran_hypergeometric (get_rng (), c, a - c, b);

// Logarithmic distribution.
extension function PDF.LOG (k >= 1, p > 0 && p <= 1)
     = gsl_ran_logarithmic_pdf (k, p);
no_opt extension function RV.LOG (p > 0 && p <= 1)
     = gsl_ran_logarithmic (get_rng (), p);

// Negative binomial distribution.
function CDF.NEGBIN (k >= 1, n == floor (n), p > 0 && p <= 1)
     = gsl_cdf_negative_binomial_P (k, p, n);
function PDF.NEGBIN (k >= 1, n == floor (n), p > 0 && p <= 1)
     = gsl_ran_negative_binomial_pdf (k, p, n);
no_opt function RV.NEGBIN (n == floor (n), p > 0 && p <= 1)
     = gsl_ran_negative_binomial (get_rng (), p, n);

// Poisson distribution.
function CDF.POISSON (k >= 0 && k == floor (k), mu > 0)
     = gsl_cdf_poisson_P (k, mu);
function PDF.POISSON (k >= 0 && k == floor (k), mu > 0)
     = gsl_ran_poisson_pdf (k, mu);
no_opt function RV.POISSON (mu > 0) = gsl_ran_poisson (get_rng (), mu);

// Weirdness.
absorb_miss boolean function MISSING (x) = x == SYSMIS || !isfinite (x);
absorb_miss boolean function SYSMIS (x) = x == SYSMIS || !isfinite (x);
no_opt boolean function SYSMIS (num_var v)
     case c;
{
  return case_num (c, v) == SYSMIS;
}
no_opt function VALUE (num_var v)
     case c;
{
  return case_num (c, v);
}
no_opt function VALUE (num_vec_elem v)
{
  return v;
}

// A numeric vector element used in a "normal" context, in which a user-missing
// value becomes system-missing.
absorb_miss no_opt operator VEC_ELEM_NUM (idx)
     vector v;
     case c;
     expression e;
     expr_node n;
{
  const struct variable *var = expr_index_vector (e, n, v, idx);
  if (var)
    {
      double d = case_num (c, var);
      if (!var_is_num_missing (var, d, MV_USER))
        return d;
    }
  return SYSMIS;
}

// A numeric vector element used as the argument to the VALUE() function, in
// which a user-missing value retains its value.
//
// All numeric vector elements are initially parsed this way.  In most contexts
// they then get coerced into numbers.
absorb_miss no_opt num_vec_elem operator VEC_ELEM_NUM_RAW (idx)
     vector v;
     case c;
     expression e;
     expr_node n;
{
  const struct variable *var = expr_index_vector (e, n, v, idx);
  return var ? case_num (c, var) : SYSMIS;
}

absorb_miss no_opt string operator VEC_ELEM_STR (idx)
     expression e;
     vector v;
     case c;
     expr_node n;
{
  const struct variable *var = expr_index_vector (e, n, v, idx);
  return (var
          ? copy_string (e, CHAR_CAST_BUG (char *, case_str (c, var)),
                         var_get_width (var))
          : empty_string);
}

// Terminals.

no_opt operator NUM_VAR ()
     case c;
     num_var v;
{
  double d = case_num (c, v);
  return !var_is_num_missing (v, d, MV_USER) ? d : SYSMIS;
}

no_opt string operator STR_VAR ()
     case c;
     expression e;
     str_var v;
{
  struct substring s = alloc_string (e, var_get_width (v));
  memcpy (s.string, case_str (c, v), var_get_width (v));
  return s;
}

no_opt perm_only function LAG (num_var v, pos_int n_before)
    dataset ds;
{
  const struct ccase *c = lagged_case (ds, n_before);
  if (c != NULL)
    {
      double x = case_num (c, v);
      return !var_is_num_missing (v, x, MV_USER) ? x : SYSMIS;
    }
  else
    return SYSMIS;
}

no_opt perm_only function LAG (num_var v)
    dataset ds;
{
  const struct ccase *c = lagged_case (ds, 1);
  if (c != NULL)
    {
      double x = case_num (c, v);
      return !var_is_num_missing (v, x, MV_USER) ? x : SYSMIS;
    }
  else
    return SYSMIS;
}

no_opt perm_only string function LAG (str_var v, pos_int n_before)
     expression e;
     dataset ds;
{
  const struct ccase *c = lagged_case (ds, n_before);
  if (c != NULL)
    return copy_string (e, CHAR_CAST_BUG (char *, case_str (c, v)),
                        var_get_width (v));
  else
    return empty_string;
}

no_opt perm_only string function LAG (str_var v)
     expression e;
     dataset ds;
{
  const struct ccase *c = lagged_case (ds, 1);
  if (c != NULL)
    return copy_string (e, CHAR_CAST_BUG (char *, case_str (c, v)),
                        var_get_width (v));
  else
    return empty_string;
}

no_opt operator NUM_SYS ()
     case c;
     num_var v;
{
  return case_num (c, v) == SYSMIS;
}

no_opt operator NUM_VAL ()
     case c;
     num_var v;
{
  return case_num (c, v);
}

no_opt operator CASENUM ()
     case_idx idx;
{
  return idx;
}