/* 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 #include #include #include #include #include #include #include #include /* 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 SHORTBITS #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)(((LONGBITS + SHORTBITS - 1) / SHORTBITS))) #define CHAR_PER_LONG ((unsigned)(((LONGBITS + CHARBITS - 1) / CHARBITS))) /* 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 // 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) delete old; return rep; } // allocate and clear IntRep* Icalloc(IntRep* old, int newlen) { IntRep* rep; if (old == 0 || newlen > old->sz) { if (old != 0) 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; 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) 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) { unsigned short src[SHORT_PER_LONG]; unsigned long u; int newsgn; if (newsgn = (x >= 0)) u = x; else u = -x; unsigned short srclen = 0; while (u != 0) { src[srclen++] = extract(u); u = down(u); } IntRep* rep; if (old == 0 || srclen > old->sz) { if (old != 0) delete old; rep = Inew(srclen); } else rep = old; rep->len = srclen; rep->sgn = newsgn; scpy(src, rep->s, srclen); return rep; } // special case for zero -- it's worth it! IntRep* Icopy_zero(IntRep* old) { IntRep* rep; if (old == 0) rep = Inew(0); else rep = old; rep->len = 0; rep->sgn = I_POSITIVE; return rep; } // special case for 1 or -1 IntRep* Icopy_one(IntRep* old, int newsgn) { IntRep* rep; if (old == 0 || 1 > old->sz) { if (old != 0) delete old; rep = Inew(1); } else rep = old; rep->sgn = newsgn; rep->len = 1; rep->s[0] = 1; return rep; } // 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) ? MAXLONG : MINLONG; 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) ? MAXLONG : MINLONG; 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 MINLONG & MAXLONG 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 = HUGE / 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 : HUGE; 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 = HUGE / 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 == HUGE || d1 == -HUGE || 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) delete yy; 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); 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); 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) 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) 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* 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 (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); } } delete t; 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); } 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; } #if defined(__GNUG__) && !defined(NO_NRV) Integer sqrt(const Integer& x) return 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; } Integer lcm(const Integer& x, const Integer& y) return 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); } #else 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; } #endif 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 = (char *) _libgxx_fmtq.alloc(fmtlen); char* f = cvtItoa(x, fmtbase, fmtlen, base, 0, width, 0, ' ', 'X', 0); return f; } ostream& operator << (ostream& s, const Integer& y) { #ifdef _OLD_STREAMS return s << Itoa(y.rep); #else int base = (s.flags() & ios::oct) ? 8 : (s.flags() & ios::hex) ? 16 : 10; int align_right = !(s.flags() & ios::left); int showpos = s.flags() & ios::showpos; int showbase = s.flags() & ios::showbase; char fillchar = s.fill(); char Xcase = (s.flags() & ios::uppercase)? 'X' : 'x'; int width = s.width(); const IntRep* x = y.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; return s; #endif } 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; } 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); } istream& operator >> (istream& s, Integer& y) { #ifdef _OLD_STREAMS if (!s.good()) return s; #else if (!s.ipfx(0)) { s.set(ios::failbit); return s; } #endif s >> ws; if (!s.good()) { s.set(_fail); return s; } #ifdef _OLD_STREAMS int know_base = 0; int base = 10; #else int know_base = s.flags() & (ios::oct | ios::hex | ios::dec); int base = (s.flags() & ios::oct) ? 8 : (s.flags() & 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; 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.set(_fail); if (sgn == '-') y.negate(); return s; } int Integer::OK() const { int v = rep != 0; // have a rep int l = rep->len; int s = rep->sgn; v &= l <= rep->sz; // length with bounds v &= s == 0 || s == 1; // legal sign Icheck(rep); // and correctly adjusted v &= rep->len == l; v &= rep->sgn == s; if (!v) error("invariant failure"); return v; } void Integer::error(const char* msg) const { (*lib_error_handler)("Integer", msg); }