lemon-project-template-glpk: 33de93886c88 deps/glpk/src/glpgmp.c

lemon-project-template-glpk

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

Import GLPK 4.47

author	Alpar Juttner <alpar@cs.elte.hu>
date	Sun, 06 Nov 2011 20:59:10 +0100
parents
children

line source

1 /* glpgmp.c */

3 /***********************************************************************

4 * This code is part of GLPK (GNU Linear Programming Kit).

5 *

7 * 2009, 2010, 2011 Andrew Makhorin, Department for Applied Informatics,

9 * E-mail: <mao@gnu.org>.

10 *

11 * GLPK is free software: you can redistribute it and/or modify it

12 * under the terms of the GNU General Public License as published by

13 * the Free Software Foundation, either version 3 of the License, or

14 * (at your option) any later version.

15 *

16 * GLPK is distributed in the hope that it will be useful, but WITHOUT

17 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY

18 * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public

19 * License for more details.

20 *

21 * You should have received a copy of the GNU General Public License

22 * along with GLPK. If not, see <http://www.gnu.org/licenses/>.

23 ***********************************************************************/

25 #define _GLPSTD_STDIO

26 #include "glpdmp.h"

27 #include "glpgmp.h"

28 #define xfault xerror

30 #ifdef HAVE_GMP /* use GNU MP bignum library */

32 int gmp_pool_count(void) { return 0; }

34 void gmp_free_mem(void) { return; }

36 #else /* use GLPK bignum module */

38 static DMP *gmp_pool = NULL;

39 static int gmp_size = 0;

40 static unsigned short *gmp_work = NULL;

42 void *gmp_get_atom(int size)

43 { if (gmp_pool == NULL)

44 gmp_pool = dmp_create_pool();

45 return dmp_get_atom(gmp_pool, size);

46 }

48 void gmp_free_atom(void *ptr, int size)

49 { xassert(gmp_pool != NULL);

50 dmp_free_atom(gmp_pool, ptr, size);

51 return;

52 }

54 int gmp_pool_count(void)

55 { if (gmp_pool == NULL)

56 return 0;

57 else

58 return dmp_in_use(gmp_pool).lo;

59 }

61 unsigned short *gmp_get_work(int size)

62 { xassert(size > 0);

63 if (gmp_size < size)

64 { if (gmp_size == 0)

65 { xassert(gmp_work == NULL);

66 gmp_size = 100;

67 }

68 else

69 { xassert(gmp_work != NULL);

70 xfree(gmp_work);

71 }

72 while (gmp_size < size) gmp_size += gmp_size;

73 gmp_work = xcalloc(gmp_size, sizeof(unsigned short));

74 }

75 return gmp_work;

76 }

78 void gmp_free_mem(void)

79 { if (gmp_pool != NULL) dmp_delete_pool(gmp_pool);

80 if (gmp_work != NULL) xfree(gmp_work);

81 gmp_pool = NULL;

82 gmp_size = 0;

83 gmp_work = NULL;

84 return;

85 }

87 /*====================================================================*/

89 mpz_t _mpz_init(void)

90 { /* initialize x, and set its value to 0 */

91 mpz_t x;

92 x = gmp_get_atom(sizeof(struct mpz));

93 x->val = 0;

94 x->ptr = NULL;

95 return x;

96 }

98 void mpz_clear(mpz_t x)

99 { /* free the space occupied by x */

100 mpz_set_si(x, 0);

101 xassert(x->ptr == NULL);

102 /* free the number descriptor */

103 gmp_free_atom(x, sizeof(struct mpz));

104 return;

105 }

106

107 void mpz_set(mpz_t z, mpz_t x)

108 { /* set the value of z from x */

109 struct mpz_seg *e, *ee, *es;

110 if (z != x)

111 { mpz_set_si(z, 0);

112 z->val = x->val;

113 xassert(z->ptr == NULL);

114 for (e = x->ptr, es = NULL; e != NULL; e = e->next)

115 { ee = gmp_get_atom(sizeof(struct mpz_seg));

116 memcpy(ee->d, e->d, 12);

117 ee->next = NULL;

118 if (z->ptr == NULL)

119 z->ptr = ee;

120 else

121 es->next = ee;

122 es = ee;

123 }

124 }

125 return;

126 }

127

128 void mpz_set_si(mpz_t x, int val)

129 { /* set the value of x to val */

130 struct mpz_seg *e;

131 /* free existing segments, if any */

132 while (x->ptr != NULL)

133 { e = x->ptr;

134 x->ptr = e->next;

135 gmp_free_atom(e, sizeof(struct mpz_seg));

136 }

137 /* assign new value */

138 if (val == 0x80000000)

139 { /* long format is needed */

140 x->val = -1;

141 x->ptr = e = gmp_get_atom(sizeof(struct mpz_seg));

142 memset(e->d, 0, 12);

143 e->d[1] = 0x8000;

144 e->next = NULL;

145 }

146 else

147 { /* short format is enough */

148 x->val = val;

149 }

150 return;

151 }

152

