1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/deps/glpk/src/glpgmp.c	Sun Nov 06 20:59:10 2011 +0100
     1.3 @@ -0,0 +1,1108 @@
     1.4 +/* glpgmp.c */
     1.5 +
     1.6 +/***********************************************************************
     1.7 +*  This code is part of GLPK (GNU Linear Programming Kit).
     1.8 +*
     1.9 +*  Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
    1.10 +*  2009, 2010, 2011 Andrew Makhorin, Department for Applied Informatics,
    1.11 +*  Moscow Aviation Institute, Moscow, Russia. All rights reserved.
    1.12 +*  E-mail: <mao@gnu.org>.
    1.13 +*
    1.14 +*  GLPK is free software: you can redistribute it and/or modify it
    1.15 +*  under the terms of the GNU General Public License as published by
    1.16 +*  the Free Software Foundation, either version 3 of the License, or
    1.17 +*  (at your option) any later version.
    1.18 +*
    1.19 +*  GLPK is distributed in the hope that it will be useful, but WITHOUT
    1.20 +*  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
    1.21 +*  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
    1.22 +*  License for more details.
    1.23 +*
    1.24 +*  You should have received a copy of the GNU General Public License
    1.25 +*  along with GLPK. If not, see <http://www.gnu.org/licenses/>.
    1.26 +***********************************************************************/
    1.27 +
    1.28 +#define _GLPSTD_STDIO
    1.29 +#include "glpdmp.h"
    1.30 +#include "glpgmp.h"
    1.31 +#define xfault xerror
    1.32 +
    1.33 +#ifdef HAVE_GMP               /* use GNU MP bignum library */
    1.34 +
    1.35 +int gmp_pool_count(void) { return 0; }
    1.36 +
    1.37 +void gmp_free_mem(void) { return; }
    1.38 +
    1.39 +#else                         /* use GLPK bignum module */
    1.40 +
    1.41 +static DMP *gmp_pool = NULL;
    1.42 +static int gmp_size = 0;
    1.43 +static unsigned short *gmp_work = NULL;
    1.44 +
    1.45 +void *gmp_get_atom(int size)
    1.46 +{     if (gmp_pool == NULL)
    1.47 +         gmp_pool = dmp_create_pool();
    1.48 +      return dmp_get_atom(gmp_pool, size);
    1.49 +}
    1.50 +
    1.51 +void gmp_free_atom(void *ptr, int size)
    1.52 +{     xassert(gmp_pool != NULL);
    1.53 +      dmp_free_atom(gmp_pool, ptr, size);
    1.54 +      return;
    1.55 +}
    1.56 +
    1.57 +int gmp_pool_count(void)
    1.58 +{     if (gmp_pool == NULL)
    1.59 +         return 0;
    1.60 +      else
    1.61 +         return dmp_in_use(gmp_pool).lo;
    1.62 +}
    1.63 +
    1.64 +unsigned short *gmp_get_work(int size)
    1.65 +{     xassert(size > 0);
    1.66 +      if (gmp_size < size)
    1.67 +      {  if (gmp_size == 0)
    1.68 +         {  xassert(gmp_work == NULL);
    1.69 +            gmp_size = 100;
    1.70 +         }
    1.71 +         else
    1.72 +         {  xassert(gmp_work != NULL);
    1.73 +            xfree(gmp_work);
    1.74 +         }
    1.75 +         while (gmp_size < size) gmp_size += gmp_size;
    1.76 +         gmp_work = xcalloc(gmp_size, sizeof(unsigned short));
    1.77 +      }
    1.78 +      return gmp_work;
    1.79 +}
    1.80 +
    1.81 +void gmp_free_mem(void)
    1.82 +{     if (gmp_pool != NULL) dmp_delete_pool(gmp_pool);
    1.83 +      if (gmp_work != NULL) xfree(gmp_work);
    1.84 +      gmp_pool = NULL;
    1.85 +      gmp_size = 0;
    1.86 +      gmp_work = NULL;
    1.87 +      return;
    1.88 +}
    1.89 +
    1.90 +/*====================================================================*/
    1.91 +
    1.92 +mpz_t _mpz_init(void)
    1.93 +{     /* initialize x, and set its value to 0 */
    1.94 +      mpz_t x;
    1.95 +      x = gmp_get_atom(sizeof(struct mpz));
    1.96 +      x->val = 0;
    1.97 +      x->ptr = NULL;
    1.98 +      return x;
    1.99 +}
   1.100 +
   1.101 +void mpz_clear(mpz_t x)
   1.102 +{     /* free the space occupied by x */
   1.103 +      mpz_set_si(x, 0);
   1.104 +      xassert(x->ptr == NULL);
   1.105 +      /* free the number descriptor */
   1.106 +      gmp_free_atom(x, sizeof(struct mpz));
   1.107 +      return;
   1.108 +}
   1.109 +
   1.110 +void mpz_set(mpz_t z, mpz_t x)
   1.111 +{     /* set the value of z from x */
   1.112 +      struct mpz_seg *e, *ee, *es;
   1.113 +      if (z != x)
   1.114 +      {  mpz_set_si(z, 0);
   1.115 +         z->val = x->val;
   1.116 +         xassert(z->ptr == NULL);
   1.117 +         for (e = x->ptr, es = NULL; e != NULL; e = e->next)
   1.118 +         {  ee = gmp_get_atom(sizeof(struct mpz_seg));
   1.119 +            memcpy(ee->d, e->d, 12);
   1.120 +            ee->next = NULL;
   1.121 +            if (z->ptr == NULL)
   1.122 +               z->ptr = ee;
   1.123 +            else
   1.124 +               es->next = ee;
   1.125 +            es = ee;
   1.126 +         }
   1.127 +      }
   1.128 +      return;
   1.129 +}
   1.130 +
   1.131 +void mpz_set_si(mpz_t x, int val)
   1.132 +{     /* set the value of x to val */
   1.133 +      struct mpz_seg *e;
   1.134 +      /* free existing segments, if any */
   1.135 +      while (x->ptr != NULL)
   1.136 +      {  e = x->ptr;
   1.137 +         x->ptr = e->next;
   1.138 +         gmp_free_atom(e, sizeof(struct mpz_seg));
   1.139 +      }
   1.140 +      /* assign new value */
   1.141 +      if (val == 0x80000000)
   1.142 +      {  /* long format is needed */
   1.143 +         x->val = -1;
   1.144 +         x->ptr = e = gmp_get_atom(sizeof(struct mpz_seg));
   1.145 +         memset(e->d, 0, 12);
   1.146 +         e->d[1] = 0x8000;
   1.147 +         e->next = NULL;
   1.148 +      }
   1.149 +      else
   1.150 +      {  /* short format is enough */
   1.151 +         x->val = val;
   1.152 +      }
   1.153 +      return;
   1.154 +}
   1.155 +
   1.156 +double mpz_get_d(mpz_t x)
   1.157 +{     /* convert x to a double, truncating if necessary */
   1.158 +      struct mpz_seg *e;
   1.159 +      int j;
   1.160 +      double val, deg;
   1.161 +      if (x->ptr == NULL)
   1.162 +         val = (double)x->val;
   1.163 +      else
   1.164 +      {  xassert(x->val != 0);
   1.165 +         val = 0.0;
   1.166 +         deg = 1.0;
   1.167 +         for (e = x->ptr; e != NULL; e = e->next)
   1.168 +         {  for (j = 0; j <= 5; j++)
   1.169 +            {  val += deg * (double)((int)e->d[j]);
   1.170 +               deg *= 65536.0;
   1.171 +            }
   1.172 +         }
   1.173 +         if (x->val < 0) val = - val;
   1.174 +      }
   1.175 +      return val;
   1.176 +}
   1.177 +
   1.178 +double mpz_get_d_2exp(int *exp, mpz_t x)
   1.179 +{     /* convert x to a double, truncating if necessary (i.e. rounding
   1.180 +         towards zero), and returning the exponent separately;
   1.181 +         the return value is in the range 0.5 <= |d| < 1 and the
   1.182 +         exponent is stored to *exp; d*2^exp is the (truncated) x value;
   1.183 +         if x is zero, the return is 0.0 and 0 is stored to *exp;
   1.184 +         this is similar to the standard C frexp function */
   1.185 +      struct mpz_seg *e;
   1.186 +      int j, n, n1;
   1.187 +      double val;
   1.188 +      if (x->ptr == NULL)
   1.189 +         val = (double)x->val, n = 0;
   1.190 +      else
   1.191 +      {  xassert(x->val != 0);
   1.192 +         val = 0.0, n = 0;
   1.193 +         for (e = x->ptr; e != NULL; e = e->next)
   1.194 +         {  for (j = 0; j <= 5; j++)
   1.195 +            {  val += (double)((int)e->d[j]);
   1.196 +               val /= 65536.0, n += 16;
   1.197 +            }
   1.198 +         }
   1.199 +         if (x->val < 0) val = - val;
   1.200 +      }
   1.201 +      val = frexp(val, &n1);
   1.202 +      *exp = n + n1;
   1.203 +      return val;
   1.204 +}
   1.205 +
   1.206 +void mpz_swap(mpz_t x, mpz_t y)
   1.207 +{     /* swap the values x and y efficiently */
   1.208 +      int val;
   1.209 +      void *ptr;
   1.210 +      val = x->val, ptr = x->ptr;
   1.211 +      x->val = y->val, x->ptr = y->ptr;
   1.212 +      y->val = val, y->ptr = ptr;
   1.213 +      return;
   1.214 +}
   1.215 +
   1.216 +static void normalize(mpz_t x)
   1.217 +{     /* normalize integer x that includes removing non-significant
   1.218 +         (leading) zeros and converting to short format, if possible */
   1.219 +      struct mpz_seg *es, *e;
   1.220 +      /* if the integer is in short format, it remains unchanged */
   1.221 +      if (x->ptr == NULL)
   1.222 +      {  xassert(x->val != 0x80000000);
   1.223 +         goto done;
   1.224 +      }
   1.225 +      xassert(x->val == +1 || x->val == -1);
   1.226 +      /* find the last (most significant) non-zero segment */
   1.227 +      es = NULL;
   1.228 +      for (e = x->ptr; e != NULL; e = e->next)
   1.229 +      {  if (e->d[0] || e->d[1] || e->d[2] ||
   1.230 +             e->d[3] || e->d[4] || e->d[5]) es = e;
   1.231 +      }
   1.232 +      /* if all segments contain zeros, the integer is zero */
   1.233 +      if (es == NULL)
   1.234 +      {  mpz_set_si(x, 0);
   1.235 +         goto done;
   1.236 +      }
   1.237 +      /* remove non-significant (leading) zero segments */
   1.238 +      while (es->next != NULL)
   1.239 +      {  e = es->next;
   1.240 +         es->next = e->next;
   1.241 +         gmp_free_atom(e, sizeof(struct mpz_seg));
   1.242 +      }
   1.243 +      /* convert the integer to short format, if possible */
   1.244 +      e = x->ptr;
   1.245 +      if (e->next == NULL && e->d[1] <= 0x7FFF &&
   1.246 +         !e->d[2] && !e->d[3] && !e->d[4] && !e->d[5])
   1.247 +      {  int val;
   1.248 +         val = (int)e->d[0] + ((int)e->d[1] << 16);
   1.249 +         if (x->val < 0) val = - val;
   1.250 +         mpz_set_si(x, val);
   1.251 +      }
   1.252 +done: return;
   1.253 +}
   1.254 +
   1.255 +void mpz_add(mpz_t z, mpz_t x, mpz_t y)
   1.256 +{     /* set z to x + y */
   1.257 +      static struct mpz_seg zero = { { 0, 0, 0, 0, 0, 0 }, NULL };
   1.258 +      struct mpz_seg dumx, dumy, *ex, *ey, *ez, *es, *ee;
   1.259 +      int k, sx, sy, sz;
   1.260 +      unsigned int t;
   1.261 +      /* if [x] = 0 then [z] = [y] */
   1.262 +      if (x->val == 0)
   1.263 +      {  xassert(x->ptr == NULL);
   1.264 +         mpz_set(z, y);
   1.265 +         goto done;
   1.266 +      }
   1.267 +      /* if [y] = 0 then [z] = [x] */
   1.268 +      if (y->val == 0)
   1.269 +      {  xassert(y->ptr == NULL);
   1.270 +         mpz_set(z, x);
   1.271 +         goto done;
   1.272 +      }
   1.273 +      /* special case when both [x] and [y] are in short format */
   1.274 +      if (x->ptr == NULL && y->ptr == NULL)
   1.275 +      {  int xval = x->val, yval = y->val, zval = x->val + y->val;
   1.276 +         xassert(xval != 0x80000000 && yval != 0x80000000);
   1.277 +         if (!(xval > 0 && yval > 0 && zval <= 0 ||
   1.278 +               xval < 0 && yval < 0 && zval >= 0))
   1.279 +         {  mpz_set_si(z, zval);
   1.280 +            goto done;
   1.281 +         }
   1.282 +      }
   1.283 +      /* convert [x] to long format, if necessary */
   1.284 +      if (x->ptr == NULL)
   1.285 +      {  xassert(x->val != 0x80000000);
   1.286 +         if (x->val >= 0)
   1.287 +         {  sx = +1;
   1.288 +            t = (unsigned int)(+ x->val);
   1.289 +         }
   1.290 +         else
   1.291 +         {  sx = -1;
   1.292 +            t = (unsigned int)(- x->val);
   1.293 +         }
   1.294 +         ex = &dumx;
   1.295 +         ex->d[0] = (unsigned short)t;
   1.296 +         ex->d[1] = (unsigned short)(t >> 16);
   1.297 +         ex->d[2] = ex->d[3] = ex->d[4] = ex->d[5] = 0;
   1.298 +         ex->next = NULL;
   1.299 +      }
   1.300 +      else
   1.301 +      {  sx = x->val;
   1.302 +         xassert(sx == +1 || sx == -1);
   1.303 +         ex = x->ptr;
   1.304 +      }
   1.305 +      /* convert [y] to long format, if necessary */
   1.306 +      if (y->ptr == NULL)
   1.307 +      {  xassert(y->val != 0x80000000);
   1.308 +         if (y->val >= 0)
   1.309 +         {  sy = +1;
   1.310 +            t = (unsigned int)(+ y->val);
   1.311 +         }
   1.312 +         else
   1.313 +         {  sy = -1;
   1.314 +            t = (unsigned int)(- y->val);
   1.315 +         }
   1.316 +         ey = &dumy;
   1.317 +         ey->d[0] = (unsigned short)t;
   1.318 +         ey->d[1] = (unsigned short)(t >> 16);
   1.319 +         ey->d[2] = ey->d[3] = ey->d[4] = ey->d[5] = 0;
   1.320 +         ey->next = NULL;
   1.321 +      }
   1.322 +      else
   1.323 +      {  sy = y->val;
   1.324 +         xassert(sy == +1 || sy == -1);
   1.325 +         ey = y->ptr;
   1.326 +      }
   1.327 +      /* main fragment */
   1.328 +      sz = sx;
   1.329 +      ez = es = NULL;
   1.330 +      if (sx > 0 && sy > 0 || sx < 0 && sy < 0)
   1.331 +      {  /* [x] and [y] have identical signs -- addition */
   1.332 +         t = 0;
   1.333 +         for (; ex || ey; ex = ex->next, ey = ey->next)
   1.334 +         {  if (ex == NULL) ex = &zero;
   1.335 +            if (ey == NULL) ey = &zero;
   1.336 +            ee = gmp_get_atom(sizeof(struct mpz_seg));
   1.337 +            for (k = 0; k <= 5; k++)
   1.338 +            {  t += (unsigned int)ex->d[k];
   1.339 +               t += (unsigned int)ey->d[k];
   1.340 +               ee->d[k] = (unsigned short)t;
   1.341 +               t >>= 16;
   1.342 +            }
   1.343 +            ee->next = NULL;
   1.344 +            if (ez == NULL)
   1.345 +               ez = ee;
   1.346 +            else
   1.347 +               es->next = ee;
   1.348 +            es = ee;
   1.349 +         }
   1.350 +         if (t)
   1.351 +         {  /* overflow -- one extra digit is needed */
   1.352 +            ee = gmp_get_atom(sizeof(struct mpz_seg));
   1.353 +            ee->d[0] = 1;
   1.354 +            ee->d[1] = ee->d[2] = ee->d[3] = ee->d[4] = ee->d[5] = 0;
   1.355 +            ee->next = NULL;
   1.356 +            xassert(es != NULL);
   1.357 +            es->next = ee;
   1.358 +         }
   1.359 +      }
   1.360 +      else
   1.361 +      {  /* [x] and [y] have different signs -- subtraction */
   1.362 +         t = 1;
   1.363 +         for (; ex || ey; ex = ex->next, ey = ey->next)
   1.364 +         {  if (ex == NULL) ex = &zero;
   1.365 +            if (ey == NULL) ey = &zero;
   1.366 +            ee = gmp_get_atom(sizeof(struct mpz_seg));
   1.367 +            for (k = 0; k <= 5; k++)
   1.368 +            {  t += (unsigned int)ex->d[k];
   1.369 +               t += (0xFFFF - (unsigned int)ey->d[k]);
   1.370 +               ee->d[k] = (unsigned short)t;
   1.371 +               t >>= 16;
   1.372 +            }
   1.373 +            ee->next = NULL;
   1.374 +            if (ez == NULL)
   1.375 +               ez = ee;
   1.376 +            else
   1.377 +               es->next = ee;
   1.378 +            es = ee;
   1.379 +         }
   1.380 +         if (!t)
   1.381 +         {  /* |[x]| < |[y]| -- result in complement coding */
   1.382 +            sz = - sz;
   1.383 +            t = 1;
   1.384 +            for (ee = ez; ee != NULL; ee = ee->next)
   1.385 +            for (k = 0; k <= 5; k++)
   1.386 +            {  t += (0xFFFF - (unsigned int)ee->d[k]);
   1.387 +               ee->d[k] = (unsigned short)t;
   1.388 +               t >>= 16;
   1.389 +            }
   1.390 +         }
   1.391 +      }
   1.392 +      /* contruct and normalize result */
   1.393 +      mpz_set_si(z, 0);
   1.394 +      z->val = sz;
   1.395 +      z->ptr = ez;
   1.396 +      normalize(z);
   1.397 +done: return;
   1.398 +}
   1.399 +
   1.400 +void mpz_sub(mpz_t z, mpz_t x, mpz_t y)
   1.401 +{     /* set z to x - y */
   1.402 +      if (x == y)
   1.403 +         mpz_set_si(z, 0);
   1.404 +      else
   1.405 +      {  y->val = - y->val;
   1.406 +         mpz_add(z, x, y);
   1.407 +         if (y != z) y->val = - y->val;
   1.408 +      }
   1.409 +      return;
   1.410 +}
   1.411 +
   1.412 +void mpz_mul(mpz_t z, mpz_t x, mpz_t y)
   1.413 +{     /* set z to x * y */
   1.414 +      struct mpz_seg dumx, dumy, *ex, *ey, *es, *e;
   1.415 +      int sx, sy, k, nx, ny, n;
   1.416 +      unsigned int t;
   1.417 +      unsigned short *work, *wx, *wy;
   1.418 +      /* if [x] = 0 then [z] = 0 */
   1.419 +      if (x->val == 0)
   1.420 +      {  xassert(x->ptr == NULL);
   1.421 +         mpz_set_si(z, 0);
   1.422 +         goto done;
   1.423 +      }
   1.424 +      /* if [y] = 0 then [z] = 0 */
   1.425 +      if (y->val == 0)
   1.426 +      {  xassert(y->ptr == NULL);
   1.427 +         mpz_set_si(z, 0);
   1.428 +         goto done;
   1.429 +      }
   1.430 +      /* special case when both [x] and [y] are in short format */
   1.431 +      if (x->ptr == NULL && y->ptr == NULL)
   1.432 +      {  int xval = x->val, yval = y->val, sz = +1;
   1.433 +         xassert(xval != 0x80000000 && yval != 0x80000000);
   1.434 +         if (xval < 0) xval = - xval, sz = - sz;
   1.435 +         if (yval < 0) yval = - yval, sz = - sz;
   1.436 +         if (xval <= 0x7FFFFFFF / yval)
   1.437 +         {  mpz_set_si(z, sz * (xval * yval));
   1.438 +            goto done;
   1.439 +         }
   1.440 +      }
   1.441 +      /* convert [x] to long format, if necessary */
   1.442 +      if (x->ptr == NULL)
   1.443 +      {  xassert(x->val != 0x80000000);
   1.444 +         if (x->val >= 0)
   1.445 +         {  sx = +1;
   1.446 +            t = (unsigned int)(+ x->val);
   1.447 +         }
   1.448 +         else
   1.449 +         {  sx = -1;
   1.450 +            t = (unsigned int)(- x->val);
   1.451 +         }
   1.452 +         ex = &dumx;
   1.453 +         ex->d[0] = (unsigned short)t;
   1.454 +         ex->d[1] = (unsigned short)(t >> 16);
   1.455 +         ex->d[2] = ex->d[3] = ex->d[4] = ex->d[5] = 0;
   1.456 +         ex->next = NULL;
   1.457 +      }
   1.458 +      else
   1.459 +      {  sx = x->val;
   1.460 +         xassert(sx == +1 || sx == -1);
   1.461 +         ex = x->ptr;
   1.462 +      }
   1.463 +      /* convert [y] to long format, if necessary */
   1.464 +      if (y->ptr == NULL)
   1.465 +      {  xassert(y->val != 0x80000000);
   1.466 +         if (y->val >= 0)
   1.467 +         {  sy = +1;
   1.468 +            t = (unsigned int)(+ y->val);
   1.469 +         }
   1.470 +         else
   1.471 +         {  sy = -1;
   1.472 +            t = (unsigned int)(- y->val);
   1.473 +         }
   1.474 +         ey = &dumy;
   1.475 +         ey->d[0] = (unsigned short)t;
   1.476 +         ey->d[1] = (unsigned short)(t >> 16);
   1.477 +         ey->d[2] = ey->d[3] = ey->d[4] = ey->d[5] = 0;
   1.478 +         ey->next = NULL;
   1.479 +      }
   1.480 +      else
   1.481 +      {  sy = y->val;
   1.482 +         xassert(sy == +1 || sy == -1);
   1.483 +         ey = y->ptr;
   1.484 +      }
   1.485 +      /* determine the number of digits of [x] */
   1.486 +      nx = n = 0;
   1.487 +      for (e = ex; e != NULL; e = e->next)
   1.488 +      for (k = 0; k <= 5; k++)
   1.489 +      {  n++;
   1.490 +         if (e->d[k]) nx = n;
   1.491 +      }
   1.492 +      xassert(nx > 0);
   1.493 +      /* determine the number of digits of [y] */
   1.494 +      ny = n = 0;
   1.495 +      for (e = ey; e != NULL; e = e->next)
   1.496 +      for (k = 0; k <= 5; k++)
   1.497 +      {  n++;
   1.498 +         if (e->d[k]) ny = n;
   1.499 +      }
   1.500 +      xassert(ny > 0);
   1.501 +      /* we need working array containing at least nx+ny+ny places */
   1.502 +      work = gmp_get_work(nx+ny+ny);
   1.503 +      /* load digits of [x] */
   1.504 +      wx = &work[0];
   1.505 +      for (n = 0; n < nx; n++) wx[ny+n] = 0;
   1.506 +      for (n = 0, e = ex; e != NULL; e = e->next)
   1.507 +         for (k = 0; k <= 5; k++, n++)
   1.508 +            if (e->d[k]) wx[ny+n] = e->d[k];
   1.509 +      /* load digits of [y] */
   1.510 +      wy = &work[nx+ny];
   1.511 +      for (n = 0; n < ny; n++) wy[n] = 0;
   1.512 +      for (n = 0, e = ey; e != NULL; e = e->next)
   1.513 +         for (k = 0; k <= 5; k++, n++)
   1.514 +            if (e->d[k]) wy[n] = e->d[k];
   1.515 +      /* compute [x] * [y] */
   1.516 +      bigmul(nx, ny, wx, wy);
   1.517 +      /* construct and normalize result */
   1.518 +      mpz_set_si(z, 0);
   1.519 +      z->val = sx * sy;
   1.520 +      es = NULL;
   1.521 +      k = 6;
   1.522 +      for (n = 0; n < nx+ny; n++)
   1.523 +      {  if (k > 5)
   1.524 +         {  e = gmp_get_atom(sizeof(struct mpz_seg));
   1.525 +            e->d[0] = e->d[1] = e->d[2] = 0;
   1.526 +            e->d[3] = e->d[4] = e->d[5] = 0;
   1.527 +            e->next = NULL;
   1.528 +            if (z->ptr == NULL)
   1.529 +               z->ptr = e;
   1.530 +            else
   1.531 +               es->next = e;
   1.532 +            es = e;
   1.533 +            k = 0;
   1.534 +         }
   1.535 +         es->d[k++] = wx[n];
   1.536 +      }
   1.537 +      normalize(z);
   1.538 +done: return;
   1.539 +}
   1.540 +
   1.541 +void mpz_neg(mpz_t z, mpz_t x)
   1.542 +{     /* set z to 0 - x */
   1.543 +      mpz_set(z, x);
   1.544 +      z->val = - z->val;
   1.545 +      return;
   1.546 +}
   1.547 +
   1.548 +void mpz_abs(mpz_t z, mpz_t x)
   1.549 +{     /* set z to the absolute value of x */
   1.550 +      mpz_set(z, x);
   1.551 +      if (z->val < 0) z->val = - z->val;
   1.552 +      return;
   1.553 +}
   1.554 +
   1.555 +void mpz_div(mpz_t q, mpz_t r, mpz_t x, mpz_t y)
   1.556 +{     /* divide x by y, forming quotient q and/or remainder r
   1.557 +         if q = NULL then quotient is not stored; if r = NULL then
   1.558 +         remainder is not stored
   1.559 +         the sign of quotient is determined as in algebra while the
   1.560 +         sign of remainder is the same as the sign of dividend:
   1.561 +         +26 : +7 = +3, remainder is +5
   1.562 +         -26 : +7 = -3, remainder is -5
   1.563 +         +26 : -7 = -3, remainder is +5
   1.564 +         -26 : -7 = +3, remainder is -5 */
   1.565 +      struct mpz_seg dumx, dumy, *ex, *ey, *es, *e;
   1.566 +      int sx, sy, k, nx, ny, n;
   1.567 +      unsigned int t;
   1.568 +      unsigned short *work, *wx, *wy;
   1.569 +      /* divide by zero is not allowed */
   1.570 +      if (y->val == 0)
   1.571 +      {  xassert(y->ptr == NULL);
   1.572 +         xfault("mpz_div: divide by zero not allowed\n");
   1.573 +      }
   1.574 +      /* if [x] = 0 then [q] = [r] = 0 */
   1.575 +      if (x->val == 0)
   1.576 +      {  xassert(x->ptr == NULL);
   1.577 +         if (q != NULL) mpz_set_si(q, 0);
   1.578 +         if (r != NULL) mpz_set_si(r, 0);
   1.579 +         goto done;
   1.580 +      }
   1.581 +      /* special case when both [x] and [y] are in short format */
   1.582 +      if (x->ptr == NULL && y->ptr == NULL)
   1.583 +      {  int xval = x->val, yval = y->val;
   1.584 +         xassert(xval != 0x80000000 && yval != 0x80000000);
   1.585 +         if (q != NULL) mpz_set_si(q, xval / yval);
   1.586 +         if (r != NULL) mpz_set_si(r, xval % yval);
   1.587 +         goto done;
   1.588 +      }
   1.589 +      /* convert [x] to long format, if necessary */
   1.590 +      if (x->ptr == NULL)
   1.591 +      {  xassert(x->val != 0x80000000);
   1.592 +         if (x->val >= 0)
   1.593 +         {  sx = +1;
   1.594 +            t = (unsigned int)(+ x->val);
   1.595 +         }
   1.596 +         else
   1.597 +         {  sx = -1;
   1.598 +            t = (unsigned int)(- x->val);
   1.599 +         }
   1.600 +         ex = &dumx;
   1.601 +         ex->d[0] = (unsigned short)t;
   1.602 +         ex->d[1] = (unsigned short)(t >> 16);
   1.603 +         ex->d[2] = ex->d[3] = ex->d[4] = ex->d[5] = 0;
   1.604 +         ex->next = NULL;
   1.605 +      }
   1.606 +      else
   1.607 +      {  sx = x->val;
   1.608 +         xassert(sx == +1 || sx == -1);
   1.609 +         ex = x->ptr;
   1.610 +      }
   1.611 +      /* convert [y] to long format, if necessary */
   1.612 +      if (y->ptr == NULL)
   1.613 +      {  xassert(y->val != 0x80000000);
   1.614 +         if (y->val >= 0)
   1.615 +         {  sy = +1;
   1.616 +            t = (unsigned int)(+ y->val);
   1.617 +         }
   1.618 +         else
   1.619 +         {  sy = -1;
   1.620 +            t = (unsigned int)(- y->val);
   1.621 +         }
   1.622 +         ey = &dumy;
   1.623 +         ey->d[0] = (unsigned short)t;
   1.624 +         ey->d[1] = (unsigned short)(t >> 16);
   1.625 +         ey->d[2] = ey->d[3] = ey->d[4] = ey->d[5] = 0;
   1.626 +         ey->next = NULL;
   1.627 +      }
   1.628 +      else
   1.629 +      {  sy = y->val;
   1.630 +         xassert(sy == +1 || sy == -1);
   1.631 +         ey = y->ptr;
   1.632 +      }
   1.633 +      /* determine the number of digits of [x] */
   1.634 +      nx = n = 0;
   1.635 +      for (e = ex; e != NULL; e = e->next)
   1.636 +      for (k = 0; k <= 5; k++)
   1.637 +      {  n++;
   1.638 +         if (e->d[k]) nx = n;
   1.639 +      }
   1.640 +      xassert(nx > 0);
   1.641 +      /* determine the number of digits of [y] */
   1.642 +      ny = n = 0;
   1.643 +      for (e = ey; e != NULL; e = e->next)
   1.644 +      for (k = 0; k <= 5; k++)
   1.645 +      {  n++;
   1.646 +         if (e->d[k]) ny = n;
   1.647 +      }
   1.648 +      xassert(ny > 0);
   1.649 +      /* if nx < ny then [q] = 0 and [r] = [x] */
   1.650 +      if (nx < ny)
   1.651 +      {  if (r != NULL) mpz_set(r, x);
   1.652 +         if (q != NULL) mpz_set_si(q, 0);
   1.653 +         goto done;
   1.654 +      }
   1.655 +      /* we need working array containing at least nx+ny+1 places */
   1.656 +      work = gmp_get_work(nx+ny+1);
   1.657 +      /* load digits of [x] */
   1.658 +      wx = &work[0];
   1.659 +      for (n = 0; n < nx; n++) wx[n] = 0;
   1.660 +      for (n = 0, e = ex; e != NULL; e = e->next)
   1.661 +         for (k = 0; k <= 5; k++, n++)
   1.662 +            if (e->d[k]) wx[n] = e->d[k];
   1.663 +      /* load digits of [y] */
   1.664 +      wy = &work[nx+1];
   1.665 +      for (n = 0; n < ny; n++) wy[n] = 0;
   1.666 +      for (n = 0, e = ey; e != NULL; e = e->next)
   1.667 +         for (k = 0; k <= 5; k++, n++)
   1.668 +            if (e->d[k]) wy[n] = e->d[k];
   1.669 +      /* compute quotient and remainder */
   1.670 +      xassert(wy[ny-1] != 0);
   1.671 +      bigdiv(nx-ny, ny, wx, wy);
   1.672 +      /* construct and normalize quotient */
   1.673 +      if (q != NULL)
   1.674 +      {  mpz_set_si(q, 0);
   1.675 +         q->val = sx * sy;
   1.676 +         es = NULL;
   1.677 +         k = 6;
   1.678 +         for (n = ny; n <= nx; n++)
   1.679 +         {  if (k > 5)
   1.680 +            {  e = gmp_get_atom(sizeof(struct mpz_seg));
   1.681 +               e->d[0] = e->d[1] = e->d[2] = 0;
   1.682 +               e->d[3] = e->d[4] = e->d[5] = 0;
   1.683 +               e->next = NULL;
   1.684 +               if (q->ptr == NULL)
   1.685 +                  q->ptr = e;
   1.686 +               else
   1.687 +                  es->next = e;
   1.688 +               es = e;
   1.689 +               k = 0;
   1.690 +            }
   1.691 +            es->d[k++] = wx[n];
   1.692 +         }
   1.693 +         normalize(q);
   1.694 +      }
   1.695 +      /* construct and normalize remainder */
   1.696 +      if (r != NULL)
   1.697 +      {  mpz_set_si(r, 0);
   1.698 +         r->val = sx;
   1.699 +         es = NULL;
   1.700 +         k = 6;
   1.701 +         for (n = 0; n < ny; n++)
   1.702 +         {  if (k > 5)
   1.703 +            {  e = gmp_get_atom(sizeof(struct mpz_seg));
   1.704 +               e->d[0] = e->d[1] = e->d[2] = 0;
   1.705 +               e->d[3] = e->d[4] = e->d[5] = 0;
   1.706 +               e->next = NULL;
   1.707 +               if (r->ptr == NULL)
   1.708 +                  r->ptr = e;
   1.709 +               else
   1.710 +                  es->next = e;
   1.711 +               es = e;
   1.712 +               k = 0;
   1.713 +            }
   1.714 +            es->d[k++] = wx[n];
   1.715 +         }
   1.716 +         normalize(r);
   1.717 +      }
   1.718 +done: return;
   1.719 +}
   1.720 +
   1.721 +void mpz_gcd(mpz_t z, mpz_t x, mpz_t y)
   1.722 +{     /* set z to the greatest common divisor of x and y */
   1.723 +      /* in case of arbitrary integers GCD(x, y) = GCD(|x|, |y|), and,
   1.724 +         in particular, GCD(0, 0) = 0 */
   1.725 +      mpz_t u, v, r;
   1.726 +      mpz_init(u);
   1.727 +      mpz_init(v);
   1.728 +      mpz_init(r);
   1.729 +      mpz_abs(u, x);
   1.730 +      mpz_abs(v, y);
   1.731 +      while (mpz_sgn(v))
   1.732 +      {  mpz_div(NULL, r, u, v);
   1.733 +         mpz_set(u, v);
   1.734 +         mpz_set(v, r);
   1.735 +      }
   1.736 +      mpz_set(z, u);
   1.737 +      mpz_clear(u);
   1.738 +      mpz_clear(v);
   1.739 +      mpz_clear(r);
   1.740 +      return;
   1.741 +}
   1.742 +
   1.743 +int mpz_cmp(mpz_t x, mpz_t y)
   1.744 +{     /* compare x and y; return a positive value if x > y, zero if
   1.745 +         x = y, or a nefative value if x < y */
   1.746 +      static struct mpz_seg zero = { { 0, 0, 0, 0, 0, 0 }, NULL };
   1.747 +      struct mpz_seg dumx, dumy, *ex, *ey;
   1.748 +      int cc, sx, sy, k;
   1.749 +      unsigned int t;
   1.750 +      if (x == y)
   1.751 +      {  cc = 0;
   1.752 +         goto done;
   1.753 +      }
   1.754 +      /* special case when both [x] and [y] are in short format */
   1.755 +      if (x->ptr == NULL && y->ptr == NULL)
   1.756 +      {  int xval = x->val, yval = y->val;
   1.757 +         xassert(xval != 0x80000000 && yval != 0x80000000);
   1.758 +         cc = (xval > yval ? +1 : xval < yval ? -1 : 0);
   1.759 +         goto done;
   1.760 +      }
   1.761 +      /* special case when [x] and [y] have different signs */
   1.762 +      if (x->val > 0 && y->val <= 0 || x->val == 0 && y->val < 0)
   1.763 +      {  cc = +1;
   1.764 +         goto done;
   1.765 +      }
   1.766 +      if (x->val < 0 && y->val >= 0 || x->val == 0 && y->val > 0)
   1.767 +      {  cc = -1;
   1.768 +         goto done;
   1.769 +      }
   1.770 +      /* convert [x] to long format, if necessary */
   1.771 +      if (x->ptr == NULL)
   1.772 +      {  xassert(x->val != 0x80000000);
   1.773 +         if (x->val >= 0)
   1.774 +         {  sx = +1;
   1.775 +            t = (unsigned int)(+ x->val);
   1.776 +         }
   1.777 +         else
   1.778 +         {  sx = -1;
   1.779 +            t = (unsigned int)(- x->val);
   1.780 +         }
   1.781 +         ex = &dumx;
   1.782 +         ex->d[0] = (unsigned short)t;
   1.783 +         ex->d[1] = (unsigned short)(t >> 16);
   1.784 +         ex->d[2] = ex->d[3] = ex->d[4] = ex->d[5] = 0;
   1.785 +         ex->next = NULL;
   1.786 +      }
   1.787 +      else
   1.788 +      {  sx = x->val;
   1.789 +         xassert(sx == +1 || sx == -1);
   1.790 +         ex = x->ptr;
   1.791 +      }
   1.792 +      /* convert [y] to long format, if necessary */
   1.793 +      if (y->ptr == NULL)
   1.794 +      {  xassert(y->val != 0x80000000);
   1.795 +         if (y->val >= 0)
   1.796 +         {  sy = +1;
   1.797 +            t = (unsigned int)(+ y->val);
   1.798 +         }
   1.799 +         else
   1.800 +         {  sy = -1;
   1.801 +            t = (unsigned int)(- y->val);
   1.802 +         }
   1.803 +         ey = &dumy;
   1.804 +         ey->d[0] = (unsigned short)t;
   1.805 +         ey->d[1] = (unsigned short)(t >> 16);
   1.806 +         ey->d[2] = ey->d[3] = ey->d[4] = ey->d[5] = 0;
   1.807 +         ey->next = NULL;
   1.808 +      }
   1.809 +      else
   1.810 +      {  sy = y->val;
   1.811 +         xassert(sy == +1 || sy == -1);
   1.812 +         ey = y->ptr;
   1.813 +      }
   1.814 +      /* main fragment */
   1.815 +      xassert(sx > 0 && sy > 0 || sx < 0 && sy < 0);
   1.816 +      cc = 0;
   1.817 +      for (; ex || ey; ex = ex->next, ey = ey->next)
   1.818 +      {  if (ex == NULL) ex = &zero;
   1.819 +         if (ey == NULL) ey = &zero;
   1.820 +         for (k = 0; k <= 5; k++)
   1.821 +         {  if (ex->d[k] > ey->d[k]) cc = +1;
   1.822 +            if (ex->d[k] < ey->d[k]) cc = -1;
   1.823 +         }
   1.824 +      }
   1.825 +      if (sx < 0) cc = - cc;
   1.826 +done: return cc;
   1.827 +}
   1.828 +
   1.829 +int mpz_sgn(mpz_t x)
   1.830 +{     /* return +1 if x > 0, 0 if x = 0, and -1 if x < 0 */
   1.831 +      int s;
   1.832 +      s = (x->val > 0 ? +1 : x->val < 0 ? -1 : 0);
   1.833 +      return s;
   1.834 +}
   1.835 +
   1.836 +int mpz_out_str(void *_fp, int base, mpz_t x)
   1.837 +{     /* output x on stream fp, as a string in given base; the base
   1.838 +         may vary from 2 to 36;
   1.839 +         return the number of bytes written, or if an error occurred,
   1.840 +         return 0 */
   1.841 +      FILE *fp = _fp;
   1.842 +      mpz_t b, y, r;
   1.843 +      int n, j, nwr = 0;
   1.844 +      unsigned char *d;
   1.845 +      static char *set = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ";
   1.846 +      if (!(2 <= base && base <= 36))
   1.847 +         xfault("mpz_out_str: base = %d; invalid base\n", base);
   1.848 +      mpz_init(b);
   1.849 +      mpz_set_si(b, base);
   1.850 +      mpz_init(y);
   1.851 +      mpz_init(r);
   1.852 +      /* determine the number of digits */
   1.853 +      mpz_abs(y, x);
   1.854 +      for (n = 0; mpz_sgn(y) != 0; n++)
   1.855 +         mpz_div(y, NULL, y, b);
   1.856 +      if (n == 0) n = 1;
   1.857 +      /* compute the digits */
   1.858 +      d = xmalloc(n);
   1.859 +      mpz_abs(y, x);
   1.860 +      for (j = 0; j < n; j++)
   1.861 +      {  mpz_div(y, r, y, b);
   1.862 +         xassert(0 <= r->val && r->val < base && r->ptr == NULL);
   1.863 +         d[j] = (unsigned char)r->val;
   1.864 +      }
   1.865 +      /* output the integer to the stream */
   1.866 +      if (fp == NULL) fp = stdout;
   1.867 +      if (mpz_sgn(x) < 0)
   1.868 +         fputc('-', fp), nwr++;
   1.869 +      for (j = n-1; j >= 0; j--)
   1.870 +         fputc(set[d[j]], fp), nwr++;
   1.871 +      if (ferror(fp)) nwr = 0;
   1.872 +      mpz_clear(b);
   1.873 +      mpz_clear(y);
   1.874 +      mpz_clear(r);
   1.875 +      xfree(d);
   1.876 +      return nwr;
   1.877 +}
   1.878 +
   1.879 +/*====================================================================*/
   1.880 +
   1.881 +mpq_t _mpq_init(void)
   1.882 +{     /* initialize x, and set its value to 0/1 */
   1.883 +      mpq_t x;
   1.884 +      x = gmp_get_atom(sizeof(struct mpq));
   1.885 +      x->p.val = 0;
   1.886 +      x->p.ptr = NULL;
   1.887 +      x->q.val = 1;
   1.888 +      x->q.ptr = NULL;
   1.889 +      return x;
   1.890 +}
   1.891 +
   1.892 +void mpq_clear(mpq_t x)
   1.893 +{     /* free the space occupied by x */
   1.894 +      mpz_set_si(&x->p, 0);
   1.895 +      xassert(x->p.ptr == NULL);
   1.896 +      mpz_set_si(&x->q, 0);
   1.897 +      xassert(x->q.ptr == NULL);
   1.898 +      /* free the number descriptor */
   1.899 +      gmp_free_atom(x, sizeof(struct mpq));
   1.900 +      return;
   1.901 +}
   1.902 +
   1.903 +void mpq_canonicalize(mpq_t x)
   1.904 +{     /* remove any factors that are common to the numerator and
   1.905 +         denominator of x, and make the denominator positive */
   1.906 +      mpz_t f;
   1.907 +      xassert(x->q.val != 0);
   1.908 +      if (x->q.val < 0)
   1.909 +      {  mpz_neg(&x->p, &x->p);
   1.910 +         mpz_neg(&x->q, &x->q);
   1.911 +      }
   1.912 +      mpz_init(f);
   1.913 +      mpz_gcd(f, &x->p, &x->q);
   1.914 +      if (!(f->val == 1 && f->ptr == NULL))
   1.915 +      {  mpz_div(&x->p, NULL, &x->p, f);
   1.916 +         mpz_div(&x->q, NULL, &x->q, f);
   1.917 +      }
   1.918 +      mpz_clear(f);
   1.919 +      return;
   1.920 +}
   1.921 +
   1.922 +void mpq_set(mpq_t z, mpq_t x)
   1.923 +{     /* set the value of z from x */
   1.924 +      if (z != x)
   1.925 +      {  mpz_set(&z->p, &x->p);
   1.926 +         mpz_set(&z->q, &x->q);
   1.927 +      }
   1.928 +      return;
   1.929 +}
   1.930 +
   1.931 +void mpq_set_si(mpq_t x, int p, unsigned int q)
   1.932 +{     /* set the value of x to p/q */
   1.933 +      if (q == 0)
   1.934 +         xfault("mpq_set_si: zero denominator not allowed\n");
   1.935 +      mpz_set_si(&x->p, p);
   1.936 +      xassert(q <= 0x7FFFFFFF);
   1.937 +      mpz_set_si(&x->q, q);
   1.938 +      return;
   1.939 +}
   1.940 +
   1.941 +double mpq_get_d(mpq_t x)
   1.942 +{     /* convert x to a double, truncating if necessary */
   1.943 +      int np, nq;
   1.944 +      double p, q;
   1.945 +      p = mpz_get_d_2exp(&np, &x->p);
   1.946 +      q = mpz_get_d_2exp(&nq, &x->q);
   1.947 +      return ldexp(p / q, np - nq);
   1.948 +}
   1.949 +
   1.950 +void mpq_set_d(mpq_t x, double val)
   1.951 +{     /* set x to val; there is no rounding, the conversion is exact */
   1.952 +      int s, n, d, j;
   1.953 +      double f;
   1.954 +      mpz_t temp;
   1.955 +      xassert(-DBL_MAX <= val && val <= +DBL_MAX);
   1.956 +      mpq_set_si(x, 0, 1);
   1.957 +      if (val > 0.0)
   1.958 +         s = +1;
   1.959 +      else if (val < 0.0)
   1.960 +         s = -1;
   1.961 +      else
   1.962 +         goto done;
   1.963 +      f = frexp(fabs(val), &n);
   1.964 +      /* |val| = f * 2^n, where 0.5 <= f < 1.0 */
   1.965 +      mpz_init(temp);
   1.966 +      while (f != 0.0)
   1.967 +      {  f *= 16.0, n -= 4;
   1.968 +         d = (int)f;
   1.969 +         xassert(0 <= d && d <= 15);
   1.970 +         f -= (double)d;
   1.971 +         /* x := 16 * x + d */
   1.972 +         mpz_set_si(temp, 16);
   1.973 +         mpz_mul(&x->p, &x->p, temp);
   1.974 +         mpz_set_si(temp, d);
   1.975 +         mpz_add(&x->p, &x->p, temp);
   1.976 +      }
   1.977 +      mpz_clear(temp);
   1.978 +      /* x := x * 2^n */
   1.979 +      if (n > 0)
   1.980 +      {  for (j = 1; j <= n; j++)
   1.981 +            mpz_add(&x->p, &x->p, &x->p);
   1.982 +      }
   1.983 +      else if (n < 0)
   1.984 +      {  for (j = 1; j <= -n; j++)
   1.985 +            mpz_add(&x->q, &x->q, &x->q);
   1.986 +         mpq_canonicalize(x);
   1.987 +      }
   1.988 +      if (s < 0) mpq_neg(x, x);
   1.989 +done: return;
   1.990 +}
   1.991 +
   1.992 +void mpq_add(mpq_t z, mpq_t x, mpq_t y)
   1.993 +{     /* set z to x + y */
   1.994 +      mpz_t p, q;
   1.995 +      mpz_init(p);
   1.996 +      mpz_init(q);
   1.997 +      mpz_mul(p, &x->p, &y->q);
   1.998 +      mpz_mul(q, &x->q, &y->p);
   1.999 +      mpz_add(p, p, q);
  1.1000 +      mpz_mul(q, &x->q, &y->q);
  1.1001 +      mpz_set(&z->p, p);
  1.1002 +      mpz_set(&z->q, q);
  1.1003 +      mpz_clear(p);
  1.1004 +      mpz_clear(q);
  1.1005 +      mpq_canonicalize(z);
  1.1006 +      return;
  1.1007 +}
  1.1008 +
  1.1009 +void mpq_sub(mpq_t z, mpq_t x, mpq_t y)
  1.1010 +{     /* set z to x - y */
  1.1011 +      mpz_t p, q;
  1.1012 +      mpz_init(p);
  1.1013 +      mpz_init(q);
  1.1014 +      mpz_mul(p, &x->p, &y->q);
  1.1015 +      mpz_mul(q, &x->q, &y->p);
  1.1016 +      mpz_sub(p, p, q);
  1.1017 +      mpz_mul(q, &x->q, &y->q);
  1.1018 +      mpz_set(&z->p, p);
  1.1019 +      mpz_set(&z->q, q);
  1.1020 +      mpz_clear(p);
  1.1021 +      mpz_clear(q);
  1.1022 +      mpq_canonicalize(z);
  1.1023 +      return;
  1.1024 +}
  1.1025 +
  1.1026 +void mpq_mul(mpq_t z, mpq_t x, mpq_t y)
  1.1027 +{     /* set z to x * y */
  1.1028 +      mpz_mul(&z->p, &x->p, &y->p);
  1.1029 +      mpz_mul(&z->q, &x->q, &y->q);
  1.1030 +      mpq_canonicalize(z);
  1.1031 +      return;
  1.1032 +}
  1.1033 +
  1.1034 +void mpq_div(mpq_t z, mpq_t x, mpq_t y)
  1.1035 +{     /* set z to x / y */
  1.1036 +      mpz_t p, q;
  1.1037 +      if (mpq_sgn(y) == 0)
  1.1038 +         xfault("mpq_div: zero divisor not allowed\n");
  1.1039 +      mpz_init(p);
  1.1040 +      mpz_init(q);
  1.1041 +      mpz_mul(p, &x->p, &y->q);
  1.1042 +      mpz_mul(q, &x->q, &y->p);
  1.1043 +      mpz_set(&z->p, p);
  1.1044 +      mpz_set(&z->q, q);
  1.1045 +      mpz_clear(p);
  1.1046 +      mpz_clear(q);
  1.1047 +      mpq_canonicalize(z);
  1.1048 +      return;
  1.1049 +}
  1.1050 +
  1.1051 +void mpq_neg(mpq_t z, mpq_t x)
  1.1052 +{     /* set z to 0 - x */
  1.1053 +      mpq_set(z, x);
  1.1054 +      mpz_neg(&z->p, &z->p);
  1.1055 +      return;
  1.1056 +}
  1.1057 +
  1.1058 +void mpq_abs(mpq_t z, mpq_t x)
  1.1059 +{     /* set z to the absolute value of x */
  1.1060 +      mpq_set(z, x);
  1.1061 +      mpz_abs(&z->p, &z->p);
  1.1062 +      xassert(mpz_sgn(&x->q) > 0);
  1.1063 +      return;
  1.1064 +}
  1.1065 +
  1.1066 +int mpq_cmp(mpq_t x, mpq_t y)
  1.1067 +{     /* compare x and y; return a positive value if x > y, zero if
  1.1068 +         x = y, or a nefative value if x < y */
  1.1069 +      mpq_t temp;
  1.1070 +      int s;
  1.1071 +      mpq_init(temp);
  1.1072 +      mpq_sub(temp, x, y);
  1.1073 +      s = mpq_sgn(temp);
  1.1074 +      mpq_clear(temp);
  1.1075 +      return s;
  1.1076 +}
  1.1077 +
  1.1078 +int mpq_sgn(mpq_t x)
  1.1079 +{     /* return +1 if x > 0, 0 if x = 0, and -1 if x < 0 */
  1.1080 +      int s;
  1.1081 +      s = mpz_sgn(&x->p);
  1.1082 +      xassert(mpz_sgn(&x->q) > 0);
  1.1083 +      return s;
  1.1084 +}
  1.1085 +
  1.1086 +int mpq_out_str(void *_fp, int base, mpq_t x)
  1.1087 +{     /* output x on stream fp, as a string in given base; the base
  1.1088 +         may vary from 2 to 36; output is in the form 'num/den' or if
  1.1089 +         the denominator is 1 then just 'num';
  1.1090 +         if the parameter fp is a null pointer, stdout is assumed;
  1.1091 +         return the number of bytes written, or if an error occurred,
  1.1092 +         return 0 */
  1.1093 +      FILE *fp = _fp;
  1.1094 +      int nwr;
  1.1095 +      if (!(2 <= base && base <= 36))
  1.1096 +         xfault("mpq_out_str: base = %d; invalid base\n", base);
  1.1097 +      if (fp == NULL) fp = stdout;
  1.1098 +      nwr = mpz_out_str(fp, base, &x->p);
  1.1099 +      if (x->q.val == 1 && x->q.ptr == NULL)
  1.1100 +         ;
  1.1101 +      else
  1.1102 +      {  fputc('/', fp), nwr++;
  1.1103 +         nwr += mpz_out_str(fp, base, &x->q);
  1.1104 +      }
  1.1105 +      if (ferror(fp)) nwr = 0;
  1.1106 +      return nwr;
  1.1107 +}
  1.1108 +
  1.1109 +#endif
  1.1110 +
  1.1111 +/* eof */