work toward better error reporting
[pspp] / src / language / stats / matrix.c
index e7d3ae8e27f380facbf3a0105d029b287a940fb4..2e7bb83dfb1c1d21d09bcc78533b1a21718b5d73 100644 (file)
@@ -17,6 +17,7 @@
 #include <config.h>
 
 #include <gsl/gsl_blas.h>
+#include <gsl/gsl_cdf.h>
 #include <gsl/gsl_eigen.h>
 #include <gsl/gsl_linalg.h>
 #include <gsl/gsl_matrix.h>
 #include "libpspp/string-array.h"
 #include "libpspp/stringi-set.h"
 #include "libpspp/u8-line.h"
+#include "math/distributions.h"
 #include "math/random.h"
 #include "output/driver.h"
 #include "output/output-item.h"
 #include "output/pivot-table.h"
 
 #include "gl/c-ctype.h"
+#include "gl/c-strcase.h"
 #include "gl/ftoastr.h"
 #include "gl/minmax.h"
 #include "gl/xsize.h"
@@ -177,87 +180,252 @@ matrix_var_set (struct matrix_var *var, gsl_matrix *value)
   var->value = value;
 }
 \f
-#define MATRIX_FUNCTIONS \
-    F(ABS, m_m) \
-    F(ALL, d_m) \
-    F(ANY, d_m) \
-    F(ARSIN, m_m) \
-    F(ARTAN, m_m) \
-    F(BLOCK, m_any) \
-    F(CHOL, m_m) \
-    F(CMIN, m_m) \
-    F(CMAX, m_m) \
-    F(COS, m_m) \
-    F(CSSQ, m_m) \
-    F(CSUM, m_m) \
-    F(DESIGN, m_m) \
-    F(DET, d_m) \
-    F(DIAG, m_m) \
-    F(EVAL, m_m) \
-    F(EXP, m_m) \
-    F(GINV, m_m) \
-    F(GRADE, m_m) \
-    F(GSCH, m_m) \
-    F(IDENT, IDENT) \
-    F(INV, m_m) \
-    F(KRONEKER, m_mm) \
-    F(LG10, m_m) \
-    F(LN, m_m) \
-    F(MAGIC, m_d) \
-    F(MAKE, m_ddd) \
-    F(MDIAG, m_v) \
-    F(MMAX, d_m) \
-    F(MMIN, d_m) \
-    F(MOD, m_md) \
-    F(MSSQ, d_m) \
-    F(MSUM, d_m) \
-    F(NCOL, d_m) \
-    F(NROW, d_m) \
-    F(RANK, d_m) \
-    F(RESHAPE, m_mdd) \
-    F(RMAX, m_m) \
-    F(RMIN, m_m) \
-    F(RND, m_m) \
-    F(RNKORDER, m_m) \
-    F(RSSQ, m_m) \
-    F(RSUM, m_m) \
-    F(SIN, m_m) \
-    F(SOLVE, m_mm) \
-    F(SQRT, m_m) \
-    F(SSCP, m_m) \
-    F(SVAL, m_m) \
-    F(SWEEP, m_md) \
-    F(T, m_m) \
-    F(TRACE, d_m) \
-    F(TRANSPOS, m_m) \
-    F(TRUNC, m_m) \
-    F(UNIFORM, m_dd)
+/* The third argument to F() is a "prototype".  For most prototypes, the first
+   letter (before the _) represents the return type and each other letter
+   (after the _) is an argument type.  The types are:
+
+     - "m": A matrix of unrestricted dimensions.
+
+     - "d": A scalar.
+
+     - "v": A row or column vector.
+
+     - "e": Primarily for the first argument, this is a matrix with
+       unrestricted dimensions treated elementwise.  Each element in the matrix
+       is passed to the implementation function separately.
+
+   The fourth argument is an optional constraints string.  For this purpose the
+   first argument is named "a", the second "b", and so on.  The following kinds
+   of constraints are supported.  For matrix arguments, the constraints are
+   applied to each value in the matrix separately:
+
+     - "a(0,1)" or "a[0,1]": 0 < a < 1 or 0 <= a <= 1, respectively.  Any
+       integer may substitute for 0 and 1.  Half-open constraints (] and [) are
+       also supported.
+
+     - "ai": Restrict a to integer values.
+
+     - "a>0", "a<0", "a>=0", "a<=0".
+
+     - "a<b", "a>b", "a<=b", "a>=b".
+*/
+#define MATRIX_FUNCTIONS                                                \
+    F(ABS,      "ABS",      m_e, NULL)                                  \
+    F(ALL,      "ALL",      d_m, NULL)                                  \
+    F(ANY,      "ANY",      d_m, NULL)                                  \
+    F(ARSIN,    "ARSIN",    m_e, "a[-1,1]")                             \
+    F(ARTAN,    "ARTAN",    m_e, NULL)                                  \
+    F(BLOCK,    "BLOCK",    m_any, NULL)                                \
+    F(CHOL,     "CHOL",     m_m, NULL)                                  \
+    F(CMIN,     "CMIN",     m_m, NULL)                                  \
+    F(CMAX,     "CMAX",     m_m, NULL)                                  \
+    F(COS,      "COS",      m_e, NULL)                                  \
+    F(CSSQ,     "CSSQ",     m_m, NULL)                                  \
+    F(CSUM,     "CSUM",     m_m, NULL)                                  \
+    F(DESIGN,   "DESIGN",   m_m, NULL)                                  \
+    F(DET,      "DET",      d_m, NULL)                                  \
+    F(DIAG,     "DIAG",     m_m, NULL)                                  \
+    F(EVAL,     "EVAL",     m_m, NULL)                                  \
+    F(EXP,      "EXP",      m_e, NULL)                                  \
+    F(GINV,     "GINV",     m_m, NULL)                                  \
+    F(GRADE,    "GRADE",    m_m, NULL)                                  \
+    F(GSCH,     "GSCH",     m_m, NULL)                                  \
+    F(IDENT,    "IDENT",    IDENT, NULL)                                \
+    F(INV,      "INV",      m_m, NULL)                                  \
+    F(KRONEKER, "KRONEKER", m_mm, NULL)                                 \
+    F(LG10,     "LG10",     m_e, "a>0")                                 \
+    F(LN,       "LN",       m_e, "a>0")                                 \
+    F(MAGIC,    "MAGIC",    m_d, "ai>=3")                               \
+    F(MAKE,     "MAKE",     m_ddd, NULL)                                \
+    F(MDIAG,    "MDIAG",    m_v, NULL)                                  \
+    F(MMAX,     "MMAX",     d_m, NULL)                                  \
+    F(MMIN,     "MMIN",     d_m, NULL)                                  \
+    F(MOD,      "MOD",      m_md, NULL)                                 \
+    F(MSSQ,     "MSSQ",     d_m, NULL)                                  \
+    F(MSUM,     "MSUM",     d_m, NULL)                                  \
+    F(NCOL,     "NCOL",     d_m, NULL)                                  \
+    F(NROW,     "NROW",     d_m, NULL)                                  \
+    F(RANK,     "RANK",     d_m, NULL)                                  \
+    F(RESHAPE,  "RESHAPE",  m_mdd, NULL)                                \
+    F(RMAX,     "RMAX",     m_m, NULL)                                  \
+    F(RMIN,     "RMIN",     m_m, NULL)                                  \
+    F(RND,      "RND",      m_e, NULL)                                  \
+    F(RNKORDER, "RNKORDER", m_m, NULL)                                  \
+    F(RSSQ,     "RSSQ",     m_m, NULL)                                  \
+    F(RSUM,     "RSUM",     m_m, NULL)                                  \
+    F(SIN,      "SIN",      m_e, NULL)                                  \
+    F(SOLVE,    "SOLVE",    m_mm, NULL)                                 \
+    F(SQRT,     "SQRT",     m_e, "a>=0")                                \
+    F(SSCP,     "SSCP",     m_m, NULL)                                  \
+    F(SVAL,     "SVAL",     m_m, NULL)                                  \
+    F(SWEEP,    "SWEEP",    m_md, NULL)                                 \
+    F(T,        "T",        m_m, NULL)                                  \
+    F(TRACE,    "TRACE",    d_m, NULL)                                  \
+    F(TRANSPOS, "TRANSPOS", m_m, NULL)                                  \
+    F(TRUNC,    "TRUNC",    m_e, NULL)                                  \
+    F(UNIFORM,  "UNIFORM",  m_dd, "ai>=0 bi>=0")                        \
+    F(PDF_BETA, "PDF.BETA", m_edd, "a[0,1] b>0 c>0")                    \
+    F(CDF_BETA, "CDF.BETA", m_edd, "a[0,1] b>0 c>0")                    \
+    F(IDF_BETA, "IDF.BETA", m_edd, "a[0,1] b>0 c>0")                    \
+    F(RV_BETA,  "RV.BETA",  d_dd, "a>0 b>0")                            \
+    F(NCDF_BETA, "NCDF.BETA", m_eddd, "a>=0 b>0 c>0 d>0")               \
+    F(NPDF_BETA, "NCDF.BETA", m_eddd, "a>=0 b>0 c>0 d>0")               \
+    F(CDF_BVNOR, "CDF.BVNOR", m_eed, "c[-1,1]")                         \
+    F(PDF_BVNOR, "PDF.BVNOR", m_eed, "c[-1,1]")                         \
+    F(CDF_CAUCHY, "CDF.CAUCHY", m_edd, "c>0")                           \
+    F(IDF_CAUCHY, "IDF.CAUCHY", m_edd, "a(0,1) c>0")                    \
+    F(PDF_CAUCHY, "PDF.CAUCHY", m_edd, "c>0")                           \
+    F(RV_CAUCHY, "RV.CAUCHY", d_dd, "b>0")                              \
+    F(CDF_CHISQ, "CDF.CHISQ", m_ed, "a>=0 b>0")                         \
+    F(CHICDF, "CHICDF", m_ed, "a>=0 b>0")                               \
+    F(IDF_CHISQ, "IDF.CHISQ", m_ed, "a[0,1) b>0")                       \
+    F(PDF_CHISQ, "PDF.CHISQ", m_ed, "a>=0 b>0")                         \
+    F(RV_CHISQ, "RV.CHISQ", d_d, "a>0")                                 \
+    F(SIG_CHISQ, "SIG.CHISQ", m_ed, "a>=0 b>0")                         \
+    F(CDF_EXP, "CDF.EXP", m_ed, "a>=0 b>=0")                            \
+    F(IDF_EXP, "IDF.EXP", m_ed, "a[0,1) b>0")                           \
+    F(PDF_EXP, "PDF.EXP", m_ed, "a>=0 b>0")                             \
+    F(RV_EXP, "RV.EXP", d_d, "a>0")                                     \
+    F(PDF_XPOWER, "PDF.XPOWER", m_edd, "b>0 c>=0")                      \
+    F(RV_XPOWER, "RV.XPOWER", d_dd, "a>0 c>=0")                         \
+    F(CDF_F, "CDF.F", m_edd, "a>=0 b>0 c>0")                            \
+    F(FCDF, "FCDF", m_edd, "a>=0 b>0 c>0")                              \
+    F(IDF_F, "IDF.F", m_edd, "a[0,1) b>0 c>0")                          \
+    F(PDF_F, "PDF.F", m_edd, "a>=0 b>0 c>0")                            \
+    F(RV_F, "RV.F", d_dd, "a>0 b>0")                                    \
+    F(SIG_F, "SIG.F", m_edd, "a>=0 b>0 c>0")                            \
+    F(CDF_GAMMA, "CDF.GAMMA", m_edd, "a>=0 b>0 c>0")                    \
+    F(IDF_GAMMA, "IDF.GAMMA", m_edd, "a[0,1] b>0 c>0")                  \
+    F(PDF_GAMMA, "PDF.GAMMA", m_edd, "a>=0 b>0 c>0")                    \
+    F(RV_GAMMA, "RV.GAMMA", d_dd, "a>0 b>0")                            \
+    F(PDF_LANDAU, "PDF.LANDAU", m_e, NULL)                              \
+    F(RV_LANDAU, "RV.LANDAU", d_none, NULL)                             \
+    F(CDF_LAPLACE, "CDF.LAPLACE", m_edd, "c>0")                         \
+    F(IDF_LAPLACE, "IDF.LAPLACE", m_edd, "a(0,1) c>0")                  \
+    F(PDF_LAPLACE, "PDF.LAPLACE", m_edd, "c>0")                         \
+    F(RV_LAPLACE, "RV.LAPLACE", d_dd, "b>0")                            \
+    F(RV_LEVY, "RV.LEVY", d_dd, "b(0,2]")                               \
+    F(RV_LVSKEW, "RV.LVSKEW", d_ddd, "b(0,2] c[-1,1]")                  \
+    F(CDF_LOGISTIC, "CDF.LOGISTIC", m_edd, "c>0")                       \
+    F(IDF_LOGISTIC, "IDF.LOGISTIC", m_edd, "a(0,1) c>0")                \
+    F(PDF_LOGISTIC, "PDF.LOGISTIC", m_edd, "c>0")                       \
+    F(RV_LOGISTIC, "RV.LOGISTIC", d_dd, "b>0")                          \
+    F(CDF_LNORMAL, "CDF.LNORMAL", m_edd, "a>=0 b>0 c>0")                \
+    F(IDF_LNORMAL, "IDF.LNORMAL", m_edd, "a[0,1) b>0 c>0")              \
+    F(PDF_LNORMAL, "PDF.LNORMAL", m_edd, "a>=0 b>0 c>0")                \
+    F(RV_LNORMAL, "RV.LNORMAL", d_dd, "a>0 b>0")                        \
+    F(CDF_NORMAL, "CDF.NORMAL", m_edd, "c>0")                           \
+    F(IDF_NORMAL, "IDF.NORMAL", m_edd, "a(0,1) c>0")                    \
+    F(PDF_NORMAL, "PDF.NORMAL", m_edd, "c>0")                           \
+    F(RV_NORMAL, "RV.NORMAL", d_dd, "b>0")                              \
+    F(CDFNORM, "CDFNORM", m_e, NULL)                                    \
+    F(PROBIT, "PROBIT", m_e, "a(0,1)")                                  \
+    F(NORMAL, "NORMAL", m_e, "a>0")                                     \
+    F(PDF_NTAIL, "PDF.NTAIL", m_edd, "b>0 c>0")                         \
+    F(RV_NTAIL, "RV.NTAIL", d_dd, "a>0 b>0")                            \
+    F(CDF_PARETO, "CDF.PARETO", m_edd, "a>=b b>0 c>0")                  \
+    F(IDF_PARETO, "IDF.PARETO", m_edd, "a[0,1) b>0 c>0")                \
+    F(PDF_PARETO, "PDF.PARETO", m_edd, "a>=b b>0 c>0")                  \
+    F(RV_PARETO, "RV.PARETO", d_dd, "a>0 b>0")                          \
+    F(CDF_RAYLEIGH, "CDF.RAYLEIGH", m_ed, "b>0")                        \
+    F(IDF_RAYLEIGH, "IDF.RAYLEIGH", m_ed, "a[0,1] b>0")                 \
+    F(PDF_RAYLEIGH, "PDF.RAYLEIGH", m_ed, "b>0")                        \
+    F(RV_RAYLEIGH, "RV.RAYLEIGH", d_d, "a>0")                           \
+    F(PDF_RTAIL, "PDF.RTAIL", m_edd, NULL)                              \
+    F(RV_RTAIL, "RV.RTAIL", d_dd, NULL)                                 \
+    F(CDF_T, "CDF.T", m_ed, "b>0")                                      \
+    F(TCDF, "TCDF", m_ed, "b>0")                                        \
+    F(IDF_T, "IDF.T", m_ed, "a(0,1) b>0")                               \
+    F(PDF_T, "PDF.T", m_ed, "b>0")                                      \
+    F(RV_T, "RV.T", d_d, "a>0")                                         \
+    F(CDF_T1G, "CDF.T1G", m_edd, NULL)                                  \
+    F(IDF_T1G, "IDF.T1G", m_edd, "a(0,1)")                              \
+    F(PDF_T1G, "PDF.T1G", m_edd, NULL)                                  \
+    F(RV_T1G, "RV.T1G", d_dd, NULL)                                     \
+    F(CDF_T2G, "CDF.T2G", m_edd, NULL)                                  \
+    F(IDF_T2G, "IDF.T2G", m_edd, "a(0,1)")                              \
+    F(PDF_T2G, "PDF.T2G", m_edd, NULL)                                  \
+    F(RV_T2G, "RV.T2G", d_dd, NULL)                                     \
+    F(CDF_UNIFORM, "CDF.UNIFORM", m_edd, "a<=c b<=c")                   \
+    F(IDF_UNIFORM, "IDF.UNIFORM", m_edd, "a[0,1] b<=c")                 \
+    F(PDF_UNIFORM, "PDF.UNIFORM", m_edd, "a<=c b<=c")                   \
+    F(RV_UNIFORM, "RV.UNIFORM", d_dd, "a<=b")                           \
+    F(CDF_WEIBULL, "CDF.WEIBULL", m_edd, "a>=0 b>0 c>0")                \
+    F(IDF_WEIBULL, "IDF.WEIBULL", m_edd, "a[0,1) b>0 c>0")              \
+    F(PDF_WEIBULL, "PDF.WEIBULL", m_edd, "a>=0 b>0 c>0")                \
+    F(RV_WEIBULL, "RV.WEIBULL", d_dd, "a>0 b>0")                        \
+    F(CDF_BERNOULLI, "CDF.BERNOULLI", m_ed, "ai[0,1] b[0,1]")           \
+    F(PDF_BERNOULLI, "PDF.BERNOULLI", m_ed, "ai[0,1] b[0,1]")           \
+    F(RV_BERNOULLI, "RV.BERNOULLI", d_d, "a[0,1]")                      \
+    F(CDF_BINOM, "CDF.BINOM", m_edd, "bi>0 c[0,1]")                     \
+    F(PDF_BINOM, "PDF.BINOM", m_edd, "ai>=0<=b bi>0 c[0,1]")            \
+    F(RV_BINOM, "RV.BINOM", d_dd, "ai>0 b[0,1]")                        \
+    F(CDF_GEOM, "CDF.GEOM", m_ed, "ai>=1 b[0,1]")                       \
+    F(PDF_GEOM, "PDF.GEOM", m_ed, "ai>=1 b[0,1]")                       \
+    F(RV_GEOM, "RV.GEOM", d_d, "a[0,1]")                                \
+    F(CDF_HYPER, "CDF.HYPER", m_eddd, "ai>=0<=d bi>0 ci>0<=b di>0<=b")  \
+    F(PDF_HYPER, "PDF.HYPER", m_eddd, "ai>=0<=d bi>0 ci>0<=b di>0<=b")  \
+    F(RV_HYPER, "RV.HYPER", d_ddd, "ai>0 bi>0<=a ci>0<=a")              \
+    F(PDF_LOG, "PDF.LOG", m_ed, "a>=1 b(0,1]")                          \
+    F(RV_LOG, "RV.LOG", d_d, "a(0,1]")                                  \
+    F(CDF_NEGBIN, "CDF.NEGBIN", m_edd, "a>=1 bi c(0,1]")                \
+    F(PDF_NEGBIN, "PDF.NEGBIN", m_edd, "a>=1 bi c(0,1]")                \
+    F(RV_NEGBIN, "RV.NEGBIN", d_dd, "ai b(0,1]")                        \
+    F(CDF_POISSON, "CDF.POISSON", m_ed, "ai>=0 b>0")                    \
+    F(PDF_POISSON, "PDF.POISSON", m_ed, "ai>=0 b>0")                    \
+    F(RV_POISSON, "RV.POISSON", d_d, "a>0")
+
+struct matrix_function_properties
+  {
+    const char *name;
+    const char *constraints;
+  };
 
