COIN-OR::LEMON - Graph Library

source: glpk-cmake/src/glpgmp.c @ 1:c445c931472f

Last change on this file since 1:c445c931472f was 1:c445c931472f, checked in by Alpar Juttner <alpar@…>, 10 years ago

Import glpk-4.45

  • Generated files and doc/notes are removed
File size: 30.9 KB
Line 
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
32int gmp_pool_count(void) { return 0; }
33
34void gmp_free_mem(void) { return; }
35
36#else                         /* use GLPK bignum module */
37
38static DMP *gmp_pool = NULL;
39static int gmp_size = 0;
40static unsigned short *gmp_work = NULL;
41
42void *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
48void gmp_free_atom(void *ptr, int size)
49{     xassert(gmp_pool != NULL);
50      dmp_free_atom(gmp_pool, ptr, size);
51      return;
52}
53
54int gmp_pool_count(void)
55{     if (gmp_pool == NULL)
56         return 0;
57      else
58         return dmp_in_use(gmp_pool).lo;
59}
60
61unsigned 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
78void 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
89mpz_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
98void 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
107void 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
128void 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
153double 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
175double 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
203void 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
213static 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      }
249done: return;
250}
251
252void 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);
394done: return;
395}
396
397void 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
409void 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);
535done: return;
536}
537
538void 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
545void 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
552void 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      }
715done: return;
716}
717
718void 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
740int 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;
823done: return cc;
824}
825
826int 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
833int 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
878mpq_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
889void 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
900void 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
919void 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
928void 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
938double 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
947void 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);
986done: return;
987}
988
989void 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
1006void 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
1023void 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
1031void 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
1048void 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
1055void 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
1063int 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
1075int 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
1083int 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 */
Note: See TracBrowser for help on using the repository browser.