help-glpk
[Top][All Lists]

## GLPK: Adding asin() and acos() to GNU MathProg

 From: Hartmut Henkel Subject: GLPK: Adding asin() and acos() to GNU MathProg Date: Sat, 24 Jul 2021 16:39:55 +0200 (CEST)

Hi Andrew,

when using gmpl/glpsol for curve fitting, the Chebyshev polynomial t(n,
x) = cos(n * acos(x)) can not easily be built due to the missing acos()
function. The patch below adds asin() and acos() to gmpl. Here is an
example for use, minimax fitting of exp(x):

# $glpsol --dual --model --tolbnd 1e-12 --toldj 1e-12 --tolpiv 1e-14 cheby.mod param x_max := 1; param tol := 0.0000001; param n := 11; param m := 1000 * n; set cols := (1..n); set rows := (1..m); param v{i in rows} := exp((i - 1) / (m - 1) * x_max); param t{i in rows, j in cols} := cos(j * acos((i - 1) / (m - 1) * x_max)); var x{j in cols}; var d; maximize slack: d; s.t. slack1: d >= 0; s.t. pass1{i in rows}: d + v[i] - sum{j in cols} (x[j] * t[i,j]) <= tol; s.t. pass2{i in rows}: d - v[i] + sum{j in cols} (x[j] * t[i,j]) <= tol; solve; printf "#d\n", d > "result"; printf{j in cols}: "%d %.10g\n", j, x[j].val >> "result"; end; Below is the simple patch, checked with the above code. I would be glad if you could include this into glpk. The documentation i could update only in the english version. Best Regards, Hartmut Hartmut Henkel, Zellhausen, Germany --- glpk-5.0/src/mpl/mpl.h 2020-12-16 10:00:00.000000000 +0100 +++ glpk-5.0-mod/src/mpl/mpl.h 2021-07-24 13:54:51.324270191 +0200 @@ -830,10 +830,18 @@ double fp_sin(MPL *mpl, double x); /* floating-point trigonometric sine */ +#define fp_asin _glp_mpl_fp_asin +double fp_asin(MPL *mpl, double x); +/* floating-point trigonometric arcsine */ + #define fp_cos _glp_mpl_fp_cos double fp_cos(MPL *mpl, double x); /* floating-point trigonometric cosine */ +#define fp_acos _glp_mpl_fp_acos +double fp_acos(MPL *mpl, double x); +/* floating-point trigonometric arccosine */ + #define fp_tan _glp_mpl_fp_tan double fp_tan(MPL *mpl, double x); /* floating-point trigonometric tangent */ @@ -2117,64 +2125,66 @@ #define O_LOG10 329 /* common (decimal) logarithm */ #define O_SQRT 330 /* square root */ #define O_SIN 331 /* trigonometric sine */ -#define O_COS 332 /* trigonometric cosine */ -#define O_TAN 333 /* trigonometric tangent */ -#define O_ATAN 334 /* trigonometric arctangent */ -#define O_ROUND 335 /* round to nearest integer */ -#define O_TRUNC 336 /* truncate to nearest integer */ -#define O_CARD 337 /* cardinality of set */ -#define O_LENGTH 338 /* length of symbolic value */ +#define O_ASIN 332 /* trigonometric arcsine */ +#define O_COS 333 /* trigonometric cosine */ +#define O_ACOS 334 /* trigonometric arccosine */ +#define O_TAN 335 /* trigonometric tangent */ +#define O_ATAN 336 /* trigonometric arctangent */ +#define O_ROUND 337 /* round to nearest integer */ +#define O_TRUNC 338 /* truncate to nearest integer */ +#define O_CARD 339 /* cardinality of set */ +#define O_LENGTH 340 /* length of symbolic value */ /* binary operations -------------------*/ -#define O_ADD 339 /* addition */ -#define O_SUB 340 /* subtraction */ -#define O_LESS 341 /* non-negative subtraction */ -#define O_MUL 342 /* multiplication */ -#define O_DIV 343 /* division */ -#define O_IDIV 344 /* quotient of exact division */ -#define O_MOD 345 /* remainder of exact division */ -#define O_POWER 346 /* exponentiation (raise to power) */ -#define O_ATAN2 347 /* trigonometric arctangent */ -#define O_ROUND2 348 /* round to n fractional digits */ -#define O_TRUNC2 349 /* truncate to n fractional digits */ -#define O_UNIFORM 350 /* pseudo-random in [a, b) */ -#define O_NORMAL 351 /* gaussian random, given mu and sigma */ -#define O_CONCAT 352 /* concatenation */ -#define O_LT 353 /* comparison on 'less than' */ -#define O_LE 354 /* comparison on 'not greater than' */ -#define O_EQ 355 /* comparison on 'equal to' */ -#define O_GE 356 /* comparison on 'not less than' */ -#define O_GT 357 /* comparison on 'greater than' */ -#define O_NE 358 /* comparison on 'not equal to' */ -#define O_AND 359 /* conjunction (logical "and") */ -#define O_OR 360 /* disjunction (logical "or") */ -#define O_UNION 361 /* union */ -#define O_DIFF 362 /* difference */ -#define O_SYMDIFF 363 /* symmetric difference */ -#define O_INTER 364 /* intersection */ -#define O_CROSS 365 /* cross (Cartesian) product */ -#define O_IN 366 /* test on 'x in Y' */ -#define O_NOTIN 367 /* test on 'x not in Y' */ -#define O_WITHIN 368 /* test on 'X within Y' */ -#define O_NOTWITHIN 369 /* test on 'X not within Y' */ -#define O_SUBSTR 370 /* substring */ -#define O_STR2TIME 371 /* convert string to time */ -#define O_TIME2STR 372 /* convert time to string */ +#define O_ADD 341 /* addition */ +#define O_SUB 342 /* subtraction */ +#define O_LESS 343 /* non-negative subtraction */ +#define O_MUL 344 /* multiplication */ +#define O_DIV 345 /* division */ +#define O_IDIV 346 /* quotient of exact division */ +#define O_MOD 347 /* remainder of exact division */ +#define O_POWER 348 /* exponentiation (raise to power) */ +#define O_ATAN2 349 /* trigonometric arctangent */ +#define O_ROUND2 350 /* round to n fractional digits */ +#define O_TRUNC2 351 /* truncate to n fractional digits */ +#define O_UNIFORM 352 /* pseudo-random in [a, b) */ +#define O_NORMAL 353 /* gaussian random, given mu and sigma */ +#define O_CONCAT 354 /* concatenation */ +#define O_LT 355 /* comparison on 'less than' */ +#define O_LE 356 /* comparison on 'not greater than' */ +#define O_EQ 357 /* comparison on 'equal to' */ +#define O_GE 358 /* comparison on 'not less than' */ +#define O_GT 359 /* comparison on 'greater than' */ +#define O_NE 360 /* comparison on 'not equal to' */ +#define O_AND 361 /* conjunction (logical "and") */ +#define O_OR 362 /* disjunction (logical "or") */ +#define O_UNION 363 /* union */ +#define O_DIFF 364 /* difference */ +#define O_SYMDIFF 365 /* symmetric difference */ +#define O_INTER 366 /* intersection */ +#define O_CROSS 367 /* cross (Cartesian) product */ +#define O_IN 368 /* test on 'x in Y' */ +#define O_NOTIN 369 /* test on 'x not in Y' */ +#define O_WITHIN 370 /* test on 'X within Y' */ +#define O_NOTWITHIN 371 /* test on 'X not within Y' */ +#define O_SUBSTR 372 /* substring */ +#define O_STR2TIME 373 /* convert string to time */ +#define O_TIME2STR 374 /* convert time to string */ /* ternary operations ------------------*/ -#define O_DOTS 373 /* build "arithmetic" set */ -#define O_FORK 374 /* if-then-else */ -#define O_SUBSTR3 375 /* substring */ +#define O_DOTS 375 /* build "arithmetic" set */ +#define O_FORK 376 /* if-then-else */ +#define O_SUBSTR3 377 /* substring */ /* n-ary operations --------------------*/ -#define O_MIN 376 /* minimal value (n-ary) */ -#define O_MAX 377 /* maximal value (n-ary) */ +#define O_MIN 378 /* minimal value (n-ary) */ +#define O_MAX 379 /* maximal value (n-ary) */ /* iterated operations -----------------*/ -#define O_SUM 378 /* summation */ -#define O_PROD 379 /* multiplication */ -#define O_MINIMUM 380 /* minimum */ -#define O_MAXIMUM 381 /* maximum */ -#define O_FORALL 382 /* conjunction (A-quantification) */ -#define O_EXISTS 383 /* disjunction (E-quantification) */ -#define O_SETOF 384 /* compute elemental set */ -#define O_BUILD 385 /* build elemental set */ +#define O_SUM 380 /* summation */ +#define O_PROD 381 /* multiplication */ +#define O_MINIMUM 382 /* minimum */ +#define O_MAXIMUM 383 /* maximum */ +#define O_FORALL 384 /* conjunction (A-quantification) */ +#define O_EXISTS 385 /* disjunction (E-quantification) */ +#define O_SETOF 386 /* compute elemental set */ +#define O_BUILD 387 /* build elemental set */ OPERANDS arg; /* operands that participate in the operation */ int type; --- glpk-5.0/src/mpl/mpl1.c 2020-12-16 10:00:00.000000000 +0100 +++ glpk-5.0-mod/src/mpl/mpl1.c 2021-07-24 14:00:35.340950329 +0200 @@ -598,7 +598,9 @@ case O_LOG10: case O_SQRT: case O_SIN: + case O_ASIN: case O_COS: + case O_ACOS: case O_TAN: case O_ATAN: case O_ROUND: @@ -1191,8 +1193,12 @@ op = O_SQRT; else if (strcmp(mpl->image, "sin") == 0) op = O_SIN; + else if (strcmp(mpl->image, "asin") == 0) + op = O_ASIN; else if (strcmp(mpl->image, "cos") == 0) op = O_COS; + else if (strcmp(mpl->image, "acos") == 0) + op = O_ACOS; else if (strcmp(mpl->image, "tan") == 0) op = O_TAN; else if (strcmp(mpl->image, "atan") == 0) --- glpk-5.0/src/mpl/mpl3.c 2020-12-16 10:00:00.000000000 +0100 +++ glpk-5.0-mod/src/mpl/mpl3.c 2021-07-24 14:36:46.492850537 +0200 @@ -214,6 +214,17 @@ } /*---------------------------------------------------------------------- +-- fp_asin - floating-point trigonometric arcsine. +-- +-- This routine computes the trigonometric arcsine asin(x). */ + +double fp_asin(MPL *mpl, double x) +{ if (!(-1 <= x && x <= +1)) + error(mpl, "asin(%.*g); argument too large", DBL_DIG, x); + return asin(x); +} + +/*---------------------------------------------------------------------- -- fp_cos - floating-point trigonometric cosine. -- -- This routine computes the trigonometric cosine cos(x). */ @@ -225,6 +236,17 @@ } /*---------------------------------------------------------------------- +-- fp_acos - floating-point trigonometric arccosine. +-- +-- This routine computes the trigonometric arccosine acos(x). */ + +double fp_acos(MPL *mpl, double x) +{ if (!(-1 <= x && x <= +1)) + error(mpl, "acos(%.*g); argument too large", DBL_DIG, x); + return acos(x); +} + +/*---------------------------------------------------------------------- -- fp_tan - floating-point trigonometric tangent. -- -- This routine computes the trigonometric tangent tan(x). */ @@ -3683,10 +3705,18 @@ /* trigonometric sine */ value = fp_sin(mpl, eval_numeric(mpl, code->arg.arg.x)); break; + case O_ASIN: + /* trigonometric arcsine */ + value = fp_asin(mpl, eval_numeric(mpl, code->arg.arg.x)); + break; case O_COS: /* trigonometric cosine */ value = fp_cos(mpl, eval_numeric(mpl, code->arg.arg.x)); break; + case O_ACOS: + /* trigonometric arccosine */ + value = fp_acos(mpl, eval_numeric(mpl, code->arg.arg.x)); + break; case O_TAN: /* trigonometric tangent */ value = fp_tan(mpl, eval_numeric(mpl, code->arg.arg.x)); @@ -4885,7 +4915,9 @@ case O_LOG10: case O_SQRT: case O_SIN: + case O_ASIN: case O_COS: + case O_ACOS: case O_TAN: case O_ATAN: case O_ROUND: --- glpk-5.0/doc/gmpl.tex 2020-12-16 10:00:00.000000000 +0100 +++ glpk-5.0-mod/doc/gmpl.tex 2021-07-24 16:13:47.320045023 +0200 @@ -582,6 +582,7 @@ {\tt ceil(}$x${\tt)}&$\lceil x\rceil$, smallest integer not less than$x$(ceiling of$x$'')\\ {\tt cos(}$x${\tt)}&$\cos x$, cosine of$x$(in radians)\\ +{\tt acos(}$x${\tt)}&$\arccos x$, arccosine (in radians) of$x$\\ {\tt exp(}$x${\tt)}&$e^x$, base-$e$exponential of$x$\\ {\tt floor(}$x${\tt)}&$\lfloor x\rfloor$, largest integer not greater than$x$(floor of$x$'')\\ @@ -599,6 +600,7 @@ {\tt round(}$x${\tt,}$n${\tt)}&rounding$x$to$n$fractional decimal digits\\ {\tt sin(}$x${\tt)}&$\sin x$, sine of$x$(in radians)\\ +{\tt asin(}$x${\tt)}&$\arcsin x$, arcsine (in radians) of$x$\\ {\tt sqrt(}$x${\tt)}&$\sqrt{x}$, non-negative square root of$x$\\ {\tt str2time(}$s${\tt,}$f${\tt)}&converting character string$s\$ to
calendar time (for details see Section \ref{str2time}, page