+enum { d_none_MIN_ARGS = 0, d_none_MAX_ARGS = 0 };
+enum { d_d_MIN_ARGS = 1, d_d_MAX_ARGS = 1 };
+enum { d_dd_MIN_ARGS = 2, d_dd_MAX_ARGS = 2 };
+enum { d_ddd_MIN_ARGS = 3, d_ddd_MAX_ARGS = 3 };
 enum { m_d_MIN_ARGS = 1, m_d_MAX_ARGS = 1 };
 enum { m_dd_MIN_ARGS = 2, m_dd_MAX_ARGS = 2 };
 enum { m_ddd_MIN_ARGS = 3, m_ddd_MAX_ARGS = 3 };
 enum { m_m_MIN_ARGS = 1, m_m_MAX_ARGS = 1 };
+enum { m_e_MIN_ARGS = 1, m_e_MAX_ARGS = 1 };
 enum { m_md_MIN_ARGS = 2, m_md_MAX_ARGS = 2 };
+enum { m_ed_MIN_ARGS = 2, m_ed_MAX_ARGS = 2 };
 enum { m_mdd_MIN_ARGS = 3, m_mdd_MAX_ARGS = 3 };
+enum { m_edd_MIN_ARGS = 3, m_edd_MAX_ARGS = 3 };
+enum { m_eddd_MIN_ARGS = 4, m_eddd_MAX_ARGS = 4 };
+enum { m_eed_MIN_ARGS = 3, m_eed_MAX_ARGS = 3 };
 enum { m_mm_MIN_ARGS = 2, m_mm_MAX_ARGS = 2 };
 enum { m_v_MIN_ARGS = 1, m_v_MAX_ARGS = 1 };
 enum { d_m_MIN_ARGS = 1, d_m_MAX_ARGS = 1 };
 enum { m_any_MIN_ARGS = 1, m_any_MAX_ARGS = INT_MAX };
 enum { IDENT_MIN_ARGS = 1, IDENT_MAX_ARGS = 2 };
 
