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