X-Git-Url: https://git.distorted.org.uk/u/mdw/catacomb/blobdiff_plain/ba6e6b64033b1f9de49feccb5c9cd438354481f7..0f00dc4c8eb47e67bc0f148c2dd109f73a451e0a:/math/mptext.c diff --git a/math/mptext.c b/math/mptext.c new file mode 100644 index 0000000..5d4640e --- /dev/null +++ b/math/mptext.c @@ -0,0 +1,849 @@ +/* -*-c-*- + * + * Textual representation of multiprecision numbers + * + * (c) 1999 Straylight/Edgeware + */ + +/*----- Licensing notice --------------------------------------------------* + * + * This file is part of Catacomb. + * + * Catacomb 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. + * + * Catacomb 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 Catacomb; if not, write to the Free + * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, + * MA 02111-1307, USA. + */ + +/*----- Header files ------------------------------------------------------*/ + +#include +#include +#include + +#include "mp.h" +#include "mptext.h" +#include "paranoia.h" + +/*----- Magical numbers ---------------------------------------------------*/ + +/* --- Maximum recursion depth --- * + * + * This is the number of bits in a @size_t@ object. Why? + * + * To see this, let %$b = \textit{MPW\_MAX} + 1$% and let %$Z$% be the + * largest @size_t@ value. Then the largest possible @mp@ is %$M - 1$% where + * %$M = b^Z$%. Let %$r$% be a radix to read or write. Since the recursion + * squares the radix at each step, the highest number reached by the + * recursion is %$d$%, where: + * + * %$r^{2^d} = b^Z$%. + * + * Solving gives that %$d = \lg \log_r b^Z$%. If %$r = 2$%, this is maximum, + * so choosing %$d = \lg \lg b^Z = \lg (Z \lg b) = \lg Z + \lg \lg b$%. + * + * Expressing %$\lg Z$% as @CHAR_BIT * sizeof(size_t)@ yields an + * overestimate, since a @size_t@ representation may contain `holes'. + * Choosing to represent %$\lg \lg b$% by 10 is almost certainly sufficient + * for `some time to come'. + */ + +#define DEPTH (CHAR_BIT * sizeof(size_t) + 10) + +/*----- Main code ---------------------------------------------------------*/ + +/* --- @mp_read@ --- * + * + * Arguments: @mp *m@ = destination multiprecision number + * @int radix@ = base to assume for data (or zero to guess) + * @const mptext_ops *ops@ = pointer to operations block + * @void *p@ = data for the operations block + * + * Returns: The integer read, or zero if it didn't work. + * + * Use: Reads an integer from some source. If the @radix@ is + * specified, the number is assumed to be given in that radix, + * with the letters `a' (either upper- or lower-case) upwards + * standing for digits greater than 9. Otherwise, base 10 is + * assumed unless the number starts with `0' (octal), `0x' (hex) + * or `nnn_' (base `nnn'). An arbitrary amount of whitespace + * before the number is ignored. + */ + +/* --- About the algorithm --- * + * + * The algorithm here is rather aggressive. I maintain an array of + * successive squarings of the radix, and a stack of partial results, each + * with a counter attached indicating which radix square to multiply by. + * Once the item at the top of the stack reaches the same counter level as + * the next item down, they are combined together and the result is given a + * counter level one higher than either of the results. + * + * Gluing the results together at the end is slightly tricky. Pay attention + * to the code. + * + * This is more complicated because of the need to handle the slightly + * bizarre syntax. + */ + +mp *mp_read(mp *m, int radix, const mptext_ops *ops, void *p) +{ + int ch; /* Current char being considered */ + unsigned f = 0; /* Flags about the current number */ + int r; /* Radix to switch over to */ + mpw rd; /* Radix as an @mp@ digit */ + mp rr; /* The @mp@ for the radix */ + unsigned nf = m ? m->f & MP_BURN : 0; /* New @mp@ flags */ + + /* --- Stacks --- */ + + mp *pow[DEPTH]; /* List of powers */ + unsigned pows; /* Next index to fill */ + struct { unsigned i; mp *m; } s[DEPTH]; /* Main stack */ + unsigned sp; /* Current stack pointer */ + + /* --- Flags --- */ + +#define f_neg 1u +#define f_ok 2u +#define f_start 4u + + /* --- Initialize the stacks --- */ + + mp_build(&rr, &rd, &rd + 1); + pow[0] = &rr; + pows = 1; + + sp = 0; + + /* --- Initialize the destination number --- */ + + if (m) + MP_DROP(m); + + /* --- Read an initial character --- */ + + ch = ops->get(p); + if (radix >= 0) { + while (isspace(ch)) + ch = ops->get(p); + } + + /* --- Handle an initial sign --- */ + + if (radix >= 0 && (ch == '-' || ch == '+')) { + if (ch == '-') + f |= f_neg; + do ch = ops->get(p); while isspace(ch); + } + + /* --- If the radix is zero, look for leading zeros --- */ + + if (radix > 0) { + assert(((void)"ascii radix must be <= 62", radix <= 62)); + rd = radix; + r = -1; + } else if (radix < 0) { + rd = -radix; + assert(((void)"binary radix must fit in a byte", rd <= UCHAR_MAX)); + r = -1; + } else if (ch != '0') { + rd = 10; + r = 0; + } else { + ch = ops->get(p); + switch (ch) { + case 'x': + rd = 16; + goto prefix; + case 'o': + rd = 8; + goto prefix; + case 'b': + rd = 2; + goto prefix; + prefix: + ch = ops->get(p); + break; + default: + rd = 8; + f |= f_ok; + } + r = -1; + } + + /* --- Use fast algorithm for binary radix --- * + * + * This is the restart point after having parsed a radix number from the + * input. We check whether the radix is binary, and if so use a fast + * algorithm which just stacks the bits up in the right order. + */ + +restart: + switch (rd) { + unsigned bit; + + case 2: bit = 1; goto bin; + case 4: bit = 2; goto bin; + case 8: bit = 3; goto bin; + case 16: bit = 4; goto bin; + case 32: bit = 5; goto bin; + case 64: bit = 6; goto bin; + case 128: bit = 7; goto bin; + default: + break; + + /* --- The fast binary algorithm --- * + * + * We stack bits up starting at the top end of a word. When one word is + * full, we write it to the integer, and start another with the left-over + * bits. When the array in the integer is full, we resize using low-level + * calls and copy the current data to the top end. Finally, we do a single + * bit-shift when we know where the end of the number is. + */ + + bin: { + mpw a = 0; + unsigned b = MPW_BITS; + size_t len, n; + mpw *v; + + m = mp_dest(MP_NEW, 1, nf); + len = n = m->sz; + n = len; + v = m->v + n; + for (;; ch = ops->get(p)) { + unsigned x; + + if (ch < 0) + break; + + /* --- Check that the character is a digit and in range --- */ + + if (radix < 0) + x = ch % rd; + else { + if (!isalnum(ch)) + break; + if (ch >= '0' && ch <= '9') + x = ch - '0'; + else { + if (rd <= 36) + ch = tolower(ch); + if (ch >= 'a' && ch <= 'z') /* ASCII dependent! */ + x = ch - 'a' + 10; + else if (ch >= 'A' && ch <= 'Z') + x = ch - 'A' + 36; + else + break; + } + } + if (x >= rd) + break; + + /* --- Feed the digit into the accumulator --- */ + + f |= f_ok; + if (!x && !(f & f_start)) + continue; + f |= f_start; + if (b > bit) { + b -= bit; + a |= MPW(x) << b; + } else { + a |= MPW(x) >> (bit - b); + b += MPW_BITS - bit; + *--v = MPW(a); + n--; + if (!n) { + n = len; + len <<= 1; + v = mpalloc(m->a, len); + memcpy(v + n, m->v, MPWS(n)); + mpfree(m->a, m->v); + m->v = v; + v = m->v + n; + } + a = (b < MPW_BITS) ? MPW(x) << b : 0; + } + } + + /* --- Finish up --- */ + + if (!(f & f_ok)) { + mp_drop(m); + m = 0; + } else { + *--v = MPW(a); + n--; + m->sz = len; + m->vl = m->v + len; + m->f &= ~MP_UNDEF; + m = mp_lsr(m, m, (unsigned long)n * MPW_BITS + b); + } + ops->unget(ch, p); + goto done; + }} + + /* --- Time to start --- */ + + for (;; ch = ops->get(p)) { + unsigned x; + + if (ch < 0) + break; + + /* --- An underscore indicates a numbered base --- */ + + if (ch == '_' && r > 0 && r <= 62) { + unsigned i; + + /* --- Clear out the stacks --- */ + + for (i = 1; i < pows; i++) + MP_DROP(pow[i]); + pows = 1; + for (i = 0; i < sp; i++) + MP_DROP(s[i].m); + sp = 0; + + /* --- Restart the search --- */ + + rd = r; + r = -1; + f &= ~f_ok; + ch = ops->get(p); + goto restart; + } + + /* --- Check that the character is a digit and in range --- */ + + if (radix < 0) + x = ch % rd; + else { + if (!isalnum(ch)) + break; + if (ch >= '0' && ch <= '9') + x = ch - '0'; + else { + if (rd <= 36) + ch = tolower(ch); + if (ch >= 'a' && ch <= 'z') /* ASCII dependent! */ + x = ch - 'a' + 10; + else if (ch >= 'A' && ch <= 'Z') + x = ch - 'A' + 36; + else + break; + } + } + + /* --- Sort out what to do with the character --- */ + + if (x >= 10 && r >= 0) + r = -1; + if (x >= rd) + break; + + if (r >= 0) + r = r * 10 + x; + + /* --- Stick the character on the end of my integer --- */ + + assert(((void)"Number is too unimaginably huge", sp < DEPTH)); + s[sp].m = m = mp_new(1, nf); + m->v[0] = x; + s[sp].i = 0; + + /* --- Now grind through the stack --- */ + + while (sp > 0 && s[sp - 1].i == s[sp].i) { + + /* --- Combine the top two items --- */ + + sp--; + m = s[sp].m; + m = mp_mul(m, m, pow[s[sp].i]); + m = mp_add(m, m, s[sp + 1].m); + s[sp].m = m; + MP_DROP(s[sp + 1].m); + s[sp].i++; + + /* --- Make a new radix power if necessary --- */ + + if (s[sp].i >= pows) { + assert(((void)"Number is too unimaginably huge", pows < DEPTH)); + pow[pows] = mp_sqr(MP_NEW, pow[pows - 1]); + pows++; + } + } + f |= f_ok; + sp++; + } + + ops->unget(ch, p); + + /* --- If we're done, compute the rest of the number --- */ + + if (f & f_ok) { + if (!sp) + return (MP_ZERO); + else { + mp *z = MP_ONE; + sp--; + + while (sp > 0) { + + /* --- Combine the top two items --- */ + + sp--; + m = s[sp].m; + z = mp_mul(z, z, pow[s[sp + 1].i]); + m = mp_mul(m, m, z); + m = mp_add(m, m, s[sp + 1].m); + s[sp].m = m; + MP_DROP(s[sp + 1].m); + + /* --- Make a new radix power if necessary --- */ + + if (s[sp].i >= pows) { + assert(((void)"Number is too unimaginably huge", pows < DEPTH)); + pow[pows] = mp_sqr(MP_NEW, pow[pows - 1]); + pows++; + } + } + MP_DROP(z); + m = s[0].m; + } + } else { + unsigned i; + for (i = 0; i < sp; i++) + MP_DROP(s[i].m); + } + + /* --- Clear the radix power list --- */ + + { + unsigned i; + for (i = 1; i < pows; i++) + MP_DROP(pow[i]); + } + + /* --- Bail out if the number was bad --- */ + +done: + if (!(f & f_ok)) + return (0); + + /* --- Set the sign and return --- */ + + if (f & f_neg) + m->f |= MP_NEG; + MP_SHRINK(m); + return (m); + +#undef f_start +#undef f_neg +#undef f_ok +} + +/* --- @mp_write@ --- * + * + * Arguments: @mp *m@ = pointer to a multi-precision integer + * @int radix@ = radix to use when writing the number out + * @const mptext_ops *ops@ = pointer to an operations block + * @void *p@ = data for the operations block + * + * Returns: Zero if it worked, nonzero otherwise. + * + * Use: Writes a large integer in textual form. + */ + +/* --- Simple case --- * + * + * Use a fixed-sized buffer and single-precision arithmetic to pick off + * low-order digits. Put each digit in a buffer, working backwards from the + * end. If the buffer becomes full, recurse to get another one. Ensure that + * there are at least @z@ digits by writing leading zeroes if there aren't + * enough real digits. + */ + +static int simple(mpw n, int radix, unsigned z, + const mptext_ops *ops, void *p) +{ + int rc = 0; + char buf[64]; + unsigned i = sizeof(buf); + int rd = radix > 0 ? radix : -radix; + + do { + int ch; + mpw x; + + x = n % rd; + n /= rd; + if (radix < 0) + ch = x; + else if (x < 10) + ch = '0' + x; + else if (x < 36) /* Ascii specific */ + ch = 'a' + x - 10; + else + ch = 'A' + x - 36; + buf[--i] = ch; + if (z) + z--; + } while (i && n); + + if (n) + rc = simple(n, radix, z, ops, p); + else { + char zbuf[32]; + memset(zbuf, (radix < 0) ? 0 : '0', sizeof(zbuf)); + while (!rc && z >= sizeof(zbuf)) { + rc = ops->put(zbuf, sizeof(zbuf), p); + z -= sizeof(zbuf); + } + if (!rc && z) + rc = ops->put(zbuf, z, p); + } + if (!rc) + rc = ops->put(buf + i, sizeof(buf) - i, p); + BURN(buf); + return (rc); +} + +/* --- Complicated case --- * + * + * If the number is small, fall back to the simple case above. Otherwise + * divide and take remainder by current large power of the radix, and emit + * each separately. Don't emit a zero quotient. Be very careful about + * leading zeroes on the remainder part, because they're deeply significant. + */ + +static int complicated(mp *m, int radix, mp **pr, unsigned i, unsigned z, + const mptext_ops *ops, void *p) +{ + int rc = 0; + mp *q = MP_NEW; + unsigned d = 1 << i; + + if (MP_LEN(m) < 2) + return (simple(MP_LEN(m) ? m->v[0] : 0, radix, z, ops, p)); + + assert(i); + mp_div(&q, &m, m, pr[i]); + if (MP_ZEROP(q)) + d = z; + else { + if (z > d) + z -= d; + else + z = 0; + rc = complicated(q, radix, pr, i - 1, z, ops, p); + } + if (!rc) + rc = complicated(m, radix, pr, i - 1, d, ops, p); + mp_drop(q); + return (rc); +} + +/* --- Binary case --- * + * + * Special case for binary output. Goes much faster. + */ + +static int binary(mp *m, int bit, int radix, const mptext_ops *ops, void *p) +{ + mpw *v; + mpw a; + int rc = 0; + unsigned b; + unsigned mask; + unsigned long n; + unsigned f = 0; + char buf[8], *q; + unsigned x; + int ch; + +#define f_out 1u + + /* --- Work out where to start --- */ + + n = mp_bits(m); + if (n % bit) + n += bit - (n % bit); + b = n % MPW_BITS; + n /= MPW_BITS; + + if (n >= MP_LEN(m)) { + n--; + b += MPW_BITS; + } + + v = m->v + n; + a = *v; + mask = (1 << bit) - 1; + q = buf; + + /* --- Main code --- */ + + for (;;) { + if (b > bit) { + b -= bit; + x = a >> b; + } else { + x = a << (bit - b); + b += MPW_BITS - bit; + if (v == m->v) + break; + a = *--v; + if (b < MPW_BITS) + x |= a >> b; + } + x &= mask; + if (!x && !(f & f_out)) + continue; + + if (radix < 0) + ch = x; + else if (x < 10) + ch = '0' + x; + else if (x < 36) + ch = 'a' + x - 10; /* Ascii specific */ + else + ch = 'A' + x - 36; + *q++ = ch; + if (q >= buf + sizeof(buf)) { + if ((rc = ops->put(buf, sizeof(buf), p)) != 0) + goto done; + q = buf; + } + f |= f_out; + } + + x &= mask; + if (radix < 0) + ch = x; + else if (x < 10) + ch = '0' + x; + else if (x < 36) + ch = 'a' + x - 10; /* Ascii specific */ + else + ch = 'A' + x - 36; + *q++ = ch; + rc = ops->put(buf, q - buf, p); + +done: + mp_drop(m); + return (rc); + +#undef f_out +} + +/* --- Main driver code --- */ + +int mp_write(mp *m, int radix, const mptext_ops *ops, void *p) +{ + int rc; + + if (MP_EQ(m, MP_ZERO)) + return (ops->put(radix > 0 ? "0" : "\0", 1, p)); + + /* --- Set various things up --- */ + + m = MP_COPY(m); + MP_SPLIT(m); + + /* --- Check the radix for sensibleness --- */ + + if (radix > 0) + assert(((void)"ascii radix must be <= 62", radix <= 62)); + else if (radix < 0) + assert(((void)"binary radix must fit in a byte", -radix <= UCHAR_MAX)); + else + assert(((void)"radix can't be zero in mp_write", 0)); + + /* --- If the number is negative, sort that out --- */ + + if (MP_NEGP(m)) { + assert(radix > 0); + if (ops->put("-", 1, p)) + return (EOF); + m->f &= ~MP_NEG; + } + + /* --- Handle binary radix --- */ + + switch (radix) { + case 2: case -2: return (binary(m, 1, radix, ops, p)); + case 4: case -4: return (binary(m, 2, radix, ops, p)); + case 8: case -8: return (binary(m, 3, radix, ops, p)); + case 16: case -16: return (binary(m, 4, radix, ops, p)); + case 32: case -32: return (binary(m, 5, radix, ops, p)); + case -64: return (binary(m, 6, radix, ops, p)); + case -128: return (binary(m, 7, radix, ops, p)); + } + + /* --- If the number is small, do it the easy way --- */ + + if (MP_LEN(m) < 2) + rc = simple(MP_LEN(m) ? m->v[0] : 0, radix, 0, ops, p); + + /* --- Use a clever algorithm --- * + * + * Square the radix repeatedly, remembering old results, until I get + * something more than half the size of the number @m@. Use this to divide + * the number: the quotient and remainder will be approximately the same + * size, and I'll have split them on a digit boundary, so I can just emit + * the quotient and remainder recursively, in order. + */ + + else { + mp *pr[DEPTH]; + size_t target = (MP_LEN(m) + 1) / 2; + unsigned i = 0; + mp *z = mp_new(1, 0); + + /* --- Set up the exponent table --- */ + + z->v[0] = (radix > 0 ? radix : -radix); + z->f = 0; + for (;;) { + assert(((void)"Number is too unimaginably huge", i < DEPTH)); + pr[i++] = z; + if (MP_LEN(z) > target) + break; + z = mp_sqr(MP_NEW, z); + } + + /* --- Write out the answer --- */ + + rc = complicated(m, radix, pr, i - 1, 0, ops, p); + + /* --- Tidy away the array --- */ + + while (i > 0) + mp_drop(pr[--i]); + } + + /* --- Tidying up code --- */ + + MP_DROP(m); + return (rc); +} + +/*----- Test rig ----------------------------------------------------------*/ + +#ifdef TEST_RIG + +#include + +static int verify(dstr *v) +{ + int ok = 1; + int ib = *(int *)v[0].buf, ob = *(int *)v[2].buf; + dstr d = DSTR_INIT; + size_t off = 0; + mp *m = mp_readdstr(MP_NEW, &v[1], &off, ib); + if (m) { + if (!ob) { + fprintf(stderr, "*** unexpected successful parse\n" + "*** input [%2i] = ", ib); + if (ib < 0) + type_hex.dump(&v[1], stderr); + else + fputs(v[1].buf, stderr); + mp_writedstr(m, &d, 10); + fprintf(stderr, "\n*** (value = %s)\n", d.buf); + ok = 0; + } else { + mp_writedstr(m, &d, ob); + if (d.len != v[3].len || memcmp(d.buf, v[3].buf, d.len) != 0) { + fprintf(stderr, "*** failed read or write\n" + "*** input [%2i] = ", ib); + if (ib < 0) + type_hex.dump(&v[1], stderr); + else + fputs(v[1].buf, stderr); + fprintf(stderr, "\n*** output [%2i] = ", ob); + if (ob < 0) + type_hex.dump(&d, stderr); + else + fputs(d.buf, stderr); + fprintf(stderr, "\n*** expected [%2i] = ", ob); + if (ob < 0) + type_hex.dump(&v[3], stderr); + else + fputs(v[3].buf, stderr); + fputc('\n', stderr); + ok = 0; + } + } + mp_drop(m); + } else { + if (ob) { + fprintf(stderr, "*** unexpected parse failure\n" + "*** input [%2i] = ", ib); + if (ib < 0) + type_hex.dump(&v[1], stderr); + else + fputs(v[1].buf, stderr); + fprintf(stderr, "\n*** expected [%2i] = ", ob); + if (ob < 0) + type_hex.dump(&v[3], stderr); + else + fputs(v[3].buf, stderr); + fputc('\n', stderr); + ok = 0; + } + } + + if (v[1].len - off != v[4].len || + memcmp(v[1].buf + off, v[4].buf, v[4].len) != 0) { + fprintf(stderr, "*** leftovers incorrect\n" + "*** input [%2i] = ", ib); + if (ib < 0) + type_hex.dump(&v[1], stderr); + else + fputs(v[1].buf, stderr); + fprintf(stderr, "\n*** expected `%s'\n" + "*** found `%s'\n", + v[4].buf, v[1].buf + off); + ok = 0; + } + + dstr_destroy(&d); + assert(mparena_count(MPARENA_GLOBAL) == 0); + return (ok); +} + +static test_chunk tests[] = { + { "mptext-ascii", verify, + { &type_int, &type_string, &type_int, &type_string, &type_string, 0 } }, + { "mptext-bin-in", verify, + { &type_int, &type_hex, &type_int, &type_string, &type_string, 0 } }, + { "mptext-bin-out", verify, + { &type_int, &type_string, &type_int, &type_hex, &type_string, 0 } }, + { 0, 0, { 0 } } +}; + +int main(int argc, char *argv[]) +{ + sub_init(); + test_run(argc, argv, tests, SRCDIR "/t/mptext"); + return (0); +} + +#endif + +/*----- That's all, folks -------------------------------------------------*/