[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
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
| [Prev in Thread] |
Current Thread |
[Next in Thread] |
- GLPK: Adding asin() and acos() to GNU MathProg,
Hartmut Henkel <=