+typedef double matrix_proto_d_none (void);
+typedef double matrix_proto_d_d (double);
+typedef double matrix_proto_d_dd (double, double);
+typedef double matrix_proto_d_dd (double, double);
+typedef double matrix_proto_d_ddd (double, double, double);
 typedef gsl_matrix *matrix_proto_m_d (double);
 typedef gsl_matrix *matrix_proto_m_dd (double, double);
 typedef gsl_matrix *matrix_proto_m_ddd (double, double, double);
 typedef gsl_matrix *matrix_proto_m_m (gsl_matrix *);
+typedef double matrix_proto_m_e (double);
 typedef gsl_matrix *matrix_proto_m_md (gsl_matrix *, double);
+typedef double matrix_proto_m_ed (double, double);
 typedef gsl_matrix *matrix_proto_m_mdd (gsl_matrix *, double, double);
+typedef double matrix_proto_m_edd (double, double, double);
+typedef double matrix_proto_m_eddd (double, double, double, double);
+typedef double matrix_proto_m_eed (double, double, double);
 typedef gsl_matrix *matrix_proto_m_mm (gsl_matrix *, gsl_matrix *);
 typedef gsl_matrix *matrix_proto_m_v (gsl_vector *);
 typedef double matrix_proto_d_m (gsl_matrix *);
 typedef gsl_matrix *matrix_proto_m_any (gsl_matrix *[], size_t n);
 typedef gsl_matrix *matrix_proto_IDENT (double, double);
 
-#define F(NAME, PROTOTYPE) static matrix_proto_##PROTOTYPE matrix_eval_##NAME;
+#define F(ENUM, STRING, PROTO, CONSTRAINTS) \
+    static matrix_proto_##PROTO matrix_eval_##ENUM;
 MATRIX_FUNCTIONS
 #undef F
 \f