153 double mpz_get_d(mpz_t x)

154 { /* convert x to a double, truncating if necessary */

155 struct mpz_seg *e;

156 int j;

157 double val, deg;

158 if (x->ptr == NULL)

159 val = (double)x->val;

160 else

161 { xassert(x->val != 0);

162 val = 0.0;

163 deg = 1.0;

164 for (e = x->ptr; e != NULL; e = e->next)

165 { for (j = 0; j <= 5; j++)

166 { val += deg * (double)((int)e->d[j]);

167 deg *= 65536.0;

168 }

169 }

170 if (x->val < 0) val = - val;

171 }

172 return val;

173 }

174

175 double mpz_get_d_2exp(int *exp, mpz_t x)

176 { /* convert x to a double, truncating if necessary (i.e. rounding

177 towards zero), and returning the exponent separately;

178 the return value is in the range 0.5 <= |d| < 1 and the

179 exponent is stored to *exp; d*2^exp is the (truncated) x value;

180 if x is zero, the return is 0.0 and 0 is stored to *exp;

181 this is similar to the standard C frexp function */

182 struct mpz_seg *e;

183 int j, n, n1;

184 double val;

185 if (x->ptr == NULL)

186 val = (double)x->val, n = 0;

187 else

188 { xassert(x->val != 0);

189 val = 0.0, n = 0;

190 for (e = x->ptr; e != NULL; e = e->next)

191 { for (j = 0; j <= 5; j++)

192 { val += (double)((int)e->d[j]);

193 val /= 65536.0, n += 16;

194 }

195 }

196 if (x->val < 0) val = - val;

197 }

198 val = frexp(val, &n1);

199 *exp = n + n1;

200 return val;

201 }

202

203 void mpz_swap(mpz_t x, mpz_t y)

204 { /* swap the values x and y efficiently */

205 int val;

206 void *ptr;

207 val = x->val, ptr = x->ptr;

208 x->val = y->val, x->ptr = y->ptr;

209 y->val = val, y->ptr = ptr;

210 return;

211 }

212

213 static void normalize(mpz_t x)

214 { /* normalize integer x that includes removing non-significant

215 (leading) zeros and converting to short format, if possible */

216 struct mpz_seg *es, *e;

217 /* if the integer is in short format, it remains unchanged */

218 if (x->ptr == NULL)

219 { xassert(x->val != 0x80000000);

220 goto done;

221 }

222 xassert(x->val == +1 || x->val == -1);

223 /* find the last (most significant) non-zero segment */

224 es = NULL;

225 for (e = x->ptr; e != NULL; e = e->next)

226 { if (e->d[0] || e->d[1] || e->d[2] ||

227 e->d[3] || e->d[4] || e->d[5]) es = e;

228 }

229 /* if all segments contain zeros, the integer is zero */

230 if (es == NULL)

231 { mpz_set_si(x, 0);

232 goto done;

233 }

234 /* remove non-significant (leading) zero segments */

235 while (es->next != NULL)

236 { e = es->next;

237 es->next = e->next;

238 gmp_free_atom(e, sizeof(struct mpz_seg));

239 }

240 /* convert the integer to short format, if possible */

241 e = x->ptr;

242 if (e->next == NULL && e->d[1] <= 0x7FFF &&

243 !e->d[2] && !e->d[3] && !e->d[4] && !e->d[5])

244 { int val;

245 val = (int)e->d[0] + ((int)e->d[1] << 16);

246 if (x->val < 0) val = - val;

247 mpz_set_si(x, val);

248 }

249 done: return;

250 }

251

252 void mpz_add(mpz_t z, mpz_t x, mpz_t y)

