1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/src/glplib03.c Mon Dec 06 13:09:21 2010 +0100
1.3 @@ -0,0 +1,692 @@
1.4 +/* glplib03.c (miscellaneous library routines) */
1.5 +
1.6 +/***********************************************************************
1.7 +* This code is part of GLPK (GNU Linear Programming Kit).
1.8 +*
1.9 +* Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
1.10 +* 2009, 2010 Andrew Makhorin, Department for Applied Informatics,
1.11 +* Moscow Aviation Institute, Moscow, Russia. All rights reserved.
1.12 +* E-mail: <mao@gnu.org>.
1.13 +*
1.14 +* GLPK is free software: you can redistribute it and/or modify it
1.15 +* under the terms of the GNU General Public License as published by
1.16 +* the Free Software Foundation, either version 3 of the License, or
1.17 +* (at your option) any later version.
1.18 +*
1.19 +* GLPK is distributed in the hope that it will be useful, but WITHOUT
1.20 +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
1.21 +* or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
1.22 +* License for more details.
1.23 +*
1.24 +* You should have received a copy of the GNU General Public License
1.25 +* along with GLPK. If not, see <http://www.gnu.org/licenses/>.
1.26 +***********************************************************************/
1.27 +
1.28 +#include "glpenv.h"
1.29 +#include "glplib.h"
1.30 +
1.31 +/***********************************************************************
1.32 +* NAME
1.33 +*
1.34 +* str2int - convert character string to value of int type
1.35 +*
1.36 +* SYNOPSIS
1.37 +*
1.38 +* #include "glplib.h"
1.39 +* int str2int(const char *str, int *val);
1.40 +*
1.41 +* DESCRIPTION
1.42 +*
1.43 +* The routine str2int converts the character string str to a value of
1.44 +* integer type and stores the value into location, which the parameter
1.45 +* val points to (in the case of error content of this location is not
1.46 +* changed).
1.47 +*
1.48 +* RETURNS
1.49 +*
1.50 +* The routine returns one of the following error codes:
1.51 +*
1.52 +* 0 - no error;
1.53 +* 1 - value out of range;
1.54 +* 2 - character string is syntactically incorrect. */
1.55 +
1.56 +int str2int(const char *str, int *_val)
1.57 +{ int d, k, s, val = 0;
1.58 + /* scan optional sign */
1.59 + if (str[0] == '+')
1.60 + s = +1, k = 1;
1.61 + else if (str[0] == '-')
1.62 + s = -1, k = 1;
1.63 + else
1.64 + s = +1, k = 0;
1.65 + /* check for the first digit */
1.66 + if (!isdigit((unsigned char)str[k])) return 2;
1.67 + /* scan digits */
1.68 + while (isdigit((unsigned char)str[k]))
1.69 + { d = str[k++] - '0';
1.70 + if (s > 0)
1.71 + { if (val > INT_MAX / 10) return 1;
1.72 + val *= 10;
1.73 + if (val > INT_MAX - d) return 1;
1.74 + val += d;
1.75 + }
1.76 + else
1.77 + { if (val < INT_MIN / 10) return 1;
1.78 + val *= 10;
1.79 + if (val < INT_MIN + d) return 1;
1.80 + val -= d;
1.81 + }
1.82 + }
1.83 + /* check for terminator */
1.84 + if (str[k] != '\0') return 2;
1.85 + /* conversion has been done */
1.86 + *_val = val;
1.87 + return 0;
1.88 +}
1.89 +
1.90 +/***********************************************************************
1.91 +* NAME
1.92 +*
1.93 +* str2num - convert character string to value of double type
1.94 +*
1.95 +* SYNOPSIS
1.96 +*
1.97 +* #include "glplib.h"
1.98 +* int str2num(const char *str, double *val);
1.99 +*
1.100 +* DESCRIPTION
1.101 +*
1.102 +* The routine str2num converts the character string str to a value of
1.103 +* double type and stores the value into location, which the parameter
1.104 +* val points to (in the case of error content of this location is not
1.105 +* changed).
1.106 +*
1.107 +* RETURNS
1.108 +*
1.109 +* The routine returns one of the following error codes:
1.110 +*
1.111 +* 0 - no error;
1.112 +* 1 - value out of range;
1.113 +* 2 - character string is syntactically incorrect. */
1.114 +
1.115 +int str2num(const char *str, double *_val)
1.116 +{ int k;
1.117 + double val;
1.118 + /* scan optional sign */
1.119 + k = (str[0] == '+' || str[0] == '-' ? 1 : 0);
1.120 + /* check for decimal point */
1.121 + if (str[k] == '.')
1.122 + { k++;
1.123 + /* a digit should follow it */
1.124 + if (!isdigit((unsigned char)str[k])) return 2;
1.125 + k++;
1.126 + goto frac;
1.127 + }
1.128 + /* integer part should start with a digit */
1.129 + if (!isdigit((unsigned char)str[k])) return 2;
1.130 + /* scan integer part */
1.131 + while (isdigit((unsigned char)str[k])) k++;
1.132 + /* check for decimal point */
1.133 + if (str[k] == '.') k++;
1.134 +frac: /* scan optional fraction part */
1.135 + while (isdigit((unsigned char)str[k])) k++;
1.136 + /* check for decimal exponent */
1.137 + if (str[k] == 'E' || str[k] == 'e')
1.138 + { k++;
1.139 + /* scan optional sign */
1.140 + if (str[k] == '+' || str[k] == '-') k++;
1.141 + /* a digit should follow E, E+ or E- */
1.142 + if (!isdigit((unsigned char)str[k])) return 2;
1.143 + }
1.144 + /* scan optional exponent part */
1.145 + while (isdigit((unsigned char)str[k])) k++;
1.146 + /* check for terminator */
1.147 + if (str[k] != '\0') return 2;
1.148 + /* perform conversion */
1.149 + { char *endptr;
1.150 + val = strtod(str, &endptr);
1.151 + if (*endptr != '\0') return 2;
1.152 + }
1.153 + /* check for overflow */
1.154 + if (!(-DBL_MAX <= val && val <= +DBL_MAX)) return 1;
1.155 + /* check for underflow */
1.156 + if (-DBL_MIN < val && val < +DBL_MIN) val = 0.0;
1.157 + /* conversion has been done */
1.158 + *_val = val;
1.159 + return 0;
1.160 +}
1.161 +
1.162 +/***********************************************************************
1.163 +* NAME
1.164 +*
1.165 +* strspx - remove all spaces from character string
1.166 +*
1.167 +* SYNOPSIS
1.168 +*
1.169 +* #include "glplib.h"
1.170 +* char *strspx(char *str);
1.171 +*
1.172 +* DESCRIPTION
1.173 +*
1.174 +* The routine strspx removes all spaces from the character string str.
1.175 +*
1.176 +* RETURNS
1.177 +*
1.178 +* The routine returns a pointer to the character string.
1.179 +*
1.180 +* EXAMPLES
1.181 +*
1.182 +* strspx(" Errare humanum est ") => "Errarehumanumest"
1.183 +*
1.184 +* strspx(" ") => "" */
1.185 +
1.186 +char *strspx(char *str)
1.187 +{ char *s, *t;
1.188 + for (s = t = str; *s; s++) if (*s != ' ') *t++ = *s;
1.189 + *t = '\0';
1.190 + return str;
1.191 +}
1.192 +
1.193 +/***********************************************************************
1.194 +* NAME
1.195 +*
1.196 +* strtrim - remove trailing spaces from character string
1.197 +*
1.198 +* SYNOPSIS
1.199 +*
1.200 +* #include "glplib.h"
1.201 +* char *strtrim(char *str);
1.202 +*
1.203 +* DESCRIPTION
1.204 +*
1.205 +* The routine strtrim removes trailing spaces from the character
1.206 +* string str.
1.207 +*
1.208 +* RETURNS
1.209 +*
1.210 +* The routine returns a pointer to the character string.
1.211 +*
1.212 +* EXAMPLES
1.213 +*
1.214 +* strtrim("Errare humanum est ") => "Errare humanum est"
1.215 +*
1.216 +* strtrim(" ") => "" */
1.217 +
1.218 +char *strtrim(char *str)
1.219 +{ char *t;
1.220 + for (t = strrchr(str, '\0') - 1; t >= str; t--)
1.221 + { if (*t != ' ') break;
1.222 + *t = '\0';
1.223 + }
1.224 + return str;
1.225 +}
1.226 +
1.227 +/***********************************************************************
1.228 +* NAME
1.229 +*
1.230 +* strrev - reverse character string
1.231 +*
1.232 +* SYNOPSIS
1.233 +*
1.234 +* #include "glplib.h"
1.235 +* char *strrev(char *s);
1.236 +*
1.237 +* DESCRIPTION
1.238 +*
1.239 +* The routine strrev changes characters in a character string s to the
1.240 +* reverse order, except the terminating null character.
1.241 +*
1.242 +* RETURNS
1.243 +*
1.244 +* The routine returns the pointer s.
1.245 +*
1.246 +* EXAMPLES
1.247 +*
1.248 +* strrev("") => ""
1.249 +*
1.250 +* strrev("Today is Monday") => "yadnoM si yadoT" */
1.251 +
1.252 +char *strrev(char *s)
1.253 +{ int i, j;
1.254 + char t;
1.255 + for (i = 0, j = strlen(s)-1; i < j; i++, j--)
1.256 + t = s[i], s[i] = s[j], s[j] = t;
1.257 + return s;
1.258 +}
1.259 +
1.260 +/***********************************************************************
1.261 +* NAME
1.262 +*
1.263 +* gcd - find greatest common divisor of two integers
1.264 +*
1.265 +* SYNOPSIS
1.266 +*
1.267 +* #include "glplib.h"
1.268 +* int gcd(int x, int y);
1.269 +*
1.270 +* RETURNS
1.271 +*
1.272 +* The routine gcd returns gcd(x, y), the greatest common divisor of
1.273 +* the two positive integers given.
1.274 +*
1.275 +* ALGORITHM
1.276 +*
1.277 +* The routine gcd is based on Euclid's algorithm.
1.278 +*
1.279 +* REFERENCES
1.280 +*
1.281 +* Don Knuth, The Art of Computer Programming, Vol.2: Seminumerical
1.282 +* Algorithms, 3rd Edition, Addison-Wesley, 1997. Section 4.5.2: The
1.283 +* Greatest Common Divisor, pp. 333-56. */
1.284 +
1.285 +int gcd(int x, int y)
1.286 +{ int r;
1.287 + xassert(x > 0 && y > 0);
1.288 + while (y > 0)
1.289 + r = x % y, x = y, y = r;
1.290 + return x;
1.291 +}
1.292 +
1.293 +/***********************************************************************
1.294 +* NAME
1.295 +*
1.296 +* gcdn - find greatest common divisor of n integers
1.297 +*
1.298 +* SYNOPSIS
1.299 +*
1.300 +* #include "glplib.h"
1.301 +* int gcdn(int n, int x[]);
1.302 +*
1.303 +* RETURNS
1.304 +*
1.305 +* The routine gcdn returns gcd(x[1], x[2], ..., x[n]), the greatest
1.306 +* common divisor of n positive integers given, n > 0.
1.307 +*
1.308 +* BACKGROUND
1.309 +*
1.310 +* The routine gcdn is based on the following identity:
1.311 +*
1.312 +* gcd(x, y, z) = gcd(gcd(x, y), z).
1.313 +*
1.314 +* REFERENCES
1.315 +*
1.316 +* Don Knuth, The Art of Computer Programming, Vol.2: Seminumerical
1.317 +* Algorithms, 3rd Edition, Addison-Wesley, 1997. Section 4.5.2: The
1.318 +* Greatest Common Divisor, pp. 333-56. */
1.319 +
1.320 +int gcdn(int n, int x[])
1.321 +{ int d, j;
1.322 + xassert(n > 0);
1.323 + for (j = 1; j <= n; j++)
1.324 + { xassert(x[j] > 0);
1.325 + if (j == 1)
1.326 + d = x[1];
1.327 + else
1.328 + d = gcd(d, x[j]);
1.329 + if (d == 1) break;
1.330 + }
1.331 + return d;
1.332 +}
1.333 +
1.334 +/***********************************************************************
1.335 +* NAME
1.336 +*
1.337 +* lcm - find least common multiple of two integers
1.338 +*
1.339 +* SYNOPSIS
1.340 +*
1.341 +* #include "glplib.h"
1.342 +* int lcm(int x, int y);
1.343 +*
1.344 +* RETURNS
1.345 +*
1.346 +* The routine lcm returns lcm(x, y), the least common multiple of the
1.347 +* two positive integers given. In case of integer overflow the routine
1.348 +* returns zero.
1.349 +*
1.350 +* BACKGROUND
1.351 +*
1.352 +* The routine lcm is based on the following identity:
1.353 +*
1.354 +* lcm(x, y) = (x * y) / gcd(x, y) = x * [y / gcd(x, y)],
1.355 +*
1.356 +* where gcd(x, y) is the greatest common divisor of x and y. */
1.357 +
1.358 +int lcm(int x, int y)
1.359 +{ xassert(x > 0);
1.360 + xassert(y > 0);
1.361 + y /= gcd(x, y);
1.362 + if (x > INT_MAX / y) return 0;
1.363 + return x * y;
1.364 +}
1.365 +
1.366 +/***********************************************************************
1.367 +* NAME
1.368 +*
1.369 +* lcmn - find least common multiple of n integers
1.370 +*
1.371 +* SYNOPSIS
1.372 +*
1.373 +* #include "glplib.h"
1.374 +* int lcmn(int n, int x[]);
1.375 +*
1.376 +* RETURNS
1.377 +*
1.378 +* The routine lcmn returns lcm(x[1], x[2], ..., x[n]), the least
1.379 +* common multiple of n positive integers given, n > 0. In case of
1.380 +* integer overflow the routine returns zero.
1.381 +*
1.382 +* BACKGROUND
1.383 +*
1.384 +* The routine lcmn is based on the following identity:
1.385 +*
1.386 +* lcmn(x, y, z) = lcm(lcm(x, y), z),
1.387 +*
1.388 +* where lcm(x, y) is the least common multiple of x and y. */
1.389 +
1.390 +int lcmn(int n, int x[])
1.391 +{ int m, j;
1.392 + xassert(n > 0);
1.393 + for (j = 1; j <= n; j++)
1.394 + { xassert(x[j] > 0);
1.395 + if (j == 1)
1.396 + m = x[1];
1.397 + else
1.398 + m = lcm(m, x[j]);
1.399 + if (m == 0) break;
1.400 + }
1.401 + return m;
1.402 +}
1.403 +
1.404 +/***********************************************************************
1.405 +* NAME
1.406 +*
1.407 +* round2n - round floating-point number to nearest power of two
1.408 +*
1.409 +* SYNOPSIS
1.410 +*
1.411 +* #include "glplib.h"
1.412 +* double round2n(double x);
1.413 +*
1.414 +* RETURNS
1.415 +*
1.416 +* Given a positive floating-point value x the routine round2n returns
1.417 +* 2^n such that |x - 2^n| is minimal.
1.418 +*
1.419 +* EXAMPLES
1.420 +*
1.421 +* round2n(10.1) = 2^3 = 8
1.422 +* round2n(15.3) = 2^4 = 16
1.423 +* round2n(0.01) = 2^(-7) = 0.0078125
1.424 +*
1.425 +* BACKGROUND
1.426 +*
1.427 +* Let x = f * 2^e, where 0.5 <= f < 1 is a normalized fractional part,
1.428 +* e is an integer exponent. Then, obviously, 0.5 * 2^e <= x < 2^e, so
1.429 +* if x - 0.5 * 2^e <= 2^e - x, we choose 0.5 * 2^e = 2^(e-1), and 2^e
1.430 +* otherwise. The latter condition can be written as 2 * x <= 1.5 * 2^e
1.431 +* or 2 * f * 2^e <= 1.5 * 2^e or, finally, f <= 0.75. */
1.432 +
1.433 +double round2n(double x)
1.434 +{ int e;
1.435 + double f;
1.436 + xassert(x > 0.0);
1.437 + f = frexp(x, &e);
1.438 + return ldexp(1.0, f <= 0.75 ? e-1 : e);
1.439 +}
1.440 +
1.441 +/***********************************************************************
1.442 +* NAME
1.443 +*
1.444 +* fp2rat - convert floating-point number to rational number
1.445 +*
1.446 +* SYNOPSIS
1.447 +*
1.448 +* #include "glplib.h"
1.449 +* int fp2rat(double x, double eps, double *p, double *q);
1.450 +*
1.451 +* DESCRIPTION
1.452 +*
1.453 +* Given a floating-point number 0 <= x < 1 the routine fp2rat finds
1.454 +* its "best" rational approximation p / q, where p >= 0 and q > 0 are
1.455 +* integer numbers, such that |x - p / q| <= eps.
1.456 +*
1.457 +* RETURNS
1.458 +*
1.459 +* The routine fp2rat returns the number of iterations used to achieve
1.460 +* the specified precision eps.
1.461 +*
1.462 +* EXAMPLES
1.463 +*
1.464 +* For x = sqrt(2) - 1 = 0.414213562373095 and eps = 1e-6 the routine
1.465 +* gives p = 408 and q = 985, where 408 / 985 = 0.414213197969543.
1.466 +*
1.467 +* BACKGROUND
1.468 +*
1.469 +* It is well known that every positive real number x can be expressed
1.470 +* as the following continued fraction:
1.471 +*
1.472 +* x = b[0] + a[1]
1.473 +* ------------------------
1.474 +* b[1] + a[2]
1.475 +* -----------------
1.476 +* b[2] + a[3]
1.477 +* ----------
1.478 +* b[3] + ...
1.479 +*
1.480 +* where:
1.481 +*
1.482 +* a[k] = 1, k = 0, 1, 2, ...
1.483 +*
1.484 +* b[k] = floor(x[k]), k = 0, 1, 2, ...
1.485 +*
1.486 +* x[0] = x,
1.487 +*
1.488 +* x[k] = 1 / frac(x[k-1]), k = 1, 2, 3, ...
1.489 +*
1.490 +* To find the "best" rational approximation of x the routine computes
1.491 +* partial fractions f[k] by dropping after k terms as follows:
1.492 +*
1.493 +* f[k] = A[k] / B[k],
1.494 +*
1.495 +* where:
1.496 +*
1.497 +* A[-1] = 1, A[0] = b[0], B[-1] = 0, B[0] = 1,
1.498 +*
1.499 +* A[k] = b[k] * A[k-1] + a[k] * A[k-2],
1.500 +*
1.501 +* B[k] = b[k] * B[k-1] + a[k] * B[k-2].
1.502 +*
1.503 +* Once the condition
1.504 +*
1.505 +* |x - f[k]| <= eps
1.506 +*
1.507 +* has been satisfied, the routine reports p = A[k] and q = B[k] as the
1.508 +* final answer.
1.509 +*
1.510 +* In the table below here is some statistics obtained for one million
1.511 +* random numbers uniformly distributed in the range [0, 1).
1.512 +*
1.513 +* eps max p mean p max q mean q max k mean k
1.514 +* -------------------------------------------------------------
1.515 +* 1e-1 8 1.6 9 3.2 3 1.4
1.516 +* 1e-2 98 6.2 99 12.4 5 2.4
1.517 +* 1e-3 997 20.7 998 41.5 8 3.4
1.518 +* 1e-4 9959 66.6 9960 133.5 10 4.4
1.519 +* 1e-5 97403 211.7 97404 424.2 13 5.3
1.520 +* 1e-6 479669 669.9 479670 1342.9 15 6.3
1.521 +* 1e-7 1579030 2127.3 3962146 4257.8 16 7.3
1.522 +* 1e-8 26188823 6749.4 26188824 13503.4 19 8.2
1.523 +*
1.524 +* REFERENCES
1.525 +*
1.526 +* W. B. Jones and W. J. Thron, "Continued Fractions: Analytic Theory
1.527 +* and Applications," Encyclopedia on Mathematics and Its Applications,
1.528 +* Addison-Wesley, 1980. */
1.529 +
1.530 +int fp2rat(double x, double eps, double *p, double *q)
1.531 +{ int k;
1.532 + double xk, Akm1, Ak, Bkm1, Bk, ak, bk, fk, temp;
1.533 + if (!(0.0 <= x && x < 1.0))
1.534 + xerror("fp2rat: x = %g; number out of range\n", x);
1.535 + for (k = 0; ; k++)
1.536 + { xassert(k <= 100);
1.537 + if (k == 0)
1.538 + { /* x[0] = x */
1.539 + xk = x;
1.540 + /* A[-1] = 1 */
1.541 + Akm1 = 1.0;
1.542 + /* A[0] = b[0] = floor(x[0]) = 0 */
1.543 + Ak = 0.0;
1.544 + /* B[-1] = 0 */
1.545 + Bkm1 = 0.0;
1.546 + /* B[0] = 1 */
1.547 + Bk = 1.0;
1.548 + }
1.549 + else
1.550 + { /* x[k] = 1 / frac(x[k-1]) */
1.551 + temp = xk - floor(xk);
1.552 + xassert(temp != 0.0);
1.553 + xk = 1.0 / temp;
1.554 + /* a[k] = 1 */
1.555 + ak = 1.0;
1.556 + /* b[k] = floor(x[k]) */
1.557 + bk = floor(xk);
1.558 + /* A[k] = b[k] * A[k-1] + a[k] * A[k-2] */
1.559 + temp = bk * Ak + ak * Akm1;
1.560 + Akm1 = Ak, Ak = temp;
1.561 + /* B[k] = b[k] * B[k-1] + a[k] * B[k-2] */
1.562 + temp = bk * Bk + ak * Bkm1;
1.563 + Bkm1 = Bk, Bk = temp;
1.564 + }
1.565 + /* f[k] = A[k] / B[k] */
1.566 + fk = Ak / Bk;
1.567 +#if 0
1.568 + print("%.*g / %.*g = %.*g", DBL_DIG, Ak, DBL_DIG, Bk, DBL_DIG,
1.569 + fk);
1.570 +#endif
1.571 + if (fabs(x - fk) <= eps) break;
1.572 + }
1.573 + *p = Ak;
1.574 + *q = Bk;
1.575 + return k;
1.576 +}
1.577 +
1.578 +/***********************************************************************
1.579 +* NAME
1.580 +*
1.581 +* jday - convert calendar date to Julian day number
1.582 +*
1.583 +* SYNOPSIS
1.584 +*
1.585 +* #include "glplib.h"
1.586 +* int jday(int d, int m, int y);
1.587 +*
1.588 +* DESCRIPTION
1.589 +*
1.590 +* The routine jday converts a calendar date, Gregorian calendar, to
1.591 +* corresponding Julian day number j.
1.592 +*
1.593 +* From the given day d, month m, and year y, the Julian day number j
1.594 +* is computed without using tables.
1.595 +*
1.596 +* The routine is valid for 1 <= y <= 4000.
1.597 +*
1.598 +* RETURNS
1.599 +*
1.600 +* The routine jday returns the Julian day number, or negative value if
1.601 +* the specified date is incorrect.
1.602 +*
1.603 +* REFERENCES
1.604 +*
1.605 +* R. G. Tantzen, Algorithm 199: conversions between calendar date and
1.606 +* Julian day number, Communications of the ACM, vol. 6, no. 8, p. 444,
1.607 +* Aug. 1963. */
1.608 +
1.609 +int jday(int d, int m, int y)
1.610 +{ int c, ya, j, dd;
1.611 + if (!(1 <= d && d <= 31 && 1 <= m && m <= 12 && 1 <= y &&
1.612 + y <= 4000))
1.613 + { j = -1;
1.614 + goto done;
1.615 + }
1.616 + if (m >= 3) m -= 3; else m += 9, y--;
1.617 + c = y / 100;
1.618 + ya = y - 100 * c;
1.619 + j = (146097 * c) / 4 + (1461 * ya) / 4 + (153 * m + 2) / 5 + d +
1.620 + 1721119;
1.621 + jdate(j, &dd, NULL, NULL);
1.622 + if (d != dd) j = -1;
1.623 +done: return j;
1.624 +}
1.625 +
1.626 +/***********************************************************************
1.627 +* NAME
1.628 +*
1.629 +* jdate - convert Julian day number to calendar date
1.630 +*
1.631 +* SYNOPSIS
1.632 +*
1.633 +* #include "glplib.h"
1.634 +* void jdate(int j, int *d, int *m, int *y);
1.635 +*
1.636 +* DESCRIPTION
1.637 +*
1.638 +* The routine jdate converts a Julian day number j to corresponding
1.639 +* calendar date, Gregorian calendar.
1.640 +*
1.641 +* The day d, month m, and year y are computed without using tables and
1.642 +* stored in corresponding locations.
1.643 +*
1.644 +* The routine is valid for 1721426 <= j <= 3182395.
1.645 +*
1.646 +* RETURNS
1.647 +*
1.648 +* If the conversion is successful, the routine returns zero, otherwise
1.649 +* non-zero.
1.650 +*
1.651 +* REFERENCES
1.652 +*
1.653 +* R. G. Tantzen, Algorithm 199: conversions between calendar date and
1.654 +* Julian day number, Communications of the ACM, vol. 6, no. 8, p. 444,
1.655 +* Aug. 1963. */
1.656 +
1.657 +int jdate(int j, int *_d, int *_m, int *_y)
1.658 +{ int d, m, y, ret = 0;
1.659 + if (!(1721426 <= j && j <= 3182395))
1.660 + { ret = 1;
1.661 + goto done;
1.662 + }
1.663 + j -= 1721119;
1.664 + y = (4 * j - 1) / 146097;
1.665 + j = (4 * j - 1) % 146097;
1.666 + d = j / 4;
1.667 + j = (4 * d + 3) / 1461;
1.668 + d = (4 * d + 3) % 1461;
1.669 + d = (d + 4) / 4;
1.670 + m = (5 * d - 3) / 153;
1.671 + d = (5 * d - 3) % 153;
1.672 + d = (d + 5) / 5;
1.673 + y = 100 * y + j;
1.674 + if (m <= 9) m += 3; else m -= 9, y++;
1.675 + if (_d != NULL) *_d = d;
1.676 + if (_m != NULL) *_m = m;
1.677 + if (_y != NULL) *_y = y;
1.678 +done: return ret;
1.679 +}
1.680 +
1.681 +#if 0
1.682 +int main(void)
1.683 +{ int jbeg, jend, j, d, m, y;
1.684 + jbeg = jday(1, 1, 1);
1.685 + jend = jday(31, 12, 4000);
1.686 + for (j = jbeg; j <= jend; j++)
1.687 + { xassert(jdate(j, &d, &m, &y) == 0);
1.688 + xassert(jday(d, m, y) == j);
1.689 + }
1.690 + xprintf("Routines jday and jdate work correctly.\n");
1.691 + return 0;
1.692 +}
1.693 +#endif
1.694 +
1.695 +/* eof */