+static double
+matrix_eval_CDF_EXP (double x, double a)
+{
+ return gsl_cdf_exponential_P (x, 1. / a);
+}
+
+static double
+matrix_eval_IDF_EXP (double P, double a)
+{
+ return gsl_cdf_exponential_Pinv (P, 1. / a);
+}
+
+static double
+matrix_eval_PDF_EXP (double x, double a)
+{
+ return gsl_ran_exponential_pdf (x, 1. / a);
+}
+
+static double
+matrix_eval_RV_EXP (double a)
+{
+ return gsl_ran_exponential (get_rng (), 1. / a);
+}
+
+static double
+matrix_eval_PDF_XPOWER (double x, double a, double b)
+{
+ return gsl_ran_exppow_pdf (x, a, b);
+}
+
+static double
+matrix_eval_RV_XPOWER (double a, double b)
+{
+ return gsl_ran_exppow (get_rng (), a, b);
+}
+
+static double
+matrix_eval_CDF_F (double x, double df1, double df2)
+{
+ return gsl_cdf_fdist_P (x, df1, df2);
+}
+
+static double
+matrix_eval_FCDF (double x, double df1, double df2)
+{
+ return matrix_eval_CDF_F (x, df1, df2);
+}
+
+static double
+matrix_eval_IDF_F (double P, double df1, double df2)
+{
+ return idf_fdist (P, df1, df2);
+}
+
+static double
+matrix_eval_RV_F (double df1, double df2)
+{
+ return gsl_ran_fdist (get_rng (), df1, df2);
+}
+
+static double
+matrix_eval_PDF_F (double x, double df1, double df2)
+{
+ return gsl_ran_fdist_pdf (x, df1, df2);
+}
+
+static double
+matrix_eval_SIG_F (double x, double df1, double df2)
+{
+ return gsl_cdf_fdist_Q (x, df1, df2);
+}
+
+static double
+matrix_eval_CDF_GAMMA (double x, double a, double b)
+{
+ return gsl_cdf_gamma_P (x, a, 1. / b);
+}
+
+static double
+matrix_eval_IDF_GAMMA (double P, double a, double b)
+{
+ return gsl_cdf_gamma_Pinv (P, a, 1. / b);
+}
+
+static double
+matrix_eval_PDF_GAMMA (double x, double a, double b)
+{
+ return gsl_ran_gamma_pdf (x, a, 1. / b);
+}
+
+static double
+matrix_eval_RV_GAMMA (double a, double b)
+{
+ return gsl_ran_gamma (get_rng (), a, 1. / b);
+}
+
+static double
+matrix_eval_PDF_LANDAU (double x)
+{
+ return gsl_ran_landau_pdf (x);
+}
+
+static double
+matrix_eval_RV_LANDAU (void)
+{
+ return gsl_ran_landau (get_rng ());
+}
+
+static double
+matrix_eval_CDF_LAPLACE (double x, double a, double b)
+{
+ return gsl_cdf_laplace_P ((x - a) / b, 1);
+}
+
+static double
+matrix_eval_IDF_LAPLACE (double P, double a, double b)
+{
+ return a + b * gsl_cdf_laplace_Pinv (P, 1);
+}
+
+static double
+matrix_eval_PDF_LAPLACE (double x, double a, double b)
+{
+ return gsl_ran_laplace_pdf ((x - a) / b, 1);
+}
+
+static double
+matrix_eval_RV_LAPLACE (double a, double b)
+{
+ return a + b * gsl_ran_laplace (get_rng (), 1);
+}
+
+static double
+matrix_eval_RV_LEVY (double c, double alpha)
+{
+ return gsl_ran_levy (get_rng (), c, alpha);
+}
+
+static double
+matrix_eval_RV_LVSKEW (double c, double alpha, double beta)
+{
+ return gsl_ran_levy_skew (get_rng (), c, alpha, beta);
+}
+
+static double
+matrix_eval_CDF_LOGISTIC (double x, double a, double b)
+{
+ return gsl_cdf_logistic_P ((x - a) / b, 1);
+}
+
+static double
+matrix_eval_IDF_LOGISTIC (double P, double a, double b)
+{
+ return a + b * gsl_cdf_logistic_Pinv (P, 1);
+}
+
+static double
+matrix_eval_PDF_LOGISTIC (double x, double a, double b)
+{
+ return gsl_ran_logistic_pdf ((x - a) / b, 1) / b;
+}
+
+static double
+matrix_eval_RV_LOGISTIC (double a, double b)
+{
+ return a + b * gsl_ran_logistic (get_rng (), 1);
+}
+
+static double
+matrix_eval_CDF_LNORMAL (double x, double m, double s)
+{
+ return gsl_cdf_lognormal_P (x, log (m), s);
+}
+
+static double
+matrix_eval_IDF_LNORMAL (double P, double m, double s)
+{
+ return gsl_cdf_lognormal_Pinv (P, log (m), s);;
+}
+
+static double
+matrix_eval_PDF_LNORMAL (double x, double m, double s)
+{
+ return gsl_ran_lognormal_pdf (x, log (m), s);
+}
+
+static double
+matrix_eval_RV_LNORMAL (double m, double s)
+{
+ return gsl_ran_lognormal (get_rng (), log (m), s);
+}
+
+static double
+matrix_eval_CDF_NORMAL (double x, double u, double s)
+{
+ return gsl_cdf_gaussian_P (x - u, s);
+}
+
+static double
+matrix_eval_IDF_NORMAL (double P, double u, double s)
+{
+ return u + gsl_cdf_gaussian_Pinv (P, s);
+}
+
+static double
+matrix_eval_PDF_NORMAL (double x, double u, double s)
+{
+ return gsl_ran_gaussian_pdf ((x - u) / s, 1) / s;
+}
+
+static double
+matrix_eval_RV_NORMAL (double u, double s)
+{
+ return u + gsl_ran_gaussian (get_rng (), s);
+}
+
+static double
+matrix_eval_CDFNORM (double x)
+{
+ return gsl_cdf_ugaussian_P (x);
+}
+
+static double
+matrix_eval_PROBIT (double P)
+{
+ return gsl_cdf_ugaussian_Pinv (P);
+}
+
+static double
+matrix_eval_NORMAL (double s)
+{
+ return gsl_ran_gaussian (get_rng (), s);
+}
+
+static double
+matrix_eval_PDF_NTAIL (double x, double a, double sigma)
+{
+ return gsl_ran_gaussian_tail_pdf (x, a, sigma);;
+}
+
+static double
+matrix_eval_RV_NTAIL (double a, double sigma)
+{
+ return gsl_ran_gaussian_tail (get_rng (), a, sigma);
+}
+
+static double
+matrix_eval_CDF_PARETO (double x, double a, double b)
+{
+ return gsl_cdf_pareto_P (x, b, a);
+}
+
+static double
+matrix_eval_IDF_PARETO (double P, double a, double b)
+{
+ return gsl_cdf_pareto_Pinv (P, b, a);
+}
+
+static double
+matrix_eval_PDF_PARETO (double x, double a, double b)
+{
+ return gsl_ran_pareto_pdf (x, b, a);
+}
+
+static double
+matrix_eval_RV_PARETO (double a, double b)
+{
+ return gsl_ran_pareto (get_rng (), b, a);
+}
+
+static double
+matrix_eval_CDF_RAYLEIGH (double x, double sigma)
+{
+ return gsl_cdf_rayleigh_P (x, sigma);
+}
+
+static double
+matrix_eval_IDF_RAYLEIGH (double P, double sigma)
+{
+ return gsl_cdf_rayleigh_Pinv (P, sigma);
+}
+
+static double
+matrix_eval_PDF_RAYLEIGH (double x, double sigma)
+{
+ return gsl_ran_rayleigh_pdf (x, sigma);
+}
+
+static double
+matrix_eval_RV_RAYLEIGH (double sigma)
+{
+ return gsl_ran_rayleigh (get_rng (), sigma);
+}
+
+static double
+matrix_eval_PDF_RTAIL (double x, double a, double sigma)
+{
+ return gsl_ran_rayleigh_tail_pdf (x, a, sigma);
+}
+
+static double
+matrix_eval_RV_RTAIL (double a, double sigma)
+{
+ return gsl_ran_rayleigh_tail (get_rng (), a, sigma);
+}
+
+static double
+matrix_eval_CDF_T (double x, double df)
+{
+ return gsl_cdf_tdist_P (x, df);
+}
+
+static double
+matrix_eval_TCDF (double x, double df)
+{
+ return matrix_eval_CDF_T (x, df);
+}
+
+static double
+matrix_eval_IDF_T (double P, double df)
+{
+ return gsl_cdf_tdist_Pinv (P, df);
+}
+
+static double
+matrix_eval_PDF_T (double x, double df)
+{
+ return gsl_ran_tdist_pdf (x, df);
+}
+
+static double
+matrix_eval_RV_T (double df)
+{
+ return gsl_ran_tdist (get_rng (), df);
+}
+
+static double
+matrix_eval_CDF_T1G (double x, double a, double b)
+{
+ return gsl_cdf_gumbel1_P (x, a, b);
+}
+
+static double
+matrix_eval_IDF_T1G (double P, double a, double b)
+{
+ return gsl_cdf_gumbel1_Pinv (P, a, b);
+}
+
+static double
+matrix_eval_PDF_T1G (double x, double a, double b)
+{
+ return gsl_ran_gumbel1_pdf (x, a, b);
+}
+
+static double
+matrix_eval_RV_T1G (double a, double b)
+{
+ return gsl_ran_gumbel1 (get_rng (), a, b);
+}
+
+static double
+matrix_eval_CDF_T2G (double x, double a, double b)
+{
+ return gsl_cdf_gumbel1_P (x, a, b);
+}
+
+static double
+matrix_eval_IDF_T2G (double P, double a, double b)
+{
+ return gsl_cdf_gumbel1_Pinv (P, a, b);
+}
+
+static double
+matrix_eval_PDF_T2G (double x, double a, double b)
+{
+ return gsl_ran_gumbel1_pdf (x, a, b);
+}
+
+static double
+matrix_eval_RV_T2G (double a, double b)
+{
+ return gsl_ran_gumbel1 (get_rng (), a, b);
+}
+
+static double
+matrix_eval_CDF_UNIFORM (double x, double a, double b)
+{
+ return gsl_cdf_flat_P (x, a, b);
+}
+
+static double
+matrix_eval_IDF_UNIFORM (double P, double a, double b)
+{
+ return gsl_cdf_flat_Pinv (P, a, b);
+}
+
+static double
+matrix_eval_PDF_UNIFORM (double x, double a, double b)
+{
+ return gsl_ran_flat_pdf (x, a, b);
+}
+
+static double
+matrix_eval_RV_UNIFORM (double a, double b)
+{
+ return gsl_ran_flat (get_rng (), a, b);
+}
+
+static double
+matrix_eval_CDF_WEIBULL (double x, double a, double b)
+{
+ return gsl_cdf_weibull_P (x, a, b);
+}
+
+static double
+matrix_eval_IDF_WEIBULL (double P, double a, double b)
+{
+ return gsl_cdf_weibull_Pinv (P, a, b);
+}
+
+static double
+matrix_eval_PDF_WEIBULL (double x, double a, double b)
+{
+ return gsl_ran_weibull_pdf (x, a, b);
+}
+
+static double
+matrix_eval_RV_WEIBULL (double a, double b)
+{
+ return gsl_ran_weibull (get_rng (), a, b);
+}
+
+static double
+matrix_eval_CDF_BERNOULLI (double k, double p)
+{
+ return k ? 1 : 1 - p;
+}
+
+static double
+matrix_eval_PDF_BERNOULLI (double k, double p)
+{
+ return gsl_ran_bernoulli_pdf (k, p);
+}
+
+static double
+matrix_eval_RV_BERNOULLI (double p)
+{
+ return gsl_ran_bernoulli (get_rng (), p);
+}
+
+static double
+matrix_eval_CDF_BINOM (double k, double n, double p)
+{
+ return gsl_cdf_binomial_P (k, p, n);
+}
+
+static double
+matrix_eval_PDF_BINOM (double k, double n, double p)
+{
+ return gsl_ran_binomial_pdf (k, p, n);
+}
+
+static double
+matrix_eval_RV_BINOM (double n, double p)
+{
+ return gsl_ran_binomial (get_rng (), p, n);
+}
+
+static double
+matrix_eval_CDF_GEOM (double k, double p)
+{
+ return gsl_cdf_geometric_P (k, p);
+}
+
+static double
+matrix_eval_PDF_GEOM (double k, double p)
+{
+ return gsl_ran_geometric_pdf (k, p);
+}
+
+static double
+matrix_eval_RV_GEOM (double p)
+{
+ return gsl_ran_geometric (get_rng (), p);
+}
+
+static double
+matrix_eval_CDF_HYPER (double k, double a, double b, double c)
+{
+ return gsl_cdf_hypergeometric_P (k, c, a - c, b);
+}
+
+static double
+matrix_eval_PDF_HYPER (double k, double a, double b, double c)
+{
+ return gsl_ran_hypergeometric_pdf (k, c, a - c, b);
+}
+
+static double
+matrix_eval_RV_HYPER (double a, double b, double c)
+{
+ return gsl_ran_hypergeometric (get_rng (), c, a - c, b);
+}
+
+static double
+matrix_eval_PDF_LOG (double k, double p)
+{
+ return gsl_ran_logarithmic_pdf (k, p);
+}
+
+static double
+matrix_eval_RV_LOG (double p)
+{
+ return gsl_ran_logarithmic (get_rng (), p);
+}
+
+static double
+matrix_eval_CDF_NEGBIN (double k, double n, double p)
+{
+ return gsl_cdf_negative_binomial_P (k, p, n);
+}
+
+static double
+matrix_eval_PDF_NEGBIN (double k, double n, double p)
+{
+ return gsl_ran_negative_binomial_pdf (k, p, n);
+}
+
+static double
+matrix_eval_RV_NEGBIN (double n, double p)
+{
+ return gsl_ran_negative_binomial (get_rng (), p, n);
+}
+
+static double
+matrix_eval_CDF_POISSON (double k, double mu)
+{
+ return gsl_cdf_poisson_P (k, mu);
+}
+
+static double
+matrix_eval_PDF_POISSON (double k, double mu)
+{
+ return gsl_ran_poisson_pdf (k, mu);
+}
+
+static double
+matrix_eval_RV_POISSON (double mu)
+{
+ return gsl_ran_poisson (get_rng (), mu);
+}
+
+struct matrix_function
+ {
+ const char *name;
+ enum matrix_op op;
+ size_t min_args, max_args;
+ };
+
+static struct matrix_expr *matrix_parse_expr (struct matrix_state *);
+
+static bool
+word_matches (const char **test, const char **name)
+{
+ size_t test_len = strcspn (*test, ".");
+ size_t name_len = strcspn (*name, ".");
+ if (test_len == name_len)
+ {
+ if (buf_compare_case (*test, *name, test_len))
+ return false;
+ }
+ else if (test_len < 3 || test_len > name_len)
+ return false;
+ else
+ {
+ if (buf_compare_case (*test, *name, test_len))
+ return false;
+ }
+
+ *test += test_len;
+ *name += name_len;
+ if (**test != **name)
+ return false;
+
+ if (**test == '.')
+ {
+ (*test)++;
+ (*name)++;
+ }
+ return true;
+}
+
+/* Returns 0 if TOKEN and FUNC do not match,
+ 1 if TOKEN is an acceptable abbreviation for FUNC,
+ 2 if TOKEN equals FUNC. */
+static int
+compare_function_names (const char *token_, const char *func_)
+{
+ const char *token = token_;
+ const char *func = func_;
+ while (*token || *func)
+ if (!word_matches (&token, &func))
+ return 0;
+ return !c_strcasecmp (token_, func_) ? 2 : 1;
+}
+
+static const struct matrix_function *
+matrix_parse_function_name (const char *token)
+{
+ static const struct matrix_function functions[] =
+ {
+#define F(ENUM, STRING, PROTO, CONSTRAINTS) \
+ { STRING, MOP_F_##ENUM, PROTO##_MIN_ARGS, PROTO##_MAX_ARGS },
+ MATRIX_FUNCTIONS
+#undef F
+ };
+ enum { N_FUNCTIONS = sizeof functions / sizeof *functions };
+
+ for (size_t i = 0; i < N_FUNCTIONS; i++)
+ {
+ if (compare_function_names (token, functions[i].name) > 0)
+ return &functions[i];
+ }
+ return NULL;
+}
+
+static struct read_file *
+read_file_create (struct matrix_state *s, struct file_handle *fh)
+{
+ for (size_t i = 0; i < s->n_read_files; i++)
+ {
+ struct read_file *rf = s->read_files[i];
+ if (rf->file == fh)
+ {
+ fh_unref (fh);
+ return rf;
+ }
+ }
+
+ struct read_file *rf = xmalloc (sizeof *rf);
+ *rf = (struct read_file) { .file = fh };
+
+ s->read_files = xrealloc (s->read_files,
+ (s->n_read_files + 1) * sizeof *s->read_files);
+ s->read_files[s->n_read_files++] = rf;
+ return rf;
+}
+
+static struct dfm_reader *
+read_file_open (struct read_file *rf)
+{
+ if (!rf->reader)
+ rf->reader = dfm_open_reader (rf->file, NULL, rf->encoding);
+ return rf->reader;
+}
+
+static void
+read_file_destroy (struct read_file *rf)
+{
+ if (rf)
+ {
+ fh_unref (rf->file);
+ dfm_close_reader (rf->reader);
+ free (rf->encoding);
+ free (rf);
+ }
+}
+
+static struct write_file *
+write_file_create (struct matrix_state *s, struct file_handle *fh)
+{
+ for (size_t i = 0; i < s->n_write_files; i++)
+ {
+ struct write_file *wf = s->write_files[i];
+ if (wf->file == fh)
+ {
+ fh_unref (fh);
+ return wf;
+ }
+ }
+
+ struct write_file *wf = xmalloc (sizeof *wf);
+ *wf = (struct write_file) { .file = fh };
+
+ s->write_files = xrealloc (s->write_files,
+ (s->n_write_files + 1) * sizeof *s->write_files);
+ s->write_files[s->n_write_files++] = wf;
+ return wf;
+}
+
+static struct dfm_writer *
+write_file_open (struct write_file *wf)
+{
+ if (!wf->writer)
+ wf->writer = dfm_open_writer (wf->file, wf->encoding);
+ return wf->writer;
+}
+
+static void
+write_file_destroy (struct write_file *wf)
+{
+ if (wf)
+ {
+ if (wf->held)
+ {
+ dfm_put_record_utf8 (wf->writer, wf->held->s.ss.string,
+ wf->held->s.ss.length);
+ u8_line_destroy (wf->held);
+ free (wf->held);
+ }
+
+ fh_unref (wf->file);
+ dfm_close_writer (wf->writer);
+ free (wf->encoding);
+ free (wf);
+ }
+}
+
+static bool
+matrix_parse_function (struct matrix_state *s, const char *token,
+ struct matrix_expr **exprp)
+{
+ *exprp = NULL;
+ if (lex_next_token (s->lexer, 1) != T_LPAREN)
+ return false;
+
+ int start_ofs = lex_ofs (s->lexer);
+ if (lex_match_id (s->lexer, "EOF"))
+ {
+ lex_get (s->lexer);
+ struct file_handle *fh = fh_parse (s->lexer, FH_REF_FILE, s->session);
+ if (!fh)
+ return true;
+
+ if (!lex_force_match (s->lexer, T_RPAREN))
+ {
+ fh_unref (fh);
+ return true;
+ }
+
+ struct read_file *rf = read_file_create (s, fh);
+
+ struct matrix_expr *e = xmalloc (sizeof *e);
+ *e = (struct matrix_expr) { .op = MOP_EOF, .eof = rf };
+ matrix_expr_wrap_location (s, start_ofs, e);
+ *exprp = e;
+ return true;
+ }
+
+ const struct matrix_function *f = matrix_parse_function_name (token);
+ if (!f)
+ return false;
+
+ struct matrix_expr *e = xmalloc (sizeof *e);
+ *e = (struct matrix_expr) { .op = f->op };
+
+ lex_get_n (s->lexer, 2);
+ if (lex_token (s->lexer) != T_RPAREN)
+ {
+ size_t allocated_subs = 0;
+ do
+ {
+ struct matrix_expr *sub = matrix_parse_expr (s);
+ if (!sub)
+ goto error;
+
+ if (e->n_subs >= allocated_subs)
+ e->subs = x2nrealloc (e->subs, &allocated_subs, sizeof *e->subs);
+ e->subs[e->n_subs++] = sub;
+ }
+ while (lex_match (s->lexer, T_COMMA));
+ }
+ if (!lex_force_match (s->lexer, T_RPAREN))
+ goto error;
+
+ if (e->n_subs < f->min_args || e->n_subs > f->max_args)
+ {
+ if (f->min_args == f->max_args)
+ msg_at (SE, e->location,
+ ngettext ("Matrix function %s requires %zu argument.",
+ "Matrix function %s requires %zu arguments.",
+ f->min_args),
+ f->name, f->min_args);
+ else if (f->min_args == 1 && f->max_args == 2)
+ msg_at (SE, e->location,
+ ngettext ("Matrix function %s requires 1 or 2 arguments, "
+ "but %zu was provided.",
+ "Matrix function %s requires 1 or 2 arguments, "
+ "but %zu were provided.",
+ e->n_subs),
+ f->name, e->n_subs);
+ else if (f->min_args == 1 && f->max_args == INT_MAX)
+ msg_at (SE, e->location,
+ _("Matrix function %s requires at least one argument."),
+ f->name);
+ else
+ NOT_REACHED ();
+
+ goto error;
+ }
+
+ matrix_expr_wrap_location (s, start_ofs, e);