1 | /* glpgmp.c */ |
---|
2 | |
---|
3 | /*********************************************************************** |
---|
4 | * This code is part of GLPK (GNU Linear Programming Kit). |
---|
5 | * |
---|
6 | * Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, |
---|
7 | * 2009, 2010 Andrew Makhorin, Department for Applied Informatics, |
---|
8 | * Moscow Aviation Institute, Moscow, Russia. All rights reserved. |
---|
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 | ***********************************************************************/ |
---|
24 | |
---|
25 | #define _GLPSTD_STDIO |
---|
26 | #include "glpdmp.h" |
---|
27 | #include "glpgmp.h" |
---|
28 | #define xfault xerror |
---|
29 | |
---|
30 | #ifdef HAVE_GMP /* use GNU MP bignum library */ |
---|
31 | |
---|
32 | int gmp_pool_count(void) { return 0; } |
---|
33 | |
---|
34 | void gmp_free_mem(void) { return; } |
---|
35 | |
---|
36 | #else /* use GLPK bignum module */ |
---|
37 | |
---|
38 | static DMP *gmp_pool = NULL; |
---|
39 | static int gmp_size = 0; |
---|
40 | static unsigned short *gmp_work = NULL; |
---|
41 | |
---|
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 | } |
---|
47 | |
---|
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 | } |
---|
53 | |
---|
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 | } |
---|
60 | |
---|
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 | } |
---|
77 | |
---|
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 | } |
---|
86 | |
---|
87 | /*====================================================================*/ |
---|
88 | |
---|
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 | } |
---|
97 | |
---|
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 */ |
---|