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 
gmp_pool_count(void)32 int gmp_pool_count(void) { return 0; }
33 
gmp_free_mem(void)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 
gmp_get_atom(int size)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 
gmp_free_atom(void * ptr,int size)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 
gmp_pool_count(void)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 
gmp_get_work(int size)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 
gmp_free_mem(void)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 
_mpz_init(void)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 
mpz_clear(mpz_t x)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 
mpz_set(mpz_t z,mpz_t x)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 
mpz_set_si(mpz_t x,int val)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 
mpz_get_d(mpz_t x)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 
mpz_get_d_2exp(int * exp,mpz_t x)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 
mpz_swap(mpz_t x,mpz_t y)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 
normalize(mpz_t x)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 
mpz_add(mpz_t z,mpz_t x,mpz_t y)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 
mpz_sub(mpz_t z,mpz_t x,mpz_t y)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 
mpz_mul(mpz_t z,mpz_t x,mpz_t y)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 
mpz_neg(mpz_t z,mpz_t x)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 
mpz_abs(mpz_t z,mpz_t x)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 
mpz_div(mpz_t q,mpz_t r,mpz_t x,mpz_t y)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 
mpz_gcd(mpz_t z,mpz_t x,mpz_t y)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 
mpz_cmp(mpz_t x,mpz_t y)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 
mpz_sgn(mpz_t x)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 
mpz_out_str(void * _fp,int base,mpz_t x)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 
_mpq_init(void)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 
mpq_clear(mpq_t x)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 
mpq_canonicalize(mpq_t x)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 
mpq_set(mpq_t z,mpq_t x)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 
mpq_set_si(mpq_t x,int p,unsigned int q)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 
mpq_get_d(mpq_t x)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 
mpq_set_d(mpq_t x,double val)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 
mpq_add(mpq_t z,mpq_t x,mpq_t y)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 
mpq_sub(mpq_t z,mpq_t x,mpq_t y)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 
mpq_mul(mpq_t z,mpq_t x,mpq_t y)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 
mpq_div(mpq_t z,mpq_t x,mpq_t y)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 
mpq_neg(mpq_t z,mpq_t x)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 
mpq_abs(mpq_t z,mpq_t x)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 
mpq_cmp(mpq_t x,mpq_t y)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 
mpq_sgn(mpq_t x)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 
mpq_out_str(void * _fp,int base,mpq_t x)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 */
1109