253 { /* set z to x + y */

254 static struct mpz_seg zero = { { 0, 0, 0, 0, 0, 0 }, NULL };

255 struct mpz_seg dumx, dumy, *ex, *ey, *ez, *es, *ee;

256 int k, sx, sy, sz;

257 unsigned int t;

258 /* if [x] = 0 then [z] = [y] */

259 if (x->val == 0)

260 { xassert(x->ptr == NULL);

261 mpz_set(z, y);

262 goto done;

263 }

264 /* if [y] = 0 then [z] = [x] */

265 if (y->val == 0)

266 { xassert(y->ptr == NULL);

267 mpz_set(z, x);

268 goto done;

269 }

270 /* special case when both [x] and [y] are in short format */

271 if (x->ptr == NULL && y->ptr == NULL)

272 { int xval = x->val, yval = y->val, zval = x->val + y->val;

273 xassert(xval != 0x80000000 && yval != 0x80000000);

274 if (!(xval > 0 && yval > 0 && zval <= 0 ||

275 xval < 0 && yval < 0 && zval >= 0))

276 { mpz_set_si(z, zval);

277 goto done;

278 }

279 }

280 /* convert [x] to long format, if necessary */

281 if (x->ptr == NULL)

282 { xassert(x->val != 0x80000000);

283 if (x->val >= 0)

284 { sx = +1;

285 t = (unsigned int)(+ x->val);

286 }

287 else

288 { sx = -1;

289 t = (unsigned int)(- x->val);

290 }

291 ex = &dumx;

292 ex->d[0] = (unsigned short)t;

293 ex->d[1] = (unsigned short)(t >> 16);

294 ex->d[2] = ex->d[3] = ex->d[4] = ex->d[5] = 0;

295 ex->next = NULL;

296 }

297 else

298 { sx = x->val;

299 xassert(sx == +1 || sx == -1);

300 ex = x->ptr;

301 }

302 /* convert [y] to long format, if necessary */

303 if (y->ptr == NULL)

304 { xassert(y->val != 0x80000000);

305 if (y->val >= 0)

306 { sy = +1;

307 t = (unsigned int)(+ y->val);

308 }

309 else

310 { sy = -1;

311 t = (unsigned int)(- y->val);

312 }

313 ey = &dumy;

314 ey->d[0] = (unsigned short)t;

315 ey->d[1] = (unsigned short)(t >> 16);

316 ey->d[2] = ey->d[3] = ey->d[4] = ey->d[5] = 0;

317 ey->next = NULL;

318 }

319 else

320 { sy = y->val;

321 xassert(sy == +1 || sy == -1);

322 ey = y->ptr;

323 }

324 /* main fragment */

325 sz = sx;

326 ez = es = NULL;

327 if (sx > 0 && sy > 0 || sx < 0 && sy < 0)

328 { /* [x] and [y] have identical signs -- addition */

329 t = 0;

330 for (; ex || ey; ex = ex->next, ey = ey->next)

331 { if (ex == NULL) ex = &zero;

332 if (ey == NULL) ey = &zero;

333 ee = gmp_get_atom(sizeof(struct mpz_seg));

334 for (k = 0; k <= 5; k++)

335 { t += (unsigned int)ex->d[k];

336 t += (unsigned int)ey->d[k];

337 ee->d[k] = (unsigned short)t;

338 t >>= 16;

339 }

340 ee->next = NULL;

341 if (ez == NULL)

342 ez = ee;

343 else

344 es->next = ee;

345 es = ee;

346 }

347 if (t)

348 { /* overflow -- one extra digit is needed */

349 ee = gmp_get_atom(sizeof(struct mpz_seg));

350 ee->d[0] = 1;

351 ee->d[1] = ee->d[2] = ee->d[3] = ee->d[4] = ee->d[5] = 0;

352 ee->next = NULL;

353 xassert(es != NULL);

354 es->next = ee;

355 }

356 }

357 else

358 { /* [x] and [y] have different signs -- subtraction */

359 t = 1;

360 for (; ex || ey; ex = ex->next, ey = ey->next)

361 { if (ex == NULL) ex = &zero;

362 if (ey == NULL) ey = &zero;

363 ee = gmp_get_atom(sizeof(struct mpz_seg));

364 for (k = 0; k <= 5; k++)

365 { t += (unsigned int)ex->d[k];

366 t += (0xFFFF - (unsigned int)ey->d[k]);

367 ee->d[k] = (unsigned short)t;

368 t >>= 16;

369 }

370 ee->next = NULL;

371 if (ez == NULL)

372 ez = ee;

373 else

374 es->next = ee;

375 es = ee;

376 }

377 if (!t)

378 { /* |[x]| < |[y]| -- result in complement coding */

379 sz = - sz;

380 t = 1;

381 for (ee = ez; ee != NULL; ee = ee->next)

382 for (k = 0; k <= 5; k++)

383 { t += (0xFFFF - (unsigned int)ee->d[k]);

384 ee->d[k] = (unsigned short)t;

385 t >>= 16;

386 }

387 }

388 }

389 /* contruct and normalize result */

390 mpz_set_si(z, 0);

391 z->val = sz;

392 z->ptr = ez;

393 normalize(z);

394 done: return;

395 }

396

397 void mpz_sub(mpz_t z, mpz_t x, mpz_t y)

398 { /* set z to x - y */

399 if (x == y)

400 mpz_set_si(z, 0);

401 else

402 { y->val = - y->val;

403 mpz_add(z, x, y);

404 if (y != z) y->val = - y->val;

405 }

406 return;

407 }

408

409 void mpz_mul(mpz_t z, mpz_t x, mpz_t y)

