glpgmp.c

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

         x = y, or a nefative value if x < y */
      static struct mpz_seg zero = { { 0, 0, 0, 0, 0, 0 }, NULL };
      struct mpz_seg dumx, dumy, *ex, *ey;
      int cc, sx, sy, k;
      unsigned int t;
      if (x == y)
      {  cc = 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);
         cc = (xval > yval ? +1 : xval < yval ? -1 : 0);
         goto done;
      }
      /* special case when [x] and [y] have different signs */
      if (x->val > 0 && y->val <= 0 || x->val == 0 && y->val < 0)
      {  cc = +1;
         goto done;
      }
      if (x->val < 0 && y->val >= 0 || x->val == 0 && y->val > 0)
      {  cc = -1;
         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 */
      xassert(sx > 0 && sy > 0 || sx < 0 && sy < 0);
      cc = 0;
      for (; ex || ey; ex = ex->next, ey = ey->next)
      {  if (ex == NULL) ex = &zero;
         if (ey == NULL) ey = &zero;
         for (k = 0; k <= 5; k++)
         {  if (ex->d[k] > ey->d[k]) cc = +1;
            if (ex->d[k] < ey->d[k]) cc = -1;
         }
      }
      if (sx < 0) cc = - cc;
done: return cc;
}

int mpz_sgn(mpz_t x)
{     /* return +1 if x > 0, 0 if x = 0, and -1 if x < 0 */
      int s;
      s = (x->val > 0 ? +1 : x->val < 0 ? -1 : 0);
      return s;
}

int mpz_out_str(void *_fp, int base, mpz_t x)
{     /* output x on stream fp, as a string in given base; the base
         may vary from 2 to 36;
         return the number of bytes written, or if an error occurred,
         return 0 */
      FILE *fp = _fp;
      mpz_t b, y, r;
      int n, j, nwr = 0;
      unsigned char *d;
      static char *set = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ";
      if (!(2 <= base && base <= 36))
         xfault("mpz_out_str: base = %d; invalid base\n", base);
      mpz_init(b);
      mpz_set_si(b, base);
      mpz_init(y);
      mpz_init(r);
      /* determine the number of digits */
      mpz_abs(y, x);
      for (n = 0; mpz_sgn(y) != 0; n++)
         mpz_div(y, NULL, y, b);
      if (n == 0) n = 1;
      /* compute the digits */
      d = xmalloc(n);
      mpz_abs(y, x);
      for (j = 0; j < n; j++)
      {  mpz_div(y, r, y, b);
         xassert(0 <= r->val && r->val < base && r->ptr == NULL);
         d[j] = (unsigned char)r->val;
      }
      /* output the integer to the stream */
      if (fp == NULL) fp = stdout;
      if (mpz_sgn(x) < 0)
         fputc('-', fp), nwr++;
      for (j = n-1; j >= 0; j--)
         fputc(set[d[j]], fp), nwr++;
      if (ferror(fp)) nwr = 0;
      mpz_clear(b);
      mpz_clear(y);
      mpz_clear(r);
      xfree(d);
      return nwr;
}

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

mpq_t _mpq_init(void)
{     /* initialize x, and set its value to 0/1 */
      mpq_t x;
      x = gmp_get_atom(sizeof(struct mpq));
      x->p.val = 0;
      x->p.ptr = NULL;
      x->q.val = 1;
      x->q.ptr = NULL;
      return x;
}

void mpq_clear(mpq_t x)
{     /* free the space occupied by x */
      mpz_set_si(&x->p, 0);
      xassert(x->p.ptr == NULL);
      mpz_set_si(&x->q, 0);
      xassert(x->q.ptr == NULL);
      /* free the number descriptor */
      gmp_free_atom(x, sizeof(struct mpq));
      return;
}

void mpq_canonicalize(mpq_t x)
{     /* remove any factors that are common to the numerator and
         denominator of x, and make the denominator positive */
      mpz_t f;
      xassert(x->q.val != 0);
      if (x->q.val < 0)
      {  mpz_neg(&x->p, &x->p);
         mpz_neg(&x->q, &x->q);
      }
      mpz_init(f);
      mpz_gcd(f, &x->p, &x->q);
      if (!(f->val == 1 && f->ptr == NULL))
      {  mpz_div(&x->p, NULL, &x->p, f);
         mpz_div(&x->q, NULL, &x->q, f);
      }
      mpz_clear(f);
      return;
}

void mpq_set(mpq_t z, mpq_t x)
{     /* set the value of z from x */
      if (z != x)
      {  mpz_set(&z->p, &x->p);
         mpz_set(&z->q, &x->q);
      }
      return;
}

void mpq_set_si(mpq_t x, int p, unsigned int q)
{     /* set the value of x to p/q */
      if (q == 0)
         xfault("mpq_set_si: zero denominator not allowed\n");
      mpz_set_si(&x->p, p);
      xassert(q <= 0x7FFFFFFF);
      mpz_set_si(&x->q, q);
      return;
}

double mpq_get_d(mpq_t x)
{     /* convert x to a double, truncating if necessary */
      int np, nq;
      double p, q;
      p = mpz_get_d_2exp(&np, &x->p);
      q = mpz_get_d_2exp(&nq, &x->q);
      return ldexp(p / q, np - nq);
}

