glpmpl03.c

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

/* glpmpl03.c */

/***********************************************************************
*  This code is part of GLPK (GNU Linear Programming Kit).
*
*  Copyright (C) 2000,01,02,03,04,05,06,07,08,2009 Andrew Makhorin,
*  Department for Applied Informatics, Moscow Aviation Institute,
*  Moscow, Russia. All rights reserved. E-mail: <[e-mail]>.
*
*  GLPK is free software: you can redistribute it and/or modify it
*  under the terms of the GNU General Public License as published by
*  the Free Software Foundation, either version 3 of the License, or
*  (at your option) any later version.
*
*  GLPK is distributed in the hope that it will be useful, but WITHOUT
*  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
*  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
*  License for more details.
*
*  You should have received a copy of the GNU General Public License
*  along with GLPK. If not, see <http://www.gnu.org/licenses/>.
***********************************************************************/

#define _GLPSTD_ERRNO
#define _GLPSTD_STDIO
#include "glplib.h"
#include "glpmpl.h"

/**********************************************************************/
/* * *                   FLOATING-POINT NUMBERS                   * * */
/**********************************************************************/

/*----------------------------------------------------------------------
-- fp_add - floating-point addition.
--
-- This routine computes the sum x + y. */

double fp_add(MPL *mpl, double x, double y)
{     if (x > 0.0 && y > 0.0 && x > + 0.999 * DBL_MAX - y ||
          x < 0.0 && y < 0.0 && x < - 0.999 * DBL_MAX - y)
         error(mpl, "%.*g + %.*g; floating-point overflow",
            DBL_DIG, x, DBL_DIG, y);
      return x + y;
}

/*----------------------------------------------------------------------
-- fp_sub - floating-point subtraction.
--
-- This routine computes the difference x - y. */

double fp_sub(MPL *mpl, double x, double y)
{     if (x > 0.0 && y < 0.0 && x > + 0.999 * DBL_MAX + y ||
          x < 0.0 && y > 0.0 && x < - 0.999 * DBL_MAX + y)
         error(mpl, "%.*g - %.*g; floating-point overflow",
            DBL_DIG, x, DBL_DIG, y);
      return x - y;
}

/*----------------------------------------------------------------------
-- fp_less - floating-point non-negative subtraction.
--
-- This routine computes the non-negative difference max(0, x - y). */

double fp_less(MPL *mpl, double x, double y)
{     if (x < y) return 0.0;
      if (x > 0.0 && y < 0.0 && x > + 0.999 * DBL_MAX + y)
         error(mpl, "%.*g less %.*g; floating-point overflow",
            DBL_DIG, x, DBL_DIG, y);
      return x - y;
}

/*----------------------------------------------------------------------
-- fp_mul - floating-point multiplication.
--
-- This routine computes the product x * y. */

double fp_mul(MPL *mpl, double x, double y)
{     if (fabs(y) > 1.0 && fabs(x) > (0.999 * DBL_MAX) / fabs(y))
         error(mpl, "%.*g * %.*g; floating-point overflow",
            DBL_DIG, x, DBL_DIG, y);
      return x * y;
}

/*----------------------------------------------------------------------
-- fp_div - floating-point division.
--
-- This routine computes the quotient x / y. */

double fp_div(MPL *mpl, double x, double y)
{     if (fabs(y) < DBL_MIN)
         error(mpl, "%.*g / %.*g; floating-point zero divide",
            DBL_DIG, x, DBL_DIG, y);
      if (fabs(y) < 1.0 && fabs(x) > (0.999 * DBL_MAX) * fabs(y))
         error(mpl, "%.*g / %.*g; floating-point overflow",
            DBL_DIG, x, DBL_DIG, y);
      return x / y;
}

/*----------------------------------------------------------------------
-- fp_idiv - floating-point quotient of exact division.
--
-- This routine computes the quotient of exact division x div y. */

double fp_idiv(MPL *mpl, double x, double y)
{     if (fabs(y) < DBL_MIN)
         error(mpl, "%.*g div %.*g; floating-point zero divide",
            DBL_DIG, x, DBL_DIG, y);
      if (fabs(y) < 1.0 && fabs(x) > (0.999 * DBL_MAX) * fabs(y))
         error(mpl, "%.*g div %.*g; floating-point overflow",
            DBL_DIG, x, DBL_DIG, y);
      x /= y;
      return x > 0.0 ? floor(x) : x < 0.0 ? ceil(x) : 0.0;
}

/*----------------------------------------------------------------------
-- fp_mod - floating-point remainder of exact division.
--
-- This routine computes the remainder of exact division x mod y.
--
-- NOTE: By definition x mod y = x - y * floor(x / y). */

