COIN-OR::LEMON - Graph Library

source: glpk-cmake/src/glplib03.c @ 1:c445c931472f

Last change on this file since 1:c445c931472f was 1:c445c931472f, checked in by Alpar Juttner <alpar@…>, 10 years ago

Import glpk-4.45

  • Generated files and doc/notes are removed
File size: 17.7 KB
Line 
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
53int 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
112int 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++;
131frac: /* 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
183char *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
215char *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
249char *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
282int 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
317int 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
355int 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
387int 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
430double 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
527int 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
606int 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;
620done: 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
654int 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;
675done: return ret;
676}
677
678#if 0
679int 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 */
Note: See TracBrowser for help on using the repository browser.