From fef0e69336ec9652c5bbf1a58c1b08c5922ae49b Mon Sep 17 00:00:00 2001 From: Arjen Baart Date: Sat, 26 Oct 2019 14:56:13 +0200 Subject: [PATCH] Ported class Integer from Gnu libg++ --- TODO | 2 +- src/Integer.cpp | 2374 +++++++++++++++++++++++++++++++++++++++++++++++++++++++ src/Integer.h | 936 ++++++++++++++++++++++ src/Makefile.am | 5 +- src/builtin.h | 144 ++++ 5 files changed, 3458 insertions(+), 3 deletions(-) create mode 100644 src/Integer.cpp create mode 100644 src/Integer.h create mode 100644 src/builtin.h diff --git a/TODO b/TODO index 7ac6bdc..fea13d1 100644 --- a/TODO +++ b/TODO @@ -10,6 +10,6 @@ Things to do: - date: Parser for stream input - date: Weekday and week number -- Include classes: UTC, complex, integer - UTC: Convert to and from time_t, struct tm +- Add classes: configuration, Integer, amount diff --git a/src/Integer.cpp b/src/Integer.cpp new file mode 100644 index 0000000..1e69029 --- /dev/null +++ b/src/Integer.cpp @@ -0,0 +1,2374 @@ +/* +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 +#include + +#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); +} + diff --git a/src/Integer.h b/src/Integer.h new file mode 100644 index 0000000..f004732 --- /dev/null +++ b/src/Integer.h @@ -0,0 +1,936 @@ +// 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 + +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); // 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; +} +*/ +#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 */ diff --git a/src/Makefile.am b/src/Makefile.am index 7e1d3cc..c9b1cf2 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -1,6 +1,7 @@ lib_LTLIBRARIES = libACL.la libACL_la_SOURCES = string.cpp regex.cpp date.cpp parsedate.c dateyacc.y datelex.c \ - hour.cpp utc.cpp + hour.cpp utc.cpp \ + Integer.cpp -include_HEADERS = String.h date.h +include_HEADERS = String.h date.h Integer.h diff --git a/src/builtin.h b/src/builtin.h new file mode 100644 index 0000000..d6a41a4 --- /dev/null +++ b/src/builtin.h @@ -0,0 +1,144 @@ +// This may look like C code, but it is really -*- C++ -*- + +/* +Copyright (C) 1988, 1992 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. +*/ + +/* + arithmetic, etc. functions on built in types +*/ + + +#ifndef _builtin_h +#ifdef __GNUG__ +#pragma interface +#endif +#define _builtin_h 1 + +#include <_G_config.h> +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#ifdef __GNUG__ +#define _VOLATILE_VOID volatile void +#else +#define _VOLATILE_VOID void +#endif + +typedef void (*one_arg_error_handler_t)(const char*); +typedef void (*two_arg_error_handler_t)(const char*, const char*); + +long gcd(long, long); +long lg(unsigned long); +double pow(double, long); +long pow(long, long); + +double start_timer(); +double return_elapsed_time(double last_time = 0.0); + +char* dtoa(double x, char cvt = 'g', int width = 0, int prec = 6); + +char* chr(char ch, int width = 0); +char* str(const char* s, int width = 0); + +unsigned int hashpjw(const char*); +unsigned int multiplicativehash(int); +unsigned int foldhash(double); + +extern _VOLATILE_VOID default_one_arg_error_handler(const char*); +extern _VOLATILE_VOID default_two_arg_error_handler(const char*, const char*); + +extern two_arg_error_handler_t lib_error_handler; + +extern two_arg_error_handler_t + set_lib_error_handler(two_arg_error_handler_t f); + + +int sign(long arg); +int sign(double arg); +long sqr(long arg); +double sqr(double arg); +int even(long arg); +int odd(long arg); +long lcm(long x, long y); +void (setbit)(long& x, long b); +void clearbit(long& x, long b); +int testbit(long x, long b); + +#if !defined(IV) + +inline int sign(long arg) +{ + return (arg == 0) ? 0 : ( (arg > 0) ? 1 : -1 ); +} + +inline int sign(double arg) +{ + return (arg == 0.0) ? 0 : ( (arg > 0.0) ? 1 : -1 ); +} + +inline long sqr(long arg) +{ + return arg * arg; +} + +#ifndef __hpux /* hpux defines this in math.h */ +inline double sqr(double arg) +{ + return arg * arg; +} +#endif + +inline int even(long arg) +{ + return !(arg & 1); +} + +inline int odd(long arg) +{ + return (arg & 1); +} + +inline long lcm(long x, long y) +{ + return x / gcd(x, y) * y; +} + +inline void (setbit)(long& x, long b) +{ + x |= (1 << b); +} + +inline void clearbit(long& x, long b) +{ + x &= ~(1 << b); +} + +inline int testbit(long x, long b) +{ + return ((x & (1 << b)) != 0); +} + +#endif +#endif -- 2.11.0