1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/src/glpgmp.c Mon Dec 06 13:09:21 2010 +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 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 */