410 { /* set z to x * y */

411 struct mpz_seg dumx, dumy, *ex, *ey, *es, *e;

412 int sx, sy, k, nx, ny, n;

413 unsigned int t;

414 unsigned short *work, *wx, *wy;

415 /* if [x] = 0 then [z] = 0 */

416 if (x->val == 0)

417 { xassert(x->ptr == NULL);

418 mpz_set_si(z, 0);

419 goto done;

420 }

421 /* if [y] = 0 then [z] = 0 */

422 if (y->val == 0)

423 { xassert(y->ptr == NULL);

424 mpz_set_si(z, 0);

425 goto done;

426 }

427 /* special case when both [x] and [y] are in short format */

428 if (x->ptr == NULL && y->ptr == NULL)

429 { int xval = x->val, yval = y->val, sz = +1;

430 xassert(xval != 0x80000000 && yval != 0x80000000);

431 if (xval < 0) xval = - xval, sz = - sz;

432 if (yval < 0) yval = - yval, sz = - sz;

433 if (xval <= 0x7FFFFFFF / yval)

434 { mpz_set_si(z, sz * (xval * yval));

435 goto done;

436 }

437 }

438 /* convert [x] to long format, if necessary */

439 if (x->ptr == NULL)

440 { xassert(x->val != 0x80000000);

441 if (x->val >= 0)

442 { sx = +1;

443 t = (unsigned int)(+ x->val);

444 }

445 else

446 { sx = -1;

447 t = (unsigned int)(- x->val);

448 }

449 ex = &dumx;

450 ex->d[0] = (unsigned short)t;

451 ex->d[1] = (unsigned short)(t >> 16);

452 ex->d[2] = ex->d[3] = ex->d[4] = ex->d[5] = 0;

453 ex->next = NULL;

454 }

455 else

456 { sx = x->val;

457 xassert(sx == +1 || sx == -1);

458 ex = x->ptr;

459 }

460 /* convert [y] to long format, if necessary */

461 if (y->ptr == NULL)

462 { xassert(y->val != 0x80000000);

463 if (y->val >= 0)

464 { sy = +1;

465 t = (unsigned int)(+ y->val);

466 }

467 else

468 { sy = -1;

469 t = (unsigned int)(- y->val);

470 }

471 ey = &dumy;

472 ey->d[0] = (unsigned short)t;

473 ey->d[1] = (unsigned short)(t >> 16);

474 ey->d[2] = ey->d[3] = ey->d[4] = ey->d[5] = 0;

475 ey->next = NULL;

476 }

477 else

478 { sy = y->val;

479 xassert(sy == +1 || sy == -1);

480 ey = y->ptr;

481 }

482 /* determine the number of digits of [x] */

483 nx = n = 0;

484 for (e = ex; e != NULL; e = e->next)

485 for (k = 0; k <= 5; k++)

486 { n++;

487 if (e->d[k]) nx = n;

488 }

489 xassert(nx > 0);

490 /* determine the number of digits of [y] */

491 ny = n = 0;

492 for (e = ey; e != NULL; e = e->next)

493 for (k = 0; k <= 5; k++)

494 { n++;

495 if (e->d[k]) ny = n;

496 }

497 xassert(ny > 0);

498 /* we need working array containing at least nx+ny+ny places */

499 work = gmp_get_work(nx+ny+ny);

500 /* load digits of [x] */

501 wx = &work[0];

502 for (n = 0; n < nx; n++) wx[ny+n] = 0;

503 for (n = 0, e = ex; e != NULL; e = e->next)

504 for (k = 0; k <= 5; k++, n++)

505 if (e->d[k]) wx[ny+n] = e->d[k];

506 /* load digits of [y] */

507 wy = &work[nx+ny];

508 for (n = 0; n < ny; n++) wy[n] = 0;

509 for (n = 0, e = ey; e != NULL; e = e->next)

510 for (k = 0; k <= 5; k++, n++)

511 if (e->d[k]) wy[n] = e->d[k];

512 /* compute [x] * [y] */

513 bigmul(nx, ny, wx, wy);

514 /* construct and normalize result */

515 mpz_set_si(z, 0);

516 z->val = sx * sy;

517 es = NULL;

518 k = 6;

519 for (n = 0; n < nx+ny; n++)

520 { if (k > 5)

521 { e = gmp_get_atom(sizeof(struct mpz_seg));

522 e->d[0] = e->d[1] = e->d[2] = 0;

523 e->d[3] = e->d[4] = e->d[5] = 0;

524 e->next = NULL;

525 if (z->ptr == NULL)

526 z->ptr = e;

527 else

528 es->next = e;

529 es = e;

530 k = 0;

531 }

532 es->d[k++] = wx[n];

533 }

534 normalize(z);

535 done: return;

536 }

537

538 void mpz_neg(mpz_t z, mpz_t x)

539 { /* set z to 0 - x */

540 mpz_set(z, x);

541 z->val = - z->val;

542 return;

543 }

544

545 void mpz_abs(mpz_t z, mpz_t x)

546 { /* set z to the absolute value of x */

547 mpz_set(z, x);

548 if (z->val < 0) z->val = - z->val;

549 return;

550 }

551

552 void mpz_div(mpz_t q, mpz_t r, mpz_t x, mpz_t y)

