glpgmp.c

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

/* glpgmp.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_STDIO
#include "glpdmp.h"
#include "glpgmp.h"
#define xfault xerror

#ifdef HAVE_GMP               /* use GNU MP bignum library */

int gmp_pool_count(void) { return 0; }

void gmp_free_mem(void) { return; }

#else                         /* use GLPK bignum module */

static DMP *gmp_pool = NULL;
static int gmp_size = 0;
static unsigned short *gmp_work = NULL;

void *gmp_get_atom(int size)
{     if (gmp_pool == NULL)
         gmp_pool = dmp_create_pool();
      return dmp_get_atom(gmp_pool, size);
}

void gmp_free_atom(void *ptr, int size)
{     xassert(gmp_pool != NULL);
      dmp_free_atom(gmp_pool, ptr, size);
      return;
}

int gmp_pool_count(void)
{     if (gmp_pool == NULL)
         return 0;
      else
         return dmp_in_use(gmp_pool).lo;
}

unsigned short *gmp_get_work(int size)
{     xassert(size > 0);
      if (gmp_size < size)
      {  if (gmp_size == 0)
         {  xassert(gmp_work == NULL);
            gmp_size = 100;
         }
         else
         {  xassert(gmp_work != NULL);
            xfree(gmp_work);
         }
         while (gmp_size < size) gmp_size += gmp_size;
         gmp_work = xcalloc(gmp_size, sizeof(unsigned short));
      }
      return gmp_work;
}

void gmp_free_mem(void)
{     if (gmp_pool != NULL) dmp_delete_pool(gmp_pool);
      if (gmp_work != NULL) xfree(gmp_work);
      gmp_pool = NULL;
      gmp_size = 0;
      gmp_work = NULL;
      return;
}

/*====================================================================*/

mpz_t _mpz_init(void)
{     /* initialize x, and set its value to 0 */
      mpz_t x;
      x = gmp_get_atom(sizeof(struct mpz));
      x->val = 0;
      x->ptr = NULL;
      return x;
}

void mpz_clear(mpz_t x)
{     /* free the space occupied by x */
      mpz_set_si(x, 0);
      xassert(x->ptr == NULL);
      /* free the number descriptor */
      gmp_free_atom(x, sizeof(struct mpz));
      return;
}

void mpz_set(mpz_t z, mpz_t x)
{     /* set the value of z from x */
      struct mpz_seg *e, *ee, *es;
      if (z != x)
      {  mpz_set_si(z, 0);
         z->val = x->val;
         xassert(z->ptr == NULL);
         for (e = x->ptr, es = NULL; e != NULL; e = e->next)
         {  ee = gmp_get_atom(sizeof(struct mpz_seg));
            memcpy(ee->d, e->d, 12);
            ee->next = NULL;
            if (z->ptr == NULL)
               z->ptr = ee;
            else
               es->next = ee;
            es = ee;
         }
      }
      return;
}

void mpz_set_si(mpz_t x, int val)
{     /* set the value of x to val */
      struct mpz_seg *e;
      /* free existing segments, if any */
      while (x->ptr != NULL)
      {  e = x->ptr;
         x->ptr = e->next;
         gmp_free_atom(e, sizeof(struct mpz_seg));
      }
      /* assign new value */
      if (val == 0x80000000)
      {  /* long format is needed */
         x->val = -1;
         x->ptr = e = gmp_get_atom(sizeof(struct mpz_seg));
         memset(e->d, 0, 12);
         e->d[1] = 0x8000;
         e->next = NULL;
      }
      else
      {  /* short format is enough */
         x->val = val;
      }
      return;
}

double mpz_get_d(mpz_t x)
{     /* convert x to a double, truncating if necessary */
      struct mpz_seg *e;
      int j;
      double val, deg;
      if (x->ptr == NULL)
         val = (double)x->val;
      else
      {  xassert(x->val != 0);
         val = 0.0;
         deg = 1.0;
         for (e = x->ptr; e != NULL; e = e->next)
         {  for (j = 0; j <= 5; j++)
            {  val += deg * (double)((int)e->d[j]);
               deg *= 65536.0;
            }
         }
         if (x->val < 0) val = - val;
      }
      return val;
}

double mpz_get_d_2exp(int *exp, mpz_t x)
{     /* convert x to a double, truncating if necessary (i.e. rounding
         towards zero), and returning the exponent separately;
         the return value is in the range 0.5 <= |d| < 1 and the
         exponent is stored to *exp; d*2^exp is the (truncated) x value;
         if x is zero, the return is 0.0 and 0 is stored to *exp;
         this is similar to the standard C frexp function */
      struct mpz_seg *e;
      int j, n, n1;
      double val;
      if (x->ptr == NULL)
         val = (double)x->val, n = 0;
      else
      {  xassert(x->val != 0);
         val = 0.0, n = 0;
         for (e = x->ptr; e != NULL; e = e->next)
         {  for (j = 0; j <= 5; j++)
            {  val += (double)((int)e->d[j]);
               val /= 65536.0, n += 16;
            }
         }
         if (x->val < 0) val = - val;
      }
      val = frexp(val, &n1);
      *exp = n + n1;
      return val;
}

void mpz_swap(mpz_t x, mpz_t y)
{     /* swap the values x and y efficiently */
      int val;
      void *ptr;
      val = x->val, ptr = x->ptr;
      x->val = y->val, x->ptr = y->ptr;
      y->val = val, y->ptr = ptr;
      return;
}

