glplib11.c

著名的大规模线性规划求解器源码GLPK.C语言版本,可以修剪.内有详细帮助文档.

*  Don Knuth, The Art of Computer Programming, Vol.2: Seminumerical
*  Algorithms, 3rd Edition, Addison-Wesley, 1997. Section 4.5.2: The
*  Greatest Common Divisor, pp. 333-56. */

int gcdn(int n, int x[])
{     int d, j;
      xassert(n > 0);
      for (j = 1; j <= n; j++)
      {  xassert(x[j] >= 0);
         if (j == 1)
            d = x[1];
         else
            d = gcd(d, x[j]);
         if (d == 1) break;
      }
      return d;
}

/***********************************************************************
*  NAME
*
*  round2n - round floating-point number to nearest power of two
*
*  SYNOPSIS
*
*  #include "glplib.h"
*  double round2n(double x);
*
*  RETURNS
*
*  Given a positive floating-point value x the routine round2n returns
*  2^n such that |x - 2^n| is minimal.
*
*  EXAMPLES
*
*  round2n(10.1) = 2^3 = 8
*  round2n(15.3) = 2^4 = 16
*  round2n(0.01) = 2^(-7) = 0.0078125
*
*  BACKGROUND
*
*  Let x = f * 2^e, where 0.5 <= f < 1 is a normalized fractional part,
*  e is an integer exponent. Then, obviously, 0.5 * 2^e <= x < 2^e, so
*  if x - 0.5 * 2^e <= 2^e - x, we choose 0.5 * 2^e = 2^(e-1), and 2^e
*  otherwise. The latter condition can be written as 2 * x <= 1.5 * 2^e
*  or 2 * f * 2^e <= 1.5 * 2^e or, finally, f <= 0.75. */

double round2n(double x)
{     int e;
      double f;
      xassert(x > 0.0);
      f = frexp(x, &e);
      return ldexp(1.0, f <= 0.75 ? e-1 : e);
}

/***********************************************************************
*  NAME
*
*  fp2rat - convert floating-point number to rational number
*
*  SYNOPSIS
*
*  #include "glplib.h"
*  int fp2rat(double x, double eps, double *p, double *q);
*
*  DESCRIPTION
*
*  Given a floating-point number 0 <= x < 1 the routine fp2rat finds
*  its "best" rational approximation p / q, where p >= 0 and q > 0 are
*  integer numbers, such that |x - p / q| <= eps.
*
*  RETURNS
*
*  The routine fp2rat returns the number of iterations used to achieve
*  the specified precision eps.
*
*  EXAMPLES
*
*  For x = sqrt(2) - 1 = 0.414213562373095 and eps = 1e-6 the routine
*  gives p = 408 and q = 985, where 408 / 985 = 0.414213197969543.
*
*  BACKGROUND
*
*  It is well known that every positive real number x can be expressed
*  as the following continued fraction:
*
*     x = b[0] + a[1]
*                ------------------------
*                b[1] + a[2]
*                       -----------------
*                       b[2] + a[3]
*                              ----------
*                              b[3] + ...
*
*  where:
*
*     a[k] = 1,                  k = 0, 1, 2, ...
*
*     b[k] = floor(x[k]),        k = 0, 1, 2, ...
*
*     x[0] = x,
*
*     x[k] = 1 / frac(x[k-1]),   k = 1, 2, 3, ...
*
*  To find the "best" rational approximation of x the routine computes
*  partial fractions f[k] by dropping after k terms as follows:
*
*     f[k] = A[k] / B[k],
*
*  where:
*
*     A[-1] = 1,   A[0] = b[0],   B[-1] = 0,   B[0] = 1,
*
*     A[k] = b[k] * A[k-1] + a[k] * A[k-2],
*
*     B[k] = b[k] * B[k-1] + a[k] * B[k-2].
*
*  Once the condition
*
*     |x - f[k]| <= eps
*
*  has been satisfied, the routine reports p = A[k] and q = B[k] as the
*  final answer.
*
*  In the table below here is some statistics obtained for one million
*  random numbers uniformly distributed in the range [0, 1).
*
*      eps      max p   mean p      max q    mean q  max k   mean k
*     -------------------------------------------------------------
*     1e-1          8      1.6          9       3.2    3      1.4
*     1e-2         98      6.2         99      12.4    5      2.4
*     1e-3        997     20.7        998      41.5    8      3.4
*     1e-4       9959     66.6       9960     133.5   10      4.4
*     1e-5      97403    211.7      97404     424.2   13      5.3
*     1e-6     479669    669.9     479670    1342.9   15      6.3
*     1e-7    1579030   2127.3    3962146    4257.8   16      7.3
*     1e-8   26188823   6749.4   26188824   13503.4   19      8.2
*
*  REFERENCES
*
*  W. B. Jones and W. J. Thron, "Continued Fractions: Analytic Theory
*  and Applications," Encyclopedia on Mathematics and Its Applications,
*  Addison-Wesley, 1980. */