double fp_mod(MPL *mpl, double x, double y)
{     double r;
      xassert(mpl == mpl);
      if (x == 0.0)
         r = 0.0;
      else if (y == 0.0)
         r = x;
      else
      {  r = fmod(fabs(x), fabs(y));
         if (r != 0.0)
         {  if (x < 0.0) r = - r;
            if (x > 0.0 && y < 0.0 || x < 0.0 && y > 0.0) r += y;
         }
      }
      return r;
}

/*----------------------------------------------------------------------
-- fp_power - floating-point exponentiation (raise to power).
--
-- This routine computes the exponentiation x ** y. */

double fp_power(MPL *mpl, double x, double y)
{     double r;
      if (x == 0.0 && y <= 0.0 || x < 0.0 && y != floor(y))
         error(mpl, "%.*g ** %.*g; result undefined",
            DBL_DIG, x, DBL_DIG, y);
      if (x == 0.0) goto eval;
      if (fabs(x) > 1.0 && y > +1.0 &&
            +log(fabs(x)) > (0.999 * log(DBL_MAX)) / y ||
          fabs(x) < 1.0 && y < -1.0 &&
            +log(fabs(x)) < (0.999 * log(DBL_MAX)) / y)
         error(mpl, "%.*g ** %.*g; floating-point overflow",
            DBL_DIG, x, DBL_DIG, y);
      if (fabs(x) > 1.0 && y < -1.0 &&
            -log(fabs(x)) < (0.999 * log(DBL_MAX)) / y ||
          fabs(x) < 1.0 && y > +1.0 &&
            -log(fabs(x)) > (0.999 * log(DBL_MAX)) / y)
         r = 0.0;
      else
eval:    r = pow(x, y);
      return r;
}

/*----------------------------------------------------------------------
-- fp_exp - floating-point base-e exponential.
--
-- This routine computes the base-e exponential e ** x. */

double fp_exp(MPL *mpl, double x)
{     if (x > 0.999 * log(DBL_MAX))
         error(mpl, "exp(%.*g); floating-point overflow", DBL_DIG, x);
      return exp(x);
}

/*----------------------------------------------------------------------
-- fp_log - floating-point natural logarithm.
--
-- This routine computes the natural logarithm log x. */

double fp_log(MPL *mpl, double x)
{     if (x <= 0.0)
         error(mpl, "log(%.*g); non-positive argument", DBL_DIG, x);
      return log(x);
}

/*----------------------------------------------------------------------
-- fp_log10 - floating-point common (decimal) logarithm.
--
-- This routine computes the common (decimal) logarithm lg x. */

double fp_log10(MPL *mpl, double x)
{     if (x <= 0.0)
         error(mpl, "log10(%.*g); non-positive argument", DBL_DIG, x);
      return log10(x);
}

/*----------------------------------------------------------------------
-- fp_sqrt - floating-point square root.
--
-- This routine computes the square root x ** 0.5. */

double fp_sqrt(MPL *mpl, double x)
{     if (x < 0.0)
         error(mpl, "sqrt(%.*g); negative argument", DBL_DIG, x);
      return sqrt(x);
}

/*----------------------------------------------------------------------
-- fp_sin - floating-point trigonometric sine.
--
-- This routine computes the trigonometric sine sin(x). */

double fp_sin(MPL *mpl, double x)
{     if (!(-1e6 <= x && x <= +1e6))
         error(mpl, "sin(%.*g); argument too large", DBL_DIG, x);
      return sin(x);
}

/*----------------------------------------------------------------------
-- fp_cos - floating-point trigonometric cosine.
--
-- This routine computes the trigonometric cosine cos(x). */

double fp_cos(MPL *mpl, double x)
{     if (!(-1e6 <= x && x <= +1e6))
         error(mpl, "cos(%.*g); argument too large", DBL_DIG, x);
      return cos(x);
}

/*----------------------------------------------------------------------
-- fp_atan - floating-point trigonometric arctangent.
--
-- This routine computes the trigonometric arctangent atan(x). */

double fp_atan(MPL *mpl, double x)
{     xassert(mpl == mpl);
      return atan(x);
}

/*----------------------------------------------------------------------
-- fp_atan2 - floating-point trigonometric arctangent.
--
-- This routine computes the trigonometric arctangent atan(y / x). */

double fp_atan2(MPL *mpl, double y, double x)
{     xassert(mpl == mpl);
      return atan2(y, x);
}

/*----------------------------------------------------------------------
-- fp_round - round floating-point value to n fractional digits.
--
-- This routine rounds given floating-point value x to n fractional
-- digits with the formula:
--
--    round(x, n) = floor(x * 10^n + 0.5) / 10^n.
--
-- The parameter n is assumed to be integer. */