553 { /* divide x by y, forming quotient q and/or remainder r

554 if q = NULL then quotient is not stored; if r = NULL then

555 remainder is not stored

556 the sign of quotient is determined as in algebra while the

557 sign of remainder is the same as the sign of dividend:

558 +26 : +7 = +3, remainder is +5

559 -26 : +7 = -3, remainder is -5

560 +26 : -7 = -3, remainder is +5

561 -26 : -7 = +3, remainder is -5 */

562 struct mpz_seg dumx, dumy, *ex, *ey, *es, *e;

563 int sx, sy, k, nx, ny, n;

564 unsigned int t;

565 unsigned short *work, *wx, *wy;

566 /* divide by zero is not allowed */

567 if (y->val == 0)

568 { xassert(y->ptr == NULL);

569 xfault("mpz_div: divide by zero not allowed\n");

570 }

571 /* if [x] = 0 then [q] = [r] = 0 */

572 if (x->val == 0)

573 { xassert(x->ptr == NULL);

574 if (q != NULL) mpz_set_si(q, 0);

575 if (r != NULL) mpz_set_si(r, 0);

576 goto done;

577 }

578 /* special case when both [x] and [y] are in short format */

579 if (x->ptr == NULL && y->ptr == NULL)

580 { int xval = x->val, yval = y->val;

581 xassert(xval != 0x80000000 && yval != 0x80000000);

582 if (q != NULL) mpz_set_si(q, xval / yval);

583 if (r != NULL) mpz_set_si(r, xval % yval);

584 goto done;

585 }

586 /* convert [x] to long format, if necessary */

587 if (x->ptr == NULL)

588 { xassert(x->val != 0x80000000);

589 if (x->val >= 0)

590 { sx = +1;

591 t = (unsigned int)(+ x->val);

592 }

593 else

594 { sx = -1;

595 t = (unsigned int)(- x->val);

596 }

597 ex = &dumx;

598 ex->d[0] = (unsigned short)t;

599 ex->d[1] = (unsigned short)(t >> 16);

600 ex->d[2] = ex->d[3] = ex->d[4] = ex->d[5] = 0;

601 ex->next = NULL;

602 }

603 else

604 { sx = x->val;

605 xassert(sx == +1 || sx == -1);

606 ex = x->ptr;

607 }

608 /* convert [y] to long format, if necessary */

609 if (y->ptr == NULL)

610 { xassert(y->val != 0x80000000);

611 if (y->val >= 0)

612 { sy = +1;

613 t = (unsigned int)(+ y->val);

614 }

615 else

616 { sy = -1;

617 t = (unsigned int)(- y->val);

618 }

619 ey = &dumy;

620 ey->d[0] = (unsigned short)t;

621 ey->d[1] = (unsigned short)(t >> 16);

622 ey->d[2] = ey->d[3] = ey->d[4] = ey->d[5] = 0;

623 ey->next = NULL;

624 }

625 else

626 { sy = y->val;

627 xassert(sy == +1 || sy == -1);

628 ey = y->ptr;

629 }

630 /* determine the number of digits of [x] */

631 nx = n = 0;

632 for (e = ex; e != NULL; e = e->next)

633 for (k = 0; k <= 5; k++)

634 { n++;

635 if (e->d[k]) nx = n;

636 }

637 xassert(nx > 0);

638 /* determine the number of digits of [y] */

639 ny = n = 0;

640 for (e = ey; e != NULL; e = e->next)

641 for (k = 0; k <= 5; k++)

642 { n++;

643 if (e->d[k]) ny = n;

644 }

645 xassert(ny > 0);

646 /* if nx < ny then [q] = 0 and [r] = [x] */

647 if (nx < ny)

648 { if (r != NULL) mpz_set(r, x);

649 if (q != NULL) mpz_set_si(q, 0);

650 goto done;

651 }

652 /* we need working array containing at least nx+ny+1 places */

653 work = gmp_get_work(nx+ny+1);

654 /* load digits of [x] */

655 wx = &work[0];

656 for (n = 0; n < nx; n++) wx[n] = 0;

657 for (n = 0, e = ex; e != NULL; e = e->next)

658 for (k = 0; k <= 5; k++, n++)

659 if (e->d[k]) wx[n] = e->d[k];

660 /* load digits of [y] */

661 wy = &work[nx+1];

662 for (n = 0; n < ny; n++) wy[n] = 0;

663 for (n = 0, e = ey; e != NULL; e = e->next)

664 for (k = 0; k <= 5; k++, n++)

665 if (e->d[k]) wy[n] = e->d[k];

666 /* compute quotient and remainder */

667 xassert(wy[ny-1] != 0);

668 bigdiv(nx-ny, ny, wx, wy);

669 /* construct and normalize quotient */

670 if (q != NULL)

671 { mpz_set_si(q, 0);

672 q->val = sx * sy;

673 es = NULL;

674 k = 6;

675 for (n = ny; n <= nx; n++)

676 { if (k > 5)

677 { e = gmp_get_atom(sizeof(struct mpz_seg));

678 e->d[0] = e->d[1] = e->d[2] = 0;

679 e->d[3] = e->d[4] = e->d[5] = 0;

680 e->next = NULL;

681 if (q->ptr == NULL)

682 q->ptr = e;

683 else

684 es->next = e;

685 es = e;

686 k = 0;

687 }

688 es->d[k++] = wx[n];

689 }

690 normalize(q);

691 }

692 /* construct and normalize remainder */

693 if (r != NULL)

694 { mpz_set_si(r, 0);

695 r->val = sx;

696 es = NULL;

697 k = 6;

698 for (n = 0; n < ny; n++)

699 { if (k > 5)

700 { e = gmp_get_atom(sizeof(struct mpz_seg));

701 e->d[0] = e->d[1] = e->d[2] = 0;

702 e->d[3] = e->d[4] = e->d[5] = 0;

703 e->next = NULL;

704 if (r->ptr == NULL)

705 r->ptr = e;

706 else

707 es->next = e;

708 es = e;

709 k = 0;

710 }

711 es->d[k++] = wx[n];

712 }

713 normalize(r);

714 }

715 done: return;

716 }