int fp2rat(double x, double eps, double *p, double *q)
{     int k;
      double xk, Akm1, Ak, Bkm1, Bk, ak, bk, fk, temp;
      if (!(0.0 <= x && x < 1.0))
         xerror("fp2rat: x = %g; number out of range\n", x);
      for (k = 0; ; k++)
      {  xassert(k <= 100);
         if (k == 0)
         {  /* x[0] = x */
            xk = x;
            /* A[-1] = 1 */
            Akm1 = 1.0;
            /* A[0] = b[0] = floor(x[0]) = 0 */
            Ak = 0.0;
            /* B[-1] = 0 */
            Bkm1 = 0.0;
            /* B[0] = 1 */
            Bk = 1.0;
         }
         else
         {  /* x[k] = 1 / frac(x[k-1]) */
            temp = xk - floor(xk);
            xassert(temp != 0.0);
            xk = 1.0 / temp;
            /* a[k] = 1 */
            ak = 1.0;
            /* b[k] = floor(x[k]) */
            bk = floor(xk);
            /* A[k] = b[k] * A[k-1] + a[k] * A[k-2] */
            temp = bk * Ak + ak * Akm1;
            Akm1 = Ak, Ak = temp;
            /* B[k] = b[k] * B[k-1] + a[k] * B[k-2] */
            temp = bk * Bk + ak * Bkm1;
            Bkm1 = Bk, Bk = temp;
         }
         /* f[k] = A[k] / B[k] */
         fk = Ak / Bk;
#if 0
         print("%.*g / %.*g = %.*g", DBL_DIG, Ak, DBL_DIG, Bk, DBL_DIG,
            fk);
#endif
         if (fabs(x - fk) <= eps) break;
      }
      *p = Ak;
      *q = Bk;
      return k;
}

/***********************************************************************
*  NAME
*
*  jday - convert calendar date to Julian day number
*
*  SYNOPSIS
*
*  #include "glplib.h"
*  int jday(int d, int m, int y);
*
*  DESCRIPTION
*
*  The routine jday converts a calendar date, Gregorian calendar, to
*  corresponding Julian day number j.
*
*  From the given day d, month m, and year y, the Julian day number j
*  is computed without using tables.
*
*  The routine is valid for 1 <= y <= 4000.
*
*  RETURNS
*
*  The routine jday returns the Julian day number, or negative value if
*  the specified date is incorrect.
*
*  REFERENCES
*
*  R. G. Tantzen, Algorithm 199: conversions between calendar date and
*  Julian day number, Communications of the ACM, vol. 6, no. 8, p. 444,
*  Aug. 1963. */

int jday(int d, int m, int y)
{     int c, ya, j, dd;
      if (!(1 <= d && d <= 31 && 1 <= m && m <= 12 && 1 <= y &&
            y <= 4000))
      {  j = -1;
         goto done;
      }
      if (m >= 3) m -= 3; else m += 9, y--;
      c = y / 100;
      ya = y - 100 * c;
      j = (146097 * c) / 4 + (1461 * ya) / 4 + (153 * m + 2) / 5 + d +
         1721119;
      jdate(j, &dd, NULL, NULL);
      if (d != dd) j = -1;
done: return j;
}

/***********************************************************************
*  NAME
*
*  jdate - convert Julian day number to calendar date
*
*  SYNOPSIS
*
*  #include "glplib.h"
*  void jdate(int j, int *d, int *m, int *y);
*
*  DESCRIPTION
*
*  The routine jdate converts a Julian day number j to corresponding
*  calendar date, Gregorian calendar.
*
*  The day d, month m, and year y are computed without using tables and
*  stored in corresponding locations.
*
*  The routine is valid for 1721426 <= j <= 3182395.
*
*  RETURNS
*
*  If the conversion is successful, the routine returns zero, otherwise
*  non-zero.
*
*  REFERENCES
*
*  R. G. Tantzen, Algorithm 199: conversions between calendar date and
*  Julian day number, Communications of the ACM, vol. 6, no. 8, p. 444,
*  Aug. 1963. */

int jdate(int j, int *_d, int *_m, int *_y)
{     int d, m, y, ret = 0;
      if (!(1721426 <= j && j <= 3182395))
      {  ret = 1;
         goto done;
      }
      j -= 1721119;
      y = (4 * j - 1) / 146097;
      j = (4 * j - 1) % 146097;
      d = j / 4;
      j = (4 * d + 3) / 1461;
      d = (4 * d + 3) % 1461;
      d = (d + 4) / 4;
      m = (5 * d - 3) / 153;
      d = (5 * d - 3) % 153;
      d = (d + 5) / 5;
      y = 100 * y + j;
      if (m <= 9) m += 3; else m -= 9, y++;
      if (_d != NULL) *_d = d;
      if (_m != NULL) *_m = m;
      if (_y != NULL) *_y = y;
done: return ret;
}

#if 0
int main(void)
{     int jbeg, jend, j, d, m, y;
      jbeg = jday(1, 1, 1);
      jend = jday(31, 12, 4000);
      for (j = jbeg; j <= jend; j++)
      {  xassert(jdate(j, &d, &m, &y) == 0);
         xassert(jday(d, m, y) == j);
      }
      xprintf("Routines jday and jdate work correctly.\n");
      return 0;
}
#endif

/* eof */