static void normalize(mpz_t x)
{     /* normalize integer x that includes removing non-significant
         (leading) zeros and converting to short format, if possible */
      struct mpz_seg *es, *e;
      /* if the integer is in short format, it remains unchanged */
      if (x->ptr == NULL)
      {  xassert(x->val != 0x80000000);
         goto done;
      }
      xassert(x->val == +1 || x->val == -1);
      /* find the last (most significant) non-zero segment */
      es = NULL;
      for (e = x->ptr; e != NULL; e = e->next)
      {  if (e->d[0] || e->d[1] || e->d[2] ||
             e->d[3] || e->d[4] || e->d[5]) es = e;
      }
      /* if all segments contain zeros, the integer is zero */
      if (es == NULL)
      {  mpz_set_si(x, 0);
         goto done;
      }
      /* remove non-significant (leading) zero segments */
      while (es->next != NULL)
      {  e = es->next;
         es->next = e->next;
         gmp_free_atom(e, sizeof(struct mpz_seg));
      }
      /* convert the integer to short format, if possible */
      e = x->ptr;
      if (e->next == NULL && e->d[1] <= 0x7FFF &&
         !e->d[2] && !e->d[3] && !e->d[4] && !e->d[5])
      {  int val;
         val = (int)e->d[0] + ((int)e->d[1] << 16);
         if (x->val < 0) val = - val;
         mpz_set_si(x, val);
      }
done: return;
}

void mpz_add(mpz_t z, mpz_t x, mpz_t y)
{     /* set z to x + y */
      static struct mpz_seg zero = { { 0, 0, 0, 0, 0, 0 }, NULL };
      struct mpz_seg dumx, dumy, *ex, *ey, *ez, *es, *ee;
      int k, sx, sy, sz;
      unsigned int t;
      /* if [x] = 0 then [z] = [y] */
      if (x->val == 0)
      {  xassert(x->ptr == NULL);
         mpz_set(z, y);
         goto done;
      }
      /* if [y] = 0 then [z] = [x] */
      if (y->val == 0)
      {  xassert(y->ptr == NULL);
         mpz_set(z, x);
         goto done;
      }
      /* special case when both [x] and [y] are in short format */
      if (x->ptr == NULL && y->ptr == NULL)
      {  int xval = x->val, yval = y->val, zval = x->val + y->val;
         xassert(xval != 0x80000000 && yval != 0x80000000);
         if (!(xval > 0 && yval > 0 && zval <= 0 ||
               xval < 0 && yval < 0 && zval >= 0))
         {  mpz_set_si(z, zval);
            goto done;
         }
      }
      /* convert [x] to long format, if necessary */
      if (x->ptr == NULL)
      {  xassert(x->val != 0x80000000);
         if (x->val >= 0)
         {  sx = +1;
            t = (unsigned int)(+ x->val);
         }
         else
         {  sx = -1;
            t = (unsigned int)(- x->val);
         }
         ex = &dumx;
         ex->d[0] = (unsigned short)t;
         ex->d[1] = (unsigned short)(t >> 16);
         ex->d[2] = ex->d[3] = ex->d[4] = ex->d[5] = 0;
         ex->next = NULL;
      }
      else
      {  sx = x->val;
         xassert(sx == +1 || sx == -1);
         ex = x->ptr;
      }
      /* convert [y] to long format, if necessary */
      if (y->ptr == NULL)
      {  xassert(y->val != 0x80000000);
         if (y->val >= 0)
         {  sy = +1;
            t = (unsigned int)(+ y->val);
         }
         else
         {  sy = -1;
            t = (unsigned int)(- y->val);
         }
         ey = &dumy;
         ey->d[0] = (unsigned short)t;
         ey->d[1] = (unsigned short)(t >> 16);
         ey->d[2] = ey->d[3] = ey->d[4] = ey->d[5] = 0;
         ey->next = NULL;
      }
      else
      {  sy = y->val;
         xassert(sy == +1 || sy == -1);
         ey = y->ptr;
      }
      /* main fragment */
      sz = sx;
      ez = es = NULL;
      if (sx > 0 && sy > 0 || sx < 0 && sy < 0)
      {  /* [x] and [y] have identical signs -- addition */
         t = 0;
         for (; ex || ey; ex = ex->next, ey = ey->next)
         {  if (ex == NULL) ex = &zero;
            if (ey == NULL) ey = &zero;
            ee = gmp_get_atom(sizeof(struct mpz_seg));
            for (k = 0; k <= 5; k++)
            {  t += (unsigned int)ex->d[k];
               t += (unsigned int)ey->d[k];
               ee->d[k] = (unsigned short)t;
               t >>= 16;
            }
            ee->next = NULL;
            if (ez == NULL)
               ez = ee;
            else
               es->next = ee;
            es = ee;
         }
         if (t)
         {  /* overflow -- one extra digit is needed */
            ee = gmp_get_atom(sizeof(struct mpz_seg));
            ee->d[0] = 1;
            ee->d[1] = ee->d[2] = ee->d[3] = ee->d[4] = ee->d[5] = 0;
            ee->next = NULL;
            xassert(es != NULL);
            es->next = ee;
         }
      }
      else
      {  /* [x] and [y] have different signs -- subtraction */
         t = 1;
         for (; ex || ey; ex = ex->next, ey = ey->next)
         {  if (ex == NULL) ex = &zero;
            if (ey == NULL) ey = &zero;
            ee = gmp_get_atom(sizeof(struct mpz_seg));
            for (k = 0; k <= 5; k++)
            {  t += (unsigned int)ex->d[k];
               t += (0xFFFF - (unsigned int)ey->d[k]);
               ee->d[k] = (unsigned short)t;
               t >>= 16;
            }
            ee->next = NULL;
            if (ez == NULL)