/* lrsmp.c library code for lrs extended precision arithmetic */
/* Version 4.0b, Feb. 8, 2000 */
/* Copyright: David Avis 1999, avis@cs.mcgill.ca */
/* This program 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 2 of the License, or
(at your option) any later version.
This program 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 this program; if not, write to the Free Software
Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
*/
/******************************************************************/
/* digit overflow is caught by digits_overflow at the end of this */
/* file, make sure it is either user supplied or uncomment */
/* the define below */
/******************************************************************/
#define digits_overflow() lrs_default_digits_overflow()
#include <stdio.h>
#include <string.h>
#include "lrsmp.h"
/*********************************************************/
/* Initialization and allocation procedures - must use! */
/******************************************************* */
long
lrs_mp_init (long dec_digits, FILE * fpin, FILE * fpout)
/* max number of decimal digits for the computation */
{
/* global variables lrs_ifp and lrs_ofp are file pointers for input and output */
lrs_ifp = fpin;
lrs_ofp = fpout;
lrs_record_digits = 0;
if (dec_digits <= 0)
dec_digits = DEFAULT_DIGITS;
lrs_digits = DEC2DIG (dec_digits); /* max permitted no. of digits */
if (lrs_digits > MAX_DIGITS)
{
fprintf (lrs_ofp, "\nDigits must be at most %ld\nChange MAX_DIGITS and recompile\n",
DIG2DEC (MAX_DIGITS));
lrs_digits = MAX_DIGITS;
return FALSE;
}
return TRUE;
}
lrs_mp_t
lrs_alloc_mp_t ()
/* dynamic allocation of lrs_mp number */
{
lrs_mp_t p;
p=(long *)calloc (lrs_digits+1, sizeof (long));
return p;
}
lrs_mp_vector
lrs_alloc_mp_vector (long n)
/* allocate lrs_mp_vector for n+1 lrs_mp numbers */
{
lrs_mp_vector p;
long i;
p = CALLOC ((n + 1), sizeof (lrs_mp *));
for (i = 0; i <= n; i++)
p[i] = CALLOC (1, sizeof (lrs_mp));
return p;
}
void
lrs_clear_mp_vector (lrs_mp_vector p, long n)
/* free space allocated to p */
{
long i;
for (i=0; i<=n; i++)
free (p[i] );
free (p);
}
lrs_mp_matrix
lrs_alloc_mp_matrix (long m, long n)
/* allocate lrs_mp_matrix for m+1 x n+1 lrs_mp numbers */
{
lrs_mp_matrix a;
long *araw;
int mp_width, row_width;
int i, j;
mp_width = lrs_digits + 1;
row_width = (n + 1) * mp_width;
araw = calloc ((m + 1) * row_width, sizeof (long));
a = calloc ((m + 1), sizeof (lrs_mp_vector));
for (i = 0; i < m + 1; i++)
{
a[i] = calloc ((n + 1), sizeof (lrs_mp *));
for (j = 0; j < n + 1; j++)
a[i][j] = (araw + i * row_width + j * mp_width);
}
return a;
}
void
lrs_clear_mp_matrix (lrs_mp_matrix p, long m, long n)
/* free space allocated to lrs_mp_matrix p */
{
long i;
/* p[0][0] is araw, the actual matrix storage address */
free(p[0][0]);
for (i = 0; i < m + 1; i++)
free (p[i]);
free(p);
}
/*********************************************************/
/* Core library functions - depend on mp implementation */
/******************************************************* */
void
copy (lrs_mp a, lrs_mp b) /* assigns a=b */
{
long i;
for (i = 0; i <= length (b); i++)
a[i] = b[i];
}
/********************************************************/
/* Divide two multiple precision integers (c=a/b). */
/* a is destroyed and contains the remainder on return. */
/* From Knuth Vol.2 SemiNumerical Algorithms */
/* coded by J. Quinn */
/********************************************************/
void
divint (lrs_mp a, lrs_mp b, lrs_mp c) /* c=a/b, a contains remainder on return */
{
long cy, la, lb, lc, d1, s, t, sig;
long i, j, qh;
/* figure out and save sign, do everything with positive numbers */
sig = sign (a) * sign (b);
la = length (a);
lb = length (b);
lc = la - lb + 2;
if (la < lb)
{
storelength (c, TWO);
storesign (c, POS);
c[1] = 0;
normalize (c);
return;
}
for (i = 1; i < lc; i++)
c[i] = 0;
storelength (c, lc);
storesign (c, (sign (a) == sign (b)) ? POS : NEG);
/******************************/
/* division by a single word: */
/* do it directly */
/******************************/
if (lb == 2)
{
cy = 0;
t = b[1];
for (i = la - 1; i > 0; i--)
{
cy = cy * BASE + a[i];
a[i] = 0;
cy -= (c[i] = cy / t) * t;
}
a[1] = cy;
storesign (a, (cy == 0) ? POS : sign (a));
storelength (a, TWO);
/* set sign of c to sig (**mod**) */
storesign (c, sig);
normalize (c);
return;
}
else
{
/* mp's are actually DIGITS+1 in length, so if length of a or b = */
/* DIGITS, there will still be room after normalization. */
/****************************************************/
/* Step D1 - normalize numbers so b > floor(BASE/2) */
d1 = BASE / (b[lb - 1] + 1);
if (d1 > 1)
{
cy = 0;
for (i = 1; i < la; i++)
{
cy = (a[i] = a[i] * d1 + cy) / BASE;
a[i] %= BASE;
}
a[i] = cy;
cy = 0;
for (i = 1; i < lb; i++)
{
cy = (b[i] = b[i] * d1 + cy) / BASE;
b[i] %= BASE;
}
b[i] = cy;
}
else
{
a[la] = 0; /* if la or lb = DIGITS this won't work */
b[lb] = 0;
}
/*********************************************/
/* Steps D2 & D7 - start and end of the loop */
for (j = 0; j <= la - lb; j++)
{
/*************************************/
/* Step D3 - determine trial divisor */
if (a[la - j] == b[lb - 1])
qh = BASE - 1;
else
{
s = (a[la - j] * BASE + a[la - j - 1]);
qh = s / b[lb - 1];
while (qh * b[lb - 2] > (s - qh * b[lb - 1]) * BASE + a[la - j - 2])
qh--;
}
/*******************************************************/
/* Step D4 - divide through using qh as quotient digit */
cy = 0;
for (i = 1; i <= lb; i++)
{
s = qh * b[i] + cy;
a[la - j - lb + i] -= s % BASE;
cy = s / BASE;
if (a[la - j - lb + i] < 0)
{
a[la - j - lb + i] += BASE;
cy++;
}
}
/*****************************************************/
/* Step D6 - adjust previous step if qh is 1 too big */
if (cy)
{
qh--;
cy = 0;
for (i = 1; i <= lb; i++) /* add a back in */
{
a[la - j - lb + i] += b[i] + cy;
cy = a[la - j - lb + i] / BASE;
a[la - j - lb + i] %= BASE;
}
}
/***********************************************************************/
/* Step D5 - write final value of qh. Saves calculating array indices */
/* to do it here instead of before D6 */
c[la - lb - j + 1] = qh;
}
/**********************************************************************/
/* Step D8 - unnormalize a and b to get correct remainder and divisor */
for (i = lc; c[i - 1] == 0 && i > 2; i--); /* strip excess 0's from quotient */
storelength (c, i);
if (i == 2 && c[1] == 0)
storesign (c, POS);
cy = 0;
for (i = lb - 1; i >= 1; i--)
{
cy = (a[i] += cy * BASE) % d1;
a[i] /= d1;
}
for (i = la; a[i - 1] == 0 && i > 2; i--); /* strip excess 0's from quotient */
storelength (a, i);
if (i == 2 && a[1] == 0)
storesign (a, POS);
if (cy)
fprintf (lrs_ofp, "divide error");
for (i = lb - 1; i >= 1; i--)
{
cy = (b[i] += cy * BASE) % d1;
b[i] /= d1;
}
}
}
/* end of divint */
void
gcd (lrs_mp u, lrs_mp v) /*returns u=gcd(u,v) destroying v */
/*Euclid's algorithm. Knuth, II, p.320
modified to avoid copies r=u,u=v,v=r
Switches to single precision when possible for greater speed */
{
lrs_mp r;
unsigned long ul, vl;
long i;
static unsigned long maxspval = MAXD; /* Max value for the last digit to guarantee */
/* fitting into a single long integer. */
static long maxsplen; /* Maximum digits for a number that will fit */
/* into a single long integer. */
static long firstime = TRUE;
if (firstime) /* initialize constants */
{
for (maxsplen = 2; maxspval >= BASE; maxsplen++)
maxspval /= BASE;
firstime = FALSE;
}
if (greater (v, u))
goto bigv;
bigu:
if (zero (v))
return;
if ((i = length (u)) < maxsplen || (i == maxsplen && u[maxsplen - 1] < maxspval))
goto quickfinish;
divint (u, v, r);
normalize (u);
bigv:
if (zero (u))
{
copy (u, v);
return;
}
if ((i = length (v)) < maxsplen || (i == maxsplen && v[maxsplen - 1] < maxspval))
goto quickfinish;
divint (v, u, r);
normalize (v);
goto bigu;
/* Base 10000 only at the moment */
/* when u and v are small enough, transfer to single precision integers */
/* and finish with euclid's algorithm, then transfer back to lrs_mp */
quickfinish:
ul = vl = 0;
for (i = length (u) - 1; i > 0; i--)
ul = BASE * ul + u[i];
for (i = length (v) - 1; i > 0; i--)
vl = BASE * vl + v[i];
if (ul > vl)
goto qv;
qu:
if (!vl)
{
for (i = 1; ul; i++)
{
u[i] = ul % BASE;
ul = ul / BASE;
}
storelength (u, i);
return;
}
ul %= vl;
qv:
if (!ul)
{
for (i = 1; vl; i++)
{
u[i] = vl % BASE;
vl = vl / BASE;
}
storelength (u, i);
return;
}
vl %= ul;
goto qu;
}
long
compare (lrs_mp a, lrs_mp b) /* a ? b and returns -1,0,1 for <,=,> */
{
long i;
if (a[0] > b[0])
return 1L;
if (a[0] < b[0])
return -1L;
for (i = length (a) - 1; i >= 1; i--)
{
if (a[i] < b[i])
{
if (sign (a) == POS)
return -1L;
else
return 1L;
}
if (a[i] > b[i])
{
if (sign (a) == NEG)
return -1L;
else
return 1L;
}
}
return 0L;
}
long
greater (lrs_mp a, lrs_mp b) /* tests if a > b and returns (TRUE=POS) */
{
long i;
if (a[0] > b[0])
return (TRUE);
if (a[0] < b[0])
return (FALSE);
for (i = length (a) - 1; i >= 1; i--)
{
if (a[i] < b[i])
{
if (sign (a) == POS)
return (0);
else
return (1);
}
if (a[i] > b[i])
{
if (sign (a) == NEG)
return (0);
else
return (1);
}
}
return (0);
}
void
itomp (long in, lrs_mp a)
/* convert integer i to multiple precision with base BASE */
{
long i;
a[0] = 2; /* initialize to zero */
for (i = 1; i < lrs_digits; i++)
a[i] = 0;
if (in < 0)
{
storesign (a, NEG);
in = in * (-1);
}
i = 0;
while (in != 0)
{
i++;
a[i] = in - BASE * (in / BASE);
in = in / BASE;
storelength (a, i + 1);
}
} /* end of itomp */
void
linint (lrs_mp a, long ka, lrs_mp b, long kb) /*compute a*ka+b*kb --> a */
/***Handbook of Algorithms and Data Structures P.239 ***/
{
long i, la, lb;
la = length (a);
lb = length (b);
for (i = 1; i < la; i++)
a[i] *= ka;
if (sign (a) != sign (b))
kb = (-kb);
if (lb > la)
{
storelength (a, lb);
for (i = la; i < lb; i++)
a[i] = 0;
}
for (i = 1; i < lb; i++)
a[i] += kb * b[i];
normalize (a);
}
/***end of linint***/
void
mptodouble (lrs_mp a, double *x) /* convert lrs_mp to double */
{
long i, la;
double y = 1.0;
(*x) = 0;
la = length (a);
for (i = 1; i < la; i++)
{
(*x) = (*x) + y * a[i];
y = y * BASE;
}
if (negative (a))
(*x)= -(*x);
}
void
mulint (lrs_mp a, lrs_mp b, lrs_mp c) /* multiply two integers a*b --> c */
/***Handbook of Algorithms and Data Structures, p239 ***/
{
long nlength, i, j, la, lb;
/*** b and c may coincide ***/
la = length (a);
lb = length (b);
nlength = la + lb - 2;
if (nlength > lrs_digits)
digits_overflow ();
for (i = 0; i < la - 2; i++)
c[lb + i] = 0;
for (i = lb - 1; i > 0; i--)
{
for (j = 2; j < la; j++)
if ((c[i + j - 1] += b[i] * a[j]) >
MAXD - (BASE - 1) * (BASE - 1) - MAXD / BASE)
{
c[i + j - 1] -= (MAXD / BASE) * BASE;
c[i + j] += MAXD / BASE;
}
c[i] = b[i] * a[1];
}
storelength (c, nlength);
storesign (c, sign (a) == sign (b) ? POS : NEG);
normalize (c);
}
/***end of mulint ***/
void
normalize (lrs_mp a)
{
long cy, i, la;
la = length (a);
start:
cy = 0;
for (i = 1; i < la; i++)
{
cy = (a[i] += cy) / BASE;
a[i] -= cy * BASE;
if (a[i] < 0)
{
a[i] += BASE;
cy--;
}
}
while (cy > 0)
{
a[i++] = cy % BASE;
cy /= BASE;
}
if (cy < 0)
{
a[la - 1] += cy * BASE;
for (i = 1; i < la; i++)
a[i] = (-a[i]);
storesign (a, sign (a) == POS ? NEG : POS);
goto start;
}
while (a[i - 1] == 0 && i > 2)
i--;
if (i > lrs_record_digits)
{
if ((lrs_record_digits = i) > lrs_digits)
digits_overflow ();
};
storelength (a, i);
if (i == 2 && a[1] == 0)
storesign (a, POS);
} /* end of normalize */
long
mptoi (lrs_mp a) /* convert lrs_mp to long integer */
{
long len = length (a);
if (len == 2)
return sign (a) * a[1];
if (len == 3)
return sign (a) * (a[1] + BASE * a[2]);
notimpl ("mp to large for conversion to long");
return 0; /* never executed */
}
void
pmp (char name[], lrs_mp a) /*print the long precision integer a */
{
long i;
fprintf (lrs_ofp, "%s", name);
if (sign (a) == NEG)
fprintf (lrs_ofp, "-");
else
fprintf (lrs_ofp, " ");
fprintf (lrs_ofp, "%lu", a[length (a) - 1]);
for (i = length (a) - 2; i >= 1; i--)
fprintf (lrs_ofp, FORMAT, a[i]);
fprintf (lrs_ofp, " ");
}
void
prat (char name[], lrs_mp Nin, lrs_mp Din) /*reduce and print Nin/Din */
{
lrs_mp Nt, Dt;
long i;
fprintf (lrs_ofp, "%s", name);
/* reduce fraction */
copy (Nt, Nin);
copy (Dt, Din);
reduce (Nt, Dt);
/* print out */
if (sign (Nin) * sign (Din) == NEG)
fprintf (lrs_ofp, "-");
else
fprintf (lrs_ofp, " ");
fprintf (lrs_ofp, "%lu", Nt[length (Nt) - 1]);
for (i = length (Nt) - 2; i >= 1; i--)
fprintf (lrs_ofp, FORMAT, Nt[i]);
if (!(Dt[0] == 2 && Dt[1] == 1)) /* rational */
{
fprintf (lrs_ofp, "/");
fprintf (lrs_ofp, "%lu", Dt[length (Dt) - 1]);
for (i = length (Dt) - 2; i >= 1; i--)
fprintf (lrs_ofp, FORMAT, Dt[i]);
}
fprintf (lrs_ofp, " ");
}
void
readmp (lrs_mp a)
/* read an integer and convert to lrs_mp with base BASE */
{
char in[MAXINPUT];
fscanf (lrs_ifp, "%s", in);
atomp (in, a);
}
long
readrat (lrs_mp Na, lrs_mp Da)
/* read a rational or integer and convert to lrs_mp with base BASE */
/* returns true if denominator is not one */
{
char in[MAXINPUT], num[MAXINPUT], den[MAXINPUT];
fscanf (lrs_ifp, "%s", in);
atoaa (in, num, den); /*convert rational to num/dem strings */
atomp (num, Na);
if (den[0] == '\0')
{
itomp (1L, Da);
return (FALSE);
}
atomp (den, Da);
return (TRUE);
}
void
addint (lrs_mp a, lrs_mp b, lrs_mp c) /* compute c=a+b */
{
copy (c, a);
linint (c, 1, b, 1);
}
void
atomp (char s[], lrs_mp a) /*convert string to lrs_mp integer */
{
lrs_mp mpone;
long diff, ten, i, sig;
itomp (1L, mpone);
ten = 10L;
for (i = 0; s[i] == ' ' || s[i] == '\n' || s[i] == '\t'; i++);
/*skip white space */
sig = POS;
if (s[i] == '+' || s[i] == '-') /* sign */
sig = (s[i++] == '+') ? POS : NEG;
itomp (0L, a);
while (s[i] >= '0' && s[i] <= '9')
{
diff = s[i] - '0';
linint (a, ten, mpone, diff);
i++;
}
storesign (a, sig);
if (s[i])
{
fprintf (stderr, "\nIllegal character in number: '%s'\n", s + i);
exit (1);
}
} /* end of atomp */
void
subint (lrs_mp a, lrs_mp b, lrs_mp c) /* compute c=a-b */
{
copy (c, a);
linint (a, 1, b, -1);
}
void
decint (lrs_mp a, lrs_mp b) /* compute a=a-b */
{
linint (a, 1, b, -1);
}
long
myrandom (long num, long nrange)
/* return a random number in range 0..nrange-1 */
{
long i;
i = (num * 401 + 673) % nrange;
return (i);
}
long
atos (char s[]) /* convert s to integer */
{
long i, j;
j = 0;
for (i = 0; s[i] >= '0' && s[i] <= '9'; ++i)
j = 10 * j + s[i] - '0';
return (j);
}
void
stringcpy (char *s, char *t) /*copy t to s pointer version */
{
while (((*s++) = (*t++)) != '\0');
}
void
rattodouble (lrs_mp a, lrs_mp b, double *x) /* convert lrs_mp rational to double */
{
double y;
mptodouble (a, &y);
mptodouble (b, x);
*x = y / (*x);
}
void
atoaa (char in[], char num[], char den[])
/* convert rational string in to num/den strings */
{
long i, j;
for (i = 0; in[i] != '\0' && in[i] != '/'; i++)
num[i] = in[i];
num[i] = '\0';
den[0] = '\0';
if (in[i] == '/')
{
for (j = 0; in[j + i + 1] != '\0'; j++)
den[j] = in[i + j + 1];
den[j] = '\0';
}
} /* end of atoaa */
void
lcm (lrs_mp a, lrs_mp b)
/* a = least common multiple of a, b; b is preserved */
{
lrs_mp u, v;
copy (u, a);
copy (v, b);
gcd (u, v);
exactdivint (a, u, v); /* v=a/u no remainder*/
mulint (v, b, a);
} /* end of lcm */
void
reducearray (lrs_mp_vector p, long n)
/* find largest gcd of p[0]..p[n-1] and divide through */
{
lrs_mp divisor;
lrs_mp Temp;
long i = 0L;
while ((i < n) && zero (p[i]))
i++;
if (i == n)
return;
copy (divisor, p[i]);
storesign (divisor, POS);
i++;
while (i < n)
{
if (!zero (p[i]))
{
copy (Temp, p[i]);
storesign (Temp, POS);
gcd (divisor, Temp);
}
i++;
}
/* reduce by divisor */
for (i = 0; i < n; i++)
if (!zero (p[i]))
reduceint (p[i], divisor);
} /* end of reducearray */
void
reduceint (lrs_mp Na, lrs_mp Da) /* divide Na by Da and return */
{
lrs_mp Temp;
copy (Temp, Na);
exactdivint (Temp, Da, Na);
}
void
reduce (lrs_mp Na, lrs_mp Da) /* reduces Na Da by gcd(Na,Da) */
{
lrs_mp Nb, Db, Nc, Dc;
copy (Nb, Na);
copy (Db, Da);
storesign (Nb, POS);
storesign (Db, POS);
copy (Nc, Na);
copy (Dc, Da);
gcd (Nb, Db); /* Nb is the gcd(Na,Da) */
exactdivint (Nc, Nb, Na);
exactdivint (Dc, Nb, Da);
}
long
comprod (lrs_mp Na, lrs_mp Nb, lrs_mp Nc, lrs_mp Nd) /* +1 if Na*Nb > Nc*Nd */
/* -1 if Na*Nb < Nc*Nd */
/* 0 if Na*Nb = Nc*Nd */
{
lrs_mp mc, md;
mulint (Na, Nb, mc);
mulint (Nc, Nd, md);
linint (mc, ONE, md, -ONE);
if (positive (mc))
return (1);
if (negative (mc))
return (-1);
return (0);
}
void
notimpl (char s[])
{
fflush (stdout);
fprintf (stderr, "\nAbnormal Termination %s\n", s);
exit (1);
}
void
getfactorial (lrs_mp factorial, long k) /* compute k factorial in lrs_mp */
{
lrs_mp temp;
long i;
itomp (ONE, factorial);
for (i = 2; i <= k; i++)
{
itomp (i, temp);
mulint (temp, factorial, factorial);
}
} /* end of getfactorial */
/***************************************************************/
/* Package of routines for rational arithmetic */
/***************************************************************/
void
scalerat (lrs_mp Na, lrs_mp Da, long ka) /* scales rational by ka */
{
lrs_mp Nt;
copy (Nt, Na);
itomp (ZERO, Na);
linint (Na, ZERO, Nt, ka);
reduce (Na, Da);
}
void
linrat (lrs_mp Na, lrs_mp Da, long ka, lrs_mp Nb, lrs_mp Db, long kb, lrs_mp Nc, lrs_mp Dc)
/* computes Nc/Dc = ka*Na/Da +kb* Nb/Db
and reduces answer by gcd(Nc,Dc) */
{
lrs_mp c;
mulint (Na, Db, Nc);
mulint (Da, Nb, c);
linint (Nc, ka, c, kb); /* Nc = (ka*Na*Db)+(kb*Da*Nb) */
mulint (Da, Db, Dc); /* Dc = Da*Db */
reduce (Nc, Dc);
}
void
divrat (lrs_mp Na, lrs_mp Da, lrs_mp Nb, lrs_mp Db, lrs_mp Nc, lrs_mp Dc)
/* computes Nc/Dc = (Na/Da) / ( Nb/Db )
and reduces answer by gcd(Nc,Dc) */
{
mulint (Na, Db, Nc);
mulint (Da, Nb, Dc);
reduce (Nc, Dc);
}
void
mulrat (lrs_mp Na, lrs_mp Da, lrs_mp Nb, lrs_mp Db, lrs_mp Nc, lrs_mp Dc)
/* computes Nc/Dc = Na/Da * Nb/Db and reduces by gcd(Nc,Dc) */
{
mulint (Na, Nb, Nc);
mulint (Da, Db, Dc);
reduce (Nc, Dc);
}
/* End package of routines for rational arithmetic */
/***************************************************************/
/* */
/* End of package for multiple precision arithmetic */
/* */
/***************************************************************/
void *
xcalloc (long n, long s, long l, char *f)
{
void *tmp;
tmp = calloc (n, s);
if (tmp == 0)
{
char buf[200];
sprintf (buf, "\n\nFatal error on line %ld of %s", l, f);
perror (buf);
exit (1);
}
return tmp;
}
void
lrs_getdigits (long *a, long *b)
{
/* send digit information to user */
*a = DIG2DEC (lrs_digits);
*b = DIG2DEC (lrs_record_digits);
return;
}
void
lrs_default_digits_overflow ()
{
fprintf (lrs_ofp, "\nOverflow at digits=%ld", DIG2DEC (lrs_digits));
fprintf (lrs_ofp, "\nInitialize lrs_mp_init with n > %ldL\n", DIG2DEC (lrs_digits));
exit (1);
}
/* end of lrsmp.c */
syntax highlighted by Code2HTML, v. 0.9.1