lemon-project-template-glpk
view deps/glpk/src/glplib03.c @ 9:33de93886c88
Import GLPK 4.47
| author | Alpar Juttner <alpar@cs.elte.hu> |
|---|---|
| date | Sun, 06 Nov 2011 20:59:10 +0100 |
| parents | |
| children |
line source
1 /* glplib03.c (miscellaneous library routines) */
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, 2011 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 ***********************************************************************/
25 #include "glpenv.h"
26 #include "glplib.h"
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. */
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 }
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. */
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 }
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(" ") => "" */
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 }
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(" ") => "" */
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 }
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" */
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 }
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. */
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 }
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. */
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 }
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. */
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 }
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. */
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 }
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. */
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 }
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. */
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 }
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. */
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 }
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. */
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 }
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
692 /* eof */