double fp_round(MPL *mpl, double x, double n)
{     double ten_to_n;
      if (n != floor(n))
         error(mpl, "round(%.*g, %.*g); non-integer second argument",
            DBL_DIG, x, DBL_DIG, n);
      if (n <= DBL_DIG + 2)
      {  ten_to_n = pow(10.0, n);
         if (fabs(x) < (0.999 * DBL_MAX) / ten_to_n)
         {  x = floor(x * ten_to_n + 0.5);
            if (x != 0.0) x /= ten_to_n;
         }
      }
      return x;
}

/*----------------------------------------------------------------------
-- fp_trunc - truncate floating-point value to n fractional digits.
--
-- This routine truncates given floating-point value x to n fractional
-- digits with the formula:
--
--                  ( floor(x * 10^n) / 10^n,  if x >= 0
--    trunc(x, n) = <
--                  ( ceil(x * 10^n) / 10^n,   if x < 0
--
-- The parameter n is assumed to be integer. */

double fp_trunc(MPL *mpl, double x, double n)
{     double ten_to_n;
      if (n != floor(n))
         error(mpl, "trunc(%.*g, %.*g); non-integer second argument",
            DBL_DIG, x, DBL_DIG, n);
      if (n <= DBL_DIG + 2)
      {  ten_to_n = pow(10.0, n);
         if (fabs(x) < (0.999 * DBL_MAX) / ten_to_n)
         {  x = (x >= 0.0 ? floor(x * ten_to_n) : ceil(x * ten_to_n));
            if (x != 0.0) x /= ten_to_n;
         }
      }
      return x;
}

/**********************************************************************/
/* * *              PSEUDO-RANDOM NUMBER GENERATORS               * * */
/**********************************************************************/

/*----------------------------------------------------------------------
-- fp_irand224 - pseudo-random integer in the range [0, 2^24).
--
-- This routine returns a next pseudo-random integer (converted to
-- floating-point) which is uniformly distributed between 0 and 2^24-1,
-- inclusive. */

#define two_to_the_24 0x1000000

double fp_irand224(MPL *mpl)
{     return
         (double)rng_unif_rand(mpl->rand, two_to_the_24);
}

/*----------------------------------------------------------------------
-- fp_uniform01 - pseudo-random number in the range [0, 1).
--
-- This routine returns a next pseudo-random number which is uniformly
-- distributed in the range [0, 1). */

#define two_to_the_31 ((unsigned int)0x80000000)

double fp_uniform01(MPL *mpl)
{     return
         (double)rng_next_rand(mpl->rand) / (double)two_to_the_31;
}

/*----------------------------------------------------------------------
-- fp_uniform - pseudo-random number in the range [a, b).
--
-- This routine returns a next pseudo-random number which is uniformly
-- distributed in the range [a, b). */

double fp_uniform(MPL *mpl, double a, double b)
{     double x;
      if (a >= b)
         error(mpl, "Uniform(%.*g, %.*g); invalid range",
            DBL_DIG, a, DBL_DIG, b);
      x = fp_uniform01(mpl);
#if 0
      x = a * (1.0 - x) + b * x;
#else
      x = fp_add(mpl, a * (1.0 - x), b * x);
#endif
      return x;
}

/*----------------------------------------------------------------------
-- fp_normal01 - Gaussian random variate with mu = 0 and sigma = 1.
--
-- This routine returns a Gaussian random variate with zero mean and
-- unit standard deviation. The polar (Box-Mueller) method is used.
--
-- This code is a modified version of the routine gsl_ran_gaussian from
-- the GNU Scientific Library Version 1.0. */

double fp_normal01(MPL *mpl)
{     double x, y, r2;
      do
      {  /* choose x, y in uniform square (-1,-1) to (+1,+1) */
         x = -1.0 + 2.0 * fp_uniform01(mpl);
         y = -1.0 + 2.0 * fp_uniform01(mpl);
         /* see if it is in the unit circle */
         r2 = x * x + y * y;
      } while (r2 > 1.0 || r2 == 0.0);
      /* Box-Muller transform */
      return y * sqrt(-2.0 * log (r2) / r2);
}

/*----------------------------------------------------------------------
-- fp_normal - Gaussian random variate with specified mu and sigma.
--
-- This routine returns a Gaussian random variate with mean mu and
-- standard deviation sigma. */

double fp_normal(MPL *mpl, double mu, double sigma)
{     double x;
#if 0
      x = mu + sigma * fp_normal01(mpl);
#else
      x = fp_add(mpl, mu, fp_mul(mpl, sigma, fp_normal01(mpl)));
#endif
      return x;
}

/**********************************************************************/
/* * *                SEGMENTED CHARACTER STRINGS                 * * */
/**********************************************************************/