717

718 void mpz_gcd(mpz_t z, mpz_t x, mpz_t y)

719 { /* set z to the greatest common divisor of x and y */

720 /* in case of arbitrary integers GCD(x, y) = GCD(|x|, |y|), and,

721 in particular, GCD(0, 0) = 0 */

722 mpz_t u, v, r;

723 mpz_init(u);

724 mpz_init(v);

725 mpz_init(r);

726 mpz_abs(u, x);

727 mpz_abs(v, y);

728 while (mpz_sgn(v))

729 { mpz_div(NULL, r, u, v);

730 mpz_set(u, v);

731 mpz_set(v, r);

732 }

733 mpz_set(z, u);

734 mpz_clear(u);

735 mpz_clear(v);

736 mpz_clear(r);

737 return;

738 }

739

740 int mpz_cmp(mpz_t x, mpz_t y)

741 { /* compare x and y; return a positive value if x > y, zero if

742 x = y, or a nefative value if x < y */

743 static struct mpz_seg zero = { { 0, 0, 0, 0, 0, 0 }, NULL };

744 struct mpz_seg dumx, dumy, *ex, *ey;

745 int cc, sx, sy, k;

746 unsigned int t;

747 if (x == y)

748 { cc = 0;

749 goto done;

750 }

751 /* special case when both [x] and [y] are in short format */

752 if (x->ptr == NULL && y->ptr == NULL)

753 { int xval = x->val, yval = y->val;

754 xassert(xval != 0x80000000 && yval != 0x80000000);

755 cc = (xval > yval ? +1 : xval < yval ? -1 : 0);

756 goto done;

757 }

758 /* special case when [x] and [y] have different signs */

759 if (x->val > 0 && y->val <= 0 || x->val == 0 && y->val < 0)

760 { cc = +1;

761 goto done;

762 }

763 if (x->val < 0 && y->val >= 0 || x->val == 0 && y->val > 0)

764 { cc = -1;

765 goto done;

766 }

767 /* convert [x] to long format, if necessary */

768 if (x->ptr == NULL)

769 { xassert(x->val != 0x80000000);

770 if (x->val >= 0)

771 { sx = +1;

772 t = (unsigned int)(+ x->val);

773 }

774 else

775 { sx = -1;

776 t = (unsigned int)(- x->val);

777 }

778 ex = &dumx;

779 ex->d[0] = (unsigned short)t;

780 ex->d[1] = (unsigned short)(t >> 16);

781 ex->d[2] = ex->d[3] = ex->d[4] = ex->d[5] = 0;

782 ex->next = NULL;

783 }

784 else

785 { sx = x->val;

786 xassert(sx == +1 || sx == -1);

787 ex = x->ptr;

788 }

789 /* convert [y] to long format, if necessary */

790 if (y->ptr == NULL)

791 { xassert(y->val != 0x80000000);

792 if (y->val >= 0)

793 { sy = +1;

794 t = (unsigned int)(+ y->val);

795 }

796 else

797 { sy = -1;

798 t = (unsigned int)(- y->val);

799 }

800 ey = &dumy;

801 ey->d[0] = (unsigned short)t;

802 ey->d[1] = (unsigned short)(t >> 16);

803 ey->d[2] = ey->d[3] = ey->d[4] = ey->d[5] = 0;

804 ey->next = NULL;

805 }

806 else

807 { sy = y->val;

808 xassert(sy == +1 || sy == -1);

809 ey = y->ptr;

810 }

811 /* main fragment */

812 xassert(sx > 0 && sy > 0 || sx < 0 && sy < 0);

813 cc = 0;

814 for (; ex || ey; ex = ex->next, ey = ey->next)

815 { if (ex == NULL) ex = &zero;

816 if (ey == NULL) ey = &zero;

817 for (k = 0; k <= 5; k++)

818 { if (ex->d[k] > ey->d[k]) cc = +1;

819 if (ex->d[k] < ey->d[k]) cc = -1;

820 }

821 }

822 if (sx < 0) cc = - cc;

823 done: return cc;

824 }

825

826 int mpz_sgn(mpz_t x)

827 { /* return +1 if x > 0, 0 if x = 0, and -1 if x < 0 */

828 int s;

829 s = (x->val > 0 ? +1 : x->val < 0 ? -1 : 0);

830 return s;

831 }

832

833 int mpz_out_str(void *_fp, int base, mpz_t x)

