You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
1116 lines
31 KiB
1116 lines
31 KiB
/* glpgmp.c */
|
|
|
|
/***********************************************************************
|
|
* This code is part of GLPK (GNU Linear Programming Kit).
|
|
*
|
|
* Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
|
|
* 2009, 2010, 2011, 2013 Andrew Makhorin, Department for Applied
|
|
* Informatics, Moscow Aviation Institute, Moscow, Russia. All rights
|
|
* reserved. E-mail: <mao@gnu.org>.
|
|
*
|
|
* GLPK is free software: you can redistribute it and/or modify it
|
|
* under the terms of the GNU General Public License as published by
|
|
* the Free Software Foundation, either version 3 of the License, or
|
|
* (at your option) any later version.
|
|
*
|
|
* GLPK is distributed in the hope that it will be useful, but WITHOUT
|
|
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
|
|
* or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
|
|
* License for more details.
|
|
*
|
|
* You should have received a copy of the GNU General Public License
|
|
* along with GLPK. If not, see <http://www.gnu.org/licenses/>.
|
|
***********************************************************************/
|
|
|
|
#define _GLPSTD_STDIO
|
|
#if 1 /* 11/VI-2013 */
|
|
#include "bignum.h"
|
|
#endif
|
|
#include "dmp.h"
|
|
#include "glpgmp.h"
|
|
#include "env.h"
|
|
#define xfault xerror
|
|
|
|
#ifdef HAVE_GMP /* use GNU MP bignum library */
|
|
|
|
int gmp_pool_count(void) { return 0; }
|
|
|
|
void gmp_free_mem(void) { return; }
|
|
|
|
#else /* use GLPK bignum module */
|
|
|
|
static DMP *gmp_pool = NULL;
|
|
static int gmp_size = 0;
|
|
static unsigned short *gmp_work = NULL;
|
|
|
|
void *gmp_get_atom(int size)
|
|
{ if (gmp_pool == NULL)
|
|
gmp_pool = dmp_create_pool();
|
|
return dmp_get_atom(gmp_pool, size);
|
|
}
|
|
|
|
void gmp_free_atom(void *ptr, int size)
|
|
{ xassert(gmp_pool != NULL);
|
|
dmp_free_atom(gmp_pool, ptr, size);
|
|
return;
|
|
}
|
|
|
|
int gmp_pool_count(void)
|
|
{ if (gmp_pool == NULL)
|
|
return 0;
|
|
else
|
|
#if 0 /* 10/VI-2013 */
|
|
return dmp_in_use(gmp_pool).lo;
|
|
#else
|
|
return dmp_in_use(gmp_pool);
|
|
#endif
|
|
}
|
|
|
|
unsigned short *gmp_get_work(int size)
|
|
{ xassert(size > 0);
|
|
if (gmp_size < size)
|
|
{ if (gmp_size == 0)
|
|
{ xassert(gmp_work == NULL);
|
|
gmp_size = 100;
|
|
}
|
|
else
|
|
{ xassert(gmp_work != NULL);
|
|
xfree(gmp_work);
|
|
}
|
|
while (gmp_size < size) gmp_size += gmp_size;
|
|
gmp_work = xcalloc(gmp_size, sizeof(unsigned short));
|
|
}
|
|
return gmp_work;
|
|
}
|
|
|
|
void gmp_free_mem(void)
|
|
{ if (gmp_pool != NULL) dmp_delete_pool(gmp_pool);
|
|
if (gmp_work != NULL) xfree(gmp_work);
|
|
gmp_pool = NULL;
|
|
gmp_size = 0;
|
|
gmp_work = NULL;
|
|
return;
|
|
}
|
|
|
|
/*====================================================================*/
|
|
|
|
mpz_t _mpz_init(void)
|
|
{ /* initialize x, and set its value to 0 */
|
|
mpz_t x;
|
|
x = gmp_get_atom(sizeof(struct mpz));
|
|
x->val = 0;
|
|
x->ptr = NULL;
|
|
return x;
|
|
}
|
|
|
|
void mpz_clear(mpz_t x)
|
|
{ /* free the space occupied by x */
|
|
mpz_set_si(x, 0);
|
|
xassert(x->ptr == NULL);
|
|
/* free the number descriptor */
|
|
gmp_free_atom(x, sizeof(struct mpz));
|
|
return;
|
|
}
|
|
|
|
void mpz_set(mpz_t z, mpz_t x)
|
|
{ /* set the value of z from x */
|
|
struct mpz_seg *e, *ee, *es;
|
|
if (z != x)
|
|
{ mpz_set_si(z, 0);
|
|
z->val = x->val;
|
|
xassert(z->ptr == NULL);
|
|
for (e = x->ptr, es = NULL; e != NULL; e = e->next)
|
|
{ ee = gmp_get_atom(sizeof(struct mpz_seg));
|
|
memcpy(ee->d, e->d, 12);
|
|
ee->next = NULL;
|
|
if (z->ptr == NULL)
|
|
z->ptr = ee;
|
|
else
|
|
es->next = ee;
|
|
es = ee;
|
|
}
|
|
}
|
|
return;
|
|
}
|
|
|
|
void mpz_set_si(mpz_t x, int val)
|
|
{ /* set the value of x to val */
|
|
struct mpz_seg *e;
|
|
/* free existing segments, if any */
|
|
while (x->ptr != NULL)
|
|
{ e = x->ptr;
|
|
x->ptr = e->next;
|
|
gmp_free_atom(e, sizeof(struct mpz_seg));
|
|
}
|
|
/* assign new value */
|
|
if (val == 0x80000000)
|
|
{ /* long format is needed */
|
|
x->val = -1;
|
|
x->ptr = e = gmp_get_atom(sizeof(struct mpz_seg));
|
|
memset(e->d, 0, 12);
|
|
e->d[1] = 0x8000;
|
|
e->next = NULL;
|
|
}
|
|
else
|
|
{ /* short format is enough */
|
|
x->val = val;
|
|
}
|
|
return;
|
|
}
|
|
|
|
double mpz_get_d(mpz_t x)
|
|
{ /* convert x to a double, truncating if necessary */
|
|
struct mpz_seg *e;
|
|
int j;
|
|
double val, deg;
|
|
if (x->ptr == NULL)
|
|
val = (double)x->val;
|
|
else
|
|
{ xassert(x->val != 0);
|
|
val = 0.0;
|
|
deg = 1.0;
|
|
for (e = x->ptr; e != NULL; e = e->next)
|
|
{ for (j = 0; j <= 5; j++)
|
|
{ val += deg * (double)((int)e->d[j]);
|
|
deg *= 65536.0;
|
|
}
|
|
}
|
|
if (x->val < 0) val = - val;
|
|
}
|
|
return val;
|
|
}
|
|
|
|
double mpz_get_d_2exp(int *exp, mpz_t x)
|
|
{ /* convert x to a double, truncating if necessary (i.e. rounding
|
|
towards zero), and returning the exponent separately;
|
|
the return value is in the range 0.5 <= |d| < 1 and the
|
|
exponent is stored to *exp; d*2^exp is the (truncated) x value;
|
|
if x is zero, the return is 0.0 and 0 is stored to *exp;
|
|
this is similar to the standard C frexp function */
|
|
struct mpz_seg *e;
|
|
int j, n, n1;
|
|
double val;
|
|
if (x->ptr == NULL)
|
|
val = (double)x->val, n = 0;
|
|
else
|
|
{ xassert(x->val != 0);
|
|
val = 0.0, n = 0;
|
|
for (e = x->ptr; e != NULL; e = e->next)
|
|
{ for (j = 0; j <= 5; j++)
|
|
{ val += (double)((int)e->d[j]);
|
|
val /= 65536.0, n += 16;
|
|
}
|
|
}
|
|
if (x->val < 0) val = - val;
|
|
}
|
|
val = frexp(val, &n1);
|
|
*exp = n + n1;
|
|
return val;
|
|
}
|
|
|
|
void mpz_swap(mpz_t x, mpz_t y)
|
|
{ /* swap the values x and y efficiently */
|
|
int val;
|
|
void *ptr;
|
|
val = x->val, ptr = x->ptr;
|
|
x->val = y->val, x->ptr = y->ptr;
|
|
y->val = val, y->ptr = ptr;
|
|
return;
|
|
}
|
|
|
|
static void normalize(mpz_t x)
|
|
{ /* normalize integer x that includes removing non-significant
|
|
(leading) zeros and converting to short format, if possible */
|
|
struct mpz_seg *es, *e;
|
|
/* if the integer is in short format, it remains unchanged */
|
|
if (x->ptr == NULL)
|
|
{ xassert(x->val != 0x80000000);
|
|
goto done;
|
|
}
|
|
xassert(x->val == +1 || x->val == -1);
|
|
/* find the last (most significant) non-zero segment */
|
|
es = NULL;
|
|
for (e = x->ptr; e != NULL; e = e->next)
|
|
{ if (e->d[0] || e->d[1] || e->d[2] ||
|
|
e->d[3] || e->d[4] || e->d[5]) es = e;
|
|
}
|
|
/* if all segments contain zeros, the integer is zero */
|
|
if (es == NULL)
|
|
{ mpz_set_si(x, 0);
|
|
goto done;
|
|
}
|
|
/* remove non-significant (leading) zero segments */
|
|
while (es->next != NULL)
|
|
{ e = es->next;
|
|
es->next = e->next;
|
|
gmp_free_atom(e, sizeof(struct mpz_seg));
|
|
}
|
|
/* convert the integer to short format, if possible */
|
|
e = x->ptr;
|
|
if (e->next == NULL && e->d[1] <= 0x7FFF &&
|
|
!e->d[2] && !e->d[3] && !e->d[4] && !e->d[5])
|
|
{ int val;
|
|
val = (int)e->d[0] + ((int)e->d[1] << 16);
|
|
if (x->val < 0) val = - val;
|
|
mpz_set_si(x, val);
|
|
}
|
|
done: return;
|
|
}
|
|
|
|
void mpz_add(mpz_t z, mpz_t x, mpz_t y)
|
|
{ /* set z to x + y */
|
|
static struct mpz_seg zero = { { 0, 0, 0, 0, 0, 0 }, NULL };
|
|
struct mpz_seg dumx, dumy, *ex, *ey, *ez, *es, *ee;
|
|
int k, sx, sy, sz;
|
|
unsigned int t;
|
|
/* if [x] = 0 then [z] = [y] */
|
|
if (x->val == 0)
|
|
{ xassert(x->ptr == NULL);
|
|
mpz_set(z, y);
|
|
goto done;
|
|
}
|
|
/* if [y] = 0 then [z] = [x] */
|
|
if (y->val == 0)
|
|
{ xassert(y->ptr == NULL);
|
|
mpz_set(z, x);
|
|
goto done;
|
|
}
|
|
/* special case when both [x] and [y] are in short format */
|
|
if (x->ptr == NULL && y->ptr == NULL)
|
|
{ int xval = x->val, yval = y->val, zval = x->val + y->val;
|
|
xassert(xval != 0x80000000 && yval != 0x80000000);
|
|
if (!(xval > 0 && yval > 0 && zval <= 0 ||
|
|
xval < 0 && yval < 0 && zval >= 0))
|
|
{ mpz_set_si(z, zval);
|
|
goto done;
|
|
}
|
|
}
|
|
/* convert [x] to long format, if necessary */
|
|
if (x->ptr == NULL)
|
|
{ xassert(x->val != 0x80000000);
|
|
if (x->val >= 0)
|
|
{ sx = +1;
|
|
t = (unsigned int)(+ x->val);
|
|
}
|
|
else
|
|
{ sx = -1;
|
|
t = (unsigned int)(- x->val);
|
|
}
|
|
ex = &dumx;
|
|
ex->d[0] = (unsigned short)t;
|
|
ex->d[1] = (unsigned short)(t >> 16);
|
|
ex->d[2] = ex->d[3] = ex->d[4] = ex->d[5] = 0;
|
|
ex->next = NULL;
|
|
}
|
|
else
|
|
{ sx = x->val;
|
|
xassert(sx == +1 || sx == -1);
|
|
ex = x->ptr;
|
|
}
|
|
/* convert [y] to long format, if necessary */
|
|
if (y->ptr == NULL)
|
|
{ xassert(y->val != 0x80000000);
|
|
if (y->val >= 0)
|
|
{ sy = +1;
|
|
t = (unsigned int)(+ y->val);
|
|
}
|
|
else
|
|
{ sy = -1;
|
|
t = (unsigned int)(- y->val);
|
|
}
|
|
ey = &dumy;
|
|
ey->d[0] = (unsigned short)t;
|
|
ey->d[1] = (unsigned short)(t >> 16);
|
|
ey->d[2] = ey->d[3] = ey->d[4] = ey->d[5] = 0;
|
|
ey->next = NULL;
|
|
}
|
|
else
|
|
{ sy = y->val;
|
|
xassert(sy == +1 || sy == -1);
|
|
ey = y->ptr;
|
|
}
|
|
/* main fragment */
|
|
sz = sx;
|
|
ez = es = NULL;
|
|
if (sx > 0 && sy > 0 || sx < 0 && sy < 0)
|
|
{ /* [x] and [y] have identical signs -- addition */
|
|
t = 0;
|
|
for (; ex || ey; ex = ex->next, ey = ey->next)
|
|
{ if (ex == NULL) ex = &zero;
|
|
if (ey == NULL) ey = &zero;
|
|
ee = gmp_get_atom(sizeof(struct mpz_seg));
|
|
for (k = 0; k <= 5; k++)
|
|
{ t += (unsigned int)ex->d[k];
|
|
t += (unsigned int)ey->d[k];
|
|
ee->d[k] = (unsigned short)t;
|
|
t >>= 16;
|
|
}
|
|
ee->next = NULL;
|
|
if (ez == NULL)
|
|
ez = ee;
|
|
else
|
|
es->next = ee;
|
|
es = ee;
|
|
}
|
|
if (t)
|
|
{ /* overflow -- one extra digit is needed */
|
|
ee = gmp_get_atom(sizeof(struct mpz_seg));
|
|
ee->d[0] = 1;
|
|
ee->d[1] = ee->d[2] = ee->d[3] = ee->d[4] = ee->d[5] = 0;
|
|
ee->next = NULL;
|
|
xassert(es != NULL);
|
|
es->next = ee;
|
|
}
|
|
}
|
|
else
|
|
{ /* [x] and [y] have different signs -- subtraction */
|
|
t = 1;
|
|
for (; ex || ey; ex = ex->next, ey = ey->next)
|
|
{ if (ex == NULL) ex = &zero;
|
|
if (ey == NULL) ey = &zero;
|
|
ee = gmp_get_atom(sizeof(struct mpz_seg));
|
|
for (k = 0; k <= 5; k++)
|
|
{ t += (unsigned int)ex->d[k];
|
|
t += (0xFFFF - (unsigned int)ey->d[k]);
|
|
ee->d[k] = (unsigned short)t;
|
|
t >>= 16;
|
|
}
|
|
ee->next = NULL;
|
|
if (ez == NULL)
|
|
ez = ee;
|
|
else
|
|
es->next = ee;
|
|
es = ee;
|
|
}
|
|
if (!t)
|
|
{ /* |[x]| < |[y]| -- result in complement coding */
|
|
sz = - sz;
|
|
t = 1;
|
|
for (ee = ez; ee != NULL; ee = ee->next)
|
|
for (k = 0; k <= 5; k++)
|
|
{ t += (0xFFFF - (unsigned int)ee->d[k]);
|
|
ee->d[k] = (unsigned short)t;
|
|
t >>= 16;
|
|
}
|
|
}
|
|
}
|
|
/* contruct and normalize result */
|
|
mpz_set_si(z, 0);
|
|
z->val = sz;
|
|
z->ptr = ez;
|
|
normalize(z);
|
|
done: return;
|
|
}
|
|
|
|
void mpz_sub(mpz_t z, mpz_t x, mpz_t y)
|
|
{ /* set z to x - y */
|
|
if (x == y)
|
|
mpz_set_si(z, 0);
|
|
else
|
|
{ y->val = - y->val;
|
|
mpz_add(z, x, y);
|
|
if (y != z) y->val = - y->val;
|
|
}
|
|
return;
|
|
}
|
|
|
|
void mpz_mul(mpz_t z, mpz_t x, mpz_t y)
|
|
{ /* set z to x * y */
|
|
struct mpz_seg dumx, dumy, *ex, *ey, *es, *e;
|
|
int sx, sy, k, nx, ny, n;
|
|
unsigned int t;
|
|
unsigned short *work, *wx, *wy;
|
|
/* if [x] = 0 then [z] = 0 */
|
|
if (x->val == 0)
|
|
{ xassert(x->ptr == NULL);
|
|
mpz_set_si(z, 0);
|
|
goto done;
|
|
}
|
|
/* if [y] = 0 then [z] = 0 */
|
|
if (y->val == 0)
|
|
{ xassert(y->ptr == NULL);
|
|
mpz_set_si(z, 0);
|
|
goto done;
|
|
}
|
|
/* special case when both [x] and [y] are in short format */
|
|
if (x->ptr == NULL && y->ptr == NULL)
|
|
{ int xval = x->val, yval = y->val, sz = +1;
|
|
xassert(xval != 0x80000000 && yval != 0x80000000);
|
|
if (xval < 0) xval = - xval, sz = - sz;
|
|
if (yval < 0) yval = - yval, sz = - sz;
|
|
if (xval <= 0x7FFFFFFF / yval)
|
|
{ mpz_set_si(z, sz * (xval * yval));
|
|
goto done;
|
|
}
|
|
}
|
|
/* convert [x] to long format, if necessary */
|
|
if (x->ptr == NULL)
|
|
{ xassert(x->val != 0x80000000);
|
|
if (x->val >= 0)
|
|
{ sx = +1;
|
|
t = (unsigned int)(+ x->val);
|
|
}
|
|
else
|
|
{ sx = -1;
|
|
t = (unsigned int)(- x->val);
|
|
}
|
|
ex = &dumx;
|
|
ex->d[0] = (unsigned short)t;
|
|
ex->d[1] = (unsigned short)(t >> 16);
|
|
ex->d[2] = ex->d[3] = ex->d[4] = ex->d[5] = 0;
|
|
ex->next = NULL;
|
|
}
|
|
else
|
|
{ sx = x->val;
|
|
xassert(sx == +1 || sx == -1);
|
|
ex = x->ptr;
|
|
}
|
|
/* convert [y] to long format, if necessary */
|
|
if (y->ptr == NULL)
|
|
{ xassert(y->val != 0x80000000);
|
|
if (y->val >= 0)
|
|
{ sy = +1;
|
|
t = (unsigned int)(+ y->val);
|
|
}
|
|
else
|
|
{ sy = -1;
|
|
t = (unsigned int)(- y->val);
|
|
}
|
|
ey = &dumy;
|
|
ey->d[0] = (unsigned short)t;
|
|
ey->d[1] = (unsigned short)(t >> 16);
|
|
ey->d[2] = ey->d[3] = ey->d[4] = ey->d[5] = 0;
|
|
ey->next = NULL;
|
|
}
|
|
else
|
|
{ sy = y->val;
|
|
xassert(sy == +1 || sy == -1);
|
|
ey = y->ptr;
|
|
}
|
|
/* determine the number of digits of [x] */
|
|
nx = n = 0;
|
|
for (e = ex; e != NULL; e = e->next)
|
|
for (k = 0; k <= 5; k++)
|
|
{ n++;
|
|
if (e->d[k]) nx = n;
|
|
}
|
|
xassert(nx > 0);
|
|
/* determine the number of digits of [y] */
|
|
ny = n = 0;
|
|
for (e = ey; e != NULL; e = e->next)
|
|
for (k = 0; k <= 5; k++)
|
|
{ n++;
|
|
if (e->d[k]) ny = n;
|
|
}
|
|
xassert(ny > 0);
|
|
/* we need working array containing at least nx+ny+ny places */
|
|
work = gmp_get_work(nx+ny+ny);
|
|
/* load digits of [x] */
|
|
wx = &work[0];
|
|
for (n = 0; n < nx; n++) wx[ny+n] = 0;
|
|
for (n = 0, e = ex; e != NULL; e = e->next)
|
|
for (k = 0; k <= 5; k++, n++)
|
|
if (e->d[k]) wx[ny+n] = e->d[k];
|
|
/* load digits of [y] */
|
|
wy = &work[nx+ny];
|
|
for (n = 0; n < ny; n++) wy[n] = 0;
|
|
for (n = 0, e = ey; e != NULL; e = e->next)
|
|
for (k = 0; k <= 5; k++, n++)
|
|
if (e->d[k]) wy[n] = e->d[k];
|
|
/* compute [x] * [y] */
|
|
bigmul(nx, ny, wx, wy);
|
|
/* construct and normalize result */
|
|
mpz_set_si(z, 0);
|
|
z->val = sx * sy;
|
|
es = NULL;
|
|
k = 6;
|
|
for (n = 0; n < nx+ny; n++)
|
|
{ if (k > 5)
|
|
{ e = gmp_get_atom(sizeof(struct mpz_seg));
|
|
e->d[0] = e->d[1] = e->d[2] = 0;
|
|
e->d[3] = e->d[4] = e->d[5] = 0;
|
|
e->next = NULL;
|
|
if (z->ptr == NULL)
|
|
z->ptr = e;
|
|
else
|
|
es->next = e;
|
|
es = e;
|
|
k = 0;
|
|
}
|
|
es->d[k++] = wx[n];
|
|
}
|
|
normalize(z);
|
|
done: return;
|
|
}
|
|
|
|
void mpz_neg(mpz_t z, mpz_t x)
|
|
{ /* set z to 0 - x */
|
|
mpz_set(z, x);
|
|
z->val = - z->val;
|
|
return;
|
|
}
|
|
|
|
void mpz_abs(mpz_t z, mpz_t x)
|
|
{ /* set z to the absolute value of x */
|
|
mpz_set(z, x);
|
|
if (z->val < 0) z->val = - z->val;
|
|
return;
|
|
}
|
|
|
|
void mpz_div(mpz_t q, mpz_t r, mpz_t x, mpz_t y)
|
|
{ /* divide x by y, forming quotient q and/or remainder r
|
|
if q = NULL then quotient is not stored; if r = NULL then
|
|
remainder is not stored
|
|
the sign of quotient is determined as in algebra while the
|
|
sign of remainder is the same as the sign of dividend:
|
|
+26 : +7 = +3, remainder is +5
|
|
-26 : +7 = -3, remainder is -5
|
|
+26 : -7 = -3, remainder is +5
|
|
-26 : -7 = +3, remainder is -5 */
|
|
struct mpz_seg dumx, dumy, *ex, *ey, *es, *e;
|
|
int sx, sy, k, nx, ny, n;
|
|
unsigned int t;
|
|
unsigned short *work, *wx, *wy;
|
|
/* divide by zero is not allowed */
|
|
if (y->val == 0)
|
|
{ xassert(y->ptr == NULL);
|
|
xfault("mpz_div: divide by zero not allowed\n");
|
|
}
|
|
/* if [x] = 0 then [q] = [r] = 0 */
|
|
if (x->val == 0)
|
|
{ xassert(x->ptr == NULL);
|
|
if (q != NULL) mpz_set_si(q, 0);
|
|
if (r != NULL) mpz_set_si(r, 0);
|
|
goto done;
|
|
}
|
|
/* special case when both [x] and [y] are in short format */
|
|
if (x->ptr == NULL && y->ptr == NULL)
|
|
{ int xval = x->val, yval = y->val;
|
|
xassert(xval != 0x80000000 && yval != 0x80000000);
|
|
if (q != NULL) mpz_set_si(q, xval / yval);
|
|
if (r != NULL) mpz_set_si(r, xval % yval);
|
|
goto done;
|
|
}
|
|
/* convert [x] to long format, if necessary */
|
|
if (x->ptr == NULL)
|
|
{ xassert(x->val != 0x80000000);
|
|
if (x->val >= 0)
|
|
{ sx = +1;
|
|
t = (unsigned int)(+ x->val);
|
|
}
|
|
else
|
|
{ sx = -1;
|
|
t = (unsigned int)(- x->val);
|
|
}
|
|
ex = &dumx;
|
|
ex->d[0] = (unsigned short)t;
|
|
ex->d[1] = (unsigned short)(t >> 16);
|
|
ex->d[2] = ex->d[3] = ex->d[4] = ex->d[5] = 0;
|
|
ex->next = NULL;
|
|
}
|
|
else
|
|
{ sx = x->val;
|
|
xassert(sx == +1 || sx == -1);
|
|
ex = x->ptr;
|
|
}
|
|
/* convert [y] to long format, if necessary */
|
|
if (y->ptr == NULL)
|
|
{ xassert(y->val != 0x80000000);
|
|
if (y->val >= 0)
|
|
{ sy = +1;
|
|
t = (unsigned int)(+ y->val);
|
|
}
|
|
else
|
|
{ sy = -1;
|
|
t = (unsigned int)(- y->val);
|
|
}
|
|
ey = &dumy;
|
|
ey->d[0] = (unsigned short)t;
|
|
ey->d[1] = (unsigned short)(t >> 16);
|
|
ey->d[2] = ey->d[3] = ey->d[4] = ey->d[5] = 0;
|
|
ey->next = NULL;
|
|
}
|
|
else
|
|
{ sy = y->val;
|
|
xassert(sy == +1 || sy == -1);
|
|
ey = y->ptr;
|
|
}
|
|
/* determine the number of digits of [x] */
|
|
nx = n = 0;
|
|
for (e = ex; e != NULL; e = e->next)
|
|
for (k = 0; k <= 5; k++)
|
|
{ n++;
|
|
if (e->d[k]) nx = n;
|
|
}
|
|
xassert(nx > 0);
|
|
/* determine the number of digits of [y] */
|
|
ny = n = 0;
|
|
for (e = ey; e != NULL; e = e->next)
|
|
for (k = 0; k <= 5; k++)
|
|
{ n++;
|
|
if (e->d[k]) ny = n;
|
|
}
|
|
xassert(ny > 0);
|
|
/* if nx < ny then [q] = 0 and [r] = [x] */
|
|
if (nx < ny)
|
|
{ if (r != NULL) mpz_set(r, x);
|
|
if (q != NULL) mpz_set_si(q, 0);
|
|
goto done;
|
|
}
|
|
/* we need working array containing at least nx+ny+1 places */
|
|
work = gmp_get_work(nx+ny+1);
|
|
/* load digits of [x] */
|
|
wx = &work[0];
|
|
for (n = 0; n < nx; n++) wx[n] = 0;
|
|
for (n = 0, e = ex; e != NULL; e = e->next)
|
|
for (k = 0; k <= 5; k++, n++)
|
|
if (e->d[k]) wx[n] = e->d[k];
|
|
/* load digits of [y] */
|
|
wy = &work[nx+1];
|
|
for (n = 0; n < ny; n++) wy[n] = 0;
|
|
for (n = 0, e = ey; e != NULL; e = e->next)
|
|
for (k = 0; k <= 5; k++, n++)
|
|
if (e->d[k]) wy[n] = e->d[k];
|
|
/* compute quotient and remainder */
|
|
xassert(wy[ny-1] != 0);
|
|
bigdiv(nx-ny, ny, wx, wy);
|
|
/* construct and normalize quotient */
|
|
if (q != NULL)
|
|
{ mpz_set_si(q, 0);
|
|
q->val = sx * sy;
|
|
es = NULL;
|
|
k = 6;
|
|
for (n = ny; n <= nx; n++)
|
|
{ if (k > 5)
|
|
{ e = gmp_get_atom(sizeof(struct mpz_seg));
|
|
e->d[0] = e->d[1] = e->d[2] = 0;
|
|
e->d[3] = e->d[4] = e->d[5] = 0;
|
|
e->next = NULL;
|
|
if (q->ptr == NULL)
|
|
q->ptr = e;
|
|
else
|
|
es->next = e;
|
|
es = e;
|
|
k = 0;
|
|
}
|
|
es->d[k++] = wx[n];
|
|
}
|
|
normalize(q);
|
|
}
|
|
/* construct and normalize remainder */
|
|
if (r != NULL)
|
|
{ mpz_set_si(r, 0);
|
|
r->val = sx;
|
|
es = NULL;
|
|
k = 6;
|
|
for (n = 0; n < ny; n++)
|
|
{ if (k > 5)
|
|
{ e = gmp_get_atom(sizeof(struct mpz_seg));
|
|
e->d[0] = e->d[1] = e->d[2] = 0;
|
|
e->d[3] = e->d[4] = e->d[5] = 0;
|
|
e->next = NULL;
|
|
if (r->ptr == NULL)
|
|
r->ptr = e;
|
|
else
|
|
es->next = e;
|
|
es = e;
|
|
k = 0;
|
|
}
|
|
es->d[k++] = wx[n];
|
|
}
|
|
normalize(r);
|
|
}
|
|
done: return;
|
|
}
|
|
|
|
void mpz_gcd(mpz_t z, mpz_t x, mpz_t y)
|
|
{ /* set z to the greatest common divisor of x and y */
|
|
/* in case of arbitrary integers GCD(x, y) = GCD(|x|, |y|), and,
|
|
in particular, GCD(0, 0) = 0 */
|
|
mpz_t u, v, r;
|
|
mpz_init(u);
|
|
mpz_init(v);
|
|
mpz_init(r);
|
|
mpz_abs(u, x);
|
|
mpz_abs(v, y);
|
|
while (mpz_sgn(v))
|
|
{ mpz_div(NULL, r, u, v);
|
|
mpz_set(u, v);
|
|
mpz_set(v, r);
|
|
}
|
|
mpz_set(z, u);
|
|
mpz_clear(u);
|
|
mpz_clear(v);
|
|
mpz_clear(r);
|
|
return;
|
|
}
|
|
|
|
int mpz_cmp(mpz_t x, mpz_t y)
|
|
{ /* compare x and y; return a positive value if x > y, zero if
|
|
x = y, or a nefative value if x < y */
|
|
static struct mpz_seg zero = { { 0, 0, 0, 0, 0, 0 }, NULL };
|
|
struct mpz_seg dumx, dumy, *ex, *ey;
|
|
int cc, sx, sy, k;
|
|
unsigned int t;
|
|
if (x == y)
|
|
{ cc = 0;
|
|
goto done;
|
|
}
|
|
/* special case when both [x] and [y] are in short format */
|
|
if (x->ptr == NULL && y->ptr == NULL)
|
|
{ int xval = x->val, yval = y->val;
|
|
xassert(xval != 0x80000000 && yval != 0x80000000);
|
|
cc = (xval > yval ? +1 : xval < yval ? -1 : 0);
|
|
goto done;
|
|
}
|
|
/* special case when [x] and [y] have different signs */
|
|
if (x->val > 0 && y->val <= 0 || x->val == 0 && y->val < 0)
|
|
{ cc = +1;
|
|
goto done;
|
|
}
|
|
if (x->val < 0 && y->val >= 0 || x->val == 0 && y->val > 0)
|
|
{ cc = -1;
|
|
goto done;
|
|
}
|
|
/* convert [x] to long format, if necessary */
|
|
if (x->ptr == NULL)
|
|
{ xassert(x->val != 0x80000000);
|
|
if (x->val >= 0)
|
|
{ sx = +1;
|
|
t = (unsigned int)(+ x->val);
|
|
}
|
|
else
|
|
{ sx = -1;
|
|
t = (unsigned int)(- x->val);
|
|
}
|
|
ex = &dumx;
|
|
ex->d[0] = (unsigned short)t;
|
|
ex->d[1] = (unsigned short)(t >> 16);
|
|
ex->d[2] = ex->d[3] = ex->d[4] = ex->d[5] = 0;
|
|
ex->next = NULL;
|
|
}
|
|
else
|
|
{ sx = x->val;
|
|
xassert(sx == +1 || sx == -1);
|
|
ex = x->ptr;
|
|
}
|
|
/* convert [y] to long format, if necessary */
|
|
if (y->ptr == NULL)
|
|
{ xassert(y->val != 0x80000000);
|
|
if (y->val >= 0)
|
|
{ sy = +1;
|
|
t = (unsigned int)(+ y->val);
|
|
}
|
|
else
|
|
{ sy = -1;
|
|
t = (unsigned int)(- y->val);
|
|
}
|
|
ey = &dumy;
|
|
ey->d[0] = (unsigned short)t;
|
|
ey->d[1] = (unsigned short)(t >> 16);
|
|
ey->d[2] = ey->d[3] = ey->d[4] = ey->d[5] = 0;
|
|
ey->next = NULL;
|
|
}
|
|
else
|
|
{ sy = y->val;
|
|
xassert(sy == +1 || sy == -1);
|
|
ey = y->ptr;
|
|
}
|
|
/* main fragment */
|
|
xassert(sx > 0 && sy > 0 || sx < 0 && sy < 0);
|
|
cc = 0;
|
|
for (; ex || ey; ex = ex->next, ey = ey->next)
|
|
{ if (ex == NULL) ex = &zero;
|
|
if (ey == NULL) ey = &zero;
|
|
for (k = 0; k <= 5; k++)
|
|
{ if (ex->d[k] > ey->d[k]) cc = +1;
|
|
if (ex->d[k] < ey->d[k]) cc = -1;
|
|
}
|
|
}
|
|
if (sx < 0) cc = - cc;
|
|
done: return cc;
|
|
}
|
|
|
|
int mpz_sgn(mpz_t x)
|
|
{ /* return +1 if x > 0, 0 if x = 0, and -1 if x < 0 */
|
|
int s;
|
|
s = (x->val > 0 ? +1 : x->val < 0 ? -1 : 0);
|
|
return s;
|
|
}
|
|
|
|
int mpz_out_str(void *_fp, int base, mpz_t x)
|
|
{ /* output x on stream fp, as a string in given base; the base
|
|
may vary from 2 to 36;
|
|
return the number of bytes written, or if an error occurred,
|
|
return 0 */
|
|
FILE *fp = _fp;
|
|
mpz_t b, y, r;
|
|
int n, j, nwr = 0;
|
|
unsigned char *d;
|
|
static char *set = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ";
|
|
if (!(2 <= base && base <= 36))
|
|
xfault("mpz_out_str: base = %d; invalid base\n", base);
|
|
mpz_init(b);
|
|
mpz_set_si(b, base);
|
|
mpz_init(y);
|
|
mpz_init(r);
|
|
/* determine the number of digits */
|
|
mpz_abs(y, x);
|
|
for (n = 0; mpz_sgn(y) != 0; n++)
|
|
mpz_div(y, NULL, y, b);
|
|
if (n == 0) n = 1;
|
|
/* compute the digits */
|
|
d = xmalloc(n);
|
|
mpz_abs(y, x);
|
|
for (j = 0; j < n; j++)
|
|
{ mpz_div(y, r, y, b);
|
|
xassert(0 <= r->val && r->val < base && r->ptr == NULL);
|
|
d[j] = (unsigned char)r->val;
|
|
}
|
|
/* output the integer to the stream */
|
|
if (fp == NULL) fp = stdout;
|
|
if (mpz_sgn(x) < 0)
|
|
fputc('-', fp), nwr++;
|
|
for (j = n-1; j >= 0; j--)
|
|
fputc(set[d[j]], fp), nwr++;
|
|
if (ferror(fp)) nwr = 0;
|
|
mpz_clear(b);
|
|
mpz_clear(y);
|
|
mpz_clear(r);
|
|
xfree(d);
|
|
return nwr;
|
|
}
|
|
|
|
/*====================================================================*/
|
|
|
|
mpq_t _mpq_init(void)
|
|
{ /* initialize x, and set its value to 0/1 */
|
|
mpq_t x;
|
|
x = gmp_get_atom(sizeof(struct mpq));
|
|
x->p.val = 0;
|
|
x->p.ptr = NULL;
|
|
x->q.val = 1;
|
|
x->q.ptr = NULL;
|
|
return x;
|
|
}
|
|
|
|
void mpq_clear(mpq_t x)
|
|
{ /* free the space occupied by x */
|
|
mpz_set_si(&x->p, 0);
|
|
xassert(x->p.ptr == NULL);
|
|
mpz_set_si(&x->q, 0);
|
|
xassert(x->q.ptr == NULL);
|
|
/* free the number descriptor */
|
|
gmp_free_atom(x, sizeof(struct mpq));
|
|
return;
|
|
}
|
|
|
|
void mpq_canonicalize(mpq_t x)
|
|
{ /* remove any factors that are common to the numerator and
|
|
denominator of x, and make the denominator positive */
|
|
mpz_t f;
|
|
xassert(x->q.val != 0);
|
|
if (x->q.val < 0)
|
|
{ mpz_neg(&x->p, &x->p);
|
|
mpz_neg(&x->q, &x->q);
|
|
}
|
|
mpz_init(f);
|
|
mpz_gcd(f, &x->p, &x->q);
|
|
if (!(f->val == 1 && f->ptr == NULL))
|
|
{ mpz_div(&x->p, NULL, &x->p, f);
|
|
mpz_div(&x->q, NULL, &x->q, f);
|
|
}
|
|
mpz_clear(f);
|
|
return;
|
|
}
|
|
|
|
void mpq_set(mpq_t z, mpq_t x)
|
|
{ /* set the value of z from x */
|
|
if (z != x)
|
|
{ mpz_set(&z->p, &x->p);
|
|
mpz_set(&z->q, &x->q);
|
|
}
|
|
return;
|
|
}
|
|
|
|
void mpq_set_si(mpq_t x, int p, unsigned int q)
|
|
{ /* set the value of x to p/q */
|
|
if (q == 0)
|
|
xfault("mpq_set_si: zero denominator not allowed\n");
|
|
mpz_set_si(&x->p, p);
|
|
xassert(q <= 0x7FFFFFFF);
|
|
mpz_set_si(&x->q, q);
|
|
return;
|
|
}
|
|
|
|
double mpq_get_d(mpq_t x)
|
|
{ /* convert x to a double, truncating if necessary */
|
|
int np, nq;
|
|
double p, q;
|
|
p = mpz_get_d_2exp(&np, &x->p);
|
|
q = mpz_get_d_2exp(&nq, &x->q);
|
|
return ldexp(p / q, np - nq);
|
|
}
|
|
|
|
void mpq_set_d(mpq_t x, double val)
|
|
{ /* set x to val; there is no rounding, the conversion is exact */
|
|
int s, n, d, j;
|
|
double f;
|
|
mpz_t temp;
|
|
xassert(-DBL_MAX <= val && val <= +DBL_MAX);
|
|
mpq_set_si(x, 0, 1);
|
|
if (val > 0.0)
|
|
s = +1;
|
|
else if (val < 0.0)
|
|
s = -1;
|
|
else
|
|
goto done;
|
|
f = frexp(fabs(val), &n);
|
|
/* |val| = f * 2^n, where 0.5 <= f < 1.0 */
|
|
mpz_init(temp);
|
|
while (f != 0.0)
|
|
{ f *= 16.0, n -= 4;
|
|
d = (int)f;
|
|
xassert(0 <= d && d <= 15);
|
|
f -= (double)d;
|
|
/* x := 16 * x + d */
|
|
mpz_set_si(temp, 16);
|
|
mpz_mul(&x->p, &x->p, temp);
|
|
mpz_set_si(temp, d);
|
|
mpz_add(&x->p, &x->p, temp);
|
|
}
|
|
mpz_clear(temp);
|
|
/* x := x * 2^n */
|
|
if (n > 0)
|
|
{ for (j = 1; j <= n; j++)
|
|
mpz_add(&x->p, &x->p, &x->p);
|
|
}
|
|
else if (n < 0)
|
|
{ for (j = 1; j <= -n; j++)
|
|
mpz_add(&x->q, &x->q, &x->q);
|
|
mpq_canonicalize(x);
|
|
}
|
|
if (s < 0) mpq_neg(x, x);
|
|
done: return;
|
|
}
|
|
|
|
void mpq_add(mpq_t z, mpq_t x, mpq_t y)
|
|
{ /* set z to x + y */
|
|
mpz_t p, q;
|
|
mpz_init(p);
|
|
mpz_init(q);
|
|
mpz_mul(p, &x->p, &y->q);
|
|
mpz_mul(q, &x->q, &y->p);
|
|
mpz_add(p, p, q);
|
|
mpz_mul(q, &x->q, &y->q);
|
|
mpz_set(&z->p, p);
|
|
mpz_set(&z->q, q);
|
|
mpz_clear(p);
|
|
mpz_clear(q);
|
|
mpq_canonicalize(z);
|
|
return;
|
|
}
|
|
|
|
void mpq_sub(mpq_t z, mpq_t x, mpq_t y)
|
|
{ /* set z to x - y */
|
|
mpz_t p, q;
|
|
mpz_init(p);
|
|
mpz_init(q);
|
|
mpz_mul(p, &x->p, &y->q);
|
|
mpz_mul(q, &x->q, &y->p);
|
|
mpz_sub(p, p, q);
|
|
mpz_mul(q, &x->q, &y->q);
|
|
mpz_set(&z->p, p);
|
|
mpz_set(&z->q, q);
|
|
mpz_clear(p);
|
|
mpz_clear(q);
|
|
mpq_canonicalize(z);
|
|
return;
|
|
}
|
|
|
|
void mpq_mul(mpq_t z, mpq_t x, mpq_t y)
|
|
{ /* set z to x * y */
|
|
mpz_mul(&z->p, &x->p, &y->p);
|
|
mpz_mul(&z->q, &x->q, &y->q);
|
|
mpq_canonicalize(z);
|
|
return;
|
|
}
|
|
|
|
void mpq_div(mpq_t z, mpq_t x, mpq_t y)
|
|
{ /* set z to x / y */
|
|
mpz_t p, q;
|
|
if (mpq_sgn(y) == 0)
|
|
xfault("mpq_div: zero divisor not allowed\n");
|
|
mpz_init(p);
|
|
mpz_init(q);
|
|
mpz_mul(p, &x->p, &y->q);
|
|
mpz_mul(q, &x->q, &y->p);
|
|
mpz_set(&z->p, p);
|
|
mpz_set(&z->q, q);
|
|
mpz_clear(p);
|
|
mpz_clear(q);
|
|
mpq_canonicalize(z);
|
|
return;
|
|
}
|
|
|
|
void mpq_neg(mpq_t z, mpq_t x)
|
|
{ /* set z to 0 - x */
|
|
mpq_set(z, x);
|
|
mpz_neg(&z->p, &z->p);
|
|
return;
|
|
}
|
|
|
|
void mpq_abs(mpq_t z, mpq_t x)
|
|
{ /* set z to the absolute value of x */
|
|
mpq_set(z, x);
|
|
mpz_abs(&z->p, &z->p);
|
|
xassert(mpz_sgn(&x->q) > 0);
|
|
return;
|
|
}
|
|
|
|
int mpq_cmp(mpq_t x, mpq_t y)
|
|
{ /* compare x and y; return a positive value if x > y, zero if
|
|
x = y, or a nefative value if x < y */
|
|
mpq_t temp;
|
|
int s;
|
|
mpq_init(temp);
|
|
mpq_sub(temp, x, y);
|
|
s = mpq_sgn(temp);
|
|
mpq_clear(temp);
|
|
return s;
|
|
}
|
|
|
|
int mpq_sgn(mpq_t x)
|
|
{ /* return +1 if x > 0, 0 if x = 0, and -1 if x < 0 */
|
|
int s;
|
|
s = mpz_sgn(&x->p);
|
|
xassert(mpz_sgn(&x->q) > 0);
|
|
return s;
|
|
}
|
|
|
|
int mpq_out_str(void *_fp, int base, mpq_t x)
|
|
{ /* output x on stream fp, as a string in given base; the base
|
|
may vary from 2 to 36; output is in the form 'num/den' or if
|
|
the denominator is 1 then just 'num';
|
|
if the parameter fp is a null pointer, stdout is assumed;
|
|
return the number of bytes written, or if an error occurred,
|
|
return 0 */
|
|
FILE *fp = _fp;
|
|
int nwr;
|
|
if (!(2 <= base && base <= 36))
|
|
xfault("mpq_out_str: base = %d; invalid base\n", base);
|
|
if (fp == NULL) fp = stdout;
|
|
nwr = mpz_out_str(fp, base, &x->p);
|
|
if (x->q.val == 1 && x->q.ptr == NULL)
|
|
;
|
|
else
|
|
{ fputc('/', fp), nwr++;
|
|
nwr += mpz_out_str(fp, base, &x->q);
|
|
}
|
|
if (ferror(fp)) nwr = 0;
|
|
return nwr;
|
|
}
|
|
|
|
#endif
|
|
|
|
/* eof */
|