@@ -268,7 +436,7 @@ struct matrix_expr
     enum matrix_op
       {
         /* Functions. */
-#define F(NAME, PROTOTYPE) MOP_F_##NAME,
+#define F(ENUM, STRING, PROTO, CONSTRAINTS) MOP_F_##ENUM,
         MATRIX_FUNCTIONS
 #undef F
 
@@ -333,6 +501,8 @@ struct matrix_expr
         struct matrix_var *variable;
         struct read_file *eof;
       };
+
+    struct msg_location *location;
   };
 
 static void
@@ -343,7 +513,7 @@ matrix_expr_destroy (struct matrix_expr *e)
 
   switch (e->op)
     {
-#define F(NAME, PROTOTYPE) case MOP_F_##NAME:
+#define F(ENUM, STRING, PROTO, CONSTRAINTS) case MOP_F_##ENUM:
 MATRIX_FUNCTIONS
 #undef F
     case MOP_NEGATE:
@@ -384,6 +554,7 @@ MATRIX_FUNCTIONS
     case MOP_EOF:
       break;
     }
+  msg_location_destroy (e->location);
   free (e);
 }
 
@@ -397,6 +568,9 @@ matrix_expr_create_subs (enum matrix_op op, struct matrix_expr **subs,
     .subs = xmemdup (subs, n_subs * sizeof *subs),
     .n_subs = n_subs
   };
+
+  for (size_t i = 0; i < n_subs; i++)
+    msg_location_merge (&e->location, subs[i]->location);
   return e;
 }
 
@@ -502,12 +676,10 @@ to_vector (gsl_matrix *m)
 }
 
 
-static gsl_matrix *
-matrix_eval_ABS (gsl_matrix *m)
+static double
+matrix_eval_ABS (double d)
 {
-  MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
-    *d = fabs (*d);
-  return m;
+  return fabs (d);
 }
 
 static double
@@ -528,20 +700,16 @@ matrix_eval_ANY (gsl_matrix *m)
   return 0.0;
 }
 
-static gsl_matrix *
-matrix_eval_ARSIN (gsl_matrix *m)
+static double
+matrix_eval_ARSIN (double d)
 {
-  MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
-    *d = asin (*d);
-  return m;
+  return asin (d);
 }
 
-static gsl_matrix *
-matrix_eval_ARTAN (gsl_matrix *m)
+static double
+matrix_eval_ARTAN (double d)
 {
-  MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
-    *d = atan (*d);
-  return m;
+  return atan (d);
 }
 
 static gsl_matrix *
@@ -623,12 +791,10 @@ matrix_eval_CMIN (gsl_matrix *m)
   return matrix_eval_col_extremum (m, true);
 }
 
-static gsl_matrix *
-matrix_eval_COS (gsl_matrix *m)
+static double
+matrix_eval_COS (double d)
 {
-  MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
-    *d = cos (*d);
-  return m;
+  return cos (d);
 }
 
 static gsl_matrix *
@@ -788,12 +954,10 @@ matrix_eval_EVAL (gsl_matrix *m)
   return eval;
 }
 
