X-Git-Url: https://pintos-os.org/cgi-bin/gitweb.cgi?a=blobdiff_plain;f=src%2Flanguage%2Fstats%2Fmatrix.c;h=076d86ae5849495dbaff87589c2cc82fdf773d4a;hb=09e02e5fc7f463a4243de3daaa9e8f2c6fb2d619;hp=a33e7693cd55dfb9558345ba6b33ddb1e603b1d7;hpb=d5c2af4608820668efdf4d40da5f5e597cebc05c;p=pspp diff --git a/src/language/stats/matrix.c b/src/language/stats/matrix.c index a33e7693cd..076d86ae58 100644 --- a/src/language/stats/matrix.c +++ b/src/language/stats/matrix.c @@ -17,6 +17,7 @@ #include #include +#include #include #include #include @@ -48,17 +49,20 @@ #include "libpspp/compiler.h" #include "libpspp/hmap.h" #include "libpspp/i18n.h" +#include "libpspp/intern.h" #include "libpspp/misc.h" #include "libpspp/str.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 +181,267 @@ matrix_var_set (struct matrix_var *var, gsl_matrix *value) var->value = value; } -#define MATRIX_FUNCTIONS \ - F(ABS, "ABS", m_m) \ - F(ALL, "ALL", d_m) \ - F(ANY, "ANY", d_m) \ - F(ARSIN, "ARSIN", m_m) \ - F(ARTAN, "ARTAN", m_m) \ - F(BLOCK, "BLOCK", m_any) \ - F(CHOL, "CHOL", m_m) \ - F(CMIN, "CMIN", m_m) \ - F(CMAX, "CMAX", m_m) \ - F(COS, "COS", m_m) \ - F(CSSQ, "CSSQ", m_m) \ - F(CSUM, "CSUM", m_m) \ - F(DESIGN, "DESIGN", m_m) \ - F(DET, "DET", d_m) \ - F(DIAG, "DIAG", m_m) \ - F(EVAL, "EVAL", m_m) \ - F(EXP, "EXP", m_m) \ - F(GINV, "GINV", m_m) \ - F(GRADE, "GRADE", m_m) \ - F(GSCH, "GSCH", m_m) \ - F(IDENT, "IDENT", IDENT) \ - F(INV, "INV", m_m) \ - F(KRONEKER, "KRONEKER", m_mm) \ - F(LG10, "LG10", m_m) \ - F(LN, "LN", m_m) \ - F(MAGIC, "MAGIC", m_d) \ - F(MAKE, "MAKE", m_ddd) \ - F(MDIAG, "MDIAG", m_v) \ - F(MMAX, "MMAX", d_m) \ - F(MMIN, "MMIN", d_m) \ - F(MOD, "MOD", m_md) \ - F(MSSQ, "MSSQ", d_m) \ - F(MSUM, "MSUM", d_m) \ - F(NCOL, "NCOL", d_m) \ - F(NROW, "NROW", d_m) \ - F(RANK, "RANK", d_m) \ - F(RESHAPE, "RESHAPE", m_mdd) \ - F(RMAX, "RMAX", m_m) \ - F(RMIN, "RMIN", m_m) \ - F(RND, "RND", m_m) \ - F(RNKORDER, "RNKORDER", m_m) \ - F(RSSQ, "RSSQ", m_m) \ - F(RSUM, "RSUM", m_m) \ - F(SIN, "SIN", m_m) \ - F(SOLVE, "SOLVE", m_mm) \ - F(SQRT, "SQRT", m_m) \ - F(SSCP, "SSCP", m_m) \ - F(SVAL, "SVAL", m_m) \ - F(SWEEP, "SWEEP", m_md) \ - F(T, "T", m_m) \ - F(TRACE, "TRACE", d_m) \ - F(TRANSPOS, "TRANSPOS", m_m) \ - F(TRUNC, "TRUNC", m_m) \ - F(UNIFORM, "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. + + - "n": This gets passed the "const struct matrix_expr *" that represents + the expression. This allows the evaluation function to grab the source + location of arguments so that it can report accurate error locations. + + 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!=0". + + - "ab", "a<=b", "a>=b", "b!=0". +*/ +#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_mn, 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_mn, NULL) \ + F(DET, "DET", d_m, NULL) \ + F(DIAG, "DIAG", m_m, NULL) \ + F(EVAL, "EVAL", m_mn, NULL) \ + F(EXP, "EXP", m_e, NULL) \ + F(GINV, "GINV", m_m, NULL) \ + F(GRADE, "GRADE", m_m, NULL) \ + F(GSCH, "GSCH", m_mn, NULL) \ + F(IDENT, "IDENT", IDENT, "ai>=0 bi>=0") \ + 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, "ai>=0 bi>=0") \ + F(MDIAG, "MDIAG", m_v, NULL) \ + F(MMAX, "MMAX", d_m, NULL) \ + F(MMIN, "MMIN", d_m, NULL) \ + F(MOD, "MOD", m_md, "b!=0") \ + 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_mddn, 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_mmn, NULL) \ + F(SQRT, "SQRT", m_e, "a>=0") \ + F(SSCP, "SSCP", m_m, NULL) \ + F(SVAL, "SVAL", m_m, NULL) \ + F(SWEEP, "SWEEP", m_mdn, 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_ddn, "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 { IDENT_MIN_ARGS = 1, IDENT_MAX_ARGS = 2 }; +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 { d_m_MIN_ARGS = 1, d_m_MAX_ARGS = 1 }; +enum { d_none_MIN_ARGS = 0, d_none_MAX_ARGS = 0 }; +enum { m_any_MIN_ARGS = 1, m_any_MAX_ARGS = INT_MAX }; 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_ddn_MIN_ARGS = 2, m_ddn_MAX_ARGS = 2 }; +enum { m_e_MIN_ARGS = 1, m_e_MAX_ARGS = 1 }; +enum { m_ed_MIN_ARGS = 2, m_ed_MAX_ARGS = 2 }; +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_m_MIN_ARGS = 1, m_m_MAX_ARGS = 1 }; enum { m_md_MIN_ARGS = 2, m_md_MAX_ARGS = 2 }; -enum { m_mdd_MIN_ARGS = 3, m_mdd_MAX_ARGS = 3 }; +enum { m_mddn_MIN_ARGS = 3, m_mddn_MAX_ARGS = 3 }; +enum { m_mdn_MIN_ARGS = 2, m_mdn_MAX_ARGS = 2 }; enum { m_mm_MIN_ARGS = 2, m_mm_MAX_ARGS = 2 }; +enum { m_mmn_MIN_ARGS = 2, m_mmn_MAX_ARGS = 2 }; +enum { m_mn_MIN_ARGS = 1, m_mn_MAX_ARGS = 1 }; 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_ddn (double, double, + const struct matrix_expr *); typedef gsl_matrix *matrix_proto_m_m (gsl_matrix *); +typedef gsl_matrix *matrix_proto_m_mn (gsl_matrix *, + const struct matrix_expr *); +typedef double matrix_proto_m_e (double); typedef gsl_matrix *matrix_proto_m_md (gsl_matrix *, double); -typedef gsl_matrix *matrix_proto_m_mdd (gsl_matrix *, double, double); +typedef gsl_matrix *matrix_proto_m_mdn (gsl_matrix *, double, + const struct matrix_expr *); +typedef double matrix_proto_m_ed (double, double); +typedef gsl_matrix *matrix_proto_m_mddn (gsl_matrix *, double, double, + const struct matrix_expr *); +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_mmn (gsl_matrix *, gsl_matrix *, + const struct matrix_expr *); 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(ENUM, STRING, PROTO) static matrix_proto_##PROTO matrix_eval_##ENUM; +#define F(ENUM, STRING, PROTO, CONSTRAINTS) \ + static matrix_proto_##PROTO matrix_eval_##ENUM; MATRIX_FUNCTIONS #undef F @@ -268,7 +452,7 @@ struct matrix_expr enum matrix_op { /* Functions. */ -#define F(ENUM, STRING, PROTO) MOP_F_##ENUM, +#define F(ENUM, STRING, PROTO, CONSTRAINTS) MOP_F_##ENUM, MATRIX_FUNCTIONS #undef F @@ -333,8 +517,69 @@ struct matrix_expr struct matrix_var *variable; struct read_file *eof; }; + + struct msg_location *location; }; +static void +matrix_expr_location__ (const struct matrix_expr *e, + const struct msg_location **minp, + const struct msg_location **maxp) +{ + struct msg_location *loc = e->location; + if (loc) + { + if (loc->p[0].line + && (!minp + || loc->p[0].line < (*minp)->p[0].line + || (loc->p[0].line == (*minp)->p[0].line + && loc->p[0].column < (*minp)->p[0].column))) + *minp = loc; + + if (loc->p[1].line + && (!maxp + || loc->p[1].line > (*maxp)->p[1].line + || (loc->p[1].line == (*maxp)->p[1].line + && loc->p[1].column > (*maxp)->p[1].column))) + *maxp = loc; + return; + } + + assert (e->op != MOP_NUMBER && e->op != MOP_VARIABLE && e->op != MOP_EOF); + for (size_t i = 0; i < e->n_subs; i++) + matrix_expr_location__ (e->subs[i], minp, maxp); +} + +static struct msg_location * +matrix_expr_location (const struct matrix_expr *e_) +{ + struct matrix_expr *e = CONST_CAST (struct matrix_expr *, e_); + + if (!e->location) + { + const struct msg_location *min = NULL; + const struct msg_location *max = NULL; + matrix_expr_location__ (e, &min, &max); + if (min && max) + { + e->location = msg_location_dup (min); + e->location->p[1] = max->p[1]; + } + } + + return e->location; +} + +static struct matrix_expr * +matrix_expr_wrap_location (struct matrix_state *s, int start_ofs, + struct matrix_expr *e) +{ + if (e && !e->location) + e->location = lex_ofs_location (s->lexer, start_ofs, + lex_ofs (s->lexer) - 1); + return e; +} + static void matrix_expr_destroy (struct matrix_expr *e) { @@ -343,7 +588,7 @@ matrix_expr_destroy (struct matrix_expr *e) switch (e->op) { -#define F(ENUM, STRING, PROTO) case MOP_F_##ENUM: +#define F(ENUM, STRING, PROTO, CONSTRAINTS) case MOP_F_##ENUM: MATRIX_FUNCTIONS #undef F case MOP_NEGATE: @@ -384,6 +629,7 @@ MATRIX_FUNCTIONS case MOP_EOF: break; } + msg_location_destroy (e->location); free (e); } @@ -397,6 +643,7 @@ matrix_expr_create_subs (enum matrix_op op, struct matrix_expr **subs, .subs = xmemdup (subs, n_subs * sizeof *subs), .n_subs = n_subs }; + return e; } @@ -430,25 +677,29 @@ matrix_expr_create_3 (enum matrix_op op, struct matrix_expr *sub0, } static struct matrix_expr * -matrix_expr_create_number (double number) +matrix_expr_create_number (struct lexer *lexer, double number) { struct matrix_expr *e = xmalloc (sizeof *e); - *e = (struct matrix_expr) { .op = MOP_NUMBER, .number = number }; + *e = (struct matrix_expr) { + .op = MOP_NUMBER, + .number = number, + }; + lex_get (lexer); return e; } -static struct matrix_expr *matrix_parse_or_xor (struct matrix_state *); +static struct matrix_expr *matrix_parse_expr (struct matrix_state *); static struct matrix_expr * matrix_parse_curly_comma (struct matrix_state *s) { - struct matrix_expr *lhs = matrix_parse_or_xor (s); + struct matrix_expr *lhs = matrix_parse_expr (s); if (!lhs) return NULL; while (lex_match (s->lexer, T_COMMA)) { - struct matrix_expr *rhs = matrix_parse_or_xor (s); + struct matrix_expr *rhs = matrix_parse_expr (s); if (!rhs) { matrix_expr_destroy (lhs); @@ -502,12 +753,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 +777,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 * @@ -568,7 +813,7 @@ matrix_eval_BLOCK (gsl_matrix *m[], size_t n) } static gsl_matrix * -matrix_eval_CHOL (gsl_matrix *m) +matrix_eval_CHOL (gsl_matrix *m, const struct matrix_expr *e) { if (!gsl_linalg_cholesky_decomp1 (m)) { @@ -583,7 +828,8 @@ matrix_eval_CHOL (gsl_matrix *m) } else { - msg (SE, _("Input to CHOL function is not positive-definite.")); + msg_at (SE, e->subs[0]->location, + _("Input to CHOL function is not positive-definite.")); return NULL; } } @@ -623,12 +869,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 * @@ -674,7 +918,7 @@ compare_double_3way (const void *a_, const void *b_) } static gsl_matrix * -matrix_eval_DESIGN (gsl_matrix *m) +matrix_eval_DESIGN (gsl_matrix *m, const struct matrix_expr *e) { double *tmp = xmalloc (m->size1 * m->size2 * sizeof *tmp); gsl_matrix m2 = gsl_matrix_view_array (tmp, m->size2, m->size1).matrix; @@ -699,7 +943,8 @@ matrix_eval_DESIGN (gsl_matrix *m) } if (n[i] <= 1) - msg (MW, _("Column %zu in DESIGN argument has constant value."), i + 1); + msg_at (MW, e->subs[0]->location, + _("Column %zu in DESIGN argument has constant value."), i + 1); else n_total += n[i]; } @@ -768,11 +1013,12 @@ compare_double_desc (const void *a_, const void *b_) } static gsl_matrix * -matrix_eval_EVAL (gsl_matrix *m) +matrix_eval_EVAL (gsl_matrix *m, const struct matrix_expr *e) { if (!is_symmetric (m)) { - msg (SE, _("Argument of EVAL must be symmetric.")); + msg_at (SE, e->subs[0]->location, + _("Argument of EVAL must be symmetric.")); return NULL; } @@ -788,12 +1034,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 @@ -935,13 +1179,14 @@ norm (gsl_vector *v) } static gsl_matrix * -matrix_eval_GSCH (gsl_matrix *v) +matrix_eval_GSCH (gsl_matrix *v, const struct matrix_expr *e) { if (v->size2 < v->size1) { - msg (SE, _("GSCH requires its argument to have at least as many columns " - "as rows, but it has dimensions %zu×%zu."), - v->size1, v->size2); + msg_at (SE, e->subs[0]->location, + _("GSCH requires its argument to have at least as many columns " + "as rows, but it has dimensions %zu×%zu."), + v->size1, v->size2); return NULL; } if (!v->size1 || !v->size2) @@ -975,9 +1220,10 @@ matrix_eval_GSCH (gsl_matrix *v) if (ux < v->size1) { - msg (SE, _("%zu×%zu argument to GSCH contains only " - "%zu linearly independent columns."), - v->size1, v->size2, ux); + msg_at (SE, e->subs[0]->location, + _("%zu×%zu argument to GSCH contains only " + "%zu linearly independent columns."), + v->size1, v->size2, ux); gsl_matrix_free (u); return NULL; } @@ -989,12 +1235,6 @@ matrix_eval_GSCH (gsl_matrix *v) static gsl_matrix * matrix_eval_IDENT (double s1, double s2) { - if (s1 < 0 || s1 > SIZE_MAX || s2 < 0 || s2 > SIZE_MAX) - { - msg (SE, _("Arguments to IDENT must be non-negative integers.")); - return NULL; - } - gsl_matrix *m = gsl_matrix_alloc (s1, s2); MATRIX_FOR_ALL_ELEMENTS (d, y, x, m) *d = x == y; @@ -1039,20 +1279,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 +1440,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); @@ -1224,12 +1455,6 @@ matrix_eval_MAGIC (double n_) static gsl_matrix * matrix_eval_MAKE (double r, double c, double value) { - if (r < 0 || r >= SIZE_MAX || c < 0 || c >= SIZE_MAX) - { - msg (SE, _("Size arguments to MAKE must be integers.")); - return NULL; - } - gsl_matrix *m = gsl_matrix_alloc (r, c); MATRIX_FOR_ALL_ELEMENTS (d, y, x, m) *d = value; @@ -1260,12 +1485,6 @@ matrix_eval_MMIN (gsl_matrix *m) static gsl_matrix * matrix_eval_MOD (gsl_matrix *m, double divisor) { - if (divisor == 0.0) - { - msg (SE, _("Divisor argument to MOD function must be nonzero.")); - return NULL; - } - MATRIX_FOR_ALL_ELEMENTS (d, y, x, m) *d = fmod (*d, divisor); return m; @@ -1312,19 +1531,30 @@ matrix_eval_RANK (gsl_matrix *m) } static gsl_matrix * -matrix_eval_RESHAPE (gsl_matrix *m, double r_, double c_) +matrix_eval_RESHAPE (gsl_matrix *m, double r_, double c_, + const struct matrix_expr *e) { - if (r_ < 0 || r_ >= SIZE_MAX || c_ < 0 || c_ >= SIZE_MAX) + bool r_ok = r_ >= 0 && r_ < SIZE_MAX; + bool c_ok = c_ >= 0 && c_ < SIZE_MAX; + if (!r_ok || !c_ok) { - msg (SE, _("Arguments 2 and 3 to RESHAPE must be integers.")); + msg_at (SE, + !r_ok ? e->subs[1]->location : e->subs[2]->location, + _("Arguments 2 and 3 to RESHAPE must be integers.")); return NULL; } size_t r = r_; size_t c = c_; if (size_overflow_p (xtimes (r, xmax (c, 1))) || c * r != m->size1 * m->size2) { - msg (SE, _("Product of RESHAPE arguments 2 and 3 differ from " - "product of matrix dimensions.")); + struct msg_location *loc = msg_location_dup (e->subs[1]->location); + loc->p[1] = e->subs[2]->location->p[1]; + msg_at (SE, loc, _("Product of RESHAPE size arguments (%zu×%zu = %zu) " + "differs from product of matrix dimensions " + "(%zu×%zu = %zu)."), + r, c, r * c, + m->size1, m->size2, m->size1 * m->size2); + msg_location_destroy (loc); return NULL; } @@ -1378,12 +1608,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,24 +1692,37 @@ 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 * -matrix_eval_SOLVE (gsl_matrix *m1, gsl_matrix *m2) +matrix_eval_SOLVE (gsl_matrix *m1, gsl_matrix *m2, const struct matrix_expr *e) { if (m1->size1 != m2->size1) { - msg (SE, _("SOLVE requires its arguments to have the same number of " - "rows, but the first argument has dimensions %zu×%zu and " - "the second %zu×%zu."), - m1->size1, m1->size2, - m2->size1, m2->size2); + struct msg_location *loc = msg_location_dup (e->subs[0]->location); + loc->p[1] = e->subs[1]->location->p[1]; + +#if 0 + msg_at (SE, loc, + _("SOLVE requires its arguments to have the same number of " + "rows, but the first argument has dimensions %zu×%zu and " + "the second %zu×%zu."), + m1->size1, m1->size2, + m2->size1, m2->size2); +#else + msg_at (SE, e->location, _("SOLVE requires its arguments to have the same " + "number of rows.")); + msg_at (SN, e->subs[0]->location, _("Argument 1 has dimensions %zu×%zu."), + m1->size1, m1->size2); + msg_at (SN, e->subs[1]->location, _("Argument 2 has dimensions %zu×%zu."), + m2->size1, m2->size2); +#endif + + msg_location_destroy (loc); return NULL; } @@ -1499,19 +1740,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 * @@ -1552,19 +1784,21 @@ matrix_eval_SVAL (gsl_matrix *m) } static gsl_matrix * -matrix_eval_SWEEP (gsl_matrix *m, double d) +matrix_eval_SWEEP (gsl_matrix *m, double d, const struct matrix_expr *e) { if (d < 1 || d > SIZE_MAX) { - msg (SE, _("Scalar argument to SWEEP must be integer.")); + msg_at (SE, e->subs[1]->location, + _("Scalar argument to SWEEP must be integer.")); return NULL; } size_t k = d - 1; if (k >= MIN (m->size1, m->size2)) { - msg (SE, _("Scalar argument to SWEEP must be integer less than or " - "equal to the smaller of the matrix argument's rows and " - "columns.")); + msg_at (SE, e->subs[1]->location, + _("Scalar argument to SWEEP must be integer less than or " + "equal to the smaller of the matrix argument's rows and " + "columns.")); return NULL; } @@ -1627,27 +1861,26 @@ 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_) +matrix_eval_UNIFORM (double r_, double c_, const struct matrix_expr *e) { - 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)))) { - msg (SE, _("Product of arguments to UNIFORM exceeds memory size.")); + struct msg_location *loc = msg_location_dup (e->subs[0]->location); + loc->p[1] = e->subs[1]->location->p[1]; + + msg_at (SE, loc, + _("Product of arguments to UNIFORM exceeds memory size.")); + + msg_location_destroy (loc); return NULL; } @@ -1657,153 +1890,860 @@ 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 double +matrix_eval_PDF_BETA (double x, double a, double b) +{ + return gsl_ran_beta_pdf (x, a, b); +} -static struct matrix_expr *matrix_parse_expr (struct matrix_state *); +static double +matrix_eval_CDF_BETA (double x, double a, double b) +{ + return gsl_cdf_beta_P (x, a, b); +} -static const struct matrix_function * -matrix_parse_function_name (const char *token) +static double +matrix_eval_IDF_BETA (double P, double a, double b) { - static const struct matrix_function functions[] = - { -#define F(ENUM, STRING, PROTO) \ - { #ENUM, MOP_F_##ENUM, PROTO##_MIN_ARGS, PROTO##_MAX_ARGS }, - MATRIX_FUNCTIONS -#undef F - }; - enum { N_FUNCTIONS = sizeof functions / sizeof *functions }; + return gsl_cdf_beta_Pinv (P, 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_RV_BETA (double a, double b) +{ + return gsl_ran_beta (get_rng (), a, b); } -static struct read_file * -read_file_create (struct matrix_state *s, struct file_handle *fh) +static double +matrix_eval_NCDF_BETA (double x, double a, double b, double lambda) { - 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 ncdf_beta (x, a, b, lambda); +} - struct read_file *rf = xmalloc (sizeof *rf); - *rf = (struct read_file) { .file = fh }; +static double +matrix_eval_NPDF_BETA (double x, double a, double b, double lambda) +{ + return npdf_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_CDF_BVNOR (double x0, double x1, double r) +{ + return cdf_bvnor (x0, x1, r); } -static struct dfm_reader * -read_file_open (struct read_file *rf) +static double +matrix_eval_PDF_BVNOR (double x0, double x1, double r) { - if (!rf->reader) - rf->reader = dfm_open_reader (rf->file, NULL, rf->encoding); - return rf->reader; + return gsl_ran_bivariate_gaussian_pdf (x0, x1, 1, 1, r); } -static void -read_file_destroy (struct read_file *rf) +static double +matrix_eval_CDF_CAUCHY (double x, double a, double b) { - if (rf) - { - fh_unref (rf->file); - dfm_close_reader (rf->reader); - free (rf->encoding); - free (rf); - } + return gsl_cdf_cauchy_P ((x - a) / b, 1); } -static struct write_file * -write_file_create (struct matrix_state *s, struct file_handle *fh) +static double +matrix_eval_IDF_CAUCHY (double P, 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 a + b * gsl_cdf_cauchy_Pinv (P, 1); +} - struct write_file *wf = xmalloc (sizeof *wf); - *wf = (struct write_file) { .file = fh }; +static double +matrix_eval_PDF_CAUCHY (double x, double a, double b) +{ + return gsl_ran_cauchy_pdf ((x - a) / b, 1) / b; +} - 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_RV_CAUCHY (double a, double b) +{ + return a + b * gsl_ran_cauchy (get_rng (), 1); } -static struct dfm_writer * -write_file_open (struct write_file *wf) +static double +matrix_eval_CDF_CHISQ (double x, double df) { - if (!wf->writer) - wf->writer = dfm_open_writer (wf->file, wf->encoding); - return wf->writer; + return gsl_cdf_chisq_P (x, df); } -static void -write_file_destroy (struct write_file *wf) +static double +matrix_eval_CHICDF (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 matrix_eval_CDF_CHISQ (x, df); +} - fh_unref (wf->file); - dfm_close_writer (wf->writer); - free (wf->encoding); - free (wf); - } +static double +matrix_eval_IDF_CHISQ (double P, double df) +{ + return gsl_cdf_chisq_Pinv (P, df); } -static bool -matrix_parse_function (struct matrix_state *s, const char *token, - struct matrix_expr **exprp) +static double +matrix_eval_PDF_CHISQ (double x, double df) { - *exprp = NULL; - if (lex_next_token (s->lexer, 1) != T_LPAREN) - return false; + return gsl_ran_chisq_pdf (x, 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_RV_CHISQ (double df) +{ + return gsl_ran_chisq (get_rng (), df); +} - if (!lex_force_match (s->lexer, T_RPAREN)) - { - fh_unref (fh); - return true; - } +static double +matrix_eval_SIG_CHISQ (double x, double df) +{ + return gsl_cdf_chisq_Q (x, df); +} - struct read_file *rf = read_file_create (s, fh); +static double +matrix_eval_CDF_EXP (double x, double a) +{ + return gsl_cdf_exponential_P (x, 1. / a); +} - struct matrix_expr *e = xmalloc (sizeof *e); - *e = (struct matrix_expr) { .op = MOP_EOF, .eof = rf }; - *exprp = e; +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; } @@ -1811,26 +2751,56 @@ matrix_parse_function (struct matrix_state *s, const char *token, if (!f) return false; - lex_get_n (s->lexer, 2); - struct matrix_expr *e = xmalloc (sizeof *e); - *e = (struct matrix_expr) { .op = f->op, .subs = NULL }; + *e = (struct matrix_expr) { .op = f->op }; - size_t allocated_subs = 0; - do + lex_get_n (s->lexer, 2); + if (lex_token (s->lexer) != T_RPAREN) { - struct matrix_expr *sub = matrix_parse_expr (s); - if (!sub) - goto error; + 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; + 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)); } - 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); + *exprp = e; return true; @@ -1840,27 +2810,22 @@ error: } static struct matrix_expr * -matrix_parse_primary (struct matrix_state *s) +matrix_parse_primary__ (struct matrix_state *s) { if (lex_is_number (s->lexer)) - { - double number = lex_number (s->lexer); - lex_get (s->lexer); - return matrix_expr_create_number (number); - } + return matrix_expr_create_number (s->lexer, lex_number (s->lexer)); else if (lex_is_string (s->lexer)) { char string[sizeof (double)]; buf_copy_str_rpad (string, sizeof string, lex_tokcstr (s->lexer), ' '); - lex_get (s->lexer); double number; memcpy (&number, string, sizeof number); - return matrix_expr_create_number (number); + return matrix_expr_create_number (s->lexer, number); } else if (lex_match (s->lexer, T_LPAREN)) { - struct matrix_expr *e = matrix_parse_or_xor (s); + struct matrix_expr *e = matrix_parse_expr (s); if (!e || !lex_force_match (s->lexer, T_RPAREN)) { matrix_expr_destroy (e); @@ -1908,6 +2873,13 @@ matrix_parse_primary (struct matrix_state *s) return NULL; } +static struct matrix_expr * +matrix_parse_primary (struct matrix_state *s) +{ + int start_ofs = lex_ofs (s->lexer); + return matrix_expr_wrap_location (s, start_ofs, matrix_parse_primary__ (s)); +} + static struct matrix_expr *matrix_parse_postfix (struct matrix_state *); static bool @@ -1968,15 +2940,28 @@ matrix_parse_postfix (struct matrix_state *s) static struct matrix_expr * matrix_parse_unary (struct matrix_state *s) { + int start_ofs = lex_ofs (s->lexer); + + struct matrix_expr *e; if (lex_match (s->lexer, T_DASH)) { - struct matrix_expr *lhs = matrix_parse_unary (s); - return lhs ? matrix_expr_create_1 (MOP_NEGATE, lhs) : NULL; + struct matrix_expr *sub = matrix_parse_unary (s); + if (!sub) + return NULL; + e = matrix_expr_create_1 (MOP_NEGATE, sub); } else if (lex_match (s->lexer, T_PLUS)) - return matrix_parse_unary (s); + { + e = matrix_parse_unary (s); + if (!e) + return NULL; + } else return matrix_parse_postfix (s); + + matrix_expr_wrap_location (s, start_ofs, e); + e->location->p[0] = lex_ofs_start_point (s->lexer, start_ofs); + return e; } static struct matrix_expr * @@ -2133,10 +3118,16 @@ matrix_parse_relational (struct matrix_state *s) static struct matrix_expr * matrix_parse_not (struct matrix_state *s) { + int start_ofs = lex_ofs (s->lexer); if (lex_match (s->lexer, T_NOT)) { - struct matrix_expr *lhs = matrix_parse_not (s); - return lhs ? matrix_expr_create_1 (MOP_NOT, lhs) : NULL; + struct matrix_expr *sub = matrix_parse_not (s); + if (!sub) + return NULL; + + matrix_expr_wrap_location (s, start_ofs, sub); + sub->location->p[0] = lex_ofs_start_point (s->lexer, start_ofs); + return sub; } else return matrix_parse_relational (s); @@ -2163,7 +3154,7 @@ matrix_parse_and (struct matrix_state *s) } static struct matrix_expr * -matrix_parse_or_xor (struct matrix_state *s) +matrix_parse_expr__ (struct matrix_state *s) { struct matrix_expr *lhs = matrix_parse_and (s); if (!lhs) @@ -2177,7 +3168,9 @@ matrix_parse_or_xor (struct matrix_state *s) else if (lex_match_id (s->lexer, "XOR")) op = MOP_XOR; else - return lhs; + { + return lhs; + } struct matrix_expr *rhs = matrix_parse_and (s); if (!rhs) @@ -2192,7 +3185,8 @@ matrix_parse_or_xor (struct matrix_state *s) static struct matrix_expr * matrix_parse_expr (struct matrix_state *s) { - return matrix_parse_or_xor (s); + int start_ofs = lex_ofs (s->lexer); + return matrix_expr_wrap_location (s, start_ofs, matrix_parse_expr__ (s));; } /* Expression evaluation. */ @@ -2217,7 +3211,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(ENUM, STRING, PROTO) case MOP_F_##ENUM: +#define F(ENUM, STRING, PROTO, CONSTRAINTS) case MOP_F_##ENUM: MATRIX_FUNCTIONS #undef F case MOP_NEGATE: @@ -2262,7 +3256,7 @@ matrix_op_name (enum matrix_op op) case MOP_OR: return "OR"; case MOP_XOR: return "XOR"; -#define F(ENUM, STRING, PROTO) case MOP_F_##ENUM: +#define F(ENUM, STRING, PROTO, CONSTRAINTS) case MOP_F_##ENUM: MATRIX_FUNCTIONS #undef F case MOP_NEGATE: @@ -2301,7 +3295,8 @@ to_scalar (const gsl_matrix *m) } static gsl_matrix * -matrix_expr_evaluate_elementwise (enum matrix_op op, +matrix_expr_evaluate_elementwise (const struct matrix_expr *e, + enum matrix_op op, gsl_matrix *a, gsl_matrix *b) { if (is_scalar (b)) @@ -2339,24 +3334,27 @@ matrix_expr_evaluate_elementwise (enum matrix_op op, } else { - msg (SE, _("Operands to %s must have the same dimensions or one " - "must be a scalar, not %zu×%zu and %zu×%zu matrices."), + msg_at (SE, e->location, + _("Operands to %s must have the same dimensions or one " + "must be a scalar, not %zu×%zu and %zu×%zu matrices."), matrix_op_name (op), a->size1, a->size2, b->size1, b->size2); return NULL; } } static gsl_matrix * -matrix_expr_evaluate_mul_mat (gsl_matrix *a, gsl_matrix *b) +matrix_expr_evaluate_mul_mat (const struct matrix_expr *e, + gsl_matrix *a, gsl_matrix *b) { if (is_scalar (a) || is_scalar (b)) - return matrix_expr_evaluate_elementwise (MOP_MUL_ELEMS, a, b); + return matrix_expr_evaluate_elementwise (e, MOP_MUL_ELEMS, a, b); if (a->size2 != b->size1) { - msg (SE, _("Matrices with dimensions %zu×%zu and %zu×%zu are " - "not conformable for multiplication."), - a->size1, a->size2, b->size1, b->size2); + msg_at (SE, e->location, + _("Matrices with dimensions %zu×%zu and %zu×%zu are " + "not conformable for multiplication."), + a->size1, a->size2, b->size1, b->size2); return NULL; } @@ -2388,28 +3386,32 @@ square_matrix (gsl_matrix **x, gsl_matrix **tmp) } static gsl_matrix * -matrix_expr_evaluate_exp_mat (gsl_matrix *x_, gsl_matrix *b) +matrix_expr_evaluate_exp_mat (const struct matrix_expr *e, + gsl_matrix *x_, gsl_matrix *b) { gsl_matrix *x = x_; if (x->size1 != x->size2) { - msg (SE, _("Matrix exponentation with ** requires a square matrix on " - "the left-hand size, not one with dimensions %zu×%zu."), - x->size1, x->size2); + msg_at (SE, matrix_expr_location (e->subs[0]), + _("Matrix exponentation with ** requires a square matrix on " + "the left-hand size, not one with dimensions %zu×%zu."), + x->size1, x->size2); return NULL; } if (!is_scalar (b)) { - msg (SE, _("Matrix exponentiation with ** requires a scalar on the " - "right-hand side, not a matrix with dimensions %zu×%zu."), - b->size1, b->size2); + msg_at (SE, matrix_expr_location (e->subs[1]), + _("Matrix exponentiation with ** requires a scalar on the " + "right-hand side, not a matrix with dimensions %zu×%zu."), + b->size1, b->size2); return NULL; } double bf = to_scalar (b); if (bf != floor (bf) || bf <= LONG_MIN || bf > LONG_MAX) { - msg (SE, _("Exponent %.1f in matrix multiplication is non-integer " - "or outside the valid range."), bf); + msg_at (SE, matrix_expr_location (e->subs[1]), + _("Exponent %.1f in matrix multiplication is non-integer " + "or outside the valid range."), bf); return NULL; } long int bl = bf; @@ -2448,13 +3450,28 @@ matrix_expr_evaluate_exp_mat (gsl_matrix *x_, gsl_matrix *b) return y; } +static void +note_nonscalar (gsl_matrix *m, const struct matrix_expr *e) +{ + if (!is_scalar (m)) + msg_at (SN, matrix_expr_location (e), + _("This operand is a %zu×%zu matrix."), m->size1, m->size2); +} + static gsl_matrix * -matrix_expr_evaluate_seq (gsl_matrix *start_, gsl_matrix *end_, +matrix_expr_evaluate_seq (const struct matrix_expr *e, + gsl_matrix *start_, gsl_matrix *end_, gsl_matrix *by_) { if (!is_scalar (start_) || !is_scalar (end_) || (by_ && !is_scalar (by_))) { - msg (SE, _("All operands of : operator must be scalars.")); + msg_at (SE, matrix_expr_location (e), + _("All operands of : operator must be scalars.")); + + note_nonscalar (start_, e->subs[0]); + note_nonscalar (end_, e->subs[1]); + if (by_) + note_nonscalar (by_, e->subs[2]); return NULL; } @@ -2720,10 +3737,10 @@ matrix_expr_evaluate_mat_index (gsl_matrix *sm, gsl_matrix *im0, return dm; } -#define F(ENUM, STRING, PROTO) \ - static gsl_matrix *matrix_expr_evaluate_##PROTO ( \ - const char *name, gsl_matrix *[], size_t, \ - matrix_proto_##PROTO *); +#define F(ENUM, STRING, PROTO, CONSTRAINTS) \ + static gsl_matrix *matrix_expr_evaluate_##PROTO ( \ + const struct matrix_function_properties *, gsl_matrix *[], \ + const struct matrix_expr *, matrix_proto_##PROTO *); MATRIX_FUNCTIONS #undef F @@ -2765,97 +3782,611 @@ 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)); +} + +enum matrix_argument_relop + { + MRR_GT, /* > */ + MRR_GE, /* >= */ + MRR_LT, /* < */ + MRR_LE, /* <= */ + MRR_NE, /* <> */ + }; + +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, + enum matrix_argument_relop relop) +{ + struct string s = DS_EMPTY_INITIALIZER; + + argument_invalid (props, a_index, a, y, x, &s); + ds_put_cstr (&s, " "); + switch (relop) + { + case MRR_GE: + ds_put_format (&s, _("This argument must be greater than or " + "equal to argument %zu, but that has value %g."), + b_index, b); + break; + + case MRR_GT: + ds_put_format (&s, _("This argument must be greater than argument %zu, " + "but that has value %g."), + b_index, b); + break; + + case MRR_LE: + ds_put_format (&s, _("This argument must be less than or " + "equal to argument %zu, but that has value %g."), + b_index, b); + break; + + case MRR_LT: + ds_put_format (&s, _("This argument must be less than argument %zu, " + "but that has value %g."), + b_index, b); + break; + + case MRR_NE: + ds_put_format (&s, + _("This argument must not be the same as argument %zu."), + b_index); + break; + } + + 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, + enum matrix_argument_relop relop) +{ + struct string s = DS_EMPTY_INITIALIZER; + argument_invalid (props, a_index, a, y, x, &s); + ds_put_cstr (&s, " "); + switch (relop) + { + case MRR_GE: + ds_put_format (&s, + _("This argument must be greater than or equal to %g."), + b); + break; + + case MRR_GT: + ds_put_format (&s, _("This argument must be greater than %g."), b); + break; + + case MRR_LE: + ds_put_format (&s, _("This argument must be less than or equal to %g."), + b); + break; + + case MRR_LT: + ds_put_format (&s, _("This argument must be less than %g."), b); + break; + + case MRR_NE: + ds_put_format (&s, _("This argument may not be %g."), b); + break; + } + + msg (ME, ds_cstr (&s)); + ds_destroy (&s); +} + +static bool +matrix_argument_relop_is_satisfied (double a, double b, + enum matrix_argument_relop relop) +{ + switch (relop) + { + case MRR_GE: return a >= b; + case MRR_GT: return a > b; + case MRR_LE: return a <= b; + case MRR_LT: return a < b; + case MRR_NE: return a != b; + } + + NOT_REACHED (); +} + +static enum matrix_argument_relop +matrix_argument_relop_flip (enum matrix_argument_relop relop) +{ + switch (relop) + { + case MRR_GE: return MRR_LE; + case MRR_GT: return MRR_LT; + case MRR_LE: return MRR_GE; + case MRR_LT: return MRR_GT; + case MRR_NE: return MRR_NE; + } + + NOT_REACHED (); +} + +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 == '<' + || *constraints == '!') + { + enum matrix_argument_relop relop; + switch (*constraints++) + { + case '>': + if (*constraints == '=') + { + constraints++; + relop = MRR_GE; + } + else + relop = MRR_GT; + break; + + case '<': + if (*constraints == '=') + { + constraints++; + relop = MRR_LE; + } + else + relop = MRR_LT; + break; + + case '!': + assert (*constraints == '='); + constraints++; + relop = MRR_NE; + } + + 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; + relop = matrix_argument_relop_flip (relop); + } + + double b = to_scalar (args[b_index]); + MATRIX_FOR_ALL_ELEMENTS (a, y, x, args[a_index]) + if (!matrix_argument_relop_is_satisfied (*a, b, relop)) + { + argument_inequality_error (props, + a_index + 1, args[a_index], y, x, + b_index + 1, b, + relop); + return false; + } + } + else + { + int comparand = parse_constraint_value (&constraints); + + MATRIX_FOR_ALL_ELEMENTS (d, y, x, args[arg_index]) + if (!matrix_argument_relop_is_satisfied (*d, comparand, relop)) + { + argument_value_error (props, + arg_index + 1, args[arg_index], y, x, + comparand, + relop); + return false; + } + } + } + else + { + assert (*constraints == ' '); + constraints++; + arg_index = SIZE_MAX; + } + } + return true; +} + static gsl_matrix * -matrix_expr_evaluate_m_d (const char *name, - gsl_matrix *subs[], size_t n_subs, - matrix_proto_m_d *f) +matrix_expr_evaluate_d_none ( + const struct matrix_function_properties *props UNUSED, + gsl_matrix *subs[] UNUSED, const struct matrix_expr *e, + matrix_proto_d_none *f) { - assert (n_subs == 1); + assert (e->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[], const struct matrix_expr *e, + matrix_proto_d_d *f) +{ + assert (e->n_subs == 1); double d; - return to_scalar_args (name, subs, n_subs, &d) ? f(d) : NULL; + if (!to_scalar_args (props->name, subs, e->n_subs, &d)) + return NULL; + + if (!check_constraints (props, subs, e->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_m_dd (const char *name, - gsl_matrix *subs[], size_t n_subs, - matrix_proto_m_dd *f) +matrix_expr_evaluate_d_dd (const struct matrix_function_properties *props, + gsl_matrix *subs[], const struct matrix_expr *e, + matrix_proto_d_dd *f) { - assert (n_subs == 2); + assert (e->n_subs == 2); double d[2]; - return to_scalar_args (name, subs, n_subs, d) ? f(d[0], d[1]) : NULL; + if (!to_scalar_args (props->name, subs, e->n_subs, d)) + return NULL; + + if (!check_constraints (props, subs, e->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[], const struct matrix_expr *e, + matrix_proto_d_ddd *f) +{ + assert (e->n_subs == 3); + + double d[3]; + if (!to_scalar_args (props->name, subs, e->n_subs, d)) + return NULL; + + if (!check_constraints (props, subs, e->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 struct matrix_function_properties *props, + gsl_matrix *subs[], const struct matrix_expr *e, + matrix_proto_m_d *f) +{ + assert (e->n_subs == 1); + + double d; + return to_scalar_args (props->name, subs, e->n_subs, &d) ? f (d) : NULL; } static gsl_matrix * -matrix_expr_evaluate_m_ddd (const char *name, - gsl_matrix *subs[], size_t n_subs, +matrix_expr_evaluate_m_ddd (const struct matrix_function_properties *props, + gsl_matrix *subs[], const struct matrix_expr *e, matrix_proto_m_ddd *f) { - assert (n_subs == 3); + assert (e->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, e->n_subs, d) ? f(d[0], d[1], d[2]) : NULL; +} + +static gsl_matrix * +matrix_expr_evaluate_m_ddn (const struct matrix_function_properties *props, + gsl_matrix *subs[], const struct matrix_expr *e, + matrix_proto_m_ddn *f) +{ + assert (e->n_subs == 2); + + double d[2]; + return (to_scalar_args (props->name, subs, e->n_subs, d) + ? f(d[0], d[1], e) + : NULL); } static gsl_matrix * -matrix_expr_evaluate_m_m (const char *name UNUSED, - gsl_matrix *subs[], size_t n_subs, +matrix_expr_evaluate_m_m (const struct matrix_function_properties *props UNUSED, + gsl_matrix *subs[], const struct matrix_expr *e, matrix_proto_m_m *f) { - assert (n_subs == 1); + assert (e->n_subs == 1); return f (subs[0]); } static gsl_matrix * -matrix_expr_evaluate_m_md (const char *name UNUSED, - gsl_matrix *subs[], size_t n_subs, +matrix_expr_evaluate_m_mn (const struct matrix_function_properties *props UNUSED, + gsl_matrix *subs[], const struct matrix_expr *e, + matrix_proto_m_mn *f) +{ + assert (e->n_subs == 1); + return f (subs[0], e); +} + +static gsl_matrix * +matrix_expr_evaluate_m_e (const struct matrix_function_properties *props, + gsl_matrix *subs[], const struct matrix_expr *e, + matrix_proto_m_e *f) +{ + assert (e->n_subs == 1); + + if (!check_constraints (props, subs, e->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[], const struct matrix_expr *e, matrix_proto_m_md *f) { - assert (n_subs == 2); - if (!check_scalar_arg (name, subs, 1)) + assert (e->n_subs == 2); + 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, - gsl_matrix *subs[], size_t n_subs, - matrix_proto_m_mdd *f) +matrix_expr_evaluate_m_mdn (const struct matrix_function_properties *props UNUSED, + gsl_matrix *subs[], const struct matrix_expr *e, + matrix_proto_m_mdn *f) +{ + assert (e->n_subs == 2); + if (!check_scalar_arg (props->name, subs, 1)) + return NULL; + return f (subs[0], to_scalar (subs[1]), e); +} + +static gsl_matrix * +matrix_expr_evaluate_m_ed (const struct matrix_function_properties *props, + gsl_matrix *subs[], const struct matrix_expr *e, + matrix_proto_m_ed *f) +{ + assert (e->n_subs == 2); + if (!check_scalar_arg (props->name, subs, 1)) + return NULL; + + if (!check_constraints (props, subs, e->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_mddn (const struct matrix_function_properties *props UNUSED, + gsl_matrix *subs[], const struct matrix_expr *e, + matrix_proto_m_mddn *f) +{ + assert (e->n_subs == 3); + 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]), e); +} + +static gsl_matrix * +matrix_expr_evaluate_m_edd (const struct matrix_function_properties *props, + gsl_matrix *subs[], const struct matrix_expr *e, + matrix_proto_m_edd *f) +{ + assert (e->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, e->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[], const struct matrix_expr *e, + matrix_proto_m_eddd *f) +{ + assert (e->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, e->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[], const struct matrix_expr *e, + matrix_proto_m_eed *f) { - assert (n_subs == 3); - if (!check_scalar_arg (name, subs, 1) || !check_scalar_arg (name, subs, 2)) + assert (e->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)) + { + struct msg_location *loc = msg_location_dup (e->subs[0]->location); + loc->p[1] = e->subs[1]->location->p[1]; + + msg_at (ME, loc, + _("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); + + msg_location_destroy (loc); + return NULL; + } + + if (!check_constraints (props, subs, e->n_subs)) return NULL; - return f (subs[0], to_scalar (subs[1]), to_scalar (subs[2])); + + 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 char *name UNUSED, - gsl_matrix *subs[], size_t n_subs, +matrix_expr_evaluate_m_mm (const struct matrix_function_properties *props UNUSED, + gsl_matrix *subs[], const struct matrix_expr *e, matrix_proto_m_mm *f) { - assert (n_subs == 2); + assert (e->n_subs == 2); return f (subs[0], subs[1]); } static gsl_matrix * -matrix_expr_evaluate_m_v (const char *name, - gsl_matrix *subs[], size_t n_subs, +matrix_expr_evaluate_m_mmn (const struct matrix_function_properties *props UNUSED, + gsl_matrix *subs[], const struct matrix_expr *e, + matrix_proto_m_mmn *f) +{ + assert (e->n_subs == 2); + return f (subs[0], subs[1], e); +} + +static gsl_matrix * +matrix_expr_evaluate_m_v (const struct matrix_function_properties *props, + gsl_matrix *subs[], const struct matrix_expr *e, matrix_proto_m_v *f) { - assert (n_subs == 1); - if (!check_vector_arg (name, subs, 0)) + assert (e->n_subs == 1); + 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, - gsl_matrix *subs[], size_t n_subs, +matrix_expr_evaluate_d_m (const struct matrix_function_properties *props UNUSED, + gsl_matrix *subs[], const struct matrix_expr *e, matrix_proto_d_m *f) { - assert (n_subs == 1); + assert (e->n_subs == 1); gsl_matrix *m = gsl_matrix_alloc (1, 1); gsl_matrix_set (m, 0, 0, f (subs[0])); @@ -2863,25 +4394,25 @@ matrix_expr_evaluate_d_m (const char *name UNUSED, } static gsl_matrix * -matrix_expr_evaluate_m_any (const char *name UNUSED, - gsl_matrix *subs[], size_t n_subs, - matrix_proto_m_any *f) +matrix_expr_evaluate_m_any (const struct matrix_function_properties *props UNUSED, + gsl_matrix *subs[], const struct matrix_expr *e, + matrix_proto_m_any *f) { - return f (subs, n_subs); + return f (subs, e->n_subs); } static gsl_matrix * -matrix_expr_evaluate_IDENT (const char *name, - gsl_matrix *subs[], size_t n_subs, +matrix_expr_evaluate_IDENT (const struct matrix_function_properties *props, + gsl_matrix *subs[], const struct matrix_expr *e, matrix_proto_IDENT *f) { - assert (n_subs <= 2); + assert (e->n_subs <= 2); double d[2]; - if (!to_scalar_args (name, subs, n_subs, d)) + if (!to_scalar_args (props->name, subs, e->n_subs, d)) return NULL; - return f (d[0], d[n_subs - 1]); + return f (d[0], d[e->n_subs - 1]); } static gsl_matrix * @@ -2898,8 +4429,9 @@ matrix_expr_evaluate (const struct matrix_expr *e) const gsl_matrix *src = e->variable->value; if (!src) { - msg (SE, _("Uninitialized variable %s used in expression."), - e->variable->name); + msg_at (SE, e->location, + _("Uninitialized variable %s used in expression."), + e->variable->name); return NULL; } @@ -2937,10 +4469,16 @@ matrix_expr_evaluate (const struct matrix_expr *e) gsl_matrix *result = NULL; switch (e->op) { -#define F(ENUM, STRING, PROTO) \ +#define F(ENUM, STRING, PROTO, CONSTRAINTS) \ case MOP_F_##ENUM: \ - result = matrix_expr_evaluate_##PROTO (STRING, subs, e->n_subs, \ - matrix_eval_##ENUM); \ + { \ + static const struct matrix_function_properties props = { \ + .name = STRING, \ + .constraints = CONSTRAINTS, \ + }; \ + result = matrix_expr_evaluate_##PROTO (&props, subs, e, \ + matrix_eval_##ENUM); \ + } \ break; MATRIX_FUNCTIONS #undef F @@ -2964,7 +4502,7 @@ matrix_expr_evaluate (const struct matrix_expr *e) case MOP_AND: case MOP_OR: case MOP_XOR: - result = matrix_expr_evaluate_elementwise (e->op, subs[0], subs[1]); + result = matrix_expr_evaluate_elementwise (e, e->op, subs[0], subs[1]); break; case MOP_NOT: @@ -2972,19 +4510,19 @@ matrix_expr_evaluate (const struct matrix_expr *e) break; case MOP_SEQ: - result = matrix_expr_evaluate_seq (subs[0], subs[1], NULL); + result = matrix_expr_evaluate_seq (e, subs[0], subs[1], NULL); break; case MOP_SEQ_BY: - result = matrix_expr_evaluate_seq (subs[0], subs[1], subs[2]); + result = matrix_expr_evaluate_seq (e, subs[0], subs[1], subs[2]); break; case MOP_MUL_MAT: - result = matrix_expr_evaluate_mul_mat (subs[0], subs[1]); + result = matrix_expr_evaluate_mul_mat (e, subs[0], subs[1]); break; case MOP_EXP_MAT: - result = matrix_expr_evaluate_exp_mat (subs[0], subs[1]); + result = matrix_expr_evaluate_exp_mat (e, subs[0], subs[1]); break; case MOP_PASTE_HORZ: @@ -3078,6 +4616,9 @@ struct matrix_lvalue struct matrix_var *var; struct matrix_expr *indexes[2]; size_t n_indexes; + + struct msg_location *var_location; /* Variable name. */ + struct msg_location *index_location; /* Variable name plus indexing. */ }; static void @@ -3085,6 +4626,8 @@ matrix_lvalue_destroy (struct matrix_lvalue *lvalue) { if (lvalue) { + msg_location_destroy (lvalue->var_location); + msg_location_destroy (lvalue->index_location); for (size_t i = 0; i < lvalue->n_indexes; i++) matrix_expr_destroy (lvalue->indexes[i]); free (lvalue); @@ -3097,14 +4640,15 @@ matrix_lvalue_parse (struct matrix_state *s) struct matrix_lvalue *lvalue = xzalloc (sizeof *lvalue); if (!lex_force_id (s->lexer)) goto error; + int start_ofs = lex_ofs (s->lexer); + 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) { if (!lvalue->var) { msg (SE, _("Undefined variable %s."), lex_tokcstr (s->lexer)); - free (lvalue); - return NULL; + goto error; } lex_get_n (s->lexer, 2); @@ -3121,6 +4665,9 @@ matrix_lvalue_parse (struct matrix_state *s) } if (!lex_force_match (s->lexer, T_RPAREN)) goto error; + + lvalue->index_location = lex_ofs_location (s->lexer, start_ofs, + lex_ofs (s->lexer) - 1); } else { @@ -3257,21 +4804,24 @@ 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->index_location, + _("Cannot index %zu×%zu matrix %s."), + dm->size1, dm->size2, lvalue->var->name); 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->index_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], @@ -4904,11 +6454,9 @@ parse_error (const struct dfm_reader *reader, enum fmt_type format, int line_number = dfm_get_line_number (reader); struct msg_location *location = xmalloc (sizeof *location); *location = (struct msg_location) { - .file_name = xstrdup (dfm_get_file_name (reader)), - .first_line = line_number, - .last_line = line_number + 1, - .first_column = first_column, - .last_column = last_column, + .file_name = intern_new (dfm_get_file_name (reader)), + .p[0] = { .line = line_number, .column = first_column }, + .p[1] = { .line = line_number, .column = last_column }, }; struct msg *m = xmalloc (sizeof *m); *m = (struct msg) { @@ -4954,7 +6502,7 @@ matrix_read_set_field (struct read_command *read, struct dfm_reader *reader, if (error) { int c1 = utf8_count_columns (line_start, p.string - line_start) + 1; - int c2 = c1 + ss_utf8_count_columns (p) - 1; + int c2 = c1 + ss_utf8_count_columns (p); parse_error (reader, read->format, p, y, x, c1, c2, error); } else