834 { /* output x on stream fp, as a string in given base; the base

835 may vary from 2 to 36;

836 return the number of bytes written, or if an error occurred,

837 return 0 */

838 FILE *fp = _fp;

839 mpz_t b, y, r;

840 int n, j, nwr = 0;

841 unsigned char *d;

842 static char *set = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ";

843 if (!(2 <= base && base <= 36))

844 xfault("mpz_out_str: base = %d; invalid base\n", base);

845 mpz_init(b);

846 mpz_set_si(b, base);

847 mpz_init(y);

848 mpz_init(r);

849 /* determine the number of digits */

850 mpz_abs(y, x);

851 for (n = 0; mpz_sgn(y) != 0; n++)

852 mpz_div(y, NULL, y, b);

853 if (n == 0) n = 1;

854 /* compute the digits */

855 d = xmalloc(n);

856 mpz_abs(y, x);

857 for (j = 0; j < n; j++)

858 { mpz_div(y, r, y, b);

859 xassert(0 <= r->val && r->val < base && r->ptr == NULL);

860 d[j] = (unsigned char)r->val;

861 }

862 /* output the integer to the stream */

863 if (fp == NULL) fp = stdout;

864 if (mpz_sgn(x) < 0)

865 fputc('-', fp), nwr++;

866 for (j = n-1; j >= 0; j--)

867 fputc(set[d[j]], fp), nwr++;

868 if (ferror(fp)) nwr = 0;

869 mpz_clear(b);

870 mpz_clear(y);

871 mpz_clear(r);

872 xfree(d);

873 return nwr;

874 }

875

876 /*====================================================================*/

877

878 mpq_t _mpq_init(void)

879 { /* initialize x, and set its value to 0/1 */

880 mpq_t x;

881 x = gmp_get_atom(sizeof(struct mpq));

882 x->p.val = 0;

883 x->p.ptr = NULL;

884 x->q.val = 1;

885 x->q.ptr = NULL;

886 return x;

887 }

888

889 void mpq_clear(mpq_t x)

890 { /* free the space occupied by x */

891 mpz_set_si(&x->p, 0);

892 xassert(x->p.ptr == NULL);

893 mpz_set_si(&x->q, 0);

894 xassert(x->q.ptr == NULL);

895 /* free the number descriptor */

896 gmp_free_atom(x, sizeof(struct mpq));

897 return;

898 }

899

900 void mpq_canonicalize(mpq_t x)

901 { /* remove any factors that are common to the numerator and

902 denominator of x, and make the denominator positive */

903 mpz_t f;

904 xassert(x->q.val != 0);

905 if (x->q.val < 0)

906 { mpz_neg(&x->p, &x->p);

907 mpz_neg(&x->q, &x->q);

908 }

909 mpz_init(f);

910 mpz_gcd(f, &x->p, &x->q);

911 if (!(f->val == 1 && f->ptr == NULL))

912 { mpz_div(&x->p, NULL, &x->p, f);

913 mpz_div(&x->q, NULL, &x->q, f);

914 }

915 mpz_clear(f);

916 return;

917 }

918

919 void mpq_set(mpq_t z, mpq_t x)

920 { /* set the value of z from x */

921 if (z != x)

922 { mpz_set(&z->p, &x->p);

923 mpz_set(&z->q, &x->q);

924 }

925 return;

926 }

927

928 void mpq_set_si(mpq_t x, int p, unsigned int q)

929 { /* set the value of x to p/q */

930 if (q == 0)

931 xfault("mpq_set_si: zero denominator not allowed\n");

932 mpz_set_si(&x->p, p);

933 xassert(q <= 0x7FFFFFFF);

934 mpz_set_si(&x->q, q);

935 return;

936 }

937

938 double mpq_get_d(mpq_t x)

939 { /* convert x to a double, truncating if necessary */

940 int np, nq;

941 double p, q;

942 p = mpz_get_d_2exp(&np, &x->p);

943 q = mpz_get_d_2exp(&nq, &x->q);

944 return ldexp(p / q, np - nq);

945 }

946

947 void mpq_set_d(mpq_t x, double val)

948 { /* set x to val; there is no rounding, the conversion is exact */

949 int s, n, d, j;

950 double f;

951 mpz_t temp;

952 xassert(-DBL_MAX <= val && val <= +DBL_MAX);

953 mpq_set_si(x, 0, 1);

954 if (val > 0.0)

955 s = +1;

956 else if (val < 0.0)

957 s = -1;

958 else

959 goto done;

960 f = frexp(fabs(val), &n);

961 /* |val| = f * 2^n, where 0.5 <= f < 1.0 */

962 mpz_init(temp);

963 while (f != 0.0)

964 { f *= 16.0, n -= 4;

965 d = (int)f;

966 xassert(0 <= d && d <= 15);

967 f -= (double)d;

968 /* x := 16 * x + d */

969 mpz_set_si(temp, 16);

970 mpz_mul(&x->p, &x->p, temp);

971 mpz_set_si(temp, d);

972 mpz_add(&x->p, &x->p, temp);

973 }

974 mpz_clear(temp);

975 /* x := x * 2^n */

976 if (n > 0)

977 { for (j = 1; j <= n; j++)

978 mpz_add(&x->p, &x->p, &x->p);

979 }

980 else if (n < 0)

981 { for (j = 1; j <= -n; j++)

982 mpz_add(&x->q, &x->q, &x->q);

983 mpq_canonicalize(x);

984 }

985 if (s < 0) mpq_neg(x, x);

986 done: return;

987 }