-static gsl_matrix *
-matrix_eval_EXP (gsl_matrix *m)
+static double
+matrix_eval_EXP (double d)
 {
-  MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
-    *d = exp (*d);
-  return m;
+  return exp (d);
 }
 
 /* From https://gist.github.com/turingbirds/5e99656e08dbe1324c99, where it was
@@ -1039,20 +1203,16 @@ matrix_eval_KRONEKER (gsl_matrix *a, gsl_matrix *b)
   return k;
 }
 
-static gsl_matrix *
-matrix_eval_LG10 (gsl_matrix *m)
+static double
+matrix_eval_LG10 (double d)
 {
-  MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
-    *d = log10 (*d);
-  return m;
+  return log10 (d);
 }
 
-static gsl_matrix *
-matrix_eval_LN (gsl_matrix *m)
+static double
+matrix_eval_LN (double d)
 {
-  MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
-    *d = log (*d);
-  return m;
+  return log (d);
 }
 
 static void
@@ -1204,11 +1364,6 @@ matrix_eval_MAGIC_singly_even (gsl_matrix *m, size_t n)
 static gsl_matrix *
 matrix_eval_MAGIC (double n_)
 {
-  if (n_ < 3 || n_ >= sqrt (SIZE_MAX))
-    {
-      msg (SE, _("MAGIC argument must be an integer 3 or greater."));
-      return NULL;
-    }
   size_t n = n_;
 
   gsl_matrix *m = gsl_matrix_calloc (n, n);
@@ -1378,12 +1533,10 @@ matrix_eval_RMIN (gsl_matrix *m)
   return matrix_eval_row_extremum (m, true);
 }
 
-static gsl_matrix *
-matrix_eval_RND (gsl_matrix *m)
+static double
+matrix_eval_RND (double d)
 {
-  MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
-    *d = rint (*d);
-  return m;
+  return rint (d);
 }
 
 struct rank
@@ -1464,12 +1617,10 @@ matrix_eval_RSUM (gsl_matrix *m)
   return matrix_eval_row_sum (m, false);
 }
 
-static gsl_matrix *
-matrix_eval_SIN (gsl_matrix *m)
+static double
+matrix_eval_SIN (double d)
 {
-  MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
-    *d = sin (*d);
-  return m;
+  return sin (d);
 }
 
 static gsl_matrix *
@@ -1499,19 +1650,10 @@ matrix_eval_SOLVE (gsl_matrix *m1, gsl_matrix *m2)
   return x;
 }
 
-static gsl_matrix *
-matrix_eval_SQRT (gsl_matrix *m)
+static double
+matrix_eval_SQRT (double d)
 {
-  MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
-    {
-      if (*d < 0)
-        {
-          msg (SE, _("Argument to SQRT must be nonnegative."));
-          return NULL;
-        }
-      *d = sqrt (*d);
-    }
-  return m;
+  return sqrt (d);
 }
 
 static gsl_matrix *
@@ -1627,22 +1769,15 @@ matrix_eval_TRANSPOS (gsl_matrix *m)
     }
 }
 
-static gsl_matrix *
-matrix_eval_TRUNC (gsl_matrix *m)
+static double
+matrix_eval_TRUNC (double d)
 {
-  MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
-    *d = trunc (*d);
-  return m;
+  return trunc (d);
 }
 
 static gsl_matrix *
 matrix_eval_UNIFORM (double r_, double c_)
 {
-  if (r_ < 0 || r_ >= SIZE_MAX || c_ < 0 || c_ >= SIZE_MAX)
-    {
-      msg (SE, _("Arguments to UNIFORM must be integers."));
-      return NULL;
-    }
   size_t r = r_;
   size_t c = c_;
   if (size_overflow_p (xtimes (r, xmax (c, 1))))
@@ -1657,189 +1792,934 @@ matrix_eval_UNIFORM (double r_, double c_)
   return m;
 }
 
-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 double
+matrix_eval_PDF_BETA (double x, double a, double b)
+{
+  return gsl_ran_beta_pdf (x, a, b);
+}
 
-static const struct matrix_function *
-matrix_parse_function_name (const char *token)
+static double
+matrix_eval_CDF_BETA (double x, double a, double b)
 {
-  static const struct matrix_function functions[] =
-    {
-#define F(NAME, PROTO) { #NAME, MOP_F_##NAME, PROTO##_MIN_ARGS, PROTO##_MAX_ARGS },
-      MATRIX_FUNCTIONS
-#undef F
-    };
-  enum { N_FUNCTIONS = sizeof functions / sizeof *functions };
+  return gsl_cdf_beta_P (x, a, b);
+}
 
-  for (size_t i = 0; i < N_FUNCTIONS; i++)
-    {
-      if (lex_id_match_n (ss_cstr (functions[i].name), ss_cstr (token), 3))
-        return &functions[i];
-    }
-  return NULL;
+static double
+matrix_eval_IDF_BETA (double P, double a, double b)
+{
+  return gsl_cdf_beta_Pinv (P, a, b);
 }
 
-static struct read_file *
-read_file_create (struct matrix_state *s, struct file_handle *fh)
+static double
+matrix_eval_RV_BETA (double a, double b)
 {
-  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;
-        }
-    }
+  return gsl_ran_beta (get_rng (), a, b);
+}
 
-  struct read_file *rf = xmalloc (sizeof *rf);
-  *rf = (struct read_file) { .file = fh };
+static double
+matrix_eval_NCDF_BETA (double x, double a, double b, double lambda)
+{
+  return ncdf_beta (x, a, b, lambda);
+}
 
-  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 double
+matrix_eval_NPDF_BETA (double x, double a, double b, double lambda)
+{
+  return npdf_beta (x, a, b, lambda);
 }
 
-static struct dfm_reader *
-read_file_open (struct read_file *rf)
+static double
+matrix_eval_CDF_BVNOR (double x0, double x1, double r)
 {
-  if (!rf->reader)
-    rf->reader = dfm_open_reader (rf->file, NULL, rf->encoding);
-  return rf->reader;
+  return cdf_bvnor (x0, x1, r);
 }
 
-static void
-read_file_destroy (struct read_file *rf)
+static double
+matrix_eval_PDF_BVNOR (double x0, double x1, double r)
 {
-  if (rf)
-    {
-      fh_unref (rf->file);
-      dfm_close_reader (rf->reader);
-      free (rf->encoding);
-      free (rf);
-    }
+  return gsl_ran_bivariate_gaussian_pdf (x0, x1, 1, 1, r);
 }
 
-static struct write_file *
-write_file_create (struct matrix_state *s, struct file_handle *fh)
+static double
+matrix_eval_CDF_CAUCHY (double x, double a, double b)
 {
-  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;
-        }
-    }
+  return gsl_cdf_cauchy_P ((x - a) / b, 1);
+}
 
-  struct write_file *wf = xmalloc (sizeof *wf);
-  *wf = (struct write_file) { .file = fh };
+static double
+matrix_eval_IDF_CAUCHY (double P, double a, double b)
+{
+  return a + b * gsl_cdf_cauchy_Pinv (P, 1);
+}
 
-  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 double
+matrix_eval_PDF_CAUCHY (double x, double a, double b)
+{
+  return gsl_ran_cauchy_pdf ((x - a) / b, 1) / b;
 }
 
-static struct dfm_writer *
-write_file_open (struct write_file *wf)
+static double
+matrix_eval_RV_CAUCHY (double a, double b)
 {
-  if (!wf->writer)
-    wf->writer = dfm_open_writer (wf->file, wf->encoding);
-  return wf->writer;
+  return a + b * gsl_ran_cauchy (get_rng (), 1);
 }
 
-static void
-write_file_destroy (struct write_file *wf)
+static double
+matrix_eval_CDF_CHISQ (double x, double df)
 {
-  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);
-        }
+  return gsl_cdf_chisq_P (x, df);
+}
 
-      fh_unref (wf->file);
-      dfm_close_writer (wf->writer);
-      free (wf->encoding);
-      free (wf);
-    }
+static double
+matrix_eval_CHICDF (double x, double df)
+{
+  return matrix_eval_CDF_CHISQ (x, df);
 }
 
-static bool
-matrix_parse_function (struct matrix_state *s, const char *token,
-                       struct matrix_expr **exprp)
+static double
+matrix_eval_IDF_CHISQ (double P, double df)
 {
-  *exprp = NULL;
-  if (lex_next_token (s->lexer, 1) != T_LPAREN)
-    return false;
+  return gsl_cdf_chisq_Pinv (P, df);
+}
 
-  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;
+static double
+matrix_eval_PDF_CHISQ (double x, double df)
+{
+  return gsl_ran_chisq_pdf (x, df);
+}
 
-      if (!lex_force_match (s->lexer, T_RPAREN))
-        {
-          fh_unref (fh);
-          return true;
-        }
+static double
+matrix_eval_RV_CHISQ (double df)
+{
+  return gsl_ran_chisq (get_rng (), df);
+}
 
-      struct read_file *rf = read_file_create (s, fh);
+static double
+matrix_eval_SIG_CHISQ (double x, double df)
+{
+  return gsl_cdf_chisq_Q (x, df);
+}
 
-      struct matrix_expr *e = xmalloc (sizeof *e);
-      *e = (struct matrix_expr) { .op = MOP_EOF, .eof = rf };
-      *exprp = e;
-      return true;
-    }
+static double
+matrix_eval_CDF_EXP (double x, double a)
+{
+  return gsl_cdf_exponential_P (x, 1. / a);
+}
 
-  const struct matrix_function *f = matrix_parse_function_name (token);
-  if (!f)
-    return false;
+static double
+matrix_eval_IDF_EXP (double P, double a)
+{
+  return gsl_cdf_exponential_Pinv (P, 1. / a);
+}
 
-  lex_get_n (s->lexer, 2);
+static double
+matrix_eval_PDF_EXP (double x, double a)
+{
+  return gsl_ran_exponential_pdf (x, 1. / a);
+}
 
-  struct matrix_expr *e = xmalloc (sizeof *e);
-  *e = (struct matrix_expr) { .op = f->op, .subs = NULL };
+static double
+matrix_eval_RV_EXP (double a)
+{
+  return gsl_ran_exponential (get_rng (), 1. / a);
+}
 
-  size_t allocated_subs = 0;
-  do
-    {
-      struct matrix_expr *sub = matrix_parse_expr (s);
-      if (!sub)
-        goto error;
+static double
+matrix_eval_PDF_XPOWER (double x, double a, double b)
+{
+  return gsl_ran_exppow_pdf (x, a, b);
+}
 
-      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;
+static double
+matrix_eval_RV_XPOWER (double a, double b)
+{
+  return gsl_ran_exppow (get_rng (), a, b);
+}
 
-  *exprp = e;
-  return true;
+static double
+matrix_eval_CDF_F (double x, double df1, double df2)
+{
+  return gsl_cdf_fdist_P (x, df1, df2);
+}
 
-error:
-  matrix_expr_destroy (e);
-  return true;
+static double
+matrix_eval_FCDF (double x, double df1, double df2)
+{
+  return matrix_eval_CDF_F (x, df1, df2);
 }
 
-static struct matrix_expr *
-matrix_parse_primary (struct matrix_state *s)
+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;
+
+  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 };
+      *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,
+    .location = lex_get_location (s->lexer, 0, 0)
+  };
+
+  lex_get_n (s->lexer, 2);
+  if (lex_token (s->lexer) != T_RPAREN)
+    {
+      size_t allocated_subs = 0;
+      do
+        {
+          struct msg_location *arg_location = lex_get_location (s->lexer, 0, 0);
+          struct matrix_expr *sub = matrix_parse_expr (s);
+          if (!sub)
+            {
+              msg_location_destroy (arg_location);
+              goto error;
+            }
+          if (!sub->location)
+            {
+              lex_extend_location (s->lexer, 0, arg_location);
+              sub->location = arg_location;
+            }
+          else
+            msg_location_destroy (arg_location);
+
+          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 (SE, 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 (SE, 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 (SE, _("Matrix function %s requires at least one argument."),
+             f->name);
+      else
+        NOT_REACHED ();
+
+      goto error;
+    }
+
+  *exprp = e;
+  return true;
+
+error:
+  matrix_expr_destroy (e);
+  return true;
+}
+
+static struct matrix_expr *
+matrix_parse_primary (struct matrix_state *s)
 {
   if (lex_is_number (s->lexer))
     {
@@ -2216,7 +3096,7 @@ matrix_op_eval (enum matrix_op op, double a, double b)
     case MOP_OR: return (a > 0) || (b > 0);
     case MOP_XOR: return (a > 0) != (b > 0);
 
-#define F(NAME, PROTOTYPE) case MOP_F_##NAME:
+#define F(ENUM, STRING, PROTO, CONSTRAINTS) case MOP_F_##ENUM:
       MATRIX_FUNCTIONS
 #undef F
     case MOP_NEGATE:
@@ -2261,7 +3141,7 @@ matrix_op_name (enum matrix_op op)
     case MOP_OR: return "OR";
     case MOP_XOR: return "XOR";
 
-#define F(NAME, PROTOTYPE) case MOP_F_##NAME:
+#define F(ENUM, STRING, PROTO, CONSTRAINTS) case MOP_F_##ENUM:
       MATRIX_FUNCTIONS
 #undef F
     case MOP_NEGATE:
@@ -2719,10 +3599,10 @@ matrix_expr_evaluate_mat_index (gsl_matrix *sm, gsl_matrix *im0,
   return dm;
 }
 
-#define F(NAME, PROTOTYPE)                              \
-  static gsl_matrix *matrix_expr_evaluate_##PROTOTYPE ( \
-    const char *name, gsl_matrix *[], size_t,           \
-    matrix_proto_##PROTOTYPE *);
+#define F(ENUM, STRING, PROTO, CONSTRAINTS)                             \
+  static gsl_matrix *matrix_expr_evaluate_##PROTO (                     \
+    const struct matrix_function_properties *, gsl_matrix *[], size_t,  \
+    matrix_proto_##PROTO *);
 MATRIX_FUNCTIONS
 #undef F
 
@@ -2764,41 +3644,311 @@ to_scalar_args (const char *name, gsl_matrix *subs[], size_t n_subs, double d[])
   return true;
 }
 
+static int
+parse_constraint_value (const char **constraintsp)
+{
+  char *tail;
+  long retval = strtol (*constraintsp, &tail, 10);
+  assert (tail > *constraintsp);
+  *constraintsp = tail;
+  return retval;
+}
+
+static void
+argument_invalid (const struct matrix_function_properties *props,
+                  size_t arg_index, gsl_matrix *a, size_t y, size_t x,
+                  struct string *out)
+{
+  if (is_scalar (a))
+    ds_put_format (out, _("Argument %zu to matrix function %s "
+                          "has invalid value %g."),
+                   arg_index, props->name, gsl_matrix_get (a, y, x));
+  else
+    ds_put_format (out, _("Row %zu, column %zu of argument %zu "
+                          "to matrix function %s has "
+                          "invalid value %g."),
+                   y + 1, x + 1, arg_index, props->name,
+                   gsl_matrix_get (a, y, x));
+}
+
+static void
+argument_inequality_error (const struct matrix_function_properties *props,
+                           size_t a_index, gsl_matrix *a, size_t y, size_t x,
+                           size_t b_index, double b,
+                           bool greater, bool equal)
+{
+  struct string s = DS_EMPTY_INITIALIZER;
+
+  argument_invalid (props, a_index, a, y, x, &s);
+  ds_put_cstr (&s, "  ");
+  if (greater && equal)
+    ds_put_format (&s, _("This argument must be greater than or "
+                         "equal to argument %zu, but that has value %g."),
+                   b_index, b);
+  else if (greater && !equal)
+    ds_put_format (&s, _("This argument must be greater than argument %zu, "
+                         "but that has value %g."),
+                   b_index, b);
+  else if (equal)
+    ds_put_format (&s, _("This argument must be less than or "
+                         "equal to argument %zu, but that has value %g."),
+                   b_index, b);
+  else
+    ds_put_format (&s, _("This argument must be less than argument %zu, "
+                         "but that has value %g."),
+                   b_index, b);
+  msg (ME, ds_cstr (&s));
+  ds_destroy (&s);
+}
+
+static void
+argument_value_error (const struct matrix_function_properties *props,
+                      size_t a_index, gsl_matrix *a, size_t y, size_t x,
+                      double b,
+                      bool greater, bool equal)
+{
+  struct string s = DS_EMPTY_INITIALIZER;
+  argument_invalid (props, a_index, a, y, x, &s);
+  ds_put_cstr (&s, "  ");
+  if (greater && equal)
+    ds_put_format (&s, _("This argument must be greater than or equal to %g."),
+                   b);
+  else if (greater && !equal)
+    ds_put_format (&s, _("This argument must be greater than %g."), b);
+  else if (equal)
+    ds_put_format (&s, _("This argument must be less than or equal to %g."), b);
+  else
+    ds_put_format (&s, _("This argument must be less than %g."), b);
+  msg (ME, ds_cstr (&s));
+  ds_destroy (&s);
+}
+
+static bool
+check_constraints (const struct matrix_function_properties *props,
+                   gsl_matrix *args[], size_t n_args)
+{
+  const char *constraints = props->constraints;
+  if (!constraints)
+    return true;
+
+  size_t arg_index = SIZE_MAX;
+  while (*constraints)
+    {
+      if (*constraints >= 'a' && *constraints <= 'd')
+        {
+          arg_index = *constraints++ - 'a';
+          assert (arg_index < n_args);
+        }
+      else if (*constraints == '[' || *constraints == '(')
+        {
+          assert (arg_index < n_args);
+          bool open_lower = *constraints++ == '(';
+          int minimum = parse_constraint_value (&constraints);
+          assert (*constraints == ',');
+          constraints++;
+          int maximum = parse_constraint_value (&constraints);
+          assert (*constraints == ']' || *constraints == ')');
+          bool open_upper = *constraints++ == ')';
+
+          MATRIX_FOR_ALL_ELEMENTS (d, y, x, args[arg_index])
+            if ((open_lower ? *d <= minimum : *d < minimum)
+                || (open_upper ? *d >= maximum : *d > maximum))
+              {
+                if (!is_scalar (args[arg_index]))
+                  msg (ME, _("Row %zu, column %zu of argument %zu to matrix "
+                             "function %s has value %g, which is outside "
+                             "the valid range %c%d,%d%c."),
+                       y + 1, x + 1, arg_index + 1, props->name, *d,
+                       open_lower ? '(' : '[',
+                       minimum, maximum,
+                       open_upper ? ')' : ']');
+                else
+                  msg (ME, _("Argument %zu to matrix function %s has value %g, "
+                             "which is outside the valid range %c%d,%d%c."),
+                       arg_index + 1, props->name, *d,
+                       open_lower ? '(' : '[',
+                       minimum, maximum,
+                       open_upper ? ')' : ']');
+                return false;
+              }
+        }
+      else if (*constraints == '>' || *constraints == '<')
+        {
+          bool greater = *constraints++ == '>';
+          bool equal;
+          if (*constraints == '=')
+            {
+              equal = true;
+              constraints++;
+            }
+          else
+            equal = false;
+
+          if (*constraints >= 'a' && *constraints <= 'd')
+            {
+              size_t a_index = arg_index;
+              size_t b_index = *constraints - 'a';
+              assert (a_index < n_args);
+              assert (b_index < n_args);
+
+              /* We only support one of the two arguments being non-scalar.
+                 It's easier to support only the first one being non-scalar, so
+                 flip things around if it's the other way. */
+              if (!is_scalar (args[b_index]))
+                {
+                  assert (is_scalar (args[a_index]));
+                  size_t tmp_index = a_index;
+                  a_index = b_index;
+                  b_index = tmp_index;
+
+                  greater = !greater;
+                }
+
+              double b = to_scalar (args[b_index]);
+              MATRIX_FOR_ALL_ELEMENTS (a, y, x, args[a_index])
+                if (greater
+                    ? (equal ? !(*a >= b) : !(*a > b))
+                    : (equal ? !(*a <= b) : !(*a < b)))
+                  {
+                    argument_inequality_error (props,
+                                               a_index + 1, args[a_index], y, x,
+                                               b_index + 1, b,
+                                               greater, equal);
+                    return false;
+                  }
+            }
+          else
+            {
+              int comparand = parse_constraint_value (&constraints);
+
+              MATRIX_FOR_ALL_ELEMENTS (d, y, x, args[arg_index])
+                if (greater
+                    ? (equal ? !(*d >= comparand) : !(*d > comparand))
+                    : (equal ? !(*d <= comparand) : !(*d < comparand)))
+                  {
+                    argument_value_error (props,
+                                          arg_index + 1, args[arg_index], y, x,
+                                          comparand,
+                                          greater, equal);
+                    return false;
+                  }
+            }
+        }
+      else
+        {
+          assert (*constraints == ' ');
+          constraints++;
+          arg_index = SIZE_MAX;
+        }
+    }
+  return true;
+}
+
+static gsl_matrix *
+matrix_expr_evaluate_d_none (
+  const struct matrix_function_properties *props UNUSED,
+  gsl_matrix *subs[] UNUSED, size_t n_subs,
+  matrix_proto_d_none *f)
+{
+  assert (n_subs == 0);
+
+  gsl_matrix *m = gsl_matrix_alloc (1, 1);
+  gsl_matrix_set (m, 0, 0, f ());
+  return m;
+}
+
+static gsl_matrix *
+matrix_expr_evaluate_d_d (const struct matrix_function_properties *props,
+                          gsl_matrix *subs[], size_t n_subs,
+                          matrix_proto_d_d *f)
+{
+  assert (n_subs == 1);
+
+  double d;
+  if (!to_scalar_args (props->name, subs, n_subs, &d))
+    return NULL;
+
+  if (!check_constraints (props, subs, n_subs))
+    return NULL;
+
+  gsl_matrix *m = gsl_matrix_alloc (1, 1);
+  gsl_matrix_set (m, 0, 0, f (d));
+  return m;
+}
+
+static gsl_matrix *
+matrix_expr_evaluate_d_dd (const struct matrix_function_properties *props,
+                           gsl_matrix *subs[], size_t n_subs,
+                           matrix_proto_d_dd *f)
+{
+  assert (n_subs == 2);
+
+  double d[2];
+  if (!to_scalar_args (props->name, subs, n_subs, d))
+    return NULL;
+
+  if (!check_constraints (props, subs, n_subs))
+    return NULL;
+
+  gsl_matrix *m = gsl_matrix_alloc (1, 1);
+  gsl_matrix_set (m, 0, 0, f (d[0], d[1]));
+  return m;
+}
+
+static gsl_matrix *
+matrix_expr_evaluate_d_ddd (const struct matrix_function_properties *props,
+                            gsl_matrix *subs[], size_t n_subs,
+                            matrix_proto_d_ddd *f)
+{
+  assert (n_subs == 3);
+
+  double d[3];
+  if (!to_scalar_args (props->name, subs, n_subs, d))
+    return NULL;
+
+  if (!check_constraints (props, subs, n_subs))
+    return NULL;
+
+  gsl_matrix *m = gsl_matrix_alloc (1, 1);
+  gsl_matrix_set (m, 0, 0, f (d[0], d[1], d[2]));
+  return m;
+}
+
 static gsl_matrix *
