/* glplib03.c (miscellaneous library routines) */ /*********************************************************************** * This code is part of GLPK (GNU Linear Programming Kit). * * Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, * 2009, 2010, 2011 Andrew Makhorin, Department for Applied Informatics, * Moscow Aviation Institute, Moscow, Russia. All rights reserved. * E-mail: . * * GLPK is free software: you can redistribute it and/or modify it * under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * GLPK is distributed in the hope that it will be useful, but WITHOUT * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public * License for more details. * * You should have received a copy of the GNU General Public License * along with GLPK. If not, see . ***********************************************************************/ #include "glpenv.h" #include "glplib.h" /*********************************************************************** * NAME * * str2int - convert character string to value of int type * * SYNOPSIS * * #include "glplib.h" * int str2int(const char *str, int *val); * * DESCRIPTION * * The routine str2int converts the character string str to a value of * integer type and stores the value into location, which the parameter * val points to (in the case of error content of this location is not * changed). * * RETURNS * * The routine returns one of the following error codes: * * 0 - no error; * 1 - value out of range; * 2 - character string is syntactically incorrect. */ int str2int(const char *str, int *_val) { int d, k, s, val = 0; /* scan optional sign */ if (str[0] == '+') s = +1, k = 1; else if (str[0] == '-') s = -1, k = 1; else s = +1, k = 0; /* check for the first digit */ if (!isdigit((unsigned char)str[k])) return 2; /* scan digits */ while (isdigit((unsigned char)str[k])) { d = str[k++] - '0'; if (s > 0) { if (val > INT_MAX / 10) return 1; val *= 10; if (val > INT_MAX - d) return 1; val += d; } else { if (val < INT_MIN / 10) return 1; val *= 10; if (val < INT_MIN + d) return 1; val -= d; } } /* check for terminator */ if (str[k] != '\0') return 2; /* conversion has been done */ *_val = val; return 0; } /*********************************************************************** * NAME * * str2num - convert character string to value of double type * * SYNOPSIS * * #include "glplib.h" * int str2num(const char *str, double *val); * * DESCRIPTION * * The routine str2num converts the character string str to a value of * double type and stores the value into location, which the parameter * val points to (in the case of error content of this location is not * changed). * * RETURNS * * The routine returns one of the following error codes: * * 0 - no error; * 1 - value out of range; * 2 - character string is syntactically incorrect. */ int str2num(const char *str, double *_val) { int k; double val; /* scan optional sign */ k = (str[0] == '+' || str[0] == '-' ? 1 : 0); /* check for decimal point */ if (str[k] == '.') { k++; /* a digit should follow it */ if (!isdigit((unsigned char)str[k])) return 2; k++; goto frac; } /* integer part should start with a digit */ if (!isdigit((unsigned char)str[k])) return 2; /* scan integer part */ while (isdigit((unsigned char)str[k])) k++; /* check for decimal point */ if (str[k] == '.') k++; frac: /* scan optional fraction part */ while (isdigit((unsigned char)str[k])) k++; /* check for decimal exponent */ if (str[k] == 'E' || str[k] == 'e') { k++; /* scan optional sign */ if (str[k] == '+' || str[k] == '-') k++; /* a digit should follow E, E+ or E- */ if (!isdigit((unsigned char)str[k])) return 2; } /* scan optional exponent part */ while (isdigit((unsigned char)str[k])) k++; /* check for terminator */ if (str[k] != '\0') return 2; /* perform conversion */ { char *endptr; val = strtod(str, &endptr); if (*endptr != '\0') return 2; } /* check for overflow */ if (!(-DBL_MAX <= val && val <= +DBL_MAX)) return 1; /* check for underflow */ if (-DBL_MIN < val && val < +DBL_MIN) val = 0.0; /* conversion has been done */ *_val = val; return 0; } /*********************************************************************** * NAME * * strspx - remove all spaces from character string * * SYNOPSIS * * #include "glplib.h" * char *strspx(char *str); * * DESCRIPTION * * The routine strspx removes all spaces from the character string str. * * RETURNS * * The routine returns a pointer to the character string. * * EXAMPLES * * strspx(" Errare humanum est ") => "Errarehumanumest" * * strspx(" ") => "" */ char *strspx(char *str) { char *s, *t; for (s = t = str; *s; s++) if (*s != ' ') *t++ = *s; *t = '\0'; return str; } /*********************************************************************** * NAME * * strtrim - remove trailing spaces from character string * * SYNOPSIS * * #include "glplib.h" * char *strtrim(char *str); * * DESCRIPTION * * The routine strtrim removes trailing spaces from the character * string str. * * RETURNS * * The routine returns a pointer to the character string. * * EXAMPLES * * strtrim("Errare humanum est ") => "Errare humanum est" * * strtrim(" ") => "" */ char *strtrim(char *str) { char *t; for (t = strrchr(str, '\0') - 1; t >= str; t--) { if (*t != ' ') break; *t = '\0'; } return str; } /*********************************************************************** * NAME * * strrev - reverse character string * * SYNOPSIS * * #include "glplib.h" * char *strrev(char *s); * * DESCRIPTION * * The routine strrev changes characters in a character string s to the * reverse order, except the terminating null character. * * RETURNS * * The routine returns the pointer s. * * EXAMPLES * * strrev("") => "" * * strrev("Today is Monday") => "yadnoM si yadoT" */ char *strrev(char *s) { int i, j; char t; for (i = 0, j = strlen(s)-1; i < j; i++, j--) t = s[i], s[i] = s[j], s[j] = t; return s; } /*********************************************************************** * NAME * * gcd - find greatest common divisor of two integers * * SYNOPSIS * * #include "glplib.h" * int gcd(int x, int y); * * RETURNS * * The routine gcd returns gcd(x, y), the greatest common divisor of * the two positive integers given. * * ALGORITHM * * The routine gcd is based on Euclid's algorithm. * * REFERENCES * * Don Knuth, The Art of Computer Programming, Vol.2: Seminumerical * Algorithms, 3rd Edition, Addison-Wesley, 1997. Section 4.5.2: The * Greatest Common Divisor, pp. 333-56. */ int gcd(int x, int y) { int r; xassert(x > 0 && y > 0); while (y > 0) r = x % y, x = y, y = r; return x; } /*********************************************************************** * NAME * * gcdn - find greatest common divisor of n integers * * SYNOPSIS * * #include "glplib.h" * int gcdn(int n, int x[]); * * RETURNS * * The routine gcdn returns gcd(x[1], x[2], ..., x[n]), the greatest * common divisor of n positive integers given, n > 0. * * BACKGROUND * * The routine gcdn is based on the following identity: * * gcd(x, y, z) = gcd(gcd(x, y), z). * * REFERENCES * * Don Knuth, The Art of Computer Programming, Vol.2: Seminumerical * Algorithms, 3rd Edition, Addison-Wesley, 1997. Section 4.5.2: The * Greatest Common Divisor, pp. 333-56. */ int gcdn(int n, int x[]) { int d, j; xassert(n > 0); for (j = 1; j <= n; j++) { xassert(x[j] > 0); if (j == 1) d = x[1]; else d = gcd(d, x[j]); if (d == 1) break; } return d; } /*********************************************************************** * NAME * * lcm - find least common multiple of two integers * * SYNOPSIS * * #include "glplib.h" * int lcm(int x, int y); * * RETURNS * * The routine lcm returns lcm(x, y), the least common multiple of the * two positive integers given. In case of integer overflow the routine * returns zero. * * BACKGROUND * * The routine lcm is based on the following identity: * * lcm(x, y) = (x * y) / gcd(x, y) = x * [y / gcd(x, y)], * * where gcd(x, y) is the greatest common divisor of x and y. */ int lcm(int x, int y) { xassert(x > 0); xassert(y > 0); y /= gcd(x, y); if (x > INT_MAX / y) return 0; return x * y; } /*********************************************************************** * NAME * * lcmn - find least common multiple of n integers * * SYNOPSIS * * #include "glplib.h" * int lcmn(int n, int x[]); * * RETURNS * * The routine lcmn returns lcm(x[1], x[2], ..., x[n]), the least * common multiple of n positive integers given, n > 0. In case of * integer overflow the routine returns zero. * * BACKGROUND * * The routine lcmn is based on the following identity: * * lcmn(x, y, z) = lcm(lcm(x, y), z), * * where lcm(x, y) is the least common multiple of x and y. */ int lcmn(int n, int x[]) { int m, j; xassert(n > 0); for (j = 1; j <= n; j++) { xassert(x[j] > 0); if (j == 1) m = x[1]; else m = lcm(m, x[j]); if (m == 0) break; } return m; } /*********************************************************************** * NAME * * round2n - round floating-point number to nearest power of two * * SYNOPSIS * * #include "glplib.h" * double round2n(double x); * * RETURNS * * Given a positive floating-point value x the routine round2n returns * 2^n such that |x - 2^n| is minimal. * * EXAMPLES * * round2n(10.1) = 2^3 = 8 * round2n(15.3) = 2^4 = 16 * round2n(0.01) = 2^(-7) = 0.0078125 * * BACKGROUND * * Let x = f * 2^e, where 0.5 <= f < 1 is a normalized fractional part, * e is an integer exponent. Then, obviously, 0.5 * 2^e <= x < 2^e, so * if x - 0.5 * 2^e <= 2^e - x, we choose 0.5 * 2^e = 2^(e-1), and 2^e * otherwise. The latter condition can be written as 2 * x <= 1.5 * 2^e * or 2 * f * 2^e <= 1.5 * 2^e or, finally, f <= 0.75. */ double round2n(double x) { int e; double f; xassert(x > 0.0); f = frexp(x, &e); return ldexp(1.0, f <= 0.75 ? e-1 : e); } /*********************************************************************** * NAME * * fp2rat - convert floating-point number to rational number * * SYNOPSIS * * #include "glplib.h" * int fp2rat(double x, double eps, double *p, double *q); * * DESCRIPTION * * Given a floating-point number 0 <= x < 1 the routine fp2rat finds * its "best" rational approximation p / q, where p >= 0 and q > 0 are * integer numbers, such that |x - p / q| <= eps. * * RETURNS * * The routine fp2rat returns the number of iterations used to achieve * the specified precision eps. * * EXAMPLES * * For x = sqrt(2) - 1 = 0.414213562373095 and eps = 1e-6 the routine * gives p = 408 and q = 985, where 408 / 985 = 0.414213197969543. * * BACKGROUND * * It is well known that every positive real number x can be expressed * as the following continued fraction: * * x = b[0] + a[1] * ------------------------ * b[1] + a[2] * ----------------- * b[2] + a[3] * ---------- * b[3] + ... * * where: * * a[k] = 1, k = 0, 1, 2, ... * * b[k] = floor(x[k]), k = 0, 1, 2, ... * * x[0] = x, * * x[k] = 1 / frac(x[k-1]), k = 1, 2, 3, ... * * To find the "best" rational approximation of x the routine computes * partial fractions f[k] by dropping after k terms as follows: * * f[k] = A[k] / B[k], * * where: * * A[-1] = 1, A[0] = b[0], B[-1] = 0, B[0] = 1, * * A[k] = b[k] * A[k-1] + a[k] * A[k-2], * * B[k] = b[k] * B[k-1] + a[k] * B[k-2]. * * Once the condition * * |x - f[k]| <= eps * * has been satisfied, the routine reports p = A[k] and q = B[k] as the * final answer. * * In the table below here is some statistics obtained for one million * random numbers uniformly distributed in the range [0, 1). * * eps max p mean p max q mean q max k mean k * ------------------------------------------------------------- * 1e-1 8 1.6 9 3.2 3 1.4 * 1e-2 98 6.2 99 12.4 5 2.4 * 1e-3 997 20.7 998 41.5 8 3.4 * 1e-4 9959 66.6 9960 133.5 10 4.4 * 1e-5 97403 211.7 97404 424.2 13 5.3 * 1e-6 479669 669.9 479670 1342.9 15 6.3 * 1e-7 1579030 2127.3 3962146 4257.8 16 7.3 * 1e-8 26188823 6749.4 26188824 13503.4 19 8.2 * * REFERENCES * * W. B. Jones and W. J. Thron, "Continued Fractions: Analytic Theory * and Applications," Encyclopedia on Mathematics and Its Applications, * Addison-Wesley, 1980. */ int fp2rat(double x, double eps, double *p, double *q) { int k; double xk, Akm1, Ak, Bkm1, Bk, ak, bk, fk, temp; if (!(0.0 <= x && x < 1.0)) xerror("fp2rat: x = %g; number out of range\n", x); for (k = 0; ; k++) { xassert(k <= 100); if (k == 0) { /* x[0] = x */ xk = x; /* A[-1] = 1 */ Akm1 = 1.0; /* A[0] = b[0] = floor(x[0]) = 0 */ Ak = 0.0; /* B[-1] = 0 */ Bkm1 = 0.0; /* B[0] = 1 */ Bk = 1.0; } else { /* x[k] = 1 / frac(x[k-1]) */ temp = xk - floor(xk); xassert(temp != 0.0); xk = 1.0 / temp; /* a[k] = 1 */ ak = 1.0; /* b[k] = floor(x[k]) */ bk = floor(xk); /* A[k] = b[k] * A[k-1] + a[k] * A[k-2] */ temp = bk * Ak + ak * Akm1; Akm1 = Ak, Ak = temp; /* B[k] = b[k] * B[k-1] + a[k] * B[k-2] */ temp = bk * Bk + ak * Bkm1; Bkm1 = Bk, Bk = temp; } /* f[k] = A[k] / B[k] */ fk = Ak / Bk; #if 0 print("%.*g / %.*g = %.*g", DBL_DIG, Ak, DBL_DIG, Bk, DBL_DIG, fk); #endif if (fabs(x - fk) <= eps) break; } *p = Ak; *q = Bk; return k; } /*********************************************************************** * NAME * * jday - convert calendar date to Julian day number * * SYNOPSIS * * #include "glplib.h" * int jday(int d, int m, int y); * * DESCRIPTION * * The routine jday converts a calendar date, Gregorian calendar, to * corresponding Julian day number j. * * From the given day d, month m, and year y, the Julian day number j * is computed without using tables. * * The routine is valid for 1 <= y <= 4000. * * RETURNS * * The routine jday returns the Julian day number, or negative value if * the specified date is incorrect. * * REFERENCES * * R. G. Tantzen, Algorithm 199: conversions between calendar date and * Julian day number, Communications of the ACM, vol. 6, no. 8, p. 444, * Aug. 1963. */ int jday(int d, int m, int y) { int c, ya, j, dd; if (!(1 <= d && d <= 31 && 1 <= m && m <= 12 && 1 <= y && y <= 4000)) { j = -1; goto done; } if (m >= 3) m -= 3; else m += 9, y--; c = y / 100; ya = y - 100 * c; j = (146097 * c) / 4 + (1461 * ya) / 4 + (153 * m + 2) / 5 + d + 1721119; jdate(j, &dd, NULL, NULL); if (d != dd) j = -1; done: return j; } /*********************************************************************** * NAME * * jdate - convert Julian day number to calendar date * * SYNOPSIS * * #include "glplib.h" * void jdate(int j, int *d, int *m, int *y); * * DESCRIPTION * * The routine jdate converts a Julian day number j to corresponding * calendar date, Gregorian calendar. * * The day d, month m, and year y are computed without using tables and * stored in corresponding locations. * * The routine is valid for 1721426 <= j <= 3182395. * * RETURNS * * If the conversion is successful, the routine returns zero, otherwise * non-zero. * * REFERENCES * * R. G. Tantzen, Algorithm 199: conversions between calendar date and * Julian day number, Communications of the ACM, vol. 6, no. 8, p. 444, * Aug. 1963. */ int jdate(int j, int *_d, int *_m, int *_y) { int d, m, y, ret = 0; if (!(1721426 <= j && j <= 3182395)) { ret = 1; goto done; } j -= 1721119; y = (4 * j - 1) / 146097; j = (4 * j - 1) % 146097; d = j / 4; j = (4 * d + 3) / 1461; d = (4 * d + 3) % 1461; d = (d + 4) / 4; m = (5 * d - 3) / 153; d = (5 * d - 3) % 153; d = (d + 5) / 5; y = 100 * y + j; if (m <= 9) m += 3; else m -= 9, y++; if (_d != NULL) *_d = d; if (_m != NULL) *_m = m; if (_y != NULL) *_y = y; done: return ret; } #if 0 int main(void) { int jbeg, jend, j, d, m, y; jbeg = jday(1, 1, 1); jend = jday(31, 12, 4000); for (j = jbeg; j <= jend; j++) { xassert(jdate(j, &d, &m, &y) == 0); xassert(jday(d, m, y) == j); } xprintf("Routines jday and jdate work correctly.\n"); return 0; } #endif /* eof */