glpgmp.c

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

               ez = ee;
            else
               es->next = ee;
            es = ee;
         }
         if (!t)
         {  /* |[x]| < |[y]| -- result in complement coding */
            sz = - sz;
            t = 1;
            for (ee = ez; ee != NULL; ee = ee->next)
            for (k = 0; k <= 5; k++)
            {  t += (0xFFFF - (unsigned int)ee->d[k]);
               ee->d[k] = (unsigned short)t;
               t >>= 16;
            }
         }
      }
      /* contruct and normalize result */
      mpz_set_si(z, 0);
      z->val = sz;
      z->ptr = ez;
      normalize(z);
done: return;
}

void mpz_sub(mpz_t z, mpz_t x, mpz_t y)
{     /* set z to x - y */
      if (x == y)
         mpz_set_si(z, 0);
      else
      {  y->val = - y->val;
         mpz_add(z, x, y);
         if (y != z) y->val = - y->val;
      }
      return;
}

void mpz_mul(mpz_t z, mpz_t x, mpz_t y)
{     /* set z to x * y */
      struct mpz_seg dumx, dumy, *ex, *ey, *es, *e;
      int sx, sy, k, nx, ny, n;
      unsigned int t;
      unsigned short *work, *wx, *wy;
      /* if [x] = 0 then [z] = 0 */
      if (x->val == 0)
      {  xassert(x->ptr == NULL);
         mpz_set_si(z, 0);
         goto done;
      }
      /* if [y] = 0 then [z] = 0 */
      if (y->val == 0)
      {  xassert(y->ptr == NULL);
         mpz_set_si(z, 0);
         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, sz = +1;
         xassert(xval != 0x80000000 && yval != 0x80000000);
         if (xval < 0) xval = - xval, sz = - sz;
         if (yval < 0) yval = - yval, sz = - sz;
         if (xval <= 0x7FFFFFFF / yval)
         {  mpz_set_si(z, sz * (xval * yval));
            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;
      }
      /* determine the number of digits of [x] */
      nx = n = 0;
      for (e = ex; e != NULL; e = e->next)
      for (k = 0; k <= 5; k++)
      {  n++;
         if (e->d[k]) nx = n;
      }
      xassert(nx > 0);
      /* determine the number of digits of [y] */
      ny = n = 0;
      for (e = ey; e != NULL; e = e->next)
      for (k = 0; k <= 5; k++)
      {  n++;
         if (e->d[k]) ny = n;
      }
      xassert(ny > 0);
      /* we need working array containing at least nx+ny+ny places */
      work = gmp_get_work(nx+ny+ny);
      /* load digits of [x] */
      wx = &work[0];
      for (n = 0; n < nx; n++) wx[ny+n] = 0;
      for (n = 0, e = ex; e != NULL; e = e->next)
         for (k = 0; k <= 5; k++, n++)
            if (e->d[k]) wx[ny+n] = e->d[k];
      /* load digits of [y] */
      wy = &work[nx+ny];
      for (n = 0; n < ny; n++) wy[n] = 0;
      for (n = 0, e = ey; e != NULL; e = e->next)
         for (k = 0; k <= 5; k++, n++)
            if (e->d[k]) wy[n] = e->d[k];
      /* compute [x] * [y] */
      bigmul(nx, ny, wx, wy);
      /* construct and normalize result */
      mpz_set_si(z, 0);
      z->val = sx * sy;
      es = NULL;
      k = 6;
      for (n = 0; n < nx+ny; n++)
      {  if (k > 5)
         {  e = gmp_get_atom(sizeof(struct mpz_seg));
            e->d[0] = e->d[1] = e->d[2] = 0;
            e->d[3] = e->d[4] = e->d[5] = 0;
            e->next = NULL;
            if (z->ptr == NULL)
               z->ptr = e;
            else
               es->next = e;
            es = e;
            k = 0;
         }
         es->d[k++] = wx[n];
      }
      normalize(z);
done: return;
}

void mpz_neg(mpz_t z, mpz_t x)
{     /* set z to 0 - x */
      mpz_set(z, x);
      z->val = - z->val;
      return;
}

void mpz_abs(mpz_t z, mpz_t x)
{     /* set z to the absolute value of x */
      mpz_set(z, x);
      if (z->val < 0) z->val = - z->val;
      return;
}

void mpz_div(mpz_t q, mpz_t r, mpz_t x, mpz_t y)
{     /* divide x by y, forming quotient q and/or remainder r
         if q = NULL then quotient is not stored; if r = NULL then
         remainder is not stored
         the sign of quotient is determined as in algebra while the
         sign of remainder is the same as the sign of dividend:
         +26 : +7 = +3, remainder is +5
         -26 : +7 = -3, remainder is -5
         +26 : -7 = -3, remainder is +5
         -26 : -7 = +3, remainder is -5 */
      struct mpz_seg dumx, dumy, *ex, *ey, *es, *e;
      int sx, sy, k, nx, ny, n;
      unsigned int t;
      unsigned short *work, *wx, *wy;
      /* divide by zero is not allowed */
      if (y->val == 0)
      {  xassert(y->ptr == NULL);
         xfault("mpz_div: divide by zero not allowed\n");
      }
      /* if [x] = 0 then [q] = [r] = 0 */
      if (x->val == 0)
      {  xassert(x->ptr == NULL);
         if (q != NULL) mpz_set_si(q, 0);
         if (r != NULL) mpz_set_si(r, 0);
         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;
         xassert(xval != 0x80000000 && yval != 0x80000000);
         if (q != NULL) mpz_set_si(q, xval / yval);
         if (r != NULL) mpz_set_si(r, xval % yval);
         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;
      }
      /* determine the number of digits of [x] */
      nx = n = 0;
      for (e = ex; e != NULL; e = e->next)
      for (k = 0; k <= 5; k++)
      {  n++;
         if (e->d[k]) nx = n;
      }
      xassert(nx > 0);
      /* determine the number of digits of [y] */
      ny = n = 0;
      for (e = ey; e != NULL; e = e->next)
      for (k = 0; k <= 5; k++)
      {  n++;
         if (e->d[k]) ny = n;
      }
      xassert(ny > 0);
      /* if nx < ny then [q] = 0 and [r] = [x] */
      if (nx < ny)
      {  if (r != NULL) mpz_set(r, x);
         if (q != NULL) mpz_set_si(q, 0);
         goto done;
      }
      /* we need working array containing at least nx+ny+1 places */
      work = gmp_get_work(nx+ny+1);
      /* load digits of [x] */
      wx = &work[0];
      for (n = 0; n < nx; n++) wx[n] = 0;
      for (n = 0, e = ex; e != NULL; e = e->next)
         for (k = 0; k <= 5; k++, n++)
            if (e->d[k]) wx[n] = e->d[k];
      /* load digits of [y] */
      wy = &work[nx+1];
      for (n = 0; n < ny; n++) wy[n] = 0;
      for (n = 0, e = ey; e != NULL; e = e->next)
         for (k = 0; k <= 5; k++, n++)
            if (e->d[k]) wy[n] = e->d[k];
      /* compute quotient and remainder */
      xassert(wy[ny-1] != 0);
      bigdiv(nx-ny, ny, wx, wy);
      /* construct and normalize quotient */
      if (q != NULL)
      {  mpz_set_si(q, 0);
         q->val = sx * sy;
         es = NULL;
         k = 6;
         for (n = ny; n <= nx; n++)
         {  if (k > 5)
            {  e = gmp_get_atom(sizeof(struct mpz_seg));
               e->d[0] = e->d[1] = e->d[2] = 0;
               e->d[3] = e->d[4] = e->d[5] = 0;
               e->next = NULL;
               if (q->ptr == NULL)
                  q->ptr = e;
               else
                  es->next = e;
               es = e;
               k = 0;
            }
            es->d[k++] = wx[n];
         }
         normalize(q);
      }
      /* construct and normalize remainder */
      if (r != NULL)
      {  mpz_set_si(r, 0);
         r->val = sx;
         es = NULL;
         k = 6;
         for (n = 0; n < ny; n++)
         {  if (k > 5)
            {  e = gmp_get_atom(sizeof(struct mpz_seg));
               e->d[0] = e->d[1] = e->d[2] = 0;
               e->d[3] = e->d[4] = e->d[5] = 0;
               e->next = NULL;
               if (r->ptr == NULL)
                  r->ptr = e;
               else
                  es->next = e;
               es = e;
               k = 0;
            }
            es->d[k++] = wx[n];
         }
         normalize(r);
      }
done: return;
}

void mpz_gcd(mpz_t z, mpz_t x, mpz_t y)
{     /* set z to the greatest common divisor of x and y */
      /* in case of arbitrary integers GCD(x, y) = GCD(|x|, |y|), and,
         in particular, GCD(0, 0) = 0 */
      mpz_t u, v, r;
      mpz_init(u);
      mpz_init(v);
      mpz_init(r);
      mpz_abs(u, x);
      mpz_abs(v, y);
      while (mpz_sgn(v))
      {  mpz_div(NULL, r, u, v);
         mpz_set(u, v);
         mpz_set(v, r);
      }
      mpz_set(z, u);
      mpz_clear(u);
      mpz_clear(v);
      mpz_clear(r);
      return;
}

int mpz_cmp(mpz_t x, mpz_t y)
{     /* compare x and y; return a positive value if x > y, zero if