--- /dev/null
+/*
+Copyright (C) 1988 Free Software Foundation
+ written by Doug Lea (dl@rocky.oswego.edu)
+
+This file is part of the GNU C++ Library. This library is free
+software; you can redistribute it and/or modify it under the terms of
+the GNU Library General Public License as published by the Free
+Software Foundation; either version 2 of the License, or (at your
+option) any later version. This library 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 Library General Public License for more details.
+You should have received a copy of the GNU Library General Public
+License along with this library; if not, write to the Free Software
+Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
+*/
+
+/*
+ Some of the following algorithms are very loosely based on those from
+ MIT C-Scheme bignum.c, which is
+ Copyright (c) 1987 Massachusetts Institute of Technology
+
+ with other guidance from Knuth, vol. 2
+
+ Thanks to the creators of the algorithms.
+*/
+
+#ifdef __GNUG__
+#pragma implementation
+#endif
+#include <Integer.h>
+//#include <std.h>
+#include <ctype.h>
+#include <math.h>
+#include <float.h>
+#include <limits.h>
+//#include <Obstack.h>
+//#include <AllocRing.h>
+//#include <new.h>
+#include <builtin.h>
+
+#ifndef HUGE_VAL
+#ifdef HUGE
+#define HUGE_VAL HUGE
+#else
+#define HUGE_VAL DBL_MAX
+#endif
+#endif
+
+/*
+ Sizes of shifts for multiple-precision arithmetic.
+ These should not be changed unless Integer representation
+ as unsigned shorts is changed in the implementation files.
+*/
+
+#define I_SHIFT (sizeof(short) * CHAR_BIT)
+#define I_RADIX ((unsigned long)(1L << I_SHIFT))
+#define I_MAXNUM ((unsigned long)((I_RADIX - 1)))
+#define I_MINNUM ((unsigned long)(I_RADIX >> 1))
+#define I_POSITIVE 1
+#define I_NEGATIVE 0
+
+/* All routines assume SHORT_PER_LONG > 1 */
+#define SHORT_PER_LONG ((unsigned)(((sizeof(long) + sizeof(short) - 1) / sizeof(short))))
+#define CHAR_PER_LONG ((unsigned)sizeof(long))
+
+/*
+ minimum and maximum sizes for an IntRep
+*/
+
+#define MINIntRep_SIZE 16
+#define MAXIntRep_SIZE I_MAXNUM
+
+#ifndef MALLOC_MIN_OVERHEAD
+#define MALLOC_MIN_OVERHEAD 4
+#endif
+
+IntRep _ZeroRep = {1, 0, 1, {0}};
+IntRep _OneRep = {1, 0, 1, {1}};
+IntRep _MinusOneRep = {1, 0, 0, {1}};
+
+
+// utilities to extract and transfer bits
+
+// get low bits
+
+inline static unsigned short extract(unsigned long x)
+{
+ return x & I_MAXNUM;
+}
+
+// transfer high bits to low
+
+inline static unsigned long down(unsigned long x)
+{
+ return (x >> I_SHIFT) & I_MAXNUM;
+}
+
+// transfer low bits to high
+
+inline static unsigned long up(unsigned long x)
+{
+ return x << I_SHIFT;
+}
+
+// compare two equal-length reps
+
+inline static int docmp(const unsigned short* x, const unsigned short* y, int l)
+{
+ int diff = 0;
+ const unsigned short* xs = &(x[l]);
+ const unsigned short* ys = &(y[l]);
+ while (l-- > 0 && (diff = (*--xs) - (*--ys)) == 0);
+ return diff;
+}
+
+// figure out max length of result of +, -, etc.
+
+inline static int calc_len(int len1, int len2, int pad)
+{
+ return (len1 >= len2)? len1 + pad : len2 + pad;
+}
+
+// ensure len & sgn are correct
+
+inline static void Icheck(IntRep* rep)
+{
+ int l = rep->len;
+ const unsigned short* p = &(rep->s[l]);
+ while (l > 0 && *--p == 0) --l;
+ if ((rep->len = l) == 0) rep->sgn = I_POSITIVE;
+}
+
+
+// zero out the end of a rep
+
+inline static void Iclear_from(IntRep* rep, int p)
+{
+ unsigned short* cp = &(rep->s[p]);
+ const unsigned short* cf = &(rep->s[rep->len]);
+ while(cp < cf) *cp++ = 0;
+}
+
+// copy parts of a rep
+
+static inline void scpy(const unsigned short* src, unsigned short* dest,int nb)
+{
+ while (--nb >= 0) *dest++ = *src++;
+}
+
+// make sure an argument is valid
+
+static inline void nonnil(const IntRep* rep)
+{
+ if (rep == 0)
+ (*lib_error_handler)("Integer", "operation on uninitialized Integer");
+}
+
+// allocate a new Irep. Pad to something close to a power of two.
+
+inline static IntRep* Inew(int newlen)
+{
+ unsigned int siz = sizeof(IntRep) + newlen * sizeof(short) +
+ MALLOC_MIN_OVERHEAD;
+ unsigned int allocsiz = MINIntRep_SIZE;
+ while (allocsiz < siz) allocsiz <<= 1; // find a power of 2
+ allocsiz -= MALLOC_MIN_OVERHEAD;
+ if (allocsiz >= MAXIntRep_SIZE * sizeof(short))
+ (*lib_error_handler)("Integer", "Requested length out of range");
+
+ IntRep* rep = (IntRep *) new char[allocsiz];
+ rep->sz = (allocsiz - sizeof(IntRep) + sizeof(short)) / sizeof(short);
+ return rep;
+}
+
+// allocate: use the bits in src if non-null, clear the rest
+
+IntRep* Ialloc(IntRep* old, const unsigned short* src, int srclen, int newsgn,
+ int newlen)
+{
+ IntRep* rep;
+ if (old == 0 || newlen > old->sz)
+ rep = Inew(newlen);
+ else
+ rep = old;
+
+ rep->len = newlen;
+ rep->sgn = newsgn;
+
+ scpy(src, rep->s, srclen);
+ Iclear_from(rep, srclen);
+
+ if (old != rep && old != 0 && !STATIC_IntRep(old)) delete old;
+ return rep;
+}
+
+// allocate and clear
+
+IntRep* Icalloc(IntRep* old, int newlen)
+{
+ IntRep* rep;
+ if (old == 0 || newlen > old->sz)
+ {
+ if (old != 0 && !STATIC_IntRep(old)) delete old;
+ rep = Inew(newlen);
+ }
+ else
+ rep = old;
+
+ rep->len = newlen;
+ rep->sgn = I_POSITIVE;
+ Iclear_from(rep, 0);
+
+ return rep;
+}
+
+// reallocate
+
+IntRep* Iresize(IntRep* old, int newlen)
+{
+ IntRep* rep;
+ unsigned short oldlen;
+ if (old == 0)
+ {
+ oldlen = 0;
+ rep = Inew(newlen);
+ rep->sgn = I_POSITIVE;
+ }
+ else
+ {
+ oldlen = old->len;
+ if (newlen > old->sz)
+ {
+ rep = Inew(newlen);
+ scpy(old->s, rep->s, oldlen);
+ rep->sgn = old->sgn;
+ if (!STATIC_IntRep(old)) delete old;
+ }
+ else
+ rep = old;
+ }
+
+ rep->len = newlen;
+ Iclear_from(rep, oldlen);
+
+ return rep;
+}
+
+
+// same, for straight copy
+
+IntRep* Icopy(IntRep* old, const IntRep* src)
+{
+ if (old == src) return old;
+ IntRep* rep;
+ if (src == 0)
+ {
+ if (old == 0)
+ rep = Inew(0);
+ else
+ {
+ rep = old;
+ Iclear_from(rep, 0);
+ }
+ rep->len = 0;
+ rep->sgn = I_POSITIVE;
+ }
+ else
+ {
+ int newlen = src->len;
+ if (old == 0 || newlen > old->sz)
+ {
+ if (old != 0 && !STATIC_IntRep(old)) delete old;
+ rep = Inew(newlen);
+ }
+ else
+ rep = old;
+
+ rep->len = newlen;
+ rep->sgn = src->sgn;
+
+ scpy(src->s, rep->s, newlen);
+ }
+
+ return rep;
+}
+
+// allocate & copy space for a long
+
+IntRep* Icopy_long(IntRep* old, long x)
+{
+ int newsgn = (x >= 0);
+ IntRep* rep = Icopy_ulong(old, newsgn ? x : -x);
+ rep->sgn = newsgn;
+ return rep;
+}
+
+IntRep* Icopy_ulong(IntRep* old, unsigned long x)
+{
+ unsigned short src[SHORT_PER_LONG];
+
+ unsigned short srclen = 0;
+ while (x != 0)
+ {
+ src[srclen++] = extract(x);
+ x = down(x);
+ }
+
+ IntRep* rep;
+ if (old == 0 || srclen > old->sz)
+ {
+ if (old != 0 && !STATIC_IntRep(old)) delete old;
+ rep = Inew(srclen);
+ }
+ else
+ rep = old;
+
+ rep->len = srclen;
+ rep->sgn = I_POSITIVE;
+
+ scpy(src, rep->s, srclen);
+
+ return rep;
+}
+
+// special case for zero -- it's worth it!
+
+IntRep* Icopy_zero(IntRep* old)
+{
+ if (old == 0)
+ return &_ZeroRep;
+
+ old->len = 0;
+ old->sgn = I_POSITIVE;
+
+ return old;
+}
+
+// special case for 1 or -1
+
+IntRep* Icopy_one(IntRep* old, int newsgn)
+{
+ if (old == 0 || 1 > old->sz)
+ {
+ if (old != 0 && !STATIC_IntRep(old)) delete old;
+ return newsgn==I_NEGATIVE ? &_MinusOneRep : &_OneRep;
+ }
+
+ old->sgn = newsgn;
+ old->len = 1;
+ old->s[0] = 1;
+
+ return old;
+}
+
+// convert to a legal two's complement long if possible
+// if too big, return most negative/positive value
+
+long Itolong(const IntRep* rep)
+{
+ if ((unsigned)(rep->len) > (unsigned)(SHORT_PER_LONG))
+ return (rep->sgn == I_POSITIVE) ? LONG_MAX : LONG_MIN;
+ else if (rep->len == 0)
+ return 0;
+ else if ((unsigned)(rep->len) < (unsigned)(SHORT_PER_LONG))
+ {
+ unsigned long a = rep->s[rep->len-1];
+ if (SHORT_PER_LONG > 2) // normally optimized out
+ {
+ for (int i = rep->len - 2; i >= 0; --i)
+ a = up(a) | rep->s[i];
+ }
+ return (rep->sgn == I_POSITIVE)? a : -((long)a);
+ }
+ else
+ {
+ unsigned long a = rep->s[SHORT_PER_LONG - 1];
+ if (a >= I_MINNUM)
+ return (rep->sgn == I_POSITIVE) ? LONG_MAX : LONG_MIN;
+ else
+ {
+ a = up(a) | rep->s[SHORT_PER_LONG - 2];
+ if (SHORT_PER_LONG > 2)
+ {
+ for (int i = SHORT_PER_LONG - 3; i >= 0; --i)
+ a = up(a) | rep->s[i];
+ }
+ return (rep->sgn == I_POSITIVE)? a : -((long)a);
+ }
+ }
+}
+
+// test whether op long() will work.
+// careful about asymmetry between LONG_MIN & LONG_MAX
+
+int Iislong(const IntRep* rep)
+{
+ unsigned int l = rep->len;
+ if (l < SHORT_PER_LONG)
+ return 1;
+ else if (l > SHORT_PER_LONG)
+ return 0;
+ else if ((unsigned)(rep->s[SHORT_PER_LONG - 1]) < (unsigned)(I_MINNUM))
+ return 1;
+ else if (rep->sgn == I_NEGATIVE && rep->s[SHORT_PER_LONG - 1] == I_MINNUM)
+ {
+ for (unsigned int i = 0; i < SHORT_PER_LONG - 1; ++i)
+ if (rep->s[i] != 0)
+ return 0;
+ return 1;
+ }
+ else
+ return 0;
+}
+
+// convert to a double
+
+double Itodouble(const IntRep* rep)
+{
+ double d = 0.0;
+ double bound = DBL_MAX / 2.0;
+ for (int i = rep->len - 1; i >= 0; --i)
+ {
+ unsigned short a = I_RADIX >> 1;
+ while (a != 0)
+ {
+ if (d >= bound)
+ return (rep->sgn == I_NEGATIVE) ? -HUGE_VAL : HUGE_VAL;
+ d *= 2.0;
+ if (rep->s[i] & a)
+ d += 1.0;
+ a >>= 1;
+ }
+ }
+ if (rep->sgn == I_NEGATIVE)
+ return -d;
+ else
+ return d;
+}
+
+// see whether op double() will work-
+// have to actually try it in order to find out
+// since otherwise might trigger fp exception
+
+int Iisdouble(const IntRep* rep)
+{
+ double d = 0.0;
+ double bound = DBL_MAX / 2.0;
+ for (int i = rep->len - 1; i >= 0; --i)
+ {
+ unsigned short a = I_RADIX >> 1;
+ while (a != 0)
+ {
+ if (d > bound || (d == bound && (i > 0 || (rep->s[i] & a))))
+ return 0;
+ d *= 2.0;
+ if (rep->s[i] & a)
+ d += 1.0;
+ a >>= 1;
+ }
+ }
+ return 1;
+}
+
+// real division of num / den
+
+double ratio(const Integer& num, const Integer& den)
+{
+ Integer q, r;
+ divide(num, den, q, r);
+ double d1 = double(q);
+
+ if (d1 >= DBL_MAX || d1 <= -DBL_MAX || sign(r) == 0)
+ return d1;
+ else // use as much precision as available for fractional part
+ {
+ double d2 = 0.0;
+ double d3 = 0.0;
+ int cont = 1;
+ for (int i = den.rep->len - 1; i >= 0 && cont; --i)
+ {
+ unsigned short a = I_RADIX >> 1;
+ while (a != 0)
+ {
+ if (d2 + 1.0 == d2) // out of precision when we get here
+ {
+ cont = 0;
+ break;
+ }
+
+ d2 *= 2.0;
+ if (den.rep->s[i] & a)
+ d2 += 1.0;
+
+ if (i < r.rep->len)
+ {
+ d3 *= 2.0;
+ if (r.rep->s[i] & a)
+ d3 += 1.0;
+ }
+
+ a >>= 1;
+ }
+ }
+
+ if (sign(r) < 0)
+ d3 = -d3;
+ return d1 + d3 / d2;
+ }
+}
+
+// comparison functions
+
+int compare(const IntRep* x, const IntRep* y)
+{
+ int diff = x->sgn - y->sgn;
+ if (diff == 0)
+ {
+ diff = x->len - y->len;
+ if (diff == 0)
+ diff = docmp(x->s, y->s, x->len);
+ if (x->sgn == I_NEGATIVE)
+ diff = -diff;
+ }
+ return diff;
+}
+
+int ucompare(const IntRep* x, const IntRep* y)
+{
+ int diff = x->len - y->len;
+ if (diff == 0)
+ {
+ int l = x->len;
+ const unsigned short* xs = &(x->s[l]);
+ const unsigned short* ys = &(y->s[l]);
+ while (l-- > 0 && (diff = (*--xs) - (*--ys)) == 0);
+ }
+ return diff;
+}
+
+int compare(const IntRep* x, long y)
+{
+ int xl = x->len;
+ int xsgn = x->sgn;
+ if (y == 0)
+ {
+ if (xl == 0)
+ return 0;
+ else if (xsgn == I_NEGATIVE)
+ return -1;
+ else
+ return 1;
+ }
+ else
+ {
+ int ysgn = y >= 0;
+ unsigned long uy = (ysgn)? y : -y;
+ int diff = xsgn - ysgn;
+ if (diff == 0)
+ {
+ diff = xl - SHORT_PER_LONG;
+ if (diff <= 0)
+ {
+ unsigned short tmp[SHORT_PER_LONG];
+ int yl = 0;
+ while (uy != 0)
+ {
+ tmp[yl++] = extract(uy);
+ uy = down(uy);
+ }
+ diff = xl - yl;
+ if (diff == 0)
+ diff = docmp(x->s, tmp, xl);
+ }
+ if (xsgn == I_NEGATIVE)
+ diff = -diff;
+ }
+ return diff;
+ }
+}
+
+int ucompare(const IntRep* x, long y)
+{
+ int xl = x->len;
+ if (y == 0)
+ return xl;
+ else
+ {
+ unsigned long uy = (y >= 0)? y : -y;
+ int diff = xl - SHORT_PER_LONG;
+ if (diff <= 0)
+ {
+ unsigned short tmp[SHORT_PER_LONG];
+ int yl = 0;
+ while (uy != 0)
+ {
+ tmp[yl++] = extract(uy);
+ uy = down(uy);
+ }
+ diff = xl - yl;
+ if (diff == 0)
+ diff = docmp(x->s, tmp, xl);
+ }
+ return diff;
+ }
+}
+
+
+
+// arithmetic functions
+
+IntRep* add(const IntRep* x, int negatex,
+ const IntRep* y, int negatey, IntRep* r)
+{
+ nonnil(x);
+ nonnil(y);
+
+ int xl = x->len;
+ int yl = y->len;
+
+ int xsgn = (negatex && xl != 0) ? !x->sgn : x->sgn;
+ int ysgn = (negatey && yl != 0) ? !y->sgn : y->sgn;
+
+ int xrsame = x == r;
+ int yrsame = y == r;
+
+ if (yl == 0)
+ r = Ialloc(r, x->s, xl, xsgn, xl);
+ else if (xl == 0)
+ r = Ialloc(r, y->s, yl, ysgn, yl);
+ else if (xsgn == ysgn)
+ {
+ if (xrsame || yrsame)
+ r = Iresize(r, calc_len(xl, yl, 1));
+ else
+ r = Icalloc(r, calc_len(xl, yl, 1));
+ r->sgn = xsgn;
+ unsigned short* rs = r->s;
+ const unsigned short* as;
+ const unsigned short* bs;
+ const unsigned short* topa;
+ const unsigned short* topb;
+ if (xl >= yl)
+ {
+ as = (xrsame)? r->s : x->s;
+ topa = &(as[xl]);
+ bs = (yrsame)? r->s : y->s;
+ topb = &(bs[yl]);
+ }
+ else
+ {
+ bs = (xrsame)? r->s : x->s;
+ topb = &(bs[xl]);
+ as = (yrsame)? r->s : y->s;
+ topa = &(as[yl]);
+ }
+ unsigned long sum = 0;
+ while (bs < topb)
+ {
+ sum += (unsigned long)(*as++) + (unsigned long)(*bs++);
+ *rs++ = extract(sum);
+ sum = down(sum);
+ }
+ while (sum != 0 && as < topa)
+ {
+ sum += (unsigned long)(*as++);
+ *rs++ = extract(sum);
+ sum = down(sum);
+ }
+ if (sum != 0)
+ *rs = extract(sum);
+ else if (rs != as)
+ while (as < topa)
+ *rs++ = *as++;
+ }
+ else
+ {
+ int comp = ucompare(x, y);
+ if (comp == 0)
+ r = Icopy_zero(r);
+ else
+ {
+ if (xrsame || yrsame)
+ r = Iresize(r, calc_len(xl, yl, 0));
+ else
+ r = Icalloc(r, calc_len(xl, yl, 0));
+ unsigned short* rs = r->s;
+ const unsigned short* as;
+ const unsigned short* bs;
+ const unsigned short* topa;
+ const unsigned short* topb;
+ if (comp > 0)
+ {
+ as = (xrsame)? r->s : x->s;
+ topa = &(as[xl]);
+ bs = (yrsame)? r->s : y->s;
+ topb = &(bs[yl]);
+ r->sgn = xsgn;
+ }
+ else
+ {
+ bs = (xrsame)? r->s : x->s;
+ topb = &(bs[xl]);
+ as = (yrsame)? r->s : y->s;
+ topa = &(as[yl]);
+ r->sgn = ysgn;
+ }
+ unsigned long hi = 1;
+ while (bs < topb)
+ {
+ hi += (unsigned long)(*as++) + I_MAXNUM - (unsigned long)(*bs++);
+ *rs++ = extract(hi);
+ hi = down(hi);
+ }
+ while (hi == 0 && as < topa)
+ {
+ hi = (unsigned long)(*as++) + I_MAXNUM;
+ *rs++ = extract(hi);
+ hi = down(hi);
+ }
+ if (rs != as)
+ while (as < topa)
+ *rs++ = *as++;
+ }
+ }
+ Icheck(r);
+ return r;
+}
+
+
+IntRep* add(const IntRep* x, int negatex, long y, IntRep* r)
+{
+ nonnil(x);
+ int xl = x->len;
+ int xsgn = (negatex && xl != 0) ? !x->sgn : x->sgn;
+ int xrsame = x == r;
+
+ int ysgn = (y >= 0);
+ unsigned long uy = (ysgn)? y : -y;
+
+ if (y == 0)
+ r = Ialloc(r, x->s, xl, xsgn, xl);
+ else if (xl == 0)
+ r = Icopy_long(r, y);
+ else if (xsgn == ysgn)
+ {
+ if (xrsame)
+ r = Iresize(r, calc_len(xl, SHORT_PER_LONG, 1));
+ else
+ r = Icalloc(r, calc_len(xl, SHORT_PER_LONG, 1));
+ r->sgn = xsgn;
+ unsigned short* rs = r->s;
+ const unsigned short* as = (xrsame)? r->s : x->s;
+ const unsigned short* topa = &(as[xl]);
+ unsigned long sum = 0;
+ while (as < topa && uy != 0)
+ {
+ unsigned long u = extract(uy);
+ uy = down(uy);
+ sum += (unsigned long)(*as++) + u;
+ *rs++ = extract(sum);
+ sum = down(sum);
+ }
+ while (sum != 0 && as < topa)
+ {
+ sum += (unsigned long)(*as++);
+ *rs++ = extract(sum);
+ sum = down(sum);
+ }
+ if (sum != 0)
+ *rs = extract(sum);
+ else if (rs != as)
+ while (as < topa)
+ *rs++ = *as++;
+ }
+ else
+ {
+ unsigned short tmp[SHORT_PER_LONG];
+ int yl = 0;
+ while (uy != 0)
+ {
+ tmp[yl++] = extract(uy);
+ uy = down(uy);
+ }
+ int comp = xl - yl;
+ if (comp == 0)
+ comp = docmp(x->s, tmp, yl);
+ if (comp == 0)
+ r = Icopy_zero(r);
+ else
+ {
+ if (xrsame)
+ r = Iresize(r, calc_len(xl, yl, 0));
+ else
+ r = Icalloc(r, calc_len(xl, yl, 0));
+ unsigned short* rs = r->s;
+ const unsigned short* as;
+ const unsigned short* bs;
+ const unsigned short* topa;
+ const unsigned short* topb;
+ if (comp > 0)
+ {
+ as = (xrsame)? r->s : x->s;
+ topa = &(as[xl]);
+ bs = tmp;
+ topb = &(bs[yl]);
+ r->sgn = xsgn;
+ }
+ else
+ {
+ bs = (xrsame)? r->s : x->s;
+ topb = &(bs[xl]);
+ as = tmp;
+ topa = &(as[yl]);
+ r->sgn = ysgn;
+ }
+ unsigned long hi = 1;
+ while (bs < topb)
+ {
+ hi += (unsigned long)(*as++) + I_MAXNUM - (unsigned long)(*bs++);
+ *rs++ = extract(hi);
+ hi = down(hi);
+ }
+ while (hi == 0 && as < topa)
+ {
+ hi = (unsigned long)(*as++) + I_MAXNUM;
+ *rs++ = extract(hi);
+ hi = down(hi);
+ }
+ if (rs != as)
+ while (as < topa)
+ *rs++ = *as++;
+ }
+ }
+ Icheck(r);
+ return r;
+}
+
+
+IntRep* multiply(const IntRep* x, const IntRep* y, IntRep* r)
+{
+ nonnil(x);
+ nonnil(y);
+ int xl = x->len;
+ int yl = y->len;
+ int rl = xl + yl;
+ int rsgn = x->sgn == y->sgn;
+ int xrsame = x == r;
+ int yrsame = y == r;
+ int xysame = x == y;
+
+ if (xl == 0 || yl == 0)
+ r = Icopy_zero(r);
+ else if (xl == 1 && x->s[0] == 1)
+ r = Icopy(r, y);
+ else if (yl == 1 && y->s[0] == 1)
+ r = Icopy(r, x);
+ else if (!(xysame && xrsame))
+ {
+ if (xrsame || yrsame)
+ r = Iresize(r, rl);
+ else
+ r = Icalloc(r, rl);
+ unsigned short* rs = r->s;
+ unsigned short* topr = &(rs[rl]);
+
+ // use best inner/outer loop params given constraints
+ unsigned short* currentr;
+ const unsigned short* bota;
+ const unsigned short* as;
+ const unsigned short* botb;
+ const unsigned short* topb;
+ if (xrsame)
+ {
+ currentr = &(rs[xl-1]);
+ bota = rs;
+ as = currentr;
+ botb = y->s;
+ topb = &(botb[yl]);
+ }
+ else if (yrsame)
+ {
+ currentr = &(rs[yl-1]);
+ bota = rs;
+ as = currentr;
+ botb = x->s;
+ topb = &(botb[xl]);
+ }
+ else if (xl <= yl)
+ {
+ currentr = &(rs[xl-1]);
+ bota = x->s;
+ as = &(bota[xl-1]);
+ botb = y->s;
+ topb = &(botb[yl]);
+ }
+ else
+ {
+ currentr = &(rs[yl-1]);
+ bota = y->s;
+ as = &(bota[yl-1]);
+ botb = x->s;
+ topb = &(botb[xl]);
+ }
+
+ while (as >= bota)
+ {
+ unsigned long ai = (unsigned long)(*as--);
+ unsigned short* rs = currentr--;
+ *rs = 0;
+ if (ai != 0)
+ {
+ unsigned long sum = 0;
+ const unsigned short* bs = botb;
+ while (bs < topb)
+ {
+ sum += ai * (unsigned long)(*bs++) + (unsigned long)(*rs);
+ *rs++ = extract(sum);
+ sum = down(sum);
+ }
+ while (sum != 0 && rs < topr)
+ {
+ sum += (unsigned long)(*rs);
+ *rs++ = extract(sum);
+ sum = down(sum);
+ }
+ }
+ }
+ }
+ else // x, y, and r same; compute over diagonals
+ {
+ r = Iresize(r, rl);
+ unsigned short* botr = r->s;
+ unsigned short* topr = &(botr[rl]);
+ unsigned short* rs = &(botr[rl - 2]);
+
+ const unsigned short* bota = (xrsame)? botr : x->s;
+ const unsigned short* loa = &(bota[xl - 1]);
+ const unsigned short* hia = loa;
+
+ for (; rs >= botr; --rs)
+ {
+ const unsigned short* h = hia;
+ const unsigned short* l = loa;
+ unsigned long prod = (unsigned long)(*h) * (unsigned long)(*l);
+ *rs = 0;
+
+ for(;;)
+ {
+ unsigned short* rt = rs;
+ unsigned long sum = prod + (unsigned long)(*rt);
+ *rt++ = extract(sum);
+ sum = down(sum);
+ while (sum != 0 && rt < topr)
+ {
+ sum += (unsigned long)(*rt);
+ *rt++ = extract(sum);
+ sum = down(sum);
+ }
+ if (h > l)
+ {
+ rt = rs;
+ sum = prod + (unsigned long)(*rt);
+ *rt++ = extract(sum);
+ sum = down(sum);
+ while (sum != 0 && rt < topr)
+ {
+ sum += (unsigned long)(*rt);
+ *rt++ = extract(sum);
+ sum = down(sum);
+ }
+ if (--h >= ++l)
+ prod = (unsigned long)(*h) * (unsigned long)(*l);
+ else
+ break;
+ }
+ else
+ break;
+ }
+ if (loa > bota)
+ --loa;
+ else
+ --hia;
+ }
+ }
+ r->sgn = rsgn;
+ Icheck(r);
+ return r;
+}
+
+
+IntRep* multiply(const IntRep* x, long y, IntRep* r)
+{
+ nonnil(x);
+ int xl = x->len;
+
+ if (xl == 0 || y == 0)
+ r = Icopy_zero(r);
+ else if (y == 1)
+ r = Icopy(r, x);
+ else
+ {
+ int ysgn = y >= 0;
+ int rsgn = x->sgn == ysgn;
+ unsigned long uy = (ysgn)? y : -y;
+ unsigned short tmp[SHORT_PER_LONG];
+ int yl = 0;
+ while (uy != 0)
+ {
+ tmp[yl++] = extract(uy);
+ uy = down(uy);
+ }
+
+ int rl = xl + yl;
+ int xrsame = x == r;
+ if (xrsame)
+ r = Iresize(r, rl);
+ else
+ r = Icalloc(r, rl);
+
+ unsigned short* rs = r->s;
+ unsigned short* topr = &(rs[rl]);
+ unsigned short* currentr;
+ const unsigned short* bota;
+ const unsigned short* as;
+ const unsigned short* botb;
+ const unsigned short* topb;
+
+ if (xrsame)
+ {
+ currentr = &(rs[xl-1]);
+ bota = rs;
+ as = currentr;
+ botb = tmp;
+ topb = &(botb[yl]);
+ }
+ else if (xl <= yl)
+ {
+ currentr = &(rs[xl-1]);
+ bota = x->s;
+ as = &(bota[xl-1]);
+ botb = tmp;
+ topb = &(botb[yl]);
+ }
+ else
+ {
+ currentr = &(rs[yl-1]);
+ bota = tmp;
+ as = &(bota[yl-1]);
+ botb = x->s;
+ topb = &(botb[xl]);
+ }
+
+ while (as >= bota)
+ {
+ unsigned long ai = (unsigned long)(*as--);
+ unsigned short* rs = currentr--;
+ *rs = 0;
+ if (ai != 0)
+ {
+ unsigned long sum = 0;
+ const unsigned short* bs = botb;
+ while (bs < topb)
+ {
+ sum += ai * (unsigned long)(*bs++) + (unsigned long)(*rs);
+ *rs++ = extract(sum);
+ sum = down(sum);
+ }
+ while (sum != 0 && rs < topr)
+ {
+ sum += (unsigned long)(*rs);
+ *rs++ = extract(sum);
+ sum = down(sum);
+ }
+ }
+ }
+ r->sgn = rsgn;
+ }
+ Icheck(r);
+ return r;
+}
+
+
+// main division routine
+
+static void do_divide(unsigned short* rs,
+ const unsigned short* ys, int yl,
+ unsigned short* qs, int ql)
+{
+ const unsigned short* topy = &(ys[yl]);
+ unsigned short d1 = ys[yl - 1];
+ unsigned short d2 = ys[yl - 2];
+
+ int l = ql - 1;
+ int i = l + yl;
+
+ for (; l >= 0; --l, --i)
+ {
+ unsigned short qhat; // guess q
+ if (d1 == rs[i])
+ qhat = I_MAXNUM;
+ else
+ {
+ unsigned long lr = up((unsigned long)rs[i]) | rs[i-1];
+ qhat = lr / d1;
+ }
+
+ for(;;) // adjust q, use docmp to avoid overflow problems
+ {
+ unsigned short ts[3];
+ unsigned long prod = (unsigned long)d2 * (unsigned long)qhat;
+ ts[0] = extract(prod);
+ prod = down(prod) + (unsigned long)d1 * (unsigned long)qhat;
+ ts[1] = extract(prod);
+ ts[2] = extract(down(prod));
+ if (docmp(ts, &(rs[i-2]), 3) > 0)
+ --qhat;
+ else
+ break;
+ };
+
+ // multiply & subtract
+
+ const unsigned short* yt = ys;
+ unsigned short* rt = &(rs[l]);
+ unsigned long prod = 0;
+ unsigned long hi = 1;
+ while (yt < topy)
+ {
+ prod = (unsigned long)qhat * (unsigned long)(*yt++) + down(prod);
+ hi += (unsigned long)(*rt) + I_MAXNUM - (unsigned long)(extract(prod));
+ *rt++ = extract(hi);
+ hi = down(hi);
+ }
+ hi += (unsigned long)(*rt) + I_MAXNUM - (unsigned long)(down(prod));
+ *rt = extract(hi);
+ hi = down(hi);
+
+ // off-by-one, add back
+
+ if (hi == 0)
+ {
+ --qhat;
+ yt = ys;
+ rt = &(rs[l]);
+ hi = 0;
+ while (yt < topy)
+ {
+ hi = (unsigned long)(*rt) + (unsigned long)(*yt++) + down(hi);
+ *rt++ = extract(hi);
+ }
+ *rt = 0;
+ }
+ if (qs != 0)
+ qs[l] = qhat;
+ }
+}
+
+// divide by single digit, return remainder
+// if q != 0, then keep the result in q, else just compute rem
+
+static int unscale(const unsigned short* x, int xl, unsigned short y,
+ unsigned short* q)
+{
+ if (xl == 0 || y == 1)
+ return 0;
+ else if (q != 0)
+ {
+ unsigned short* botq = q;
+ unsigned short* qs = &(botq[xl - 1]);
+ const unsigned short* xs = &(x[xl - 1]);
+ unsigned long rem = 0;
+ while (qs >= botq)
+ {
+ rem = up(rem) | *xs--;
+ unsigned long u = rem / y;
+ *qs-- = extract(u);
+ rem -= u * y;
+ }
+ int r = extract(rem);
+ return r;
+ }
+ else // same loop, a bit faster if just need rem
+ {
+ const unsigned short* botx = x;
+ const unsigned short* xs = &(botx[xl - 1]);
+ unsigned long rem = 0;
+ while (xs >= botx)
+ {
+ rem = up(rem) | *xs--;
+ unsigned long u = rem / y;
+ rem -= u * y;
+ }
+ int r = extract(rem);
+ return r;
+ }
+}
+
+
+IntRep* div(const IntRep* x, const IntRep* y, IntRep* q)
+{
+ nonnil(x);
+ nonnil(y);
+ int xl = x->len;
+ int yl = y->len;
+ if (yl == 0) (*lib_error_handler)("Integer", "attempted division by zero");
+
+ int comp = ucompare(x, y);
+ int xsgn = x->sgn;
+ int ysgn = y->sgn;
+
+ int samesign = xsgn == ysgn;
+
+ if (comp < 0)
+ q = Icopy_zero(q);
+ else if (comp == 0)
+ q = Icopy_one(q, samesign);
+ else if (yl == 1)
+ {
+ q = Icopy(q, x);
+ unscale(q->s, q->len, y->s[0], q->s);
+ }
+ else
+ {
+ IntRep* yy = 0;
+ IntRep* r = 0;
+ unsigned short prescale = (I_RADIX / (1 + y->s[yl - 1]));
+ if (prescale != 1 || y == q)
+ {
+ yy = multiply(y, ((long)prescale & I_MAXNUM), yy);
+ r = multiply(x, ((long)prescale & I_MAXNUM), r);
+ }
+ else
+ {
+ yy = (IntRep*)y;
+ r = Icalloc(r, xl + 1);
+ scpy(x->s, r->s, xl);
+ }
+
+ int ql = xl - yl + 1;
+
+ q = Icalloc(q, ql);
+ do_divide(r->s, yy->s, yl, q->s, ql);
+
+ if (yy != y && !STATIC_IntRep(yy)) delete yy;
+ if (!STATIC_IntRep(r)) delete r;
+ }
+ q->sgn = samesign;
+ Icheck(q);
+ return q;
+}
+
+IntRep* div(const IntRep* x, long y, IntRep* q)
+{
+ nonnil(x);
+ int xl = x->len;
+ if (y == 0) (*lib_error_handler)("Integer", "attempted division by zero");
+
+ unsigned short ys[SHORT_PER_LONG];
+ unsigned long u;
+ int ysgn = y >= 0;
+ if (ysgn)
+ u = y;
+ else
+ u = -y;
+ int yl = 0;
+ while (u != 0)
+ {
+ ys[yl++] = extract(u);
+ u = down(u);
+ }
+
+ int comp = xl - yl;
+ if (comp == 0) comp = docmp(x->s, ys, xl);
+
+ int xsgn = x->sgn;
+ int samesign = xsgn == ysgn;
+
+ if (comp < 0)
+ q = Icopy_zero(q);
+ else if (comp == 0)
+ {
+ q = Icopy_one(q, samesign);
+ }
+ else if (yl == 1)
+ {
+ q = Icopy(q, x);
+ unscale(q->s, q->len, ys[0], q->s);
+ }
+ else
+ {
+ IntRep* r = 0;
+ unsigned short prescale = (I_RADIX / (1 + ys[yl - 1]));
+ if (prescale != 1)
+ {
+ unsigned long prod = (unsigned long)prescale * (unsigned long)ys[0];
+ ys[0] = extract(prod);
+ prod = down(prod) + (unsigned long)prescale * (unsigned long)ys[1];
+ ys[1] = extract(prod);
+ r = multiply(x, ((long)prescale & I_MAXNUM), r);
+ }
+ else
+ {
+ r = Icalloc(r, xl + 1);
+ scpy(x->s, r->s, xl);
+ }
+
+ int ql = xl - yl + 1;
+
+ q = Icalloc(q, ql);
+ do_divide(r->s, ys, yl, q->s, ql);
+
+ if (!STATIC_IntRep(r)) delete r;
+ }
+ q->sgn = samesign;
+ Icheck(q);
+ return q;
+}
+
+
+void divide(const Integer& Ix, long y, Integer& Iq, long& rem)
+{
+ const IntRep* x = Ix.rep;
+ nonnil(x);
+ IntRep* q = Iq.rep;
+ int xl = x->len;
+ if (y == 0) (*lib_error_handler)("Integer", "attempted division by zero");
+
+ unsigned short ys[SHORT_PER_LONG];
+ unsigned long u;
+ int ysgn = y >= 0;
+ if (ysgn)
+ u = y;
+ else
+ u = -y;
+ int yl = 0;
+ while (u != 0)
+ {
+ ys[yl++] = extract(u);
+ u = down(u);
+ }
+
+ int comp = xl - yl;
+ if (comp == 0) comp = docmp(x->s, ys, xl);
+
+ int xsgn = x->sgn;
+ int samesign = xsgn == ysgn;
+
+ if (comp < 0)
+ {
+ rem = Itolong(x);
+ q = Icopy_zero(q);
+ }
+ else if (comp == 0)
+ {
+ q = Icopy_one(q, samesign);
+ rem = 0;
+ }
+ else if (yl == 1)
+ {
+ q = Icopy(q, x);
+ rem = unscale(q->s, q->len, ys[0], q->s);
+ }
+ else
+ {
+ IntRep* r = 0;
+ unsigned short prescale = (I_RADIX / (1 + ys[yl - 1]));
+ if (prescale != 1)
+ {
+ unsigned long prod = (unsigned long)prescale * (unsigned long)ys[0];
+ ys[0] = extract(prod);
+ prod = down(prod) + (unsigned long)prescale * (unsigned long)ys[1];
+ ys[1] = extract(prod);
+ r = multiply(x, ((long)prescale & I_MAXNUM), r);
+ }
+ else
+ {
+ r = Icalloc(r, xl + 1);
+ scpy(x->s, r->s, xl);
+ }
+
+ int ql = xl - yl + 1;
+
+ q = Icalloc(q, ql);
+
+ do_divide(r->s, ys, yl, q->s, ql);
+
+ if (prescale != 1)
+ {
+ Icheck(r);
+ unscale(r->s, r->len, prescale, r->s);
+ }
+ Icheck(r);
+ rem = Itolong(r);
+ if (!STATIC_IntRep(r)) delete r;
+ }
+ rem = abs(rem);
+ if (xsgn == I_NEGATIVE) rem = -rem;
+ q->sgn = samesign;
+ Icheck(q);
+ Iq.rep = q;
+}
+
+
+void divide(const Integer& Ix, const Integer& Iy, Integer& Iq, Integer& Ir)
+{
+ const IntRep* x = Ix.rep;
+ nonnil(x);
+ const IntRep* y = Iy.rep;
+ nonnil(y);
+ IntRep* q = Iq.rep;
+ IntRep* r = Ir.rep;
+
+ int xl = x->len;
+ int yl = y->len;
+ if (yl == 0)
+ (*lib_error_handler)("Integer", "attempted division by zero");
+
+ int comp = ucompare(x, y);
+ int xsgn = x->sgn;
+ int ysgn = y->sgn;
+
+ int samesign = xsgn == ysgn;
+
+ if (comp < 0)
+ {
+ q = Icopy_zero(q);
+ r = Icopy(r, x);
+ }
+ else if (comp == 0)
+ {
+ q = Icopy_one(q, samesign);
+ r = Icopy_zero(r);
+ }
+ else if (yl == 1)
+ {
+ q = Icopy(q, x);
+ int rem = unscale(q->s, q->len, y->s[0], q->s);
+ r = Icopy_long(r, rem);
+ if (rem != 0)
+ r->sgn = xsgn;
+ }
+ else
+ {
+ IntRep* yy = 0;
+ unsigned short prescale = (I_RADIX / (1 + y->s[yl - 1]));
+ if (prescale != 1 || y == q || y == r)
+ {
+ yy = multiply(y, ((long)prescale & I_MAXNUM), yy);
+ r = multiply(x, ((long)prescale & I_MAXNUM), r);
+ }
+ else
+ {
+ yy = (IntRep*)y;
+ r = Icalloc(r, xl + 1);
+ scpy(x->s, r->s, xl);
+ }
+
+ int ql = xl - yl + 1;
+
+ q = Icalloc(q, ql);
+ do_divide(r->s, yy->s, yl, q->s, ql);
+
+ if (yy != y && !STATIC_IntRep(yy)) delete yy;
+ if (prescale != 1)
+ {
+ Icheck(r);
+ unscale(r->s, r->len, prescale, r->s);
+ }
+ }
+ q->sgn = samesign;
+ Icheck(q);
+ Iq.rep = q;
+ Icheck(r);
+ Ir.rep = r;
+}
+
+IntRep* mod(const IntRep* x, const IntRep* y, IntRep* r)
+{
+ nonnil(x);
+ nonnil(y);
+ int xl = x->len;
+ int yl = y->len;
+ if (yl == 0) (*lib_error_handler)("Integer", "attempted division by zero");
+
+ int comp = ucompare(x, y);
+ int xsgn = x->sgn;
+
+ if (comp < 0)
+ r = Icopy(r, x);
+ else if (comp == 0)
+ r = Icopy_zero(r);
+ else if (yl == 1)
+ {
+ int rem = unscale(x->s, xl, y->s[0], 0);
+ r = Icopy_long(r, rem);
+ if (rem != 0)
+ r->sgn = xsgn;
+ }
+ else
+ {
+ IntRep* yy = 0;
+ unsigned short prescale = (I_RADIX / (1 + y->s[yl - 1]));
+ if (prescale != 1 || y == r)
+ {
+ yy = multiply(y, ((long)prescale & I_MAXNUM), yy);
+ r = multiply(x, ((long)prescale & I_MAXNUM), r);
+ }
+ else
+ {
+ yy = (IntRep*)y;
+ r = Icalloc(r, xl + 1);
+ scpy(x->s, r->s, xl);
+ }
+
+ do_divide(r->s, yy->s, yl, 0, xl - yl + 1);
+
+ if (yy != y && !STATIC_IntRep(yy)) delete yy;
+
+ if (prescale != 1)
+ {
+ Icheck(r);
+ unscale(r->s, r->len, prescale, r->s);
+ }
+ }
+ Icheck(r);
+ return r;
+}
+
+IntRep* mod(const IntRep* x, long y, IntRep* r)
+{
+ nonnil(x);
+ int xl = x->len;
+ if (y == 0) (*lib_error_handler)("Integer", "attempted division by zero");
+
+ unsigned short ys[SHORT_PER_LONG];
+ unsigned long u;
+ int ysgn = y >= 0;
+ if (ysgn)
+ u = y;
+ else
+ u = -y;
+ int yl = 0;
+ while (u != 0)
+ {
+ ys[yl++] = extract(u);
+ u = down(u);
+ }
+
+ int comp = xl - yl;
+ if (comp == 0) comp = docmp(x->s, ys, xl);
+
+ int xsgn = x->sgn;
+
+ if (comp < 0)
+ r = Icopy(r, x);
+ else if (comp == 0)
+ r = Icopy_zero(r);
+ else if (yl == 1)
+ {
+ int rem = unscale(x->s, xl, ys[0], 0);
+ r = Icopy_long(r, rem);
+ if (rem != 0)
+ r->sgn = xsgn;
+ }
+ else
+ {
+ unsigned short prescale = (I_RADIX / (1 + ys[yl - 1]));
+ if (prescale != 1)
+ {
+ unsigned long prod = (unsigned long)prescale * (unsigned long)ys[0];
+ ys[0] = extract(prod);
+ prod = down(prod) + (unsigned long)prescale * (unsigned long)ys[1];
+ ys[1] = extract(prod);
+ r = multiply(x, ((long)prescale & I_MAXNUM), r);
+ }
+ else
+ {
+ r = Icalloc(r, xl + 1);
+ scpy(x->s, r->s, xl);
+ }
+
+ do_divide(r->s, ys, yl, 0, xl - yl + 1);
+
+ if (prescale != 1)
+ {
+ Icheck(r);
+ unscale(r->s, r->len, prescale, r->s);
+ }
+ }
+ Icheck(r);
+ return r;
+}
+
+IntRep* lshift(const IntRep* x, long y, IntRep* r)
+{
+ nonnil(x);
+ int xl = x->len;
+ if (xl == 0 || y == 0)
+ {
+ r = Icopy(r, x);
+ return r;
+ }
+
+ int xrsame = x == r;
+ int rsgn = x->sgn;
+
+ long ay = (y < 0)? -y : y;
+ int bw = ay / I_SHIFT;
+ int sw = ay % I_SHIFT;
+
+ if (y > 0)
+ {
+ int rl = bw + xl + 1;
+ if (xrsame)
+ r = Iresize(r, rl);
+ else
+ r = Icalloc(r, rl);
+
+ unsigned short* botr = r->s;
+ unsigned short* rs = &(botr[rl - 1]);
+ const unsigned short* botx = (xrsame)? botr : x->s;
+ const unsigned short* xs = &(botx[xl - 1]);
+ unsigned long a = 0;
+ while (xs >= botx)
+ {
+ a = up(a) | ((unsigned long)(*xs--) << sw);
+ *rs-- = extract(down(a));
+ }
+ *rs-- = extract(a);
+ while (rs >= botr)
+ *rs-- = 0;
+ }
+ else
+ {
+ int rl = xl - bw;
+ if (rl < 0)
+ r = Icopy_zero(r);
+ else
+ {
+ if (xrsame)
+ r = Iresize(r, rl);
+ else
+ r = Icalloc(r, rl);
+ int rw = I_SHIFT - sw;
+ unsigned short* rs = r->s;
+ unsigned short* topr = &(rs[rl]);
+ const unsigned short* botx = (xrsame)? rs : x->s;
+ const unsigned short* xs = &(botx[bw]);
+ const unsigned short* topx = &(botx[xl]);
+ unsigned long a = (unsigned long)(*xs++) >> sw;
+ while (xs < topx)
+ {
+ a |= (unsigned long)(*xs++) << rw;
+ *rs++ = extract(a);
+ a = down(a);
+ }
+ *rs++ = extract(a);
+ if (xrsame) topr = (unsigned short*)topx;
+ while (rs < topr)
+ *rs++ = 0;
+ }
+ }
+ r->sgn = rsgn;
+ Icheck(r);
+ return r;
+}
+
+IntRep* lshift(const IntRep* x, const IntRep* yy, int negatey, IntRep* r)
+{
+ long y = Itolong(yy);
+ if (negatey)
+ y = -y;
+
+ return lshift(x, y, r);
+}
+
+IntRep* bitop(const IntRep* x, const IntRep* y, IntRep* r, char op)
+{
+ nonnil(x);
+ nonnil(y);
+ int xl = x->len;
+ int yl = y->len;
+ int xsgn = x->sgn;
+ int xrsame = x == r;
+ int yrsame = y == r;
+ if (xrsame || yrsame)
+ r = Iresize(r, calc_len(xl, yl, 0));
+ else
+ r = Icalloc(r, calc_len(xl, yl, 0));
+ r->sgn = xsgn;
+ unsigned short* rs = r->s;
+ unsigned short* topr = &(rs[r->len]);
+ const unsigned short* as;
+ const unsigned short* bs;
+ const unsigned short* topb;
+ if (xl >= yl)
+ {
+ as = (xrsame)? rs : x->s;
+ bs = (yrsame)? rs : y->s;
+ topb = &(bs[yl]);
+ }
+ else
+ {
+ bs = (xrsame)? rs : x->s;
+ topb = &(bs[xl]);
+ as = (yrsame)? rs : y->s;
+ }
+
+ switch (op)
+ {
+ case '&':
+ while (bs < topb) *rs++ = *as++ & *bs++;
+ while (rs < topr) *rs++ = 0;
+ break;
+ case '|':
+ while (bs < topb) *rs++ = *as++ | *bs++;
+ while (rs < topr) *rs++ = *as++;
+ break;
+ case '^':
+ while (bs < topb) *rs++ = *as++ ^ *bs++;
+ while (rs < topr) *rs++ = *as++;
+ break;
+ }
+ Icheck(r);
+ return r;
+}
+
+IntRep* bitop(const IntRep* x, long y, IntRep* r, char op)
+{
+ nonnil(x);
+ unsigned short tmp[SHORT_PER_LONG];
+ unsigned long u;
+ int newsgn;
+ if (newsgn = (y >= 0))
+ u = y;
+ else
+ u = -y;
+
+ int l = 0;
+ while (u != 0)
+ {
+ tmp[l++] = extract(u);
+ u = down(u);
+ }
+
+ int xl = x->len;
+ int yl = l;
+ int xsgn = x->sgn;
+ int xrsame = x == r;
+ if (xrsame)
+ r = Iresize(r, calc_len(xl, yl, 0));
+ else
+ r = Icalloc(r, calc_len(xl, yl, 0));
+ r->sgn = xsgn;
+ unsigned short* rs = r->s;
+ unsigned short* topr = &(rs[r->len]);
+ const unsigned short* as;
+ const unsigned short* bs;
+ const unsigned short* topb;
+ if (xl >= yl)
+ {
+ as = (xrsame)? rs : x->s;
+ bs = tmp;
+ topb = &(bs[yl]);
+ }
+ else
+ {
+ bs = (xrsame)? rs : x->s;
+ topb = &(bs[xl]);
+ as = tmp;
+ }
+
+ switch (op)
+ {
+ case '&':
+ while (bs < topb) *rs++ = *as++ & *bs++;
+ while (rs < topr) *rs++ = 0;
+ break;
+ case '|':
+ while (bs < topb) *rs++ = *as++ | *bs++;
+ while (rs < topr) *rs++ = *as++;
+ break;
+ case '^':
+ while (bs < topb) *rs++ = *as++ ^ *bs++;
+ while (rs < topr) *rs++ = *as++;
+ break;
+ }
+ Icheck(r);
+ return r;
+}
+
+
+
+IntRep* I_compl(const IntRep* src, IntRep* r)
+{
+ nonnil(src);
+ r = Icopy(r, src);
+ unsigned short* s = r->s;
+ unsigned short* top = &(s[r->len - 1]);
+ while (s < top)
+ {
+ unsigned short cmp = ~(*s);
+ *s++ = cmp;
+ }
+ unsigned short a = *s;
+ unsigned short b = 0;
+ while (a != 0)
+ {
+ b <<= 1;
+ if (!(a & 1)) b |= 1;
+ a >>= 1;
+ }
+ *s = b;
+ Icheck(r);
+ return r;
+}
+
+void (setbit)(Integer& x, long b)
+{
+ if (b >= 0)
+ {
+ int bw = (unsigned long)b / I_SHIFT;
+ int sw = (unsigned long)b % I_SHIFT;
+ if (x.rep == 0)
+ x.rep = Icalloc(0, bw + 1);
+ else if (x.rep->len < bw)
+ {
+ int xl = x.rep->len;
+ x.rep = Iresize(x.rep, calc_len(xl, bw+1, 0));
+ }
+ x.rep->s[bw] |= (1 << sw);
+ Icheck(x.rep);
+ }
+}
+
+void clearbit(Integer& x, long b)
+{
+ if (b >= 0)
+ {
+ int bw = (unsigned long)b / I_SHIFT;
+ int sw = (unsigned long)b % I_SHIFT;
+ if (x.rep == 0)
+ x.rep = Icalloc(0, bw + 1);
+ else if (x.rep->len < bw)
+ {
+ int xl = x.rep->len;
+ x.rep = Iresize(x.rep, calc_len(xl, bw+1, 0));
+ }
+ x.rep->s[bw] &= ~(1 << sw);
+ Icheck(x.rep);
+ }
+}
+
+int testbit(const Integer& x, long b)
+{
+ if (x.rep != 0 && b >= 0)
+ {
+ int bw = (unsigned long)b / I_SHIFT;
+ int sw = (unsigned long)b % I_SHIFT;
+ return (bw < x.rep->len && (x.rep->s[bw] & (1 << sw)) != 0);
+ }
+ else
+ return 0;
+}
+
+// A version of knuth's algorithm B / ex. 4.5.3.34
+// A better version that doesn't bother shifting all of `t' forthcoming
+
+IntRep* gcd(const IntRep* x, const IntRep* y)
+{
+ nonnil(x);
+ nonnil(y);
+ int ul = x->len;
+ int vl = y->len;
+
+ if (vl == 0)
+ return Ialloc(0, x->s, ul, I_POSITIVE, ul);
+ else if (ul == 0)
+ return Ialloc(0, y->s, vl, I_POSITIVE, vl);
+
+ IntRep* u = Ialloc(0, x->s, ul, I_POSITIVE, ul);
+ IntRep* v = Ialloc(0, y->s, vl, I_POSITIVE, vl);
+
+// find shift so that both not even
+
+ long k = 0;
+ int l = (ul <= vl)? ul : vl;
+ int cont = 1;
+ for (int i = 0; i < l && cont; ++i)
+ {
+ unsigned long a = (i < ul)? u->s[i] : 0;
+ unsigned long b = (i < vl)? v->s[i] : 0;
+ for (int j = 0; j < I_SHIFT; ++j)
+ {
+ if ((a | b) & 1)
+ {
+ cont = 0;
+ break;
+ }
+ else
+ {
+ ++k;
+ a >>= 1;
+ b >>= 1;
+ }
+ }
+ }
+
+ if (k != 0)
+ {
+ u = lshift(u, -k, u);
+ v = lshift(v, -k, v);
+ }
+
+ IntRep* t;
+ if (u->s[0] & 01)
+ t = Ialloc(0, v->s, v->len, !v->sgn, v->len);
+ else
+ t = Ialloc(0, u->s, u->len, u->sgn, u->len);
+
+ while (t->len != 0)
+ {
+ long s = 0; // shift t until odd
+ cont = 1;
+ int tl = t->len;
+ for (int i = 0; i < tl && cont; ++i)
+ {
+ unsigned long a = t->s[i];
+ for (int j = 0; j < I_SHIFT; ++j)
+ {
+ if (a & 1)
+ {
+ cont = 0;
+ break;
+ }
+ else
+ {
+ ++s;
+ a >>= 1;
+ }
+ }
+ }
+
+ if (s != 0) t = lshift(t, -s, t);
+
+ if (t->sgn == I_POSITIVE)
+ {
+ u = Icopy(u, t);
+ t = add(t, 0, v, 1, t);
+ }
+ else
+ {
+ v = Ialloc(v, t->s, t->len, !t->sgn, t->len);
+ t = add(t, 0, u, 0, t);
+ }
+ }
+ if (!STATIC_IntRep(t)) delete t;
+ if (!STATIC_IntRep(v)) delete v;
+ if (k != 0) u = lshift(u, k, u);
+ return u;
+}
+
+
+
+long lg(const IntRep* x)
+{
+ nonnil(x);
+ int xl = x->len;
+ if (xl == 0)
+ return 0;
+
+ long l = (xl - 1) * I_SHIFT - 1;
+ unsigned short a = x->s[xl-1];
+
+ while (a != 0)
+ {
+ a = a >> 1;
+ ++l;
+ }
+ return l;
+}
+
+IntRep* power(const IntRep* x, long y, IntRep* r)
+{
+ nonnil(x);
+ int sgn;
+ if (x->sgn == I_POSITIVE || (!(y & 1)))
+ sgn = I_POSITIVE;
+ else
+ sgn = I_NEGATIVE;
+
+ int xl = x->len;
+
+ if (y == 0 || (xl == 1 && x->s[0] == 1))
+ r = Icopy_one(r, sgn);
+ else if (xl == 0 || y < 0)
+ r = Icopy_zero(r);
+ else if (y == 1 || y == -1)
+ r = Icopy(r, x);
+ else
+ {
+ int maxsize = ((lg(x) + 1) * y) / I_SHIFT + 2; // pre-allocate space
+ IntRep* b = Ialloc(0, x->s, xl, I_POSITIVE, maxsize);
+ b->len = xl;
+ r = Icalloc(r, maxsize);
+ r = Icopy_one(r, I_POSITIVE);
+ for(;;)
+ {
+ if (y & 1)
+ r = multiply(r, b, r);
+ if ((y >>= 1) == 0)
+ break;
+ else
+ b = multiply(b, b, b);
+ }
+ if (!STATIC_IntRep(b)) delete b;
+ }
+ r->sgn = sgn;
+ Icheck(r);
+ return r;
+}
+
+IntRep* abs(const IntRep* src, IntRep* dest)
+{
+ nonnil(src);
+ if (src != dest)
+ dest = Icopy(dest, src);
+ dest->sgn = I_POSITIVE;
+ return dest;
+}
+
+IntRep* negate(const IntRep* src, IntRep* dest)
+{
+ nonnil(src);
+ if (src != dest)
+ dest = Icopy(dest, src);
+ if (dest->len != 0)
+ dest->sgn = !dest->sgn;
+ return dest;
+}
+
+Integer sqrt(const Integer& x)
+{
+ Integer r(x);
+ int s = sign(x);
+ if (s < 0) x.error("Attempted square root of negative Integer");
+ if (s != 0)
+ {
+ r >>= (lg(x) / 2); // get close
+ Integer q;
+ div(x, r, q);
+ while (q < r)
+ {
+ r += q;
+ r >>= 1;
+ div(x, r, q);
+ }
+ }
+ return r;
+}
+
+Integer lcm(const Integer& x, const Integer& y)
+{
+ Integer r;
+ if (!x.initialized() || !y.initialized())
+ x.error("operation on uninitialized Integer");
+ Integer g;
+ if (sign(x) == 0 || sign(y) == 0)
+ g = 1;
+ else
+ g = gcd(x, y);
+ div(x, g, r);
+ mul(r, y, r);
+ return r;
+}
+
+
+
+IntRep* atoIntRep(const char* s, int base)
+{
+ int sl = strlen(s);
+ IntRep* r = Icalloc(0, sl * (lg(base) + 1) / I_SHIFT + 1);
+ if (s != 0)
+ {
+ char sgn;
+ while (isspace(*s)) ++s;
+ if (*s == '-')
+ {
+ sgn = I_NEGATIVE;
+ s++;
+ }
+ else if (*s == '+')
+ {
+ sgn = I_POSITIVE;
+ s++;
+ }
+ else
+ sgn = I_POSITIVE;
+ for (;;)
+ {
+ long digit;
+ if (*s >= '0' && *s <= '9') digit = *s - '0';
+ else if (*s >= 'a' && *s <= 'z') digit = *s - 'a' + 10;
+ else if (*s >= 'A' && *s <= 'Z') digit = *s - 'A' + 10;
+ else break;
+ if (digit >= base) break;
+ r = multiply(r, base, r);
+ r = add(r, 0, digit, r);
+ ++s;
+ }
+ r->sgn = sgn;
+ }
+ return r;
+}
+
+
+
+//extern AllocRing _libgxx_fmtq;
+
+char* Itoa(const IntRep* x, int base, int width)
+{
+ int fmtlen = (x->len + 1) * I_SHIFT / lg(base) + 4 + width;
+ char* fmtbase = new char[fmtlen]; //(char *) _libgxx_fmtq.alloc(fmtlen);
+ char* f = cvtItoa(x, fmtbase, fmtlen, base, 0, width, 0, ' ', 'X', 0);
+ return f;
+}
+
+std::ostream& operator << (std::ostream& s, const Integer& y)
+{
+ return s << Itoa(y.rep);
+}
+
+void Integer::printon(std::ostream& s, int base /* =10 */, int width /* =0 */) const
+{
+ int align_right = !(s.flags() & std::ios::left);
+ int showpos = s.flags() & std::ios::showpos;
+ int showbase = s.flags() & std::ios::showbase;
+ char fillchar = s.fill();
+ char Xcase = (s.flags() & std::ios::uppercase)? 'X' : 'x';
+ const IntRep* x = rep;
+ int fmtlen = (x->len + 1) * I_SHIFT / lg(base) + 4 + width;
+ char* fmtbase = new char[fmtlen];
+ char* f = cvtItoa(x, fmtbase, fmtlen, base, showbase, width, align_right,
+ fillchar, Xcase, showpos);
+ s.write(f, fmtlen);
+ delete fmtbase;
+}
+
+char* cvtItoa(const IntRep* x, char* fmt, int& fmtlen, int base, int showbase,
+ int width, int align_right, char fillchar, char Xcase,
+ int showpos)
+{
+ char* e = fmt + fmtlen - 1;
+ char* s = e;
+ *--s = 0;
+
+ if (x->len == 0)
+ *--s = '0';
+ else
+ {
+ IntRep* z = Icopy(0, x);
+
+ // split division by base into two parts:
+ // first divide by biggest power of base that fits in an unsigned short,
+ // then use straight signed div/mods from there.
+
+ // find power
+ int bpower = 1;
+ unsigned short b = base;
+ unsigned short maxb = I_MAXNUM / base;
+ while (b < maxb)
+ {
+ b *= base;
+ ++bpower;
+ }
+ for(;;)
+ {
+ int rem = unscale(z->s, z->len, b, z->s);
+ Icheck(z);
+ if (z->len == 0)
+ {
+ while (rem != 0)
+ {
+ char ch = rem % base;
+ rem /= base;
+ if (ch >= 10)
+ ch += 'a' - 10;
+ else
+ ch += '0';
+ *--s = ch;
+ }
+ if (!STATIC_IntRep(z)) delete z;
+ break;
+ }
+ else
+ {
+ for (int i = 0; i < bpower; ++i)
+ {
+ char ch = rem % base;
+ rem /= base;
+ if (ch >= 10)
+ ch += 'a' - 10;
+ else
+ ch += '0';
+ *--s = ch;
+ }
+ }
+ }
+ }
+
+ if (base == 8 && showbase)
+ *--s = '0';
+ else if (base == 16 && showbase)
+ {
+ *--s = Xcase;
+ *--s = '0';
+ }
+ if (x->sgn == I_NEGATIVE) *--s = '-';
+ else if (showpos) *--s = '+';
+ int w = e - s - 1;
+ if (!align_right || w >= width)
+ {
+ while (w++ < width)
+ *--s = fillchar;
+ fmtlen = e - s - 1;
+ return s;
+ }
+ else
+ {
+ char* p = fmt;
+ int gap = s - p;
+ for (char* t = s; *t != 0; ++t, ++p) *p = *t;
+ while (w++ < width) *p++ = fillchar;
+ *p = 0;
+ fmtlen = p - fmt;
+ return fmt;
+ }
+}
+
+char* dec(const Integer& x, int width)
+{
+ return Itoa(x, 10, width);
+}
+
+char* oct(const Integer& x, int width)
+{
+ return Itoa(x, 8, width);
+}
+
+char* hex(const Integer& x, int width)
+{
+ return Itoa(x, 16, width);
+}
+
+std::istream& operator >> (std::istream& s, Integer& y)
+{
+ if (!s.good())
+ return s;
+
+ s >> std::ws;
+ if (!s.good())
+ {
+ s.setstate(std::ios_base::failbit);
+ return s;
+ }
+
+#ifdef _OLD_STREAMS
+ int know_base = 0;
+ int base = 10;
+#else
+ int know_base = s.flags() & (std::ios::oct | std::ios::hex | std::ios::dec);
+ int base = (s.flags() & std::ios::oct) ? 8 : (s.flags() & std::ios::hex) ? 16 : 10;
+#endif
+
+ int got_one = 0;
+ char sgn = 0;
+ char ch;
+ y.rep = Icopy_zero(y.rep);
+
+ while (s.get(ch))
+ {
+ if (ch == '-')
+ {
+ if (sgn == 0)
+ sgn = '-';
+ else
+ break;
+ }
+ else if (!know_base & !got_one && ch == '0')
+ base = 8, got_one = 1;
+ else if (!know_base & !got_one && base == 8 && (ch == 'X' || ch == 'x'))
+ base = 16;
+ else if (base == 8)
+ {
+ if (ch >= '0' && ch <= '7')
+ {
+ long digit = ch - '0';
+ y <<= 3;
+ y += digit;
+ got_one = 1;
+ }
+ else
+ break;
+ }
+ else if (base == 16)
+ {
+ long digit;
+ if (ch >= '0' && ch <= '9')
+ digit = ch - '0';
+ else if (ch >= 'A' && ch <= 'F')
+ digit = ch - 'A' + 10;
+ else if (ch >= 'a' && ch <= 'f')
+ digit = ch - 'a' + 10;
+ else
+ digit = base;
+ if (digit < base)
+ {
+ y <<= 4;
+ y += digit;
+ got_one = 1;
+ }
+ else
+ break;
+ }
+ else if (base == 10)
+ {
+ if (ch >= '0' && ch <= '9')
+ {
+ long digit = ch - '0';
+ y *= 10;
+ y += digit;
+ got_one = 1;
+ }
+ else
+ break;
+ }
+ else
+ abort(); // can't happen for now
+ }
+ if (s.good())
+ s.putback(ch);
+ if (!got_one)
+ s.setstate(std::ios_base::failbit);
+
+ if (sgn == '-')
+ y.negate();
+
+ return s;
+}
+
+int Integer::OK() const
+{
+ if (rep != 0)
+ {
+ int l = rep->len;
+ int s = rep->sgn;
+ int v = l <= rep->sz || STATIC_IntRep(rep); // length within bounds
+ v &= s == 0 || s == 1; // legal sign
+ Icheck(rep); // and correctly adjusted
+ v &= rep->len == l;
+ v &= rep->sgn == s;
+ if (v)
+ return v;
+ }
+ error("invariant failure");
+ return 0;
+}
+
+void Integer::error(const char* msg) const
+{
+ (*lib_error_handler)("Integer", msg);
+}
+
--- /dev/null
+// This may look like C code, but it is really -*- C++ -*-
+
+/*
+Copyright (C) 1988 Free Software Foundation
+ written by Doug Lea (dl@rocky.oswego.edu)
+
+This file is part of the GNU C++ Library. This library is free
+software; you can redistribute it and/or modify it under the terms of
+the GNU Library General Public License as published by the Free
+Software Foundation; either version 2 of the License, or (at your
+option) any later version. This library 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 Library General Public License for more details.
+You should have received a copy of the GNU Library General Public
+License along with this library; if not, write to the Free Software
+Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
+*/
+
+#ifndef _Integer_h
+#ifdef __GNUG__
+#pragma interface
+#endif
+#define _Integer_h 1
+
+#include <iostream>
+
+struct IntRep // internal Integer representations
+{
+ unsigned short len; // current length
+ unsigned short sz; // allocated space (0 means static).
+ short sgn; // 1 means >= 0; 0 means < 0
+ unsigned short s[1]; // represented as ushort array starting here
+};
+
+// True if REP is staticly (or manually) allocated,
+// and should not be deleted by an Integer destructor.
+#define STATIC_IntRep(rep) ((rep)->sz==0)
+
+extern IntRep* Ialloc(IntRep*, const unsigned short *, int, int, int);
+extern IntRep* Icalloc(IntRep*, int);
+extern IntRep* Icopy_ulong(IntRep*, unsigned long);
+extern IntRep* Icopy_long(IntRep*, long);
+extern IntRep* Icopy(IntRep*, const IntRep*);
+extern IntRep* Iresize(IntRep*, int);
+extern IntRep* add(const IntRep*, int, const IntRep*, int, IntRep*);
+extern IntRep* add(const IntRep*, int, long, IntRep*);
+extern IntRep* multiply(const IntRep*, const IntRep*, IntRep*);
+extern IntRep* multiply(const IntRep*, long, IntRep*);
+extern IntRep* lshift(const IntRep*, long, IntRep*);
+extern IntRep* lshift(const IntRep*, const IntRep*, int, IntRep*);
+extern IntRep* bitop(const IntRep*, const IntRep*, IntRep*, char);
+extern IntRep* bitop(const IntRep*, long, IntRep*, char);
+extern IntRep* power(const IntRep*, long, IntRep*);
+extern IntRep* div(const IntRep*, const IntRep*, IntRep*);
+extern IntRep* mod(const IntRep*, const IntRep*, IntRep*);
+extern IntRep* div(const IntRep*, long, IntRep*);
+extern IntRep* mod(const IntRep*, long, IntRep*);
+extern IntRep* I_compl(const IntRep*, IntRep*);
+extern IntRep* abs(const IntRep*, IntRep*);
+extern IntRep* negate(const IntRep*, IntRep*);
+extern IntRep* pow(const IntRep*, long);
+extern IntRep* gcd(const IntRep*, const IntRep* y);
+extern int compare(const IntRep*, const IntRep*);
+extern int compare(const IntRep*, long);
+extern int ucompare(const IntRep*, const IntRep*);
+extern int ucompare(const IntRep*, long);
+extern char* Itoa(const IntRep* x, int base = 10, int width = 0);
+extern char* cvtItoa(const IntRep* x, char* fmt, int& fmtlen, int base,
+ int showbase, int width, int align_right,
+ char fillchar, char Xcase, int showpos);
+extern IntRep* atoIntRep(const char* s, int base = 10);
+extern long Itolong(const IntRep*);
+extern double Itodouble(const IntRep*);
+extern int Iislong(const IntRep*);
+extern int Iisdouble(const IntRep*);
+extern long lg(const IntRep*);
+
+extern IntRep _ZeroRep, _OneRep, _MinusOneRep;
+
+class Integer
+{
+protected:
+ IntRep* rep;
+public:
+ Integer();
+ Integer(long);
+ Integer(unsigned long);
+ Integer(IntRep*);
+ Integer(const Integer&);
+
+ ~Integer();
+
+ void operator = (const Integer&);
+ void operator = (long);
+
+// unary operations to self
+
+ void operator ++ ();
+ void operator -- ();
+ void negate(); // negate in-place
+ void abs(); // absolute-value in-place
+ void complement(); // bitwise complement in-place
+
+// assignment-based operations
+
+ void operator += (const Integer&);
+ void operator -= (const Integer&);
+ void operator *= (const Integer&);
+ void operator /= (const Integer&);
+ void operator %= (const Integer&);
+ void operator <<=(const Integer&);
+ void operator >>=(const Integer&);
+ void operator &= (const Integer&);
+ void operator |= (const Integer&);
+ void operator ^= (const Integer&);
+
+ void operator += (long);
+ void operator -= (long);
+ void operator *= (long);
+ void operator /= (long);
+ void operator %= (long);
+ void operator <<=(long);
+ void operator >>=(long);
+ void operator &= (long);
+ void operator |= (long);
+ void operator ^= (long);
+
+// (constructive binary operations are inlined below)
+
+#ifdef __GNUG__
+ //friend Integer operator <? (const Integer& x, const Integer& y); // min
+ //friend Integer operator >? (const Integer& x, const Integer& y); // max
+#endif
+
+// builtin Integer functions that must be friends
+
+ friend long lg (const Integer&); // floor log base 2 of abs(x)
+ friend double ratio(const Integer& x, const Integer& y);
+ // return x/y as a double
+
+ friend Integer gcd(const Integer&, const Integer&);
+ friend int even(const Integer&); // true if even
+ friend int odd(const Integer&); // true if odd
+ friend int sign(const Integer&); // returns -1, 0, +1
+
+ friend void (setbit)(Integer& x, long b); // set b'th bit of x
+ friend void clearbit(Integer& x, long b); // clear b'th bit
+ friend int testbit(const Integer& x, long b); // return b'th bit
+
+// procedural versions of operators
+
+ friend void abs(const Integer& x, Integer& dest);
+ friend void negate(const Integer& x, Integer& dest);
+ friend void complement(const Integer& x, Integer& dest);
+
+ friend int compare(const Integer&, const Integer&);
+ friend int ucompare(const Integer&, const Integer&);
+ friend void add(const Integer& x, const Integer& y, Integer& dest);
+ friend void sub(const Integer& x, const Integer& y, Integer& dest);
+ friend void mul(const Integer& x, const Integer& y, Integer& dest);
+ friend void div(const Integer& x, const Integer& y, Integer& dest);
+ friend void mod(const Integer& x, const Integer& y, Integer& dest);
+ friend void divide(const Integer& x, const Integer& y,
+ Integer& q, Integer& r);
+ friend void I_and(const Integer& x, const Integer& y, Integer& dest);
+ friend void I_or(const Integer& x, const Integer& y, Integer& dest);
+ friend void I_xor(const Integer& x, const Integer& y, Integer& dest);
+ friend void lshift(const Integer& x, const Integer& y, Integer& dest);
+ friend void rshift(const Integer& x, const Integer& y, Integer& dest);
+ friend void pow(const Integer& x, const Integer& y, Integer& dest);
+
+ friend int compare(const Integer&, long);
+ friend int ucompare(const Integer&, long);
+ friend void add(const Integer& x, long y, Integer& dest);
+ friend void sub(const Integer& x, long y, Integer& dest);
+ friend void mul(const Integer& x, long y, Integer& dest);
+ friend void div(const Integer& x, long y, Integer& dest);
+ friend void mod(const Integer& x, long y, Integer& dest);
+ friend void divide(const Integer& x, long y, Integer& q, long& r);
+ friend void I_and(const Integer& x, long y, Integer& dest);
+ friend void I_or(const Integer& x, long y, Integer& dest);
+ friend void I_xor(const Integer& x, long y, Integer& dest);
+ friend void lshift(const Integer& x, long y, Integer& dest);
+ friend void rshift(const Integer& x, long y, Integer& dest);
+ friend void pow(const Integer& x, long y, Integer& dest);
+
+ friend int compare(long, const Integer&);
+ friend int ucompare(long, const Integer&);
+ friend void add(long x, const Integer& y, Integer& dest);
+ friend void sub(long x, const Integer& y, Integer& dest);
+ friend void mul(long x, const Integer& y, Integer& dest);
+ friend void I_and(long x, const Integer& y, Integer& dest);
+ friend void I_or(long x, const Integer& y, Integer& dest);
+ friend void I_xor(long x, const Integer& y, Integer& dest);
+
+// coercion & conversion
+
+ int fits_in_long() const;
+ int fits_in_double() const;
+
+ operator long() const;
+ operator double() const;
+
+ friend char* Itoa(const Integer& x, int base = 10, int width = 0);
+ friend Integer atoI(const char* s, int base = 10);
+ void printon(std::ostream& s, int base = 10, int width = 0) const;
+
+ friend std::istream& operator >> (std::istream& s, Integer& y);
+ friend std::ostream& operator << (std::ostream& s, const Integer& y);
+
+// error detection
+
+ int initialized() const;
+ void error(const char* msg) const;
+ int OK() const;
+};
+
+
+// (These are declared inline)
+
+ int operator == (const Integer&, const Integer&);
+ int operator == (const Integer&, long);
+ int operator != (const Integer&, const Integer&);
+ int operator != (const Integer&, long);
+ int operator < (const Integer&, const Integer&);
+ int operator < (const Integer&, long);
+ int operator <= (const Integer&, const Integer&);
+ int operator <= (const Integer&, long);
+ int operator > (const Integer&, const Integer&);
+ int operator > (const Integer&, long);
+ int operator >= (const Integer&, const Integer&);
+ int operator >= (const Integer&, long);
+ Integer operator - (const Integer&);
+ Integer operator ~ (const Integer&);
+ Integer operator + (const Integer&, const Integer&);
+ Integer operator + (const Integer&, long);
+ Integer operator + (long, const Integer&);
+ Integer operator - (const Integer&, const Integer&);
+ Integer operator - (const Integer&, long);
+ Integer operator - (long, const Integer&);
+ Integer operator * (const Integer&, const Integer&);
+ Integer operator * (const Integer&, long);
+ Integer operator * (long, const Integer&);
+ Integer operator / (const Integer&, const Integer&);
+ Integer operator / (const Integer&, long);
+ Integer operator % (const Integer&, const Integer&);
+ Integer operator % (const Integer&, long);
+ Integer operator << (const Integer&, const Integer&);
+ Integer operator << (const Integer&, long);
+ Integer operator >> (const Integer&, const Integer&);
+ Integer operator >> (const Integer&, long);
+ Integer operator & (const Integer&, const Integer&);
+ Integer operator & (const Integer&, long);
+ Integer operator & (long, const Integer&);
+ Integer operator | (const Integer&, const Integer&);
+ Integer operator | (const Integer&, long);
+ Integer operator | (long, const Integer&);
+ Integer operator ^ (const Integer&, const Integer&);
+ Integer operator ^ (const Integer&, long);
+ Integer operator ^ (long, const Integer&);
+
+ Integer abs(const Integer&); // absolute value
+ Integer sqr(const Integer&); // square
+
+ Integer pow(const Integer& x, const Integer& y);
+ Integer pow(const Integer& x, long y);
+ Integer Ipow(long x, long y); // x to the y as Integer
+
+
+extern char* dec(const Integer& x, int width = 0);
+extern char* oct(const Integer& x, int width = 0);
+extern char* hex(const Integer& x, int width = 0);
+extern Integer sqrt(const Integer&); // floor of square root
+extern Integer lcm(const Integer& x, const Integer& y); // least common mult
+
+
+typedef Integer IntTmp; // for backward compatibility
+
+inline Integer::Integer() :rep(&_ZeroRep) {}
+
+inline Integer::Integer(IntRep* r) :rep(r) {}
+
+inline Integer::Integer(long y) :rep(Icopy_long(0, y)) {}
+
+inline Integer::Integer(unsigned long y) :rep(Icopy_ulong(0, y)) {}
+
+inline Integer::Integer(const Integer& y) :rep(Icopy(0, y.rep)) {}
+
+inline Integer::~Integer() { if (rep && !STATIC_IntRep(rep)) delete rep; }
+
+inline void Integer::operator = (const Integer& y)
+{
+ rep = Icopy(rep, y.rep);
+}
+
+inline void Integer::operator = (long y)
+{
+ rep = Icopy_long(rep, y);
+}
+
+inline Integer::operator long() const
+{
+ return Itolong(rep);
+}
+
+inline int Integer::initialized() const
+{
+ return rep != 0;
+}
+
+inline int Integer::fits_in_long() const
+{
+ return Iislong(rep);
+}
+
+inline Integer::operator double() const
+{
+ return Itodouble(rep);
+}
+
+inline int Integer::fits_in_double() const
+{
+ return Iisdouble(rep);
+}
+
+// procedural versions
+
+inline int compare(const Integer& x, const Integer& y)
+{
+ return compare(x.rep, y.rep);
+}
+
+inline int ucompare(const Integer& x, const Integer& y)
+{
+ return ucompare(x.rep, y.rep);
+}
+
+inline int compare(const Integer& x, long y)
+{
+ return compare(x.rep, y);
+}
+
+inline int ucompare(const Integer& x, long y)
+{
+ return ucompare(x.rep, y);
+}
+
+inline int compare(long x, const Integer& y)
+{
+ return -compare(y.rep, x);
+}
+
+inline int ucompare(long x, const Integer& y)
+{
+ return -ucompare(y.rep, x);
+}
+
+inline void add(const Integer& x, const Integer& y, Integer& dest)
+{
+ dest.rep = add(x.rep, 0, y.rep, 0, dest.rep);
+}
+
+inline void sub(const Integer& x, const Integer& y, Integer& dest)
+{
+ dest.rep = add(x.rep, 0, y.rep, 1, dest.rep);
+}
+
+inline void mul(const Integer& x, const Integer& y, Integer& dest)
+{
+ dest.rep = multiply(x.rep, y.rep, dest.rep);
+}
+
+inline void div(const Integer& x, const Integer& y, Integer& dest)
+{
+ dest.rep = div(x.rep, y.rep, dest.rep);
+}
+
+inline void mod(const Integer& x, const Integer& y, Integer& dest)
+{
+ dest.rep = mod(x.rep, y.rep, dest.rep);
+}
+
+inline void I_and(const Integer& x, const Integer& y, Integer& dest)
+{
+ dest.rep = bitop(x.rep, y.rep, dest.rep, '&');
+}
+
+inline void I_or(const Integer& x, const Integer& y, Integer& dest)
+{
+ dest.rep = bitop(x.rep, y.rep, dest.rep, '|');
+}
+
+inline void I_xor(const Integer& x, const Integer& y, Integer& dest)
+{
+ dest.rep = bitop(x.rep, y.rep, dest.rep, '^');
+}
+
+inline void lshift(const Integer& x, const Integer& y, Integer& dest)
+{
+ dest.rep = lshift(x.rep, y.rep, 0, dest.rep);
+}
+
+inline void rshift(const Integer& x, const Integer& y, Integer& dest)
+{
+ dest.rep = lshift(x.rep, y.rep, 1, dest.rep);
+}
+
+inline void pow(const Integer& x, const Integer& y, Integer& dest)
+{
+ dest.rep = power(x.rep, long(y), dest.rep); // not incorrect
+}
+
+inline void add(const Integer& x, long y, Integer& dest)
+{
+ dest.rep = add(x.rep, 0, y, dest.rep);
+}
+
+inline void sub(const Integer& x, long y, Integer& dest)
+{
+ dest.rep = add(x.rep, 0, -y, dest.rep);
+}
+
+inline void mul(const Integer& x, long y, Integer& dest)
+{
+ dest.rep = multiply(x.rep, y, dest.rep);
+}
+
+inline void div(const Integer& x, long y, Integer& dest)
+{
+ dest.rep = div(x.rep, y, dest.rep);
+}
+
+inline void mod(const Integer& x, long y, Integer& dest)
+{
+ dest.rep = mod(x.rep, y, dest.rep);
+}
+
+inline void I_and(const Integer& x, long y, Integer& dest)
+{
+ dest.rep = bitop(x.rep, y, dest.rep, '&');
+}
+
+inline void I_or(const Integer& x, long y, Integer& dest)
+{
+ dest.rep = bitop(x.rep, y, dest.rep, '|');
+}
+
+inline void I_xor(const Integer& x, long y, Integer& dest)
+{
+ dest.rep = bitop(x.rep, y, dest.rep, '^');
+}
+
+inline void lshift(const Integer& x, long y, Integer& dest)
+{
+ dest.rep = lshift(x.rep, y, dest.rep);
+}
+
+inline void rshift(const Integer& x, long y, Integer& dest)
+{
+ dest.rep = lshift(x.rep, -y, dest.rep);
+}
+
+inline void pow(const Integer& x, long y, Integer& dest)
+{
+ dest.rep = power(x.rep, y, dest.rep);
+}
+
+inline void abs(const Integer& x, Integer& dest)
+{
+ dest.rep = abs(x.rep, dest.rep);
+}
+
+inline void negate(const Integer& x, Integer& dest)
+{
+ dest.rep = negate(x.rep, dest.rep);
+}
+
+inline void complement(const Integer& x, Integer& dest)
+{
+ dest.rep = I_compl(x.rep, dest.rep);
+}
+
+inline void add(long x, const Integer& y, Integer& dest)
+{
+ dest.rep = add(y.rep, 0, x, dest.rep);
+}
+
+inline void sub(long x, const Integer& y, Integer& dest)
+{
+ dest.rep = add(y.rep, 1, x, dest.rep);
+}
+
+inline void mul(long x, const Integer& y, Integer& dest)
+{
+ dest.rep = multiply(y.rep, x, dest.rep);
+}
+
+inline void I_and(long x, const Integer& y, Integer& dest)
+{
+ dest.rep = bitop(y.rep, x, dest.rep, '&');
+}
+
+inline void I_or(long x, const Integer& y, Integer& dest)
+{
+ dest.rep = bitop(y.rep, x, dest.rep, '|');
+}
+
+inline void I_xor(long x, const Integer& y, Integer& dest)
+{
+ dest.rep = bitop(y.rep, x, dest.rep, '^');
+}
+
+
+// operator versions
+
+inline int operator == (const Integer& x, const Integer& y)
+{
+ return compare(x, y) == 0;
+}
+
+inline int operator == (const Integer& x, long y)
+{
+ return compare(x, y) == 0;
+}
+
+inline int operator != (const Integer& x, const Integer& y)
+{
+ return compare(x, y) != 0;
+}
+
+inline int operator != (const Integer& x, long y)
+{
+ return compare(x, y) != 0;
+}
+
+inline int operator < (const Integer& x, const Integer& y)
+{
+ return compare(x, y) < 0;
+}
+
+inline int operator < (const Integer& x, long y)
+{
+ return compare(x, y) < 0;
+}
+
+inline int operator <= (const Integer& x, const Integer& y)
+{
+ return compare(x, y) <= 0;
+}
+
+inline int operator <= (const Integer& x, long y)
+{
+ return compare(x, y) <= 0;
+}
+
+inline int operator > (const Integer& x, const Integer& y)
+{
+ return compare(x, y) > 0;
+}
+
+inline int operator > (const Integer& x, long y)
+{
+ return compare(x, y) > 0;
+}
+
+inline int operator >= (const Integer& x, const Integer& y)
+{
+ return compare(x, y) >= 0;
+}
+
+inline int operator >= (const Integer& x, long y)
+{
+ return compare(x, y) >= 0;
+}
+
+
+inline void Integer::operator += (const Integer& y)
+{
+ add(*this, y, *this);
+}
+
+inline void Integer::operator += (long y)
+{
+ add(*this, y, *this);
+}
+
+inline void Integer::operator ++ ()
+{
+ add(*this, 1, *this);
+}
+
+
+inline void Integer::operator -= (const Integer& y)
+{
+ sub(*this, y, *this);
+}
+
+inline void Integer::operator -= (long y)
+{
+ sub(*this, y, *this);
+}
+
+inline void Integer::operator -- ()
+{
+ add(*this, -1, *this);
+}
+
+
+
+inline void Integer::operator *= (const Integer& y)
+{
+ mul(*this, y, *this);
+}
+
+inline void Integer::operator *= (long y)
+{
+ mul(*this, y, *this);
+}
+
+
+inline void Integer::operator &= (const Integer& y)
+{
+ I_and(*this, y, *this);
+}
+
+inline void Integer::operator &= (long y)
+{
+ I_and(*this, y, *this);
+}
+
+inline void Integer::operator |= (const Integer& y)
+{
+ I_or(*this, y, *this);
+}
+
+inline void Integer::operator |= (long y)
+{
+ I_or(*this, y, *this);
+}
+
+
+inline void Integer::operator ^= (const Integer& y)
+{
+ I_xor(*this, y, *this);
+}
+
+inline void Integer::operator ^= (long y)
+{
+ I_xor(*this, y, *this);
+}
+
+
+
+inline void Integer::operator /= (const Integer& y)
+{
+ div(*this, y, *this);
+}
+
+inline void Integer::operator /= (long y)
+{
+ div(*this, y, *this);
+}
+
+
+inline void Integer::operator %= (const Integer& y)
+{
+ *this = *this % y; // mod(*this, y, *this) doesn't work.
+}
+
+inline void Integer::operator %= (long y)
+{
+ *this = *this % y; // mod(*this, y, *this) doesn't work.
+}
+
+
+inline void Integer::operator <<= (const Integer& y)
+{
+ lshift(*this, y, *this);
+}
+
+inline void Integer::operator <<= (long y)
+{
+ lshift(*this, y, *this);
+}
+
+
+inline void Integer::operator >>= (const Integer& y)
+{
+ rshift(*this, y, *this);
+}
+
+inline void Integer::operator >>= (long y)
+{
+ rshift(*this, y, *this);
+}
+
+#ifdef __GNUG__
+/*
+inline Integer operator <? (const Integer& x, const Integer& y)
+{
+ return (compare(x.rep, y.rep) <= 0) ? x : y;
+}
+
+inline Integer operator >? (const Integer& x, const Integer& y)
+{
+ return (compare(x.rep, y.rep) >= 0)? x : y;
+}
+*/
+#endif
+
+
+inline void Integer::abs()
+{
+ ::abs(*this, *this);
+}
+
+inline void Integer::negate()
+{
+ ::negate(*this, *this);
+}
+
+
+inline void Integer::complement()
+{
+ ::complement(*this, *this);
+}
+
+
+inline int sign(const Integer& x)
+{
+ return (x.rep->len == 0) ? 0 : ( (x.rep->sgn == 1) ? 1 : -1 );
+}
+
+inline int even(const Integer& y)
+{
+ return y.rep->len == 0 || !(y.rep->s[0] & 1);
+}
+
+inline int odd(const Integer& y)
+{
+ return y.rep->len > 0 && (y.rep->s[0] & 1);
+}
+
+inline char* Itoa(const Integer& y, int base, int width)
+{
+ return Itoa(y.rep, base, width);
+}
+
+
+
+inline long lg(const Integer& x)
+{
+ return lg(x.rep);
+}
+
+// constructive operations
+
+inline Integer operator + (const Integer& x, const Integer& y)
+{
+ Integer r; add(x, y, r); return r;
+}
+
+inline Integer operator + (const Integer& x, long y)
+{
+ Integer r; add(x, y, r); return r;
+}
+
+inline Integer operator + (long x, const Integer& y)
+{
+ Integer r; add(x, y, r); return r;
+}
+
+inline Integer operator - (const Integer& x, const Integer& y)
+{
+ Integer r; sub(x, y, r); return r;
+}
+
+inline Integer operator - (const Integer& x, long y)
+{
+ Integer r; sub(x, y, r); return r;
+}
+
+inline Integer operator - (long x, const Integer& y)
+{
+ Integer r; sub(x, y, r); return r;
+}
+
+inline Integer operator * (const Integer& x, const Integer& y)
+{
+ Integer r; mul(x, y, r); return r;
+}
+
+inline Integer operator * (const Integer& x, long y)
+{
+ Integer r; mul(x, y, r); return r;
+}
+
+inline Integer operator * (long x, const Integer& y)
+{
+ Integer r; mul(x, y, r); return r;
+}
+
+inline Integer sqr(const Integer& x)
+{
+ Integer r; mul(x, x, r); return r;
+}
+
+inline Integer operator & (const Integer& x, const Integer& y)
+{
+ Integer r; I_and(x, y, r); return r;
+}
+
+inline Integer operator & (const Integer& x, long y)
+{
+ Integer r; I_and(x, y, r); return r;
+}
+
+inline Integer operator & (long x, const Integer& y)
+{
+ Integer r; I_and(x, y, r); return r;
+}
+
+inline Integer operator | (const Integer& x, const Integer& y)
+{
+ Integer r; I_or(x, y, r); return r;
+}
+
+inline Integer operator | (const Integer& x, long y)
+{
+ Integer r; I_or(x, y, r); return r;
+}
+
+inline Integer operator | (long x, const Integer& y)
+{
+ Integer r; I_or(x, y, r); return r;
+}
+
+inline Integer operator ^ (const Integer& x, const Integer& y)
+{
+ Integer r; I_xor(x, y, r); return r;
+}
+
+inline Integer operator ^ (const Integer& x, long y)
+{
+ Integer r; I_xor(x, y, r); return r;
+}
+
+inline Integer operator ^ (long x, const Integer& y)
+{
+ Integer r; I_xor(x, y, r); return r;
+}
+
+inline Integer operator / (const Integer& x, const Integer& y)
+{
+ Integer r; div(x, y, r); return r;
+}
+
+inline Integer operator / (const Integer& x, long y)
+{
+ Integer r; div(x, y, r); return r;
+}
+
+inline Integer operator % (const Integer& x, const Integer& y)
+{
+ Integer r; mod(x, y, r); return r;
+}
+
+inline Integer operator % (const Integer& x, long y)
+{
+ Integer r; mod(x, y, r); return r;
+}
+
+inline Integer operator << (const Integer& x, const Integer& y)
+{
+ Integer r; lshift(x, y, r); return r;
+}
+
+inline Integer operator << (const Integer& x, long y)
+{
+ Integer r; lshift(x, y, r); return r;
+}
+
+inline Integer operator >> (const Integer& x, const Integer& y)
+{
+ Integer r; rshift(x, y, r); return r;
+}
+
+inline Integer operator >> (const Integer& x, long y)
+{
+ Integer r; rshift(x, y, r); return r;
+}
+
+inline Integer pow(const Integer& x, long y)
+{
+ Integer r; pow(x, y, r); return r;
+}
+
+inline Integer Ipow(long x, long y)
+{
+ Integer r(x); pow(r, y, r); return r;
+}
+
+inline Integer pow(const Integer& x, const Integer& y)
+{
+ Integer r; pow(x, y, r); return r;
+}
+
+
+
+inline Integer abs(const Integer& x)
+{
+ Integer r; abs(x, r); return r;
+}
+
+inline Integer operator - (const Integer& x)
+{
+ Integer r; negate(x, r); return r;
+}
+
+inline Integer operator ~ (const Integer& x)
+{
+ Integer r; complement(x, r); return r;
+}
+
+inline Integer atoI(const char* s, int base)
+{
+ Integer r; r.rep = atoIntRep(s, base); return r;
+}
+
+inline Integer gcd(const Integer& x, const Integer& y)
+{
+ Integer r; r.rep = gcd(x.rep, y.rep); return r;
+}
+
+#endif /* !_Integer_h */
projects / ACL.git / commitdiff