1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/deps/glpk/src/glplib03.c	Sun Nov 06 20:59:10 2011 +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, 2011 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 */