lemon-project-template-glpk

annotate deps/glpk/src/glplib03.c @ 9:33de93886c88

find changesets by author, revision, files, or words in the commit message
Import GLPK 4.47
author Alpar Juttner <alpar@cs.elte.hu>
date Sun, 06 Nov 2011 20:59:10 +0100
parents
children
rev   line source
alpar@9 1 /* glplib03.c (miscellaneous library routines) */
alpar@9 2
alpar@9 3 /***********************************************************************
alpar@9 4 * This code is part of GLPK (GNU Linear Programming Kit).
alpar@9 5 *
alpar@9 6 * Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
alpar@9 7 * 2009, 2010, 2011 Andrew Makhorin, Department for Applied Informatics,
alpar@9 8 * Moscow Aviation Institute, Moscow, Russia. All rights reserved.
alpar@9 9 * E-mail: <mao@gnu.org>.
alpar@9 10 *
alpar@9 11 * GLPK is free software: you can redistribute it and/or modify it
alpar@9 12 * under the terms of the GNU General Public License as published by
alpar@9 13 * the Free Software Foundation, either version 3 of the License, or
alpar@9 14 * (at your option) any later version.
alpar@9 15 *
alpar@9 16 * GLPK is distributed in the hope that it will be useful, but WITHOUT
alpar@9 17 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
alpar@9 18 * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
alpar@9 19 * License for more details.
alpar@9 20 *
alpar@9 21 * You should have received a copy of the GNU General Public License
alpar@9 22 * along with GLPK. If not, see <http://www.gnu.org/licenses/>.
alpar@9 23 ***********************************************************************/
alpar@9 24
alpar@9 25 #include "glpenv.h"
alpar@9 26 #include "glplib.h"
alpar@9 27
alpar@9 28 /***********************************************************************
alpar@9 29 * NAME
alpar@9 30 *
alpar@9 31 * str2int - convert character string to value of int type
alpar@9 32 *
alpar@9 33 * SYNOPSIS
alpar@9 34 *
alpar@9 35 * #include "glplib.h"
alpar@9 36 * int str2int(const char *str, int *val);
alpar@9 37 *
alpar@9 38 * DESCRIPTION
alpar@9 39 *
alpar@9 40 * The routine str2int converts the character string str to a value of
alpar@9 41 * integer type and stores the value into location, which the parameter
alpar@9 42 * val points to (in the case of error content of this location is not
alpar@9 43 * changed).
alpar@9 44 *
alpar@9 45 * RETURNS
alpar@9 46 *
alpar@9 47 * The routine returns one of the following error codes:
alpar@9 48 *
alpar@9 49 * 0 - no error;
alpar@9 50 * 1 - value out of range;
alpar@9 51 * 2 - character string is syntactically incorrect. */
alpar@9 52
alpar@9 53 int str2int(const char *str, int *_val)
alpar@9 54 { int d, k, s, val = 0;
alpar@9 55 /* scan optional sign */
alpar@9 56 if (str[0] == '+')
alpar@9 57 s = +1, k = 1;
alpar@9 58 else if (str[0] == '-')
alpar@9 59 s = -1, k = 1;
alpar@9 60 else
alpar@9 61 s = +1, k = 0;
alpar@9 62 /* check for the first digit */
alpar@9 63 if (!isdigit((unsigned char)str[k])) return 2;
alpar@9 64 /* scan digits */
alpar@9 65 while (isdigit((unsigned char)str[k]))
alpar@9 66 { d = str[k++] - '0';
alpar@9 67 if (s > 0)
alpar@9 68 { if (val > INT_MAX / 10) return 1;
alpar@9 69 val *= 10;
alpar@9 70 if (val > INT_MAX - d) return 1;
alpar@9 71 val += d;
alpar@9 72 }
alpar@9 73 else
alpar@9 74 { if (val < INT_MIN / 10) return 1;
alpar@9 75 val *= 10;
alpar@9 76 if (val < INT_MIN + d) return 1;
alpar@9 77 val -= d;
alpar@9 78 }
alpar@9 79 }
alpar@9 80 /* check for terminator */
alpar@9 81 if (str[k] != '\0') return 2;
alpar@9 82 /* conversion has been done */
alpar@9 83 *_val = val;
alpar@9 84 return 0;
alpar@9 85 }
alpar@9 86
alpar@9 87 /***********************************************************************
alpar@9 88 * NAME
alpar@9 89 *
alpar@9 90 * str2num - convert character string to value of double type
alpar@9 91 *
alpar@9 92 * SYNOPSIS
alpar@9 93 *
alpar@9 94 * #include "glplib.h"
alpar@9 95 * int str2num(const char *str, double *val);
alpar@9 96 *
alpar@9 97 * DESCRIPTION
alpar@9 98 *
alpar@9 99 * The routine str2num converts the character string str to a value of
alpar@9 100 * double type and stores the value into location, which the parameter
alpar@9 101 * val points to (in the case of error content of this location is not
alpar@9 102 * changed).
alpar@9 103 *
alpar@9 104 * RETURNS
alpar@9 105 *
alpar@9 106 * The routine returns one of the following error codes:
alpar@9 107 *
alpar@9 108 * 0 - no error;
alpar@9 109 * 1 - value out of range;
alpar@9 110 * 2 - character string is syntactically incorrect. */
alpar@9 111
alpar@9 112 int str2num(const char *str, double *_val)
alpar@9 113 { int k;
alpar@9 114 double val;
alpar@9 115 /* scan optional sign */
alpar@9 116 k = (str[0] == '+' || str[0] == '-' ? 1 : 0);
alpar@9 117 /* check for decimal point */
alpar@9 118 if (str[k] == '.')
alpar@9 119 { k++;
alpar@9 120 /* a digit should follow it */
alpar@9 121 if (!isdigit((unsigned char)str[k])) return 2;
alpar@9 122 k++;
alpar@9 123 goto frac;
alpar@9 124 }
alpar@9 125 /* integer part should start with a digit */
alpar@9 126 if (!isdigit((unsigned char)str[k])) return 2;
alpar@9 127 /* scan integer part */
alpar@9 128 while (isdigit((unsigned char)str[k])) k++;
alpar@9 129 /* check for decimal point */
alpar@9 130 if (str[k] == '.') k++;
alpar@9 131 frac: /* scan optional fraction part */
alpar@9 132 while (isdigit((unsigned char)str[k])) k++;
alpar@9 133 /* check for decimal exponent */
alpar@9 134 if (str[k] == 'E' || str[k] == 'e')
alpar@9 135 { k++;
alpar@9 136 /* scan optional sign */
alpar@9 137 if (str[k] == '+' || str[k] == '-') k++;
alpar@9 138 /* a digit should follow E, E+ or E- */
alpar@9 139 if (!isdigit((unsigned char)str[k])) return 2;
alpar@9 140 }
alpar@9 141 /* scan optional exponent part */
alpar@9 142 while (isdigit((unsigned char)str[k])) k++;
alpar@9 143 /* check for terminator */
alpar@9 144 if (str[k] != '\0') return 2;
alpar@9 145 /* perform conversion */
alpar@9 146 { char *endptr;
alpar@9 147 val = strtod(str, &endptr);
alpar@9 148 if (*endptr != '\0') return 2;
alpar@9 149 }
alpar@9 150 /* check for overflow */
alpar@9 151 if (!(-DBL_MAX <= val && val <= +DBL_MAX)) return 1;
alpar@9 152 /* check for underflow */
alpar@9 153 if (-DBL_MIN < val && val < +DBL_MIN) val = 0.0;
alpar@9 154 /* conversion has been done */
alpar@9 155 *_val = val;
alpar@9 156 return 0;
alpar@9 157 }
alpar@9 158
alpar@9 159 /***********************************************************************
alpar@9 160 * NAME
alpar@9 161 *
alpar@9 162 * strspx - remove all spaces from character string
alpar@9 163 *
alpar@9 164 * SYNOPSIS
alpar@9 165 *
alpar@9 166 * #include "glplib.h"
alpar@9 167 * char *strspx(char *str);
alpar@9 168 *
alpar@9 169 * DESCRIPTION
alpar@9 170 *
alpar@9 171 * The routine strspx removes all spaces from the character string str.
alpar@9 172 *
alpar@9 173 * RETURNS
alpar@9 174 *
alpar@9 175 * The routine returns a pointer to the character string.
alpar@9 176 *
alpar@9 177 * EXAMPLES
alpar@9 178 *
alpar@9 179 * strspx(" Errare humanum est ") => "Errarehumanumest"
alpar@9 180 *
alpar@9 181 * strspx(" ") => "" */
alpar@9 182
alpar@9 183 char *strspx(char *str)
alpar@9 184 { char *s, *t;
alpar@9 185 for (s = t = str; *s; s++) if (*s != ' ') *t++ = *s;
alpar@9 186 *t = '\0';
alpar@9 187 return str;
alpar@9 188 }
alpar@9 189
alpar@9 190 /***********************************************************************
alpar@9 191 * NAME
alpar@9 192 *
alpar@9 193 * strtrim - remove trailing spaces from character string
alpar@9 194 *
alpar@9 195 * SYNOPSIS
alpar@9 196 *
alpar@9 197 * #include "glplib.h"
alpar@9 198 * char *strtrim(char *str);
alpar@9 199 *
alpar@9 200 * DESCRIPTION
alpar@9 201 *
alpar@9 202 * The routine strtrim removes trailing spaces from the character
alpar@9 203 * string str.
alpar@9 204 *
alpar@9 205 * RETURNS
alpar@9 206 *
alpar@9 207 * The routine returns a pointer to the character string.
alpar@9 208 *
alpar@9 209 * EXAMPLES
alpar@9 210 *
alpar@9 211 * strtrim("Errare humanum est ") => "Errare humanum est"
alpar@9 212 *
alpar@9 213 * strtrim(" ") => "" */
alpar@9 214
alpar@9 215 char *strtrim(char *str)
alpar@9 216 { char *t;
alpar@9 217 for (t = strrchr(str, '\0') - 1; t >= str; t--)
alpar@9 218 { if (*t != ' ') break;
alpar@9 219 *t = '\0';
alpar@9 220 }
alpar@9 221 return str;
alpar@9 222 }
alpar@9 223
alpar@9 224 /***********************************************************************
alpar@9 225 * NAME
alpar@9 226 *
alpar@9 227 * strrev - reverse character string
alpar@9 228 *
alpar@9 229 * SYNOPSIS
alpar@9 230 *
alpar@9 231 * #include "glplib.h"
alpar@9 232 * char *strrev(char *s);
alpar@9 233 *
alpar@9 234 * DESCRIPTION
alpar@9 235 *
alpar@9 236 * The routine strrev changes characters in a character string s to the
alpar@9 237 * reverse order, except the terminating null character.
alpar@9 238 *
alpar@9 239 * RETURNS
alpar@9 240 *
alpar@9 241 * The routine returns the pointer s.
alpar@9 242 *
alpar@9 243 * EXAMPLES
alpar@9 244 *
alpar@9 245 * strrev("") => ""
alpar@9 246 *
alpar@9 247 * strrev("Today is Monday") => "yadnoM si yadoT" */
alpar@9 248
alpar@9 249 char *strrev(char *s)
alpar@9 250 { int i, j;
alpar@9 251 char t;
alpar@9 252 for (i = 0, j = strlen(s)-1; i < j; i++, j--)
alpar@9 253 t = s[i], s[i] = s[j], s[j] = t;
alpar@9 254 return s;
alpar@9 255 }
alpar@9 256
alpar@9 257 /***********************************************************************
alpar@9 258 * NAME
alpar@9 259 *
alpar@9 260 * gcd - find greatest common divisor of two integers
alpar@9 261 *
alpar@9 262 * SYNOPSIS
alpar@9 263 *
alpar@9 264 * #include "glplib.h"
alpar@9 265 * int gcd(int x, int y);
alpar@9 266 *
alpar@9 267 * RETURNS
alpar@9 268 *
alpar@9 269 * The routine gcd returns gcd(x, y), the greatest common divisor of
alpar@9 270 * the two positive integers given.
alpar@9 271 *
alpar@9 272 * ALGORITHM
alpar@9 273 *
alpar@9 274 * The routine gcd is based on Euclid's algorithm.
alpar@9 275 *
alpar@9 276 * REFERENCES
alpar@9 277 *
alpar@9 278 * Don Knuth, The Art of Computer Programming, Vol.2: Seminumerical
alpar@9 279 * Algorithms, 3rd Edition, Addison-Wesley, 1997. Section 4.5.2: The
alpar@9 280 * Greatest Common Divisor, pp. 333-56. */
alpar@9 281
alpar@9 282 int gcd(int x, int y)
alpar@9 283 { int r;
alpar@9 284 xassert(x > 0 && y > 0);
alpar@9 285 while (y > 0)
alpar@9 286 r = x % y, x = y, y = r;
alpar@9 287 return x;
alpar@9 288 }
alpar@9 289
alpar@9 290 /***********************************************************************
alpar@9 291 * NAME
alpar@9 292 *
alpar@9 293 * gcdn - find greatest common divisor of n integers
alpar@9 294 *
alpar@9 295 * SYNOPSIS
alpar@9 296 *
alpar@9 297 * #include "glplib.h"
alpar@9 298 * int gcdn(int n, int x[]);
alpar@9 299 *
alpar@9 300 * RETURNS
alpar@9 301 *
alpar@9 302 * The routine gcdn returns gcd(x[1], x[2], ..., x[n]), the greatest
alpar@9 303 * common divisor of n positive integers given, n > 0.
alpar@9 304 *
alpar@9 305 * BACKGROUND
alpar@9 306 *
alpar@9 307 * The routine gcdn is based on the following identity:
alpar@9 308 *
alpar@9 309 * gcd(x, y, z) = gcd(gcd(x, y), z).
alpar@9 310 *
alpar@9 311 * REFERENCES
alpar@9 312 *
alpar@9 313 * Don Knuth, The Art of Computer Programming, Vol.2: Seminumerical
alpar@9 314 * Algorithms, 3rd Edition, Addison-Wesley, 1997. Section 4.5.2: The
alpar@9 315 * Greatest Common Divisor, pp. 333-56. */
alpar@9 316
alpar@9 317 int gcdn(int n, int x[])
alpar@9 318 { int d, j;
alpar@9 319 xassert(n > 0);
alpar@9 320 for (j = 1; j <= n; j++)
alpar@9 321 { xassert(x[j] > 0);
alpar@9 322 if (j == 1)
alpar@9 323 d = x[1];
alpar@9 324 else
alpar@9 325 d = gcd(d, x[j]);
alpar@9 326 if (d == 1) break;
alpar@9 327 }
alpar@9 328 return d;
alpar@9 329 }
alpar@9 330
alpar@9 331 /***********************************************************************
alpar@9 332 * NAME
alpar@9 333 *
alpar@9 334 * lcm - find least common multiple of two integers
alpar@9 335 *
alpar@9 336 * SYNOPSIS
alpar@9 337 *
alpar@9 338 * #include "glplib.h"
alpar@9 339 * int lcm(int x, int y);
alpar@9 340 *
alpar@9 341 * RETURNS
alpar@9 342 *
alpar@9 343 * The routine lcm returns lcm(x, y), the least common multiple of the
alpar@9 344 * two positive integers given. In case of integer overflow the routine
alpar@9 345 * returns zero.
alpar@9 346 *
alpar@9 347 * BACKGROUND
alpar@9 348 *
alpar@9 349 * The routine lcm is based on the following identity:
alpar@9 350 *
alpar@9 351 * lcm(x, y) = (x * y) / gcd(x, y) = x * [y / gcd(x, y)],
alpar@9 352 *
alpar@9 353 * where gcd(x, y) is the greatest common divisor of x and y. */
alpar@9 354
alpar@9 355 int lcm(int x, int y)
alpar@9 356 { xassert(x > 0);
alpar@9 357 xassert(y > 0);
alpar@9 358 y /= gcd(x, y);
alpar@9 359 if (x > INT_MAX / y) return 0;
alpar@9 360 return x * y;
alpar@9 361 }
alpar@9 362
alpar@9 363 /***********************************************************************
alpar@9 364 * NAME
alpar@9 365 *
alpar@9 366 * lcmn - find least common multiple of n integers
alpar@9 367 *
alpar@9 368 * SYNOPSIS
alpar@9 369 *
alpar@9 370 * #include "glplib.h"
alpar@9 371 * int lcmn(int n, int x[]);
alpar@9 372 *
alpar@9 373 * RETURNS
alpar@9 374 *
alpar@9 375 * The routine lcmn returns lcm(x[1], x[2], ..., x[n]), the least
alpar@9 376 * common multiple of n positive integers given, n > 0. In case of
alpar@9 377 * integer overflow the routine returns zero.
alpar@9 378 *
alpar@9 379 * BACKGROUND
alpar@9 380 *
alpar@9 381 * The routine lcmn is based on the following identity:
alpar@9 382 *
alpar@9 383 * lcmn(x, y, z) = lcm(lcm(x, y), z),
alpar@9 384 *
alpar@9 385 * where lcm(x, y) is the least common multiple of x and y. */
alpar@9 386
alpar@9 387 int lcmn(int n, int x[])
alpar@9 388 { int m, j;
alpar@9 389 xassert(n > 0);
alpar@9 390 for (j = 1; j <= n; j++)
alpar@9 391 { xassert(x[j] > 0);
alpar@9 392 if (j == 1)
alpar@9 393 m = x[1];
alpar@9 394 else
alpar@9 395 m = lcm(m, x[j]);
alpar@9 396 if (m == 0) break;
alpar@9 397 }
alpar@9 398 return m;
alpar@9 399 }
alpar@9 400
alpar@9 401 /***********************************************************************
alpar@9 402 * NAME
alpar@9 403 *
alpar@9 404 * round2n - round floating-point number to nearest power of two
alpar@9 405 *
alpar@9 406 * SYNOPSIS
alpar@9 407 *
alpar@9 408 * #include "glplib.h"
alpar@9 409 * double round2n(double x);
alpar@9 410 *
alpar@9 411 * RETURNS
alpar@9 412 *
alpar@9 413 * Given a positive floating-point value x the routine round2n returns
alpar@9 414 * 2^n such that |x - 2^n| is minimal.
alpar@9 415 *
alpar@9 416 * EXAMPLES
alpar@9 417 *
alpar@9 418 * round2n(10.1) = 2^3 = 8
alpar@9 419 * round2n(15.3) = 2^4 = 16
alpar@9 420 * round2n(0.01) = 2^(-7) = 0.0078125
alpar@9 421 *
alpar@9 422 * BACKGROUND
alpar@9 423 *
alpar@9 424 * Let x = f * 2^e, where 0.5 <= f < 1 is a normalized fractional part,
alpar@9 425 * e is an integer exponent. Then, obviously, 0.5 * 2^e <= x < 2^e, so
alpar@9 426 * if x - 0.5 * 2^e <= 2^e - x, we choose 0.5 * 2^e = 2^(e-1), and 2^e
alpar@9 427 * otherwise. The latter condition can be written as 2 * x <= 1.5 * 2^e
alpar@9 428 * or 2 * f * 2^e <= 1.5 * 2^e or, finally, f <= 0.75. */
alpar@9 429
alpar@9 430 double round2n(double x)
alpar@9 431 { int e;
alpar@9 432 double f;
alpar@9 433 xassert(x > 0.0);
alpar@9 434 f = frexp(x, &e);
alpar@9 435 return ldexp(1.0, f <= 0.75 ? e-1 : e);
alpar@9 436 }
alpar@9 437
alpar@9 438 /***********************************************************************
alpar@9 439 * NAME
alpar@9 440 *
alpar@9 441 * fp2rat - convert floating-point number to rational number
alpar@9 442 *
alpar@9 443 * SYNOPSIS
alpar@9 444 *
alpar@9 445 * #include "glplib.h"
alpar@9 446 * int fp2rat(double x, double eps, double *p, double *q);
alpar@9 447 *
alpar@9 448 * DESCRIPTION
alpar@9 449 *
alpar@9 450 * Given a floating-point number 0 <= x < 1 the routine fp2rat finds
alpar@9 451 * its "best" rational approximation p / q, where p >= 0 and q > 0 are
alpar@9 452 * integer numbers, such that |x - p / q| <= eps.
alpar@9 453 *
alpar@9 454 * RETURNS
alpar@9 455 *
alpar@9 456 * The routine fp2rat returns the number of iterations used to achieve
alpar@9 457 * the specified precision eps.
alpar@9 458 *
alpar@9 459 * EXAMPLES
alpar@9 460 *
alpar@9 461 * For x = sqrt(2) - 1 = 0.414213562373095 and eps = 1e-6 the routine
alpar@9 462 * gives p = 408 and q = 985, where 408 / 985 = 0.414213197969543.
alpar@9 463 *
alpar@9 464 * BACKGROUND
alpar@9 465 *
alpar@9 466 * It is well known that every positive real number x can be expressed
alpar@9 467 * as the following continued fraction:
alpar@9 468 *
alpar@9 469 * x = b[0] + a[1]
alpar@9 470 * ------------------------
alpar@9 471 * b[1] + a[2]
alpar@9 472 * -----------------
alpar@9 473 * b[2] + a[3]
alpar@9 474 * ----------
alpar@9 475 * b[3] + ...
alpar@9 476 *
alpar@9 477 * where:
alpar@9 478 *
alpar@9 479 * a[k] = 1, k = 0, 1, 2, ...
alpar@9 480 *
alpar@9 481 * b[k] = floor(x[k]), k = 0, 1, 2, ...
alpar@9 482 *
alpar@9 483 * x[0] = x,
alpar@9 484 *
alpar@9 485 * x[k] = 1 / frac(x[k-1]), k = 1, 2, 3, ...
alpar@9 486 *
alpar@9 487 * To find the "best" rational approximation of x the routine computes
alpar@9 488 * partial fractions f[k] by dropping after k terms as follows:
alpar@9 489 *
alpar@9 490 * f[k] = A[k] / B[k],
alpar@9 491 *
alpar@9 492 * where:
alpar@9 493 *
alpar@9 494 * A[-1] = 1, A[0] = b[0], B[-1] = 0, B[0] = 1,
alpar@9 495 *
alpar@9 496 * A[k] = b[k] * A[k-1] + a[k] * A[k-2],
alpar@9 497 *
alpar@9 498 * B[k] = b[k] * B[k-1] + a[k] * B[k-2].
alpar@9 499 *
alpar@9 500 * Once the condition
alpar@9 501 *
alpar@9 502 * |x - f[k]| <= eps
alpar@9 503 *
alpar@9 504 * has been satisfied, the routine reports p = A[k] and q = B[k] as the
alpar@9 505 * final answer.
alpar@9 506 *
alpar@9 507 * In the table below here is some statistics obtained for one million
alpar@9 508 * random numbers uniformly distributed in the range [0, 1).
alpar@9 509 *
alpar@9 510 * eps max p mean p max q mean q max k mean k
alpar@9 511 * -------------------------------------------------------------
alpar@9 512 * 1e-1 8 1.6 9 3.2 3 1.4
alpar@9 513 * 1e-2 98 6.2 99 12.4 5 2.4
alpar@9 514 * 1e-3 997 20.7 998 41.5 8 3.4
alpar@9 515 * 1e-4 9959 66.6 9960 133.5 10 4.4
alpar@9 516 * 1e-5 97403 211.7 97404 424.2 13 5.3
alpar@9 517 * 1e-6 479669 669.9 479670 1342.9 15 6.3
alpar@9 518 * 1e-7 1579030 2127.3 3962146 4257.8 16 7.3
alpar@9 519 * 1e-8 26188823 6749.4 26188824 13503.4 19 8.2
alpar@9 520 *
alpar@9 521 * REFERENCES
alpar@9 522 *
alpar@9 523 * W. B. Jones and W. J. Thron, "Continued Fractions: Analytic Theory
alpar@9 524 * and Applications," Encyclopedia on Mathematics and Its Applications,
alpar@9 525 * Addison-Wesley, 1980. */
alpar@9 526
alpar@9 527 int fp2rat(double x, double eps, double *p, double *q)
alpar@9 528 { int k;
alpar@9 529 double xk, Akm1, Ak, Bkm1, Bk, ak, bk, fk, temp;
alpar@9 530 if (!(0.0 <= x && x < 1.0))
alpar@9 531 xerror("fp2rat: x = %g; number out of range\n", x);
alpar@9 532 for (k = 0; ; k++)
alpar@9 533 { xassert(k <= 100);
alpar@9 534 if (k == 0)
alpar@9 535 { /* x[0] = x */
alpar@9 536 xk = x;
alpar@9 537 /* A[-1] = 1 */
alpar@9 538 Akm1 = 1.0;
alpar@9 539 /* A[0] = b[0] = floor(x[0]) = 0 */
alpar@9 540 Ak = 0.0;
alpar@9 541 /* B[-1] = 0 */
alpar@9 542 Bkm1 = 0.0;
alpar@9 543 /* B[0] = 1 */
alpar@9 544 Bk = 1.0;
alpar@9 545 }
alpar@9 546 else
alpar@9 547 { /* x[k] = 1 / frac(x[k-1]) */
alpar@9 548 temp = xk - floor(xk);
alpar@9 549 xassert(temp != 0.0);
alpar@9 550 xk = 1.0 / temp;
alpar@9 551 /* a[k] = 1 */
alpar@9 552 ak = 1.0;
alpar@9 553 /* b[k] = floor(x[k]) */
alpar@9 554 bk = floor(xk);
alpar@9 555 /* A[k] = b[k] * A[k-1] + a[k] * A[k-2] */
alpar@9 556 temp = bk * Ak + ak * Akm1;
alpar@9 557 Akm1 = Ak, Ak = temp;
alpar@9 558 /* B[k] = b[k] * B[k-1] + a[k] * B[k-2] */
alpar@9 559 temp = bk * Bk + ak * Bkm1;
alpar@9 560 Bkm1 = Bk, Bk = temp;
alpar@9 561 }
alpar@9 562 /* f[k] = A[k] / B[k] */
alpar@9 563 fk = Ak / Bk;
alpar@9 564 #if 0
alpar@9 565 print("%.*g / %.*g = %.*g", DBL_DIG, Ak, DBL_DIG, Bk, DBL_DIG,
alpar@9 566 fk);
alpar@9 567 #endif
alpar@9 568 if (fabs(x - fk) <= eps) break;
alpar@9 569 }
alpar@9 570 *p = Ak;
alpar@9 571 *q = Bk;
alpar@9 572 return k;
alpar@9 573 }
alpar@9 574
alpar@9 575 /***********************************************************************
alpar@9 576 * NAME
alpar@9 577 *
alpar@9 578 * jday - convert calendar date to Julian day number
alpar@9 579 *
alpar@9 580 * SYNOPSIS
alpar@9 581 *
alpar@9 582 * #include "glplib.h"
alpar@9 583 * int jday(int d, int m, int y);
alpar@9 584 *
alpar@9 585 * DESCRIPTION
alpar@9 586 *
alpar@9 587 * The routine jday converts a calendar date, Gregorian calendar, to
alpar@9 588 * corresponding Julian day number j.
alpar@9 589 *
alpar@9 590 * From the given day d, month m, and year y, the Julian day number j
alpar@9 591 * is computed without using tables.
alpar@9 592 *
alpar@9 593 * The routine is valid for 1 <= y <= 4000.
alpar@9 594 *
alpar@9 595 * RETURNS
alpar@9 596 *
alpar@9 597 * The routine jday returns the Julian day number, or negative value if
alpar@9 598 * the specified date is incorrect.
alpar@9 599 *
alpar@9 600 * REFERENCES
alpar@9 601 *
alpar@9 602 * R. G. Tantzen, Algorithm 199: conversions between calendar date and
alpar@9 603 * Julian day number, Communications of the ACM, vol. 6, no. 8, p. 444,
alpar@9 604 * Aug. 1963. */
alpar@9 605
alpar@9 606 int jday(int d, int m, int y)
alpar@9 607 { int c, ya, j, dd;
alpar@9 608 if (!(1 <= d && d <= 31 && 1 <= m && m <= 12 && 1 <= y &&
alpar@9 609 y <= 4000))
alpar@9 610 { j = -1;
alpar@9 611 goto done;
alpar@9 612 }
alpar@9 613 if (m >= 3) m -= 3; else m += 9, y--;
alpar@9 614 c = y / 100;
alpar@9 615 ya = y - 100 * c;
alpar@9 616 j = (146097 * c) / 4 + (1461 * ya) / 4 + (153 * m + 2) / 5 + d +
alpar@9 617 1721119;
alpar@9 618 jdate(j, &dd, NULL, NULL);
alpar@9 619 if (d != dd) j = -1;
alpar@9 620 done: return j;
alpar@9 621 }
alpar@9 622
alpar@9 623 /***********************************************************************
alpar@9 624 * NAME
alpar@9 625 *
alpar@9 626 * jdate - convert Julian day number to calendar date
alpar@9 627 *
alpar@9 628 * SYNOPSIS
alpar@9 629 *
alpar@9 630 * #include "glplib.h"
alpar@9 631 * void jdate(int j, int *d, int *m, int *y);
alpar@9 632 *
alpar@9 633 * DESCRIPTION
alpar@9 634 *
alpar@9 635 * The routine jdate converts a Julian day number j to corresponding
alpar@9 636 * calendar date, Gregorian calendar.
alpar@9 637 *
alpar@9 638 * The day d, month m, and year y are computed without using tables and
alpar@9 639 * stored in corresponding locations.
alpar@9 640 *
alpar@9 641 * The routine is valid for 1721426 <= j <= 3182395.
alpar@9 642 *
alpar@9 643 * RETURNS
alpar@9 644 *
alpar@9 645 * If the conversion is successful, the routine returns zero, otherwise
alpar@9 646 * non-zero.
alpar@9 647 *
alpar@9 648 * REFERENCES
alpar@9 649 *
alpar@9 650 * R. G. Tantzen, Algorithm 199: conversions between calendar date and
alpar@9 651 * Julian day number, Communications of the ACM, vol. 6, no. 8, p. 444,
alpar@9 652 * Aug. 1963. */
alpar@9 653
alpar@9 654 int jdate(int j, int *_d, int *_m, int *_y)
alpar@9 655 { int d, m, y, ret = 0;
alpar@9 656 if (!(1721426 <= j && j <= 3182395))
alpar@9 657 { ret = 1;
alpar@9 658 goto done;
alpar@9 659 }
alpar@9 660 j -= 1721119;
alpar@9 661 y = (4 * j - 1) / 146097;
alpar@9 662 j = (4 * j - 1) % 146097;
alpar@9 663 d = j / 4;
alpar@9 664 j = (4 * d + 3) / 1461;
alpar@9 665 d = (4 * d + 3) % 1461;
alpar@9 666 d = (d + 4) / 4;
alpar@9 667 m = (5 * d - 3) / 153;
alpar@9 668 d = (5 * d - 3) % 153;
alpar@9 669 d = (d + 5) / 5;
alpar@9 670 y = 100 * y + j;
alpar@9 671 if (m <= 9) m += 3; else m -= 9, y++;
alpar@9 672 if (_d != NULL) *_d = d;
alpar@9 673 if (_m != NULL) *_m = m;
alpar@9 674 if (_y != NULL) *_y = y;
alpar@9 675 done: return ret;
alpar@9 676 }
alpar@9 677
alpar@9 678 #if 0
alpar@9 679 int main(void)
alpar@9 680 { int jbeg, jend, j, d, m, y;
alpar@9 681 jbeg = jday(1, 1, 1);
alpar@9 682 jend = jday(31, 12, 4000);
alpar@9 683 for (j = jbeg; j <= jend; j++)
alpar@9 684 { xassert(jdate(j, &d, &m, &y) == 0);
alpar@9 685 xassert(jday(d, m, y) == j);
alpar@9 686 }
alpar@9 687 xprintf("Routines jday and jdate work correctly.\n");
alpar@9 688 return 0;
alpar@9 689 }
alpar@9 690 #endif
alpar@9 691
alpar@9 692 /* eof */