X-Git-Url: https://git.distorted.org.uk/~mdw/catacomb/blobdiff_plain/ba6e6b64033b1f9de49feccb5c9cd438354481f7..0f00dc4c8eb47e67bc0f148c2dd109f73a451e0a:/mptext.c?ds=sidebyside diff --git a/mptext.c b/mptext.c deleted file mode 100644 index 8c00e346..00000000 --- a/mptext.c +++ /dev/null @@ -1,851 +0,0 @@ -/* -*-c-*- - * - * $Id$ - * - * 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 "/tests/mptext"); - return (0); -} - -#endif - -/*----- That's all, folks -------------------------------------------------*/