void mpq_set_d(mpq_t x, double val)
{     /* set x to val; there is no rounding, the conversion is exact */
      int s, n, d, j;
      double f;
      mpz_t temp;
      xassert(-DBL_MAX <= val && val <= +DBL_MAX);
      mpq_set_si(x, 0, 1);
      if (val > 0.0)
         s = +1;
      else if (val < 0.0)
         s = -1;
      else
         goto done;
      f = frexp(fabs(val), &n);
      /* |val| = f * 2^n, where 0.5 <= f < 1.0 */
      mpz_init(temp);
      while (f != 0.0)
      {  f *= 16.0, n -= 4;
         d = (int)f;
         xassert(0 <= d && d <= 15);
         f -= (double)d;
         /* x := 16 * x + d */
         mpz_set_si(temp, 16);
         mpz_mul(&x->p, &x->p, temp);
         mpz_set_si(temp, d);
         mpz_add(&x->p, &x->p, temp);
      }
      mpz_clear(temp);
      /* x := x * 2^n */
      if (n > 0)
      {  for (j = 1; j <= n; j++)
            mpz_add(&x->p, &x->p, &x->p);
      }
      else if (n < 0)
      {  for (j = 1; j <= -n; j++)
            mpz_add(&x->q, &x->q, &x->q);
         mpq_canonicalize(x);
      }
      if (s < 0) mpq_neg(x, x);
done: return;
}

void mpq_add(mpq_t z, mpq_t x, mpq_t y)
{     /* set z to x + y */
      mpz_t p, q;
      mpz_init(p);
      mpz_init(q);
      mpz_mul(p, &x->p, &y->q);
      mpz_mul(q, &x->q, &y->p);
      mpz_add(p, p, q);
      mpz_mul(q, &x->q, &y->q);
      mpz_set(&z->p, p);
      mpz_set(&z->q, q);
      mpz_clear(p);
      mpz_clear(q);
      mpq_canonicalize(z);
      return;
}

void mpq_sub(mpq_t z, mpq_t x, mpq_t y)
{     /* set z to x - y */
      mpz_t p, q;
      mpz_init(p);
      mpz_init(q);
      mpz_mul(p, &x->p, &y->q);
      mpz_mul(q, &x->q, &y->p);
      mpz_sub(p, p, q);
      mpz_mul(q, &x->q, &y->q);
      mpz_set(&z->p, p);
      mpz_set(&z->q, q);
      mpz_clear(p);
      mpz_clear(q);
      mpq_canonicalize(z);
      return;
}

void mpq_mul(mpq_t z, mpq_t x, mpq_t y)
{     /* set z to x * y */
      mpz_mul(&z->p, &x->p, &y->p);
      mpz_mul(&z->q, &x->q, &y->q);
      mpq_canonicalize(z);
      return;
}

void mpq_div(mpq_t z, mpq_t x, mpq_t y)
{     /* set z to x / y */
      mpz_t p, q;
      if (mpq_sgn(y) == 0)
         xfault("mpq_div: zero divisor not allowed\n");
      mpz_init(p);
      mpz_init(q);
      mpz_mul(p, &x->p, &y->q);
      mpz_mul(q, &x->q, &y->p);
      mpz_set(&z->p, p);
      mpz_set(&z->q, q);
      mpz_clear(p);
      mpz_clear(q);
      mpq_canonicalize(z);
      return;
}

void mpq_neg(mpq_t z, mpq_t x)
{     /* set z to 0 - x */
      mpq_set(z, x);
      mpz_neg(&z->p, &z->p);
      return;
}

void mpq_abs(mpq_t z, mpq_t x)
{     /* set z to the absolute value of x */
      mpq_set(z, x);
      mpz_abs(&z->p, &z->p);
      xassert(mpz_sgn(&x->q) > 0);
      return;
}

int mpq_cmp(mpq_t x, mpq_t y)
{     /* compare x and y; return a positive value if x > y, zero if
         x = y, or a nefative value if x < y */
      mpq_t temp;
      int s;
      mpq_init(temp);
      mpq_sub(temp, x, y);
      s = mpq_sgn(temp);
      mpq_clear(temp);
      return s;
}

int mpq_sgn(mpq_t x)
{     /* return +1 if x > 0, 0 if x = 0, and -1 if x < 0 */
      int s;
      s = mpz_sgn(&x->p);
      xassert(mpz_sgn(&x->q) > 0);
      return s;
}

int mpq_out_str(void *_fp, int base, mpq_t x)
{     /* output x on stream fp, as a string in given base; the base
         may vary from 2 to 36; output is in the form 'num/den' or if
         the denominator is 1 then just 'num';
         if the parameter fp is a null pointer, stdout is assumed;
         return the number of bytes written, or if an error occurred,
         return 0 */
      FILE *fp = _fp;
      int nwr;
      if (!(2 <= base && base <= 36))
         xfault("mpq_out_str: base = %d; invalid base\n", base);
      if (fp == NULL) fp = stdout;
      nwr = mpz_out_str(fp, base, &x->p);
      if (x->q.val == 1 && x->q.ptr == NULL)
         ;
      else
      {  fputc('/', fp), nwr++;
         nwr += mpz_out_str(fp, base, &x->q);
      }
      if (ferror(fp)) nwr = 0;
      return nwr;
}

#endif

/* eof */