988

989 void mpq_add(mpq_t z, mpq_t x, mpq_t y)

990 { /* set z to x + y */

991 mpz_t p, q;

992 mpz_init(p);

993 mpz_init(q);

994 mpz_mul(p, &x->p, &y->q);

995 mpz_mul(q, &x->q, &y->p);

996 mpz_add(p, p, q);

997 mpz_mul(q, &x->q, &y->q);

998 mpz_set(&z->p, p);

999 mpz_set(&z->q, q);

1000 mpz_clear(p);

1001 mpz_clear(q);

1002 mpq_canonicalize(z);

1003 return;

1004 }

1005

1006 void mpq_sub(mpq_t z, mpq_t x, mpq_t y)

1007 { /* set z to x - y */

1008 mpz_t p, q;

1009 mpz_init(p);

1010 mpz_init(q);

1011 mpz_mul(p, &x->p, &y->q);

1012 mpz_mul(q, &x->q, &y->p);

1013 mpz_sub(p, p, q);

1014 mpz_mul(q, &x->q, &y->q);

1015 mpz_set(&z->p, p);

1016 mpz_set(&z->q, q);

1017 mpz_clear(p);

1018 mpz_clear(q);

1019 mpq_canonicalize(z);

1020 return;

1021 }

1022

1023 void mpq_mul(mpq_t z, mpq_t x, mpq_t y)

1024 { /* set z to x * y */

1025 mpz_mul(&z->p, &x->p, &y->p);

1026 mpz_mul(&z->q, &x->q, &y->q);

1027 mpq_canonicalize(z);

1028 return;

1029 }

1030

1031 void mpq_div(mpq_t z, mpq_t x, mpq_t y)

1032 { /* set z to x / y */

1033 mpz_t p, q;

1034 if (mpq_sgn(y) == 0)

1035 xfault("mpq_div: zero divisor not allowed\n");

1036 mpz_init(p);

1037 mpz_init(q);

1038 mpz_mul(p, &x->p, &y->q);

1039 mpz_mul(q, &x->q, &y->p);

1040 mpz_set(&z->p, p);

1041 mpz_set(&z->q, q);

1042 mpz_clear(p);

1043 mpz_clear(q);

1044 mpq_canonicalize(z);

1045 return;

1046 }

1047

1048 void mpq_neg(mpq_t z, mpq_t x)

1049 { /* set z to 0 - x */

1050 mpq_set(z, x);

1051 mpz_neg(&z->p, &z->p);

1052 return;

1053 }

1054

1055 void mpq_abs(mpq_t z, mpq_t x)

1056 { /* set z to the absolute value of x */

1057 mpq_set(z, x);

1058 mpz_abs(&z->p, &z->p);

1059 xassert(mpz_sgn(&x->q) > 0);

1060 return;

1061 }

1062

1063 int mpq_cmp(mpq_t x, mpq_t y)

1064 { /* compare x and y; return a positive value if x > y, zero if

1065 x = y, or a nefative value if x < y */

1066 mpq_t temp;

1067 int s;

1068 mpq_init(temp);

1069 mpq_sub(temp, x, y);

1070 s = mpq_sgn(temp);

1071 mpq_clear(temp);

1072 return s;

1073 }

1074

1075 int mpq_sgn(mpq_t x)

1076 { /* return +1 if x > 0, 0 if x = 0, and -1 if x < 0 */

1077 int s;

1078 s = mpz_sgn(&x->p);

1079 xassert(mpz_sgn(&x->q) > 0);

1080 return s;

1081 }

1082

1083 int mpq_out_str(void *_fp, int base, mpq_t x)

1084 { /* output x on stream fp, as a string in given base; the base

1085 may vary from 2 to 36; output is in the form 'num/den' or if

1086 the denominator is 1 then just 'num';

1087 if the parameter fp is a null pointer, stdout is assumed;

1088 return the number of bytes written, or if an error occurred,

1089 return 0 */

1090 FILE *fp = _fp;

1091 int nwr;

1092 if (!(2 <= base && base <= 36))

1093 xfault("mpq_out_str: base = %d; invalid base\n", base);

1094 if (fp == NULL) fp = stdout;

1095 nwr = mpz_out_str(fp, base, &x->p);

1096 if (x->q.val == 1 && x->q.ptr == NULL)

1097 ;

1098 else

1099 { fputc('/', fp), nwr++;

1100 nwr += mpz_out_str(fp, base, &x->q);

1101 }

1102 if (ferror(fp)) nwr = 0;

1103 return nwr;

1104 }

1105

1106 #endif

1107

1108 /* eof */