lemon-project-template-glpk

annotate deps/glpk/src/glpgmp.c @ 9:33de93886c88

find changesets by author, revision, files, or words in the commit message
Import GLPK 4.47
author Alpar Juttner <alpar@cs.elte.hu>
date Sun, 06 Nov 2011 20:59:10 +0100
parents
children
rev   line source
alpar@9 1 /* glpgmp.c */
alpar@9 2
alpar@9 3 /***********************************************************************
alpar@9 4 * This code is part of GLPK (GNU Linear Programming Kit).
alpar@9 5 *
alpar@9 6 * Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
alpar@9 7 * 2009, 2010, 2011 Andrew Makhorin, Department for Applied Informatics,
alpar@9 8 * Moscow Aviation Institute, Moscow, Russia. All rights reserved.
alpar@9 9 * E-mail: <mao@gnu.org>.
alpar@9 10 *
alpar@9 11 * GLPK is free software: you can redistribute it and/or modify it
alpar@9 12 * under the terms of the GNU General Public License as published by
alpar@9 13 * the Free Software Foundation, either version 3 of the License, or
alpar@9 14 * (at your option) any later version.
alpar@9 15 *
alpar@9 16 * GLPK is distributed in the hope that it will be useful, but WITHOUT
alpar@9 17 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
alpar@9 18 * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
alpar@9 19 * License for more details.
alpar@9 20 *
alpar@9 21 * You should have received a copy of the GNU General Public License
alpar@9 22 * along with GLPK. If not, see <http://www.gnu.org/licenses/>.
alpar@9 23 ***********************************************************************/
alpar@9 24
alpar@9 25 #define _GLPSTD_STDIO
alpar@9 26 #include "glpdmp.h"
alpar@9 27 #include "glpgmp.h"
alpar@9 28 #define xfault xerror
alpar@9 29
alpar@9 30 #ifdef HAVE_GMP /* use GNU MP bignum library */
alpar@9 31
alpar@9 32 int gmp_pool_count(void) { return 0; }
alpar@9 33
alpar@9 34 void gmp_free_mem(void) { return; }
alpar@9 35
alpar@9 36 #else /* use GLPK bignum module */
alpar@9 37
alpar@9 38 static DMP *gmp_pool = NULL;
alpar@9 39 static int gmp_size = 0;
alpar@9 40 static unsigned short *gmp_work = NULL;
alpar@9 41
alpar@9 42 void *gmp_get_atom(int size)
alpar@9 43 { if (gmp_pool == NULL)
alpar@9 44 gmp_pool = dmp_create_pool();
alpar@9 45 return dmp_get_atom(gmp_pool, size);
alpar@9 46 }
alpar@9 47
alpar@9 48 void gmp_free_atom(void *ptr, int size)
alpar@9 49 { xassert(gmp_pool != NULL);
alpar@9 50 dmp_free_atom(gmp_pool, ptr, size);
alpar@9 51 return;
alpar@9 52 }
alpar@9 53
alpar@9 54 int gmp_pool_count(void)
alpar@9 55 { if (gmp_pool == NULL)
alpar@9 56 return 0;
alpar@9 57 else
alpar@9 58 return dmp_in_use(gmp_pool).lo;
alpar@9 59 }
alpar@9 60
alpar@9 61 unsigned short *gmp_get_work(int size)
alpar@9 62 { xassert(size > 0);
alpar@9 63 if (gmp_size < size)
alpar@9 64 { if (gmp_size == 0)
alpar@9 65 { xassert(gmp_work == NULL);
alpar@9 66 gmp_size = 100;
alpar@9 67 }
alpar@9 68 else
alpar@9 69 { xassert(gmp_work != NULL);
alpar@9 70 xfree(gmp_work);
alpar@9 71 }
alpar@9 72 while (gmp_size < size) gmp_size += gmp_size;
alpar@9 73 gmp_work = xcalloc(gmp_size, sizeof(unsigned short));
alpar@9 74 }
alpar@9 75 return gmp_work;
alpar@9 76 }
alpar@9 77
alpar@9 78 void gmp_free_mem(void)
alpar@9 79 { if (gmp_pool != NULL) dmp_delete_pool(gmp_pool);
alpar@9 80 if (gmp_work != NULL) xfree(gmp_work);
alpar@9 81 gmp_pool = NULL;
alpar@9 82 gmp_size = 0;
alpar@9 83 gmp_work = NULL;
alpar@9 84 return;
alpar@9 85 }
alpar@9 86
alpar@9 87 /*====================================================================*/
alpar@9 88
alpar@9 89 mpz_t _mpz_init(void)
alpar@9 90 { /* initialize x, and set its value to 0 */
alpar@9 91 mpz_t x;
alpar@9 92 x = gmp_get_atom(sizeof(struct mpz));
alpar@9 93 x->val = 0;
alpar@9 94 x->ptr = NULL;
alpar@9 95 return x;
alpar@9 96 }
alpar@9 97
alpar@9 98 void mpz_clear(mpz_t x)
alpar@9 99 { /* free the space occupied by x */
alpar@9 100 mpz_set_si(x, 0);
alpar@9 101 xassert(x->ptr == NULL);
alpar@9 102 /* free the number descriptor */
alpar@9 103 gmp_free_atom(x, sizeof(struct mpz));
alpar@9 104 return;
alpar@9 105 }
alpar@9 106
alpar@9 107 void mpz_set(mpz_t z, mpz_t x)
alpar@9 108 { /* set the value of z from x */
alpar@9 109 struct mpz_seg *e, *ee, *es;
alpar@9 110 if (z != x)
alpar@9 111 { mpz_set_si(z, 0);
alpar@9 112 z->val = x->val;
alpar@9 113 xassert(z->ptr == NULL);
alpar@9 114 for (e = x->ptr, es = NULL; e != NULL; e = e->next)
alpar@9 115 { ee = gmp_get_atom(sizeof(struct mpz_seg));
alpar@9 116 memcpy(ee->d, e->d, 12);
alpar@9 117 ee->next = NULL;
alpar@9 118 if (z->ptr == NULL)
alpar@9 119 z->ptr = ee;
alpar@9 120 else
alpar@9 121 es->next = ee;
alpar@9 122 es = ee;
alpar@9 123 }
alpar@9 124 }
alpar@9 125 return;
alpar@9 126 }
alpar@9 127
alpar@9 128 void mpz_set_si(mpz_t x, int val)
alpar@9 129 { /* set the value of x to val */
alpar@9 130 struct mpz_seg *e;
alpar@9 131 /* free existing segments, if any */
alpar@9 132 while (x->ptr != NULL)
alpar@9 133 { e = x->ptr;
alpar@9 134 x->ptr = e->next;
alpar@9 135 gmp_free_atom(e, sizeof(struct mpz_seg));
alpar@9 136 }
alpar@9 137 /* assign new value */
alpar@9 138 if (val == 0x80000000)
alpar@9 139 { /* long format is needed */
alpar@9 140 x->val = -1;
alpar@9 141 x->ptr = e = gmp_get_atom(sizeof(struct mpz_seg));
alpar@9 142 memset(e->d, 0, 12);
alpar@9 143 e->d[1] = 0x8000;
alpar@9 144 e->next = NULL;
alpar@9 145 }
alpar@9 146 else
alpar@9 147 { /* short format is enough */
alpar@9 148 x->val = val;
alpar@9 149 }
alpar@9 150 return;
alpar@9 151 }
alpar@9 152
alpar@9 153 double mpz_get_d(mpz_t x)
alpar@9 154 { /* convert x to a double, truncating if necessary */
alpar@9 155 struct mpz_seg *e;
alpar@9 156 int j;
alpar@9 157 double val, deg;
alpar@9 158 if (x->ptr == NULL)
alpar@9 159 val = (double)x->val;
alpar@9 160 else
alpar@9 161 { xassert(x->val != 0);
alpar@9 162 val = 0.0;
alpar@9 163 deg = 1.0;
alpar@9 164 for (e = x->ptr; e != NULL; e = e->next)
alpar@9 165 { for (j = 0; j <= 5; j++)
alpar@9 166 { val += deg * (double)((int)e->d[j]);
alpar@9 167 deg *= 65536.0;
alpar@9 168 }
alpar@9 169 }
alpar@9 170 if (x->val < 0) val = - val;
alpar@9 171 }
alpar@9 172 return val;
alpar@9 173 }
alpar@9 174
alpar@9 175 double mpz_get_d_2exp(int *exp, mpz_t x)
alpar@9 176 { /* convert x to a double, truncating if necessary (i.e. rounding
alpar@9 177 towards zero), and returning the exponent separately;
alpar@9 178 the return value is in the range 0.5 <= |d| < 1 and the
alpar@9 179 exponent is stored to *exp; d*2^exp is the (truncated) x value;
alpar@9 180 if x is zero, the return is 0.0 and 0 is stored to *exp;
alpar@9 181 this is similar to the standard C frexp function */
alpar@9 182 struct mpz_seg *e;
alpar@9 183 int j, n, n1;
alpar@9 184 double val;
alpar@9 185 if (x->ptr == NULL)
alpar@9 186 val = (double)x->val, n = 0;
alpar@9 187 else
alpar@9 188 { xassert(x->val != 0);
alpar@9 189 val = 0.0, n = 0;
alpar@9 190 for (e = x->ptr; e != NULL; e = e->next)
alpar@9 191 { for (j = 0; j <= 5; j++)
alpar@9 192 { val += (double)((int)e->d[j]);
alpar@9 193 val /= 65536.0, n += 16;
alpar@9 194 }
alpar@9 195 }
alpar@9 196 if (x->val < 0) val = - val;
alpar@9 197 }
alpar@9 198 val = frexp(val, &n1);
alpar@9 199 *exp = n + n1;
alpar@9 200 return val;
alpar@9 201 }
alpar@9 202
alpar@9 203 void mpz_swap(mpz_t x, mpz_t y)
alpar@9 204 { /* swap the values x and y efficiently */
alpar@9 205 int val;
alpar@9 206 void *ptr;
alpar@9 207 val = x->val, ptr = x->ptr;
alpar@9 208 x->val = y->val, x->ptr = y->ptr;
alpar@9 209 y->val = val, y->ptr = ptr;
alpar@9 210 return;
alpar@9 211 }
alpar@9 212
alpar@9 213 static void normalize(mpz_t x)
alpar@9 214 { /* normalize integer x that includes removing non-significant
alpar@9 215 (leading) zeros and converting to short format, if possible */
alpar@9 216 struct mpz_seg *es, *e;
alpar@9 217 /* if the integer is in short format, it remains unchanged */
alpar@9 218 if (x->ptr == NULL)
alpar@9 219 { xassert(x->val != 0x80000000);
alpar@9 220 goto done;
alpar@9 221 }
alpar@9 222 xassert(x->val == +1 || x->val == -1);
alpar@9 223 /* find the last (most significant) non-zero segment */
alpar@9 224 es = NULL;
alpar@9 225 for (e = x->ptr; e != NULL; e = e->next)
alpar@9 226 { if (e->d[0] || e->d[1] || e->d[2] ||
alpar@9 227 e->d[3] || e->d[4] || e->d[5]) es = e;
alpar@9 228 }
alpar@9 229 /* if all segments contain zeros, the integer is zero */
alpar@9 230 if (es == NULL)
alpar@9 231 { mpz_set_si(x, 0);
alpar@9 232 goto done;
alpar@9 233 }
alpar@9 234 /* remove non-significant (leading) zero segments */
alpar@9 235 while (es->next != NULL)
alpar@9 236 { e = es->next;
alpar@9 237 es->next = e->next;
alpar@9 238 gmp_free_atom(e, sizeof(struct mpz_seg));
alpar@9 239 }
alpar@9 240 /* convert the integer to short format, if possible */
alpar@9 241 e = x->ptr;
alpar@9 242 if (e->next == NULL && e->d[1] <= 0x7FFF &&
alpar@9 243 !e->d[2] && !e->d[3] && !e->d[4] && !e->d[5])
alpar@9 244 { int val;
alpar@9 245 val = (int)e->d[0] + ((int)e->d[1] << 16);
alpar@9 246 if (x->val < 0) val = - val;
alpar@9 247 mpz_set_si(x, val);
alpar@9 248 }
alpar@9 249 done: return;
alpar@9 250 }
alpar@9 251
alpar@9 252 void mpz_add(mpz_t z, mpz_t x, mpz_t y)
alpar@9 253 { /* set z to x + y */
alpar@9 254 static struct mpz_seg zero = { { 0, 0, 0, 0, 0, 0 }, NULL };
alpar@9 255 struct mpz_seg dumx, dumy, *ex, *ey, *ez, *es, *ee;
alpar@9 256 int k, sx, sy, sz;
alpar@9 257 unsigned int t;
alpar@9 258 /* if [x] = 0 then [z] = [y] */
alpar@9 259 if (x->val == 0)
alpar@9 260 { xassert(x->ptr == NULL);
alpar@9 261 mpz_set(z, y);
alpar@9 262 goto done;
alpar@9 263 }
alpar@9 264 /* if [y] = 0 then [z] = [x] */
alpar@9 265 if (y->val == 0)
alpar@9 266 { xassert(y->ptr == NULL);
alpar@9 267 mpz_set(z, x);
alpar@9 268 goto done;
alpar@9 269 }
alpar@9 270 /* special case when both [x] and [y] are in short format */
alpar@9 271 if (x->ptr == NULL && y->ptr == NULL)
alpar@9 272 { int xval = x->val, yval = y->val, zval = x->val + y->val;
alpar@9 273 xassert(xval != 0x80000000 && yval != 0x80000000);
alpar@9 274 if (!(xval > 0 && yval > 0 && zval <= 0 ||
alpar@9 275 xval < 0 && yval < 0 && zval >= 0))
alpar@9 276 { mpz_set_si(z, zval);
alpar@9 277 goto done;
alpar@9 278 }
alpar@9 279 }
alpar@9 280 /* convert [x] to long format, if necessary */
alpar@9 281 if (x->ptr == NULL)
alpar@9 282 { xassert(x->val != 0x80000000);
alpar@9 283 if (x->val >= 0)
alpar@9 284 { sx = +1;
alpar@9 285 t = (unsigned int)(+ x->val);
alpar@9 286 }
alpar@9 287 else
alpar@9 288 { sx = -1;
alpar@9 289 t = (unsigned int)(- x->val);
alpar@9 290 }
alpar@9 291 ex = &dumx;
alpar@9 292 ex->d[0] = (unsigned short)t;
alpar@9 293 ex->d[1] = (unsigned short)(t >> 16);
alpar@9 294 ex->d[2] = ex->d[3] = ex->d[4] = ex->d[5] = 0;
alpar@9 295 ex->next = NULL;
alpar@9 296 }
alpar@9 297 else
alpar@9 298 { sx = x->val;
alpar@9 299 xassert(sx == +1 || sx == -1);
alpar@9 300 ex = x->ptr;
alpar@9 301 }
alpar@9 302 /* convert [y] to long format, if necessary */
alpar@9 303 if (y->ptr == NULL)
alpar@9 304 { xassert(y->val != 0x80000000);
alpar@9 305 if (y->val >= 0)
alpar@9 306 { sy = +1;
alpar@9 307 t = (unsigned int)(+ y->val);
alpar@9 308 }
alpar@9 309 else
alpar@9 310 { sy = -1;
alpar@9 311 t = (unsigned int)(- y->val);
alpar@9 312 }
alpar@9 313 ey = &dumy;
alpar@9 314 ey->d[0] = (unsigned short)t;
alpar@9 315 ey->d[1] = (unsigned short)(t >> 16);
alpar@9 316 ey->d[2] = ey->d[3] = ey->d[4] = ey->d[5] = 0;
alpar@9 317 ey->next = NULL;
alpar@9 318 }
alpar@9 319 else
alpar@9 320 { sy = y->val;
alpar@9 321 xassert(sy == +1 || sy == -1);
alpar@9 322 ey = y->ptr;
alpar@9 323 }
alpar@9 324 /* main fragment */
alpar@9 325 sz = sx;
alpar@9 326 ez = es = NULL;
alpar@9 327 if (sx > 0 && sy > 0 || sx < 0 && sy < 0)
alpar@9 328 { /* [x] and [y] have identical signs -- addition */
alpar@9 329 t = 0;
alpar@9 330 for (; ex || ey; ex = ex->next, ey = ey->next)
alpar@9 331 { if (ex == NULL) ex = &zero;
alpar@9 332 if (ey == NULL) ey = &zero;
alpar@9 333 ee = gmp_get_atom(sizeof(struct mpz_seg));
alpar@9 334 for (k = 0; k <= 5; k++)
alpar@9 335 { t += (unsigned int)ex->d[k];
alpar@9 336 t += (unsigned int)ey->d[k];
alpar@9 337 ee->d[k] = (unsigned short)t;
alpar@9 338 t >>= 16;
alpar@9 339 }
alpar@9 340 ee->next = NULL;
alpar@9 341 if (ez == NULL)
alpar@9 342 ez = ee;
alpar@9 343 else
alpar@9 344 es->next = ee;
alpar@9 345 es = ee;
alpar@9 346 }
alpar@9 347 if (t)
alpar@9 348 { /* overflow -- one extra digit is needed */
alpar@9 349 ee = gmp_get_atom(sizeof(struct mpz_seg));
alpar@9 350 ee->d[0] = 1;
alpar@9 351 ee->d[1] = ee->d[2] = ee->d[3] = ee->d[4] = ee->d[5] = 0;
alpar@9 352 ee->next = NULL;
alpar@9 353 xassert(es != NULL);
alpar@9 354 es->next = ee;
alpar@9 355 }
alpar@9 356 }
alpar@9 357 else
alpar@9 358 { /* [x] and [y] have different signs -- subtraction */
alpar@9 359 t = 1;
alpar@9 360 for (; ex || ey; ex = ex->next, ey = ey->next)
alpar@9 361 { if (ex == NULL) ex = &zero;
alpar@9 362 if (ey == NULL) ey = &zero;
alpar@9 363 ee = gmp_get_atom(sizeof(struct mpz_seg));
alpar@9 364 for (k = 0; k <= 5; k++)
alpar@9 365 { t += (unsigned int)ex->d[k];
alpar@9 366 t += (0xFFFF - (unsigned int)ey->d[k]);
alpar@9 367 ee->d[k] = (unsigned short)t;
alpar@9 368 t >>= 16;
alpar@9 369 }
alpar@9 370 ee->next = NULL;
alpar@9 371 if (ez == NULL)
alpar@9 372 ez = ee;
alpar@9 373 else
alpar@9 374 es->next = ee;
alpar@9 375 es = ee;
alpar@9 376 }
alpar@9 377 if (!t)
alpar@9 378 { /* |[x]| < |[y]| -- result in complement coding */
alpar@9 379 sz = - sz;
alpar@9 380 t = 1;
alpar@9 381 for (ee = ez; ee != NULL; ee = ee->next)
alpar@9 382 for (k = 0; k <= 5; k++)
alpar@9 383 { t += (0xFFFF - (unsigned int)ee->d[k]);
alpar@9 384 ee->d[k] = (unsigned short)t;
alpar@9 385 t >>= 16;
alpar@9 386 }
alpar@9 387 }
alpar@9 388 }
alpar@9 389 /* contruct and normalize result */
alpar@9 390 mpz_set_si(z, 0);
alpar@9 391 z->val = sz;
alpar@9 392 z->ptr = ez;
alpar@9 393 normalize(z);
alpar@9 394 done: return;
alpar@9 395 }
alpar@9 396
alpar@9 397 void mpz_sub(mpz_t z, mpz_t x, mpz_t y)
alpar@9 398 { /* set z to x - y */
alpar@9 399 if (x == y)
alpar@9 400 mpz_set_si(z, 0);
alpar@9 401 else
alpar@9 402 { y->val = - y->val;
alpar@9 403 mpz_add(z, x, y);
alpar@9 404 if (y != z) y->val = - y->val;
alpar@9 405 }
alpar@9 406 return;
alpar@9 407 }
alpar@9 408
alpar@9 409 void mpz_mul(mpz_t z, mpz_t x, mpz_t y)
alpar@9 410 { /* set z to x * y */
alpar@9 411 struct mpz_seg dumx, dumy, *ex, *ey, *es, *e;
alpar@9 412 int sx, sy, k, nx, ny, n;
alpar@9 413 unsigned int t;
alpar@9 414 unsigned short *work, *wx, *wy;
alpar@9 415 /* if [x] = 0 then [z] = 0 */
alpar@9 416 if (x->val == 0)
alpar@9 417 { xassert(x->ptr == NULL);
alpar@9 418 mpz_set_si(z, 0);
alpar@9 419 goto done;
alpar@9 420 }
alpar@9 421 /* if [y] = 0 then [z] = 0 */
alpar@9 422 if (y->val == 0)
alpar@9 423 { xassert(y->ptr == NULL);
alpar@9 424 mpz_set_si(z, 0);
alpar@9 425 goto done;
alpar@9 426 }
alpar@9 427 /* special case when both [x] and [y] are in short format */
alpar@9 428 if (x->ptr == NULL && y->ptr == NULL)
alpar@9 429 { int xval = x->val, yval = y->val, sz = +1;
alpar@9 430 xassert(xval != 0x80000000 && yval != 0x80000000);
alpar@9 431 if (xval < 0) xval = - xval, sz = - sz;
alpar@9 432 if (yval < 0) yval = - yval, sz = - sz;
alpar@9 433 if (xval <= 0x7FFFFFFF / yval)
alpar@9 434 { mpz_set_si(z, sz * (xval * yval));
alpar@9 435 goto done;
alpar@9 436 }
alpar@9 437 }
alpar@9 438 /* convert [x] to long format, if necessary */
alpar@9 439 if (x->ptr == NULL)
alpar@9 440 { xassert(x->val != 0x80000000);
alpar@9 441 if (x->val >= 0)
alpar@9 442 { sx = +1;
alpar@9 443 t = (unsigned int)(+ x->val);
alpar@9 444 }
alpar@9 445 else
alpar@9 446 { sx = -1;
alpar@9 447 t = (unsigned int)(- x->val);
alpar@9 448 }
alpar@9 449 ex = &dumx;
alpar@9 450 ex->d[0] = (unsigned short)t;
alpar@9 451 ex->d[1] = (unsigned short)(t >> 16);
alpar@9 452 ex->d[2] = ex->d[3] = ex->d[4] = ex->d[5] = 0;
alpar@9 453 ex->next = NULL;
alpar@9 454 }
alpar@9 455 else
alpar@9 456 { sx = x->val;
alpar@9 457 xassert(sx == +1 || sx == -1);
alpar@9 458 ex = x->ptr;
alpar@9 459 }
alpar@9 460 /* convert [y] to long format, if necessary */
alpar@9 461 if (y->ptr == NULL)
alpar@9 462 { xassert(y->val != 0x80000000);
alpar@9 463 if (y->val >= 0)
alpar@9 464 { sy = +1;
alpar@9 465 t = (unsigned int)(+ y->val);
alpar@9 466 }
alpar@9 467 else
alpar@9 468 { sy = -1;
alpar@9 469 t = (unsigned int)(- y->val);
alpar@9 470 }
alpar@9 471 ey = &dumy;
alpar@9 472 ey->d[0] = (unsigned short)t;
alpar@9 473 ey->d[1] = (unsigned short)(t >> 16);
alpar@9 474 ey->d[2] = ey->d[3] = ey->d[4] = ey->d[5] = 0;
alpar@9 475 ey->next = NULL;
alpar@9 476 }
alpar@9 477 else
alpar@9 478 { sy = y->val;
alpar@9 479 xassert(sy == +1 || sy == -1);
alpar@9 480 ey = y->ptr;
alpar@9 481 }
alpar@9 482 /* determine the number of digits of [x] */
alpar@9 483 nx = n = 0;
alpar@9 484 for (e = ex; e != NULL; e = e->next)
alpar@9 485 for (k = 0; k <= 5; k++)
alpar@9 486 { n++;
alpar@9 487 if (e->d[k]) nx = n;
alpar@9 488 }
alpar@9 489 xassert(nx > 0);
alpar@9 490 /* determine the number of digits of [y] */
alpar@9 491 ny = n = 0;
alpar@9 492 for (e = ey; e != NULL; e = e->next)
alpar@9 493 for (k = 0; k <= 5; k++)
alpar@9 494 { n++;
alpar@9 495 if (e->d[k]) ny = n;
alpar@9 496 }
alpar@9 497 xassert(ny > 0);
alpar@9 498 /* we need working array containing at least nx+ny+ny places */
alpar@9 499 work = gmp_get_work(nx+ny+ny);
alpar@9 500 /* load digits of [x] */
alpar@9 501 wx = &work[0];
alpar@9 502 for (n = 0; n < nx; n++) wx[ny+n] = 0;
alpar@9 503 for (n = 0, e = ex; e != NULL; e = e->next)
alpar@9 504 for (k = 0; k <= 5; k++, n++)
alpar@9 505 if (e->d[k]) wx[ny+n] = e->d[k];
alpar@9 506 /* load digits of [y] */
alpar@9 507 wy = &work[nx+ny];
alpar@9 508 for (n = 0; n < ny; n++) wy[n] = 0;
alpar@9 509 for (n = 0, e = ey; e != NULL; e = e->next)
alpar@9 510 for (k = 0; k <= 5; k++, n++)
alpar@9 511 if (e->d[k]) wy[n] = e->d[k];
alpar@9 512 /* compute [x] * [y] */
alpar@9 513 bigmul(nx, ny, wx, wy);
alpar@9 514 /* construct and normalize result */
alpar@9 515 mpz_set_si(z, 0);
alpar@9 516 z->val = sx * sy;
alpar@9 517 es = NULL;
alpar@9 518 k = 6;
alpar@9 519 for (n = 0; n < nx+ny; n++)
alpar@9 520 { if (k > 5)
alpar@9 521 { e = gmp_get_atom(sizeof(struct mpz_seg));
alpar@9 522 e->d[0] = e->d[1] = e->d[2] = 0;
alpar@9 523 e->d[3] = e->d[4] = e->d[5] = 0;
alpar@9 524 e->next = NULL;
alpar@9 525 if (z->ptr == NULL)
alpar@9 526 z->ptr = e;
alpar@9 527 else
alpar@9 528 es->next = e;
alpar@9 529 es = e;
alpar@9 530 k = 0;
alpar@9 531 }
alpar@9 532 es->d[k++] = wx[n];
alpar@9 533 }
alpar@9 534 normalize(z);
alpar@9 535 done: return;
alpar@9 536 }
alpar@9 537
alpar@9 538 void mpz_neg(mpz_t z, mpz_t x)
alpar@9 539 { /* set z to 0 - x */
alpar@9 540 mpz_set(z, x);
alpar@9 541 z->val = - z->val;
alpar@9 542 return;
alpar@9 543 }
alpar@9 544
alpar@9 545 void mpz_abs(mpz_t z, mpz_t x)
alpar@9 546 { /* set z to the absolute value of x */
alpar@9 547 mpz_set(z, x);
alpar@9 548 if (z->val < 0) z->val = - z->val;
alpar@9 549 return;
alpar@9 550 }
alpar@9 551
alpar@9 552 void mpz_div(mpz_t q, mpz_t r, mpz_t x, mpz_t y)
alpar@9 553 { /* divide x by y, forming quotient q and/or remainder r
alpar@9 554 if q = NULL then quotient is not stored; if r = NULL then
alpar@9 555 remainder is not stored
alpar@9 556 the sign of quotient is determined as in algebra while the
alpar@9 557 sign of remainder is the same as the sign of dividend:
alpar@9 558 +26 : +7 = +3, remainder is +5
alpar@9 559 -26 : +7 = -3, remainder is -5
alpar@9 560 +26 : -7 = -3, remainder is +5
alpar@9 561 -26 : -7 = +3, remainder is -5 */
alpar@9 562 struct mpz_seg dumx, dumy, *ex, *ey, *es, *e;
alpar@9 563 int sx, sy, k, nx, ny, n;
alpar@9 564 unsigned int t;
alpar@9 565 unsigned short *work, *wx, *wy;
alpar@9 566 /* divide by zero is not allowed */
alpar@9 567 if (y->val == 0)
alpar@9 568 { xassert(y->ptr == NULL);
alpar@9 569 xfault("mpz_div: divide by zero not allowed\n");
alpar@9 570 }
alpar@9 571 /* if [x] = 0 then [q] = [r] = 0 */
alpar@9 572 if (x->val == 0)
alpar@9 573 { xassert(x->ptr == NULL);
alpar@9 574 if (q != NULL) mpz_set_si(q, 0);
alpar@9 575 if (r != NULL) mpz_set_si(r, 0);
alpar@9 576 goto done;
alpar@9 577 }
alpar@9 578 /* special case when both [x] and [y] are in short format */
alpar@9 579 if (x->ptr == NULL && y->ptr == NULL)
alpar@9 580 { int xval = x->val, yval = y->val;
alpar@9 581 xassert(xval != 0x80000000 && yval != 0x80000000);
alpar@9 582 if (q != NULL) mpz_set_si(q, xval / yval);
alpar@9 583 if (r != NULL) mpz_set_si(r, xval % yval);
alpar@9 584 goto done;
alpar@9 585 }
alpar@9 586 /* convert [x] to long format, if necessary */
alpar@9 587 if (x->ptr == NULL)
alpar@9 588 { xassert(x->val != 0x80000000);
alpar@9 589 if (x->val >= 0)
alpar@9 590 { sx = +1;
alpar@9 591 t = (unsigned int)(+ x->val);
alpar@9 592 }
alpar@9 593 else
alpar@9 594 { sx = -1;
alpar@9 595 t = (unsigned int)(- x->val);
alpar@9 596 }
alpar@9 597 ex = &dumx;
alpar@9 598 ex->d[0] = (unsigned short)t;
alpar@9 599 ex->d[1] = (unsigned short)(t >> 16);
alpar@9 600 ex->d[2] = ex->d[3] = ex->d[4] = ex->d[5] = 0;
alpar@9 601 ex->next = NULL;
alpar@9 602 }
alpar@9 603 else
alpar@9 604 { sx = x->val;
alpar@9 605 xassert(sx == +1 || sx == -1);
alpar@9 606 ex = x->ptr;
alpar@9 607 }
alpar@9 608 /* convert [y] to long format, if necessary */
alpar@9 609 if (y->ptr == NULL)
alpar@9 610 { xassert(y->val != 0x80000000);
alpar@9 611 if (y->val >= 0)
alpar@9 612 { sy = +1;
alpar@9 613 t = (unsigned int)(+ y->val);
alpar@9 614 }
alpar@9 615 else
alpar@9 616 { sy = -1;
alpar@9 617 t = (unsigned int)(- y->val);
alpar@9 618 }
alpar@9 619 ey = &dumy;
alpar@9 620 ey->d[0] = (unsigned short)t;
alpar@9 621 ey->d[1] = (unsigned short)(t >> 16);
alpar@9 622 ey->d[2] = ey->d[3] = ey->d[4] = ey->d[5] = 0;
alpar@9 623 ey->next = NULL;
alpar@9 624 }
alpar@9 625 else
alpar@9 626 { sy = y->val;
alpar@9 627 xassert(sy == +1 || sy == -1);
alpar@9 628 ey = y->ptr;
alpar@9 629 }
alpar@9 630 /* determine the number of digits of [x] */
alpar@9 631 nx = n = 0;
alpar@9 632 for (e = ex; e != NULL; e = e->next)
alpar@9 633 for (k = 0; k <= 5; k++)
alpar@9 634 { n++;
alpar@9 635 if (e->d[k]) nx = n;
alpar@9 636 }
alpar@9 637 xassert(nx > 0);
alpar@9 638 /* determine the number of digits of [y] */
alpar@9 639 ny = n = 0;
alpar@9 640 for (e = ey; e != NULL; e = e->next)
alpar@9 641 for (k = 0; k <= 5; k++)
alpar@9 642 { n++;
alpar@9 643 if (e->d[k]) ny = n;
alpar@9 644 }
alpar@9 645 xassert(ny > 0);
alpar@9 646 /* if nx < ny then [q] = 0 and [r] = [x] */
alpar@9 647 if (nx < ny)
alpar@9 648 { if (r != NULL) mpz_set(r, x);
alpar@9 649 if (q != NULL) mpz_set_si(q, 0);
alpar@9 650 goto done;
alpar@9 651 }
alpar@9 652 /* we need working array containing at least nx+ny+1 places */
alpar@9 653 work = gmp_get_work(nx+ny+1);
alpar@9 654 /* load digits of [x] */
alpar@9 655 wx = &work[0];
alpar@9 656 for (n = 0; n < nx; n++) wx[n] = 0;
alpar@9 657 for (n = 0, e = ex; e != NULL; e = e->next)
alpar@9 658 for (k = 0; k <= 5; k++, n++)
alpar@9 659 if (e->d[k]) wx[n] = e->d[k];
alpar@9 660 /* load digits of [y] */
alpar@9 661 wy = &work[nx+1];
alpar@9 662 for (n = 0; n < ny; n++) wy[n] = 0;
alpar@9 663 for (n = 0, e = ey; e != NULL; e = e->next)
alpar@9 664 for (k = 0; k <= 5; k++, n++)
alpar@9 665 if (e->d[k]) wy[n] = e->d[k];
alpar@9 666 /* compute quotient and remainder */
alpar@9 667 xassert(wy[ny-1] != 0);
alpar@9 668 bigdiv(nx-ny, ny, wx, wy);
alpar@9 669 /* construct and normalize quotient */
alpar@9 670 if (q != NULL)
alpar@9 671 { mpz_set_si(q, 0);
alpar@9 672 q->val = sx * sy;
alpar@9 673 es = NULL;
alpar@9 674 k = 6;
alpar@9 675 for (n = ny; n <= nx; n++)
alpar@9 676 { if (k > 5)
alpar@9 677 { e = gmp_get_atom(sizeof(struct mpz_seg));
alpar@9 678 e->d[0] = e->d[1] = e->d[2] = 0;
alpar@9 679 e->d[3] = e->d[4] = e->d[5] = 0;
alpar@9 680 e->next = NULL;
alpar@9 681 if (q->ptr == NULL)
alpar@9 682 q->ptr = e;
alpar@9 683 else
alpar@9 684 es->next = e;
alpar@9 685 es = e;
alpar@9 686 k = 0;
alpar@9 687 }
alpar@9 688 es->d[k++] = wx[n];
alpar@9 689 }
alpar@9 690 normalize(q);
alpar@9 691 }
alpar@9 692 /* construct and normalize remainder */
alpar@9 693 if (r != NULL)
alpar@9 694 { mpz_set_si(r, 0);
alpar@9 695 r->val = sx;
alpar@9 696 es = NULL;
alpar@9 697 k = 6;
alpar@9 698 for (n = 0; n < ny; n++)
alpar@9 699 { if (k > 5)
alpar@9 700 { e = gmp_get_atom(sizeof(struct mpz_seg));
alpar@9 701 e->d[0] = e->d[1] = e->d[2] = 0;
alpar@9 702 e->d[3] = e->d[4] = e->d[5] = 0;
alpar@9 703 e->next = NULL;
alpar@9 704 if (r->ptr == NULL)
alpar@9 705 r->ptr = e;
alpar@9 706 else
alpar@9 707 es->next = e;
alpar@9 708 es = e;
alpar@9 709 k = 0;
alpar@9 710 }
alpar@9 711 es->d[k++] = wx[n];
alpar@9 712 }
alpar@9 713 normalize(r);
alpar@9 714 }
alpar@9 715 done: return;
alpar@9 716 }
alpar@9 717
alpar@9 718 void mpz_gcd(mpz_t z, mpz_t x, mpz_t y)
alpar@9 719 { /* set z to the greatest common divisor of x and y */
alpar@9 720 /* in case of arbitrary integers GCD(x, y) = GCD(|x|, |y|), and,
alpar@9 721 in particular, GCD(0, 0) = 0 */
alpar@9 722 mpz_t u, v, r;
alpar@9 723 mpz_init(u);
alpar@9 724 mpz_init(v);
alpar@9 725 mpz_init(r);
alpar@9 726 mpz_abs(u, x);
alpar@9 727 mpz_abs(v, y);
alpar@9 728 while (mpz_sgn(v))
alpar@9 729 { mpz_div(NULL, r, u, v);
alpar@9 730 mpz_set(u, v);
alpar@9 731 mpz_set(v, r);
alpar@9 732 }
alpar@9 733 mpz_set(z, u);
alpar@9 734 mpz_clear(u);
alpar@9 735 mpz_clear(v);
alpar@9 736 mpz_clear(r);
alpar@9 737 return;
alpar@9 738 }
alpar@9 739
alpar@9 740 int mpz_cmp(mpz_t x, mpz_t y)
alpar@9 741 { /* compare x and y; return a positive value if x > y, zero if
alpar@9 742 x = y, or a nefative value if x < y */
alpar@9 743 static struct mpz_seg zero = { { 0, 0, 0, 0, 0, 0 }, NULL };
alpar@9 744 struct mpz_seg dumx, dumy, *ex, *ey;
alpar@9 745 int cc, sx, sy, k;
alpar@9 746 unsigned int t;
alpar@9 747 if (x == y)
alpar@9 748 { cc = 0;
alpar@9 749 goto done;
alpar@9 750 }
alpar@9 751 /* special case when both [x] and [y] are in short format */
alpar@9 752 if (x->ptr == NULL && y->ptr == NULL)
alpar@9 753 { int xval = x->val, yval = y->val;
alpar@9 754 xassert(xval != 0x80000000 && yval != 0x80000000);
alpar@9 755 cc = (xval > yval ? +1 : xval < yval ? -1 : 0);
alpar@9 756 goto done;
alpar@9 757 }
alpar@9 758 /* special case when [x] and [y] have different signs */
alpar@9 759 if (x->val > 0 && y->val <= 0 || x->val == 0 && y->val < 0)
alpar@9 760 { cc = +1;
alpar@9 761 goto done;
alpar@9 762 }
alpar@9 763 if (x->val < 0 && y->val >= 0 || x->val == 0 && y->val > 0)
alpar@9 764 { cc = -1;
alpar@9 765 goto done;
alpar@9 766 }
alpar@9 767 /* convert [x] to long format, if necessary */
alpar@9 768 if (x->ptr == NULL)
alpar@9 769 { xassert(x->val != 0x80000000);
alpar@9 770 if (x->val >= 0)
alpar@9 771 { sx = +1;
alpar@9 772 t = (unsigned int)(+ x->val);
alpar@9 773 }
alpar@9 774 else
alpar@9 775 { sx = -1;
alpar@9 776 t = (unsigned int)(- x->val);
alpar@9 777 }
alpar@9 778 ex = &dumx;
alpar@9 779 ex->d[0] = (unsigned short)t;
alpar@9 780 ex->d[1] = (unsigned short)(t >> 16);
alpar@9 781 ex->d[2] = ex->d[3] = ex->d[4] = ex->d[5] = 0;
alpar@9 782 ex->next = NULL;
alpar@9 783 }
alpar@9 784 else
alpar@9 785 { sx = x->val;
alpar@9 786 xassert(sx == +1 || sx == -1);
alpar@9 787 ex = x->ptr;
alpar@9 788 }
alpar@9 789 /* convert [y] to long format, if necessary */
alpar@9 790 if (y->ptr == NULL)
alpar@9 791 { xassert(y->val != 0x80000000);
alpar@9 792 if (y->val >= 0)
alpar@9 793 { sy = +1;
alpar@9 794 t = (unsigned int)(+ y->val);
alpar@9 795 }
alpar@9 796 else
alpar@9 797 { sy = -1;
alpar@9 798 t = (unsigned int)(- y->val);
alpar@9 799 }
alpar@9 800 ey = &dumy;
alpar@9 801 ey->d[0] = (unsigned short)t;
alpar@9 802 ey->d[1] = (unsigned short)(t >> 16);
alpar@9 803 ey->d[2] = ey->d[3] = ey->d[4] = ey->d[5] = 0;
alpar@9 804 ey->next = NULL;
alpar@9 805 }
alpar@9 806 else
alpar@9 807 { sy = y->val;
alpar@9 808 xassert(sy == +1 || sy == -1);
alpar@9 809 ey = y->ptr;
alpar@9 810 }
alpar@9 811 /* main fragment */
alpar@9 812 xassert(sx > 0 && sy > 0 || sx < 0 && sy < 0);
alpar@9 813 cc = 0;
alpar@9 814 for (; ex || ey; ex = ex->next, ey = ey->next)
alpar@9 815 { if (ex == NULL) ex = &zero;
alpar@9 816 if (ey == NULL) ey = &zero;
alpar@9 817 for (k = 0; k <= 5; k++)
alpar@9 818 { if (ex->d[k] > ey->d[k]) cc = +1;
alpar@9 819 if (ex->d[k] < ey->d[k]) cc = -1;
alpar@9 820 }
alpar@9 821 }
alpar@9 822 if (sx < 0) cc = - cc;
alpar@9 823 done: return cc;
alpar@9 824 }
alpar@9 825
alpar@9 826 int mpz_sgn(mpz_t x)
alpar@9 827 { /* return +1 if x > 0, 0 if x = 0, and -1 if x < 0 */
alpar@9 828 int s;
alpar@9 829 s = (x->val > 0 ? +1 : x->val < 0 ? -1 : 0);
alpar@9 830 return s;
alpar@9 831 }
alpar@9 832
alpar@9 833 int mpz_out_str(void *_fp, int base, mpz_t x)
alpar@9 834 { /* output x on stream fp, as a string in given base; the base
alpar@9 835 may vary from 2 to 36;
alpar@9 836 return the number of bytes written, or if an error occurred,
alpar@9 837 return 0 */
alpar@9 838 FILE *fp = _fp;
alpar@9 839 mpz_t b, y, r;
alpar@9 840 int n, j, nwr = 0;
alpar@9 841 unsigned char *d;
alpar@9 842 static char *set = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ";
alpar@9 843 if (!(2 <= base && base <= 36))
alpar@9 844 xfault("mpz_out_str: base = %d; invalid base\n", base);
alpar@9 845 mpz_init(b);
alpar@9 846 mpz_set_si(b, base);
alpar@9 847 mpz_init(y);
alpar@9 848 mpz_init(r);
alpar@9 849 /* determine the number of digits */
alpar@9 850 mpz_abs(y, x);
alpar@9 851 for (n = 0; mpz_sgn(y) != 0; n++)
alpar@9 852 mpz_div(y, NULL, y, b);
alpar@9 853 if (n == 0) n = 1;
alpar@9 854 /* compute the digits */
alpar@9 855 d = xmalloc(n);
alpar@9 856 mpz_abs(y, x);
alpar@9 857 for (j = 0; j < n; j++)
alpar@9 858 { mpz_div(y, r, y, b);
alpar@9 859 xassert(0 <= r->val && r->val < base && r->ptr == NULL);
alpar@9 860 d[j] = (unsigned char)r->val;
alpar@9 861 }
alpar@9 862 /* output the integer to the stream */
alpar@9 863 if (fp == NULL) fp = stdout;
alpar@9 864 if (mpz_sgn(x) < 0)
alpar@9 865 fputc('-', fp), nwr++;
alpar@9 866 for (j = n-1; j >= 0; j--)
alpar@9 867 fputc(set[d[j]], fp), nwr++;
alpar@9 868 if (ferror(fp)) nwr = 0;
alpar@9 869 mpz_clear(b);
alpar@9 870 mpz_clear(y);
alpar@9 871 mpz_clear(r);
alpar@9 872 xfree(d);
alpar@9 873 return nwr;
alpar@9 874 }
alpar@9 875
alpar@9 876 /*====================================================================*/
alpar@9 877
alpar@9 878 mpq_t _mpq_init(void)
alpar@9 879 { /* initialize x, and set its value to 0/1 */
alpar@9 880 mpq_t x;
alpar@9 881 x = gmp_get_atom(sizeof(struct mpq));
alpar@9 882 x->p.val = 0;
alpar@9 883 x->p.ptr = NULL;
alpar@9 884 x->q.val = 1;
alpar@9 885 x->q.ptr = NULL;
alpar@9 886 return x;
alpar@9 887 }
alpar@9 888
alpar@9 889 void mpq_clear(mpq_t x)
alpar@9 890 { /* free the space occupied by x */
alpar@9 891 mpz_set_si(&x->p, 0);
alpar@9 892 xassert(x->p.ptr == NULL);
alpar@9 893 mpz_set_si(&x->q, 0);
alpar@9 894 xassert(x->q.ptr == NULL);
alpar@9 895 /* free the number descriptor */
alpar@9 896 gmp_free_atom(x, sizeof(struct mpq));
alpar@9 897 return;
alpar@9 898 }
alpar@9 899
alpar@9 900 void mpq_canonicalize(mpq_t x)
alpar@9 901 { /* remove any factors that are common to the numerator and
alpar@9 902 denominator of x, and make the denominator positive */
alpar@9 903 mpz_t f;
alpar@9 904 xassert(x->q.val != 0);
alpar@9 905 if (x->q.val < 0)
alpar@9 906 { mpz_neg(&x->p, &x->p);
alpar@9 907 mpz_neg(&x->q, &x->q);
alpar@9 908 }
alpar@9 909 mpz_init(f);
alpar@9 910 mpz_gcd(f, &x->p, &x->q);
alpar@9 911 if (!(f->val == 1 && f->ptr == NULL))
alpar@9 912 { mpz_div(&x->p, NULL, &x->p, f);
alpar@9 913 mpz_div(&x->q, NULL, &x->q, f);
alpar@9 914 }
alpar@9 915 mpz_clear(f);
alpar@9 916 return;
alpar@9 917 }
alpar@9 918
alpar@9 919 void mpq_set(mpq_t z, mpq_t x)
alpar@9 920 { /* set the value of z from x */
alpar@9 921 if (z != x)
alpar@9 922 { mpz_set(&z->p, &x->p);
alpar@9 923 mpz_set(&z->q, &x->q);
alpar@9 924 }
alpar@9 925 return;
alpar@9 926 }
alpar@9 927
alpar@9 928 void mpq_set_si(mpq_t x, int p, unsigned int q)
alpar@9 929 { /* set the value of x to p/q */
alpar@9 930 if (q == 0)
alpar@9 931 xfault("mpq_set_si: zero denominator not allowed\n");
alpar@9 932 mpz_set_si(&x->p, p);
alpar@9 933 xassert(q <= 0x7FFFFFFF);
alpar@9 934 mpz_set_si(&x->q, q);
alpar@9 935 return;
alpar@9 936 }
alpar@9 937
alpar@9 938 double mpq_get_d(mpq_t x)
alpar@9 939 { /* convert x to a double, truncating if necessary */
alpar@9 940 int np, nq;
alpar@9 941 double p, q;
alpar@9 942 p = mpz_get_d_2exp(&np, &x->p);
alpar@9 943 q = mpz_get_d_2exp(&nq, &x->q);
alpar@9 944 return ldexp(p / q, np - nq);
alpar@9 945 }
alpar@9 946
alpar@9 947 void mpq_set_d(mpq_t x, double val)
alpar@9 948 { /* set x to val; there is no rounding, the conversion is exact */
alpar@9 949 int s, n, d, j;
alpar@9 950 double f;
alpar@9 951 mpz_t temp;
alpar@9 952 xassert(-DBL_MAX <= val && val <= +DBL_MAX);
alpar@9 953 mpq_set_si(x, 0, 1);
alpar@9 954 if (val > 0.0)
alpar@9 955 s = +1;
alpar@9 956 else if (val < 0.0)
alpar@9 957 s = -1;
alpar@9 958 else
alpar@9 959 goto done;
alpar@9 960 f = frexp(fabs(val), &n);
alpar@9 961 /* |val| = f * 2^n, where 0.5 <= f < 1.0 */
alpar@9 962 mpz_init(temp);
alpar@9 963 while (f != 0.0)
alpar@9 964 { f *= 16.0, n -= 4;
alpar@9 965 d = (int)f;
alpar@9 966 xassert(0 <= d && d <= 15);
alpar@9 967 f -= (double)d;
alpar@9 968 /* x := 16 * x + d */
alpar@9 969 mpz_set_si(temp, 16);
alpar@9 970 mpz_mul(&x->p, &x->p, temp);
alpar@9 971 mpz_set_si(temp, d);
alpar@9 972 mpz_add(&x->p, &x->p, temp);
alpar@9 973 }
alpar@9 974 mpz_clear(temp);
alpar@9 975 /* x := x * 2^n */
alpar@9 976 if (n > 0)
alpar@9 977 { for (j = 1; j <= n; j++)
alpar@9 978 mpz_add(&x->p, &x->p, &x->p);
alpar@9 979 }
alpar@9 980 else if (n < 0)
alpar@9 981 { for (j = 1; j <= -n; j++)
alpar@9 982 mpz_add(&x->q, &x->q, &x->q);
alpar@9 983 mpq_canonicalize(x);
alpar@9 984 }
alpar@9 985 if (s < 0) mpq_neg(x, x);
alpar@9 986 done: return;
alpar@9 987 }
alpar@9 988
alpar@9 989 void mpq_add(mpq_t z, mpq_t x, mpq_t y)
alpar@9 990 { /* set z to x + y */
alpar@9 991 mpz_t p, q;
alpar@9 992 mpz_init(p);
alpar@9 993 mpz_init(q);
alpar@9 994 mpz_mul(p, &x->p, &y->q);
alpar@9 995 mpz_mul(q, &x->q, &y->p);
alpar@9 996 mpz_add(p, p, q);
alpar@9 997 mpz_mul(q, &x->q, &y->q);
alpar@9 998 mpz_set(&z->p, p);
alpar@9 999 mpz_set(&z->q, q);
alpar@9 1000 mpz_clear(p);
alpar@9 1001 mpz_clear(q);
alpar@9 1002 mpq_canonicalize(z);
alpar@9 1003 return;
alpar@9 1004 }
alpar@9 1005
alpar@9 1006 void mpq_sub(mpq_t z, mpq_t x, mpq_t y)
alpar@9 1007 { /* set z to x - y */
alpar@9 1008 mpz_t p, q;
alpar@9 1009 mpz_init(p);
alpar@9 1010 mpz_init(q);
alpar@9 1011 mpz_mul(p, &x->p, &y->q);
alpar@9 1012 mpz_mul(q, &x->q, &y->p);
alpar@9 1013 mpz_sub(p, p, q);
alpar@9 1014 mpz_mul(q, &x->q, &y->q);
alpar@9 1015 mpz_set(&z->p, p);
alpar@9 1016 mpz_set(&z->q, q);
alpar@9 1017 mpz_clear(p);
alpar@9 1018 mpz_clear(q);
alpar@9 1019 mpq_canonicalize(z);
alpar@9 1020 return;
alpar@9 1021 }
alpar@9 1022
alpar@9 1023 void mpq_mul(mpq_t z, mpq_t x, mpq_t y)
alpar@9 1024 { /* set z to x * y */
alpar@9 1025 mpz_mul(&z->p, &x->p, &y->p);
alpar@9 1026 mpz_mul(&z->q, &x->q, &y->q);
alpar@9 1027 mpq_canonicalize(z);
alpar@9 1028 return;
alpar@9 1029 }
alpar@9 1030
alpar@9 1031 void mpq_div(mpq_t z, mpq_t x, mpq_t y)
alpar@9 1032 { /* set z to x / y */
alpar@9 1033 mpz_t p, q;
alpar@9 1034 if (mpq_sgn(y) == 0)
alpar@9 1035 xfault("mpq_div: zero divisor not allowed\n");
alpar@9 1036 mpz_init(p);
alpar@9 1037 mpz_init(q);
alpar@9 1038 mpz_mul(p, &x->p, &y->q);
alpar@9 1039 mpz_mul(q, &x->q, &y->p);
alpar@9 1040 mpz_set(&z->p, p);
alpar@9 1041 mpz_set(&z->q, q);
alpar@9 1042 mpz_clear(p);
alpar@9 1043 mpz_clear(q);
alpar@9 1044 mpq_canonicalize(z);
alpar@9 1045 return;
alpar@9 1046 }
alpar@9 1047
alpar@9 1048 void mpq_neg(mpq_t z, mpq_t x)
alpar@9 1049 { /* set z to 0 - x */
alpar@9 1050 mpq_set(z, x);
alpar@9 1051 mpz_neg(&z->p, &z->p);
alpar@9 1052 return;
alpar@9 1053 }
alpar@9 1054
alpar@9 1055 void mpq_abs(mpq_t z, mpq_t x)
alpar@9 1056 { /* set z to the absolute value of x */
alpar@9 1057 mpq_set(z, x);
alpar@9 1058 mpz_abs(&z->p, &z->p);
alpar@9 1059 xassert(mpz_sgn(&x->q) > 0);
alpar@9 1060 return;
alpar@9 1061 }
alpar@9 1062
alpar@9 1063 int mpq_cmp(mpq_t x, mpq_t y)
alpar@9 1064 { /* compare x and y; return a positive value if x > y, zero if
alpar@9 1065 x = y, or a nefative value if x < y */
alpar@9 1066 mpq_t temp;
alpar@9 1067 int s;
alpar@9 1068 mpq_init(temp);
alpar@9 1069 mpq_sub(temp, x, y);
alpar@9 1070 s = mpq_sgn(temp);
alpar@9 1071 mpq_clear(temp);
alpar@9 1072 return s;
alpar@9 1073 }
alpar@9 1074
alpar@9 1075 int mpq_sgn(mpq_t x)
alpar@9 1076 { /* return +1 if x > 0, 0 if x = 0, and -1 if x < 0 */
alpar@9 1077 int s;
alpar@9 1078 s = mpz_sgn(&x->p);
alpar@9 1079 xassert(mpz_sgn(&x->q) > 0);
alpar@9 1080 return s;
alpar@9 1081 }
alpar@9 1082
alpar@9 1083 int mpq_out_str(void *_fp, int base, mpq_t x)
alpar@9 1084 { /* output x on stream fp, as a string in given base; the base
alpar@9 1085 may vary from 2 to 36; output is in the form 'num/den' or if
alpar@9 1086 the denominator is 1 then just 'num';
alpar@9 1087 if the parameter fp is a null pointer, stdout is assumed;
alpar@9 1088 return the number of bytes written, or if an error occurred,
alpar@9 1089 return 0 */
alpar@9 1090 FILE *fp = _fp;
alpar@9 1091 int nwr;
alpar@9 1092 if (!(2 <= base && base <= 36))
alpar@9 1093 xfault("mpq_out_str: base = %d; invalid base\n", base);
alpar@9 1094 if (fp == NULL) fp = stdout;
alpar@9 1095 nwr = mpz_out_str(fp, base, &x->p);
alpar@9 1096 if (x->q.val == 1 && x->q.ptr == NULL)
alpar@9 1097 ;
alpar@9 1098 else
alpar@9 1099 { fputc('/', fp), nwr++;
alpar@9 1100 nwr += mpz_out_str(fp, base, &x->q);
alpar@9 1101 }
alpar@9 1102 if (ferror(fp)) nwr = 0;
alpar@9 1103 return nwr;
alpar@9 1104 }
alpar@9 1105
alpar@9 1106 #endif
alpar@9 1107
alpar@9 1108 /* eof */