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