-matrix_expr_evaluate_m_d (const char *name,
+matrix_expr_evaluate_m_d (const struct matrix_function_properties *props,
                           gsl_matrix *subs[], size_t n_subs,
                           matrix_proto_m_d *f)
 {
   assert (n_subs == 1);
 
   double d;
-  return to_scalar_args (name, subs, n_subs, &d) ? f(d) : NULL;
+  return to_scalar_args (props->name, subs, n_subs, &d) ? f (d) : NULL;
 }
 
 static gsl_matrix *
-matrix_expr_evaluate_m_dd (const char *name,
+matrix_expr_evaluate_m_dd (const struct matrix_function_properties *props,
                            gsl_matrix *subs[], size_t n_subs,
                            matrix_proto_m_dd *f)
 {
   assert (n_subs == 2);
 
   double d[2];
-  return to_scalar_args (name, subs, n_subs, d) ? f(d[0], d[1]) : NULL;
+  return to_scalar_args (props->name, subs, n_subs, d) ? f(d[0], d[1]) : NULL;
 }
 
 static gsl_matrix *
-matrix_expr_evaluate_m_ddd (const char *name,
+matrix_expr_evaluate_m_ddd (const struct matrix_function_properties *props,
                            gsl_matrix *subs[], size_t n_subs,
                            matrix_proto_m_ddd *f)
 {
   assert (n_subs == 3);
 
   double d[3];
-  return to_scalar_args (name, subs, n_subs, d) ? f(d[0], d[1], d[2]) : NULL;
+  return to_scalar_args (props->name, subs, n_subs, d) ? f(d[0], d[1], d[2]) : NULL;
 }
 
 static gsl_matrix *
-matrix_expr_evaluate_m_m (const char *name UNUSED,
+matrix_expr_evaluate_m_m (const struct matrix_function_properties *props UNUSED,
                           gsl_matrix *subs[], size_t n_subs,
                           matrix_proto_m_m *f)
 {
@@ -2807,29 +3957,145 @@ matrix_expr_evaluate_m_m (const char *name UNUSED,
 }
 
 static gsl_matrix *
-matrix_expr_evaluate_m_md (const char *name UNUSED,
+matrix_expr_evaluate_m_e (const struct matrix_function_properties *props,
+                            gsl_matrix *subs[], size_t n_subs,
+                            matrix_proto_m_e *f)
+{
+  assert (n_subs == 1);
+
+  if (!check_constraints (props, subs, n_subs))
+    return NULL;
+
+  MATRIX_FOR_ALL_ELEMENTS (a, y, x, subs[0])
+      *a = f (*a);
+  return subs[0];
+}
+
+static gsl_matrix *
+matrix_expr_evaluate_m_md (const struct matrix_function_properties *props UNUSED,
                            gsl_matrix *subs[], size_t n_subs,
                            matrix_proto_m_md *f)
 {
   assert (n_subs == 2);
-  if (!check_scalar_arg (name, subs, 1))
+  if (!check_scalar_arg (props->name, subs, 1))
     return NULL;
   return f (subs[0], to_scalar (subs[1]));
 }
 
 static gsl_matrix *
-matrix_expr_evaluate_m_mdd (const char *name UNUSED,
+matrix_expr_evaluate_m_ed (const struct matrix_function_properties *props,
+                           gsl_matrix *subs[], size_t n_subs,
+                           matrix_proto_m_ed *f)
+{
+  assert (n_subs == 2);
+  if (!check_scalar_arg (props->name, subs, 1))
+    return NULL;
+
+  if (!check_constraints (props, subs, n_subs))
+    return NULL;
+
+  double b = to_scalar (subs[1]);
+  MATRIX_FOR_ALL_ELEMENTS (a, y, x, subs[0])
+    *a = f (*a, b);
+  return subs[0];
+}
+
+static gsl_matrix *
+matrix_expr_evaluate_m_mdd (const struct matrix_function_properties *props UNUSED,
                             gsl_matrix *subs[], size_t n_subs,
                             matrix_proto_m_mdd *f)
 {
   assert (n_subs == 3);
-  if (!check_scalar_arg (name, subs, 1) || !check_scalar_arg (name, subs, 2))
+  if (!check_scalar_arg (props->name, subs, 1) || !check_scalar_arg (props->name, subs, 2))
     return NULL;
   return f (subs[0], to_scalar (subs[1]), to_scalar (subs[2]));
 }
 
 static gsl_matrix *
-matrix_expr_evaluate_m_mm (const char *name UNUSED,
+matrix_expr_evaluate_m_edd (const struct matrix_function_properties *props,
+                            gsl_matrix *subs[], size_t n_subs,
+                            matrix_proto_m_edd *f)
+{
+  assert (n_subs == 3);
+  if (!check_scalar_arg (props->name, subs, 1) || !check_scalar_arg (props->name, subs, 2))
+    return NULL;
+
+  if (!check_constraints (props, subs, n_subs))
+    return NULL;
+
+  double b = to_scalar (subs[1]);
+  double c = to_scalar (subs[2]);
+  MATRIX_FOR_ALL_ELEMENTS (a, y, x, subs[0])
+    *a = f (*a, b, c);
+  return subs[0];
+}
+
+static gsl_matrix *
+matrix_expr_evaluate_m_eddd (const struct matrix_function_properties *props,
+                             gsl_matrix *subs[], size_t n_subs,
+                             matrix_proto_m_eddd *f)
+{
+  assert (n_subs == 4);
+  for (size_t i = 1; i < 4; i++)
+    if (!check_scalar_arg (props->name, subs, i))
+    return NULL;
+
+  if (!check_constraints (props, subs, n_subs))
+    return NULL;
+
+  double b = to_scalar (subs[1]);
+  double c = to_scalar (subs[2]);
+  double d = to_scalar (subs[3]);
+  MATRIX_FOR_ALL_ELEMENTS (a, y, x, subs[0])
+    *a = f (*a, b, c, d);
+  return subs[0];
+}
+
+static gsl_matrix *
+matrix_expr_evaluate_m_eed (const struct matrix_function_properties *props,
+                            gsl_matrix *subs[], size_t n_subs,
+                            matrix_proto_m_eed *f)
+{
+  assert (n_subs == 3);
+  if (!check_scalar_arg (props->name, subs, 2))
+    return NULL;
+
+  if (!is_scalar (subs[0]) && !is_scalar (subs[1])
+      && (subs[0]->size1 != subs[1]->size1 || subs[0]->size2 != subs[1]->size2))
+    {
+      msg (ME, _("Arguments 1 and 2 to %s have dimensions %zu×%zu and "
+                 "%zu×%zu, but %s requires these arguments either to have "
+                 "the same dimensions or for one of them to be a scalar."),
+           props->name,
+           subs[0]->size1, subs[0]->size2,
+           subs[1]->size1, subs[1]->size2,
+           props->name);
+      return NULL;
+    }
+
+  if (!check_constraints (props, subs, n_subs))
+    return NULL;
+
+  double c = to_scalar (subs[2]);
+
+  if (is_scalar (subs[0]))
+    {
+      double a = to_scalar (subs[0]);
+      MATRIX_FOR_ALL_ELEMENTS (b, y, x, subs[1])
+        *b = f (a, *b, c);
+      return subs[1];
+    }
+  else
+    {
+      double b = to_scalar (subs[1]);
+      MATRIX_FOR_ALL_ELEMENTS (a, y, x, subs[0])
+        *a = f (*a, b, c);
+      return subs[0];
+    }
+}
+
+static gsl_matrix *
+matrix_expr_evaluate_m_mm (const struct matrix_function_properties *props UNUSED,
                            gsl_matrix *subs[], size_t n_subs,
                            matrix_proto_m_mm *f)
 {
@@ -2838,19 +4104,19 @@ matrix_expr_evaluate_m_mm (const char *name UNUSED,
 }
 
 static gsl_matrix *
-matrix_expr_evaluate_m_v (const char *name,
+matrix_expr_evaluate_m_v (const struct matrix_function_properties *props,
                           gsl_matrix *subs[], size_t n_subs,
                           matrix_proto_m_v *f)
 {
   assert (n_subs == 1);
-  if (!check_vector_arg (name, subs, 0))
+  if (!check_vector_arg (props->name, subs, 0))
     return NULL;
   gsl_vector v = to_vector (subs[0]);
   return f (&v);
 }
 
 static gsl_matrix *
-matrix_expr_evaluate_d_m (const char *name UNUSED,
+matrix_expr_evaluate_d_m (const struct matrix_function_properties *props UNUSED,
                           gsl_matrix *subs[], size_t n_subs,
                           matrix_proto_d_m *f)
 {
@@ -2862,7 +4128,7 @@ matrix_expr_evaluate_d_m (const char *name UNUSED,
 }
 
 static gsl_matrix *
-matrix_expr_evaluate_m_any (const char *name UNUSED,
+matrix_expr_evaluate_m_any (const struct matrix_function_properties *props UNUSED,
                           gsl_matrix *subs[], size_t n_subs,
                           matrix_proto_m_any *f)
 {
@@ -2870,14 +4136,14 @@ matrix_expr_evaluate_m_any (const char *name UNUSED,
 }
 
 static gsl_matrix *
-matrix_expr_evaluate_IDENT (const char *name,
+matrix_expr_evaluate_IDENT (const struct matrix_function_properties *props,
                             gsl_matrix *subs[], size_t n_subs,
                             matrix_proto_IDENT *f)
 {
   assert (n_subs <= 2);
 
   double d[2];
-  if (!to_scalar_args (name, subs, n_subs, d))
+  if (!to_scalar_args (props->name, subs, n_subs, d))
     return NULL;
 
   return f (d[0], d[n_subs - 1]);
@@ -2936,11 +4202,16 @@ matrix_expr_evaluate (const struct matrix_expr *e)
   gsl_matrix *result = NULL;
   switch (e->op)
     {
-#define F(NAME, PROTOTYPE)                                              \
-      case MOP_F_##NAME:                                                \
-        result = matrix_expr_evaluate_##PROTOTYPE (#NAME,               \
-                                                   subs, e->n_subs,     \
-                                                   matrix_eval_##NAME); \
+#define F(ENUM, STRING, PROTO, CONSTRAINTS)                             \
+      case MOP_F_##ENUM:                                                \
+        {                                                               \
+          static const struct matrix_function_properties props = {      \
+            .name = STRING,                                             \
+            .constraints = CONSTRAINTS,                                 \
+          };                                                            \
+          result = matrix_expr_evaluate_##PROTO (&props, subs, e->n_subs, \
+                                                 matrix_eval_##ENUM);   \
+        }                                                               \
       break;
       MATRIX_FUNCTIONS
 #undef F
@@ -3076,6 +4347,7 @@ matrix_expr_evaluate_integer (const struct matrix_expr *e, const char *context,
 struct matrix_lvalue
   {
     struct matrix_var *var;
+    struct msg_location *var_location;
     struct matrix_expr *indexes[2];
     size_t n_indexes;
   };
@@ -3097,6 +4369,7 @@ matrix_lvalue_parse (struct matrix_state *s)
   struct matrix_lvalue *lvalue = xzalloc (sizeof *lvalue);
   if (!lex_force_id (s->lexer))
     goto error;
+  lvalue->var_location = lex_get_location (s->lexer, 0, 0);
   lvalue->var = matrix_var_lookup (s, lex_tokss (s->lexer));
   if (lex_next_token (s->lexer, 1) == T_LPAREN)
     {
@@ -3257,21 +4530,23 @@ matrix_lvalue_evaluate (struct matrix_lvalue *lvalue,
   gsl_matrix *dm = lvalue->var->value;
   if (!dm)
     {
-      msg (SE, _("Undefined variable %s."), lvalue->var->name);
+      msg_at (SE, lvalue->var_location,
+              _("Undefined variable %s."), lvalue->var->name);
       return false;
     }
   else if (dm->size1 == 0 || dm->size2 == 0)
     {
-      msg (SE, _("Cannot index %zu×%zu matrix."),
-           dm->size1, dm->size2);
+      msg_at (SE, lvalue->var_location,
+              _("Cannot index %zu×%zu matrix."), dm->size1, dm->size2);
       return false;
     }
   else if (lvalue->n_indexes == 1)
     {
       if (!is_vector (dm))
         {
-          msg (SE, _("Can't use vector indexing on %zu×%zu matrix %s."),
-               dm->size1, dm->size2, lvalue->var->name);
+          msg_at (SE, lvalue->var_location,
+                  _("Can't use vector indexing on %zu×%zu matrix %s."),
+                  dm->size1, dm->size2, lvalue->var->name);
           return false;
         }
       return matrix_lvalue_evaluate_vector (lvalue->indexes[0],