From e360a4f2010027d1bd17d9376d7c4beb91dbb715 Mon Sep 17 00:00:00 2001 From: mdw Date: Wed, 22 Dec 1999 15:56:56 +0000 Subject: [PATCH] Use clever recursive algorithm for writing numbers out. --- mptext.c | 176 +++++++++++++++++++++++++++++++++++++++++++++++++-------------- 1 file changed, 137 insertions(+), 39 deletions(-) diff --git a/mptext.c b/mptext.c index b216fea..c5d20e4 100644 --- a/mptext.c +++ b/mptext.c @@ -1,6 +1,6 @@ /* -*-c-*- * - * $Id: mptext.c,v 1.3 1999/12/10 23:23:26 mdw Exp $ + * $Id: mptext.c,v 1.4 1999/12/22 15:56:56 mdw Exp $ * * Textual representation of multiprecision numbers * @@ -30,6 +30,9 @@ /*----- Revision history --------------------------------------------------* * * $Log: mptext.c,v $ + * Revision 1.4 1999/12/22 15:56:56 mdw + * Use clever recursive algorithm for writing numbers out. + * * Revision 1.3 1999/12/10 23:23:26 mdw * Allocate slightly less memory. * @@ -46,17 +49,9 @@ #include #include -#include - #include "mp.h" #include "mptext.h" - -/*----- Data structures ---------------------------------------------------*/ - -#ifndef CHARV - DA_DECL(charv, char); -# define CHARV -#endif +#include "paranoia.h" /*----- Main code ---------------------------------------------------------*/ @@ -206,54 +201,157 @@ mp *mp_read(mp *m, int radix, const mptext_ops *ops, void *p) * Use: Writes a large integer in textual form. */ +/* --- Simple case --- * + * + * Use a fixed-sized buffer and the simple single-precision division + * algorithm 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(mp *m, int radix, unsigned z, + const mptext_ops *ops, void *p) +{ + int rc = 0; + char buf[64]; + unsigned i = sizeof(buf); + + do { + int ch; + mpw x; + + x = mpx_udivn(m->v, m->vl, m->v, m->vl, radix); + MP_SHRINK(m); + if (x < 10) + ch = '0' + x; + else + ch = 'a' + x - 10; + buf[--i] = ch; + if (z) + z--; + } while (i && MP_LEN(m)); + + if (MP_LEN(m)) + rc = simple(m, radix, z, ops, p); + else { + static const char zero[32] = "00000000000000000000000000000000"; + while (!rc && z >= sizeof(zero)) { + rc = ops->put(zero, sizeof(zero), p); + z -= sizeof(zero); + } + if (!rc && z) + rc = ops->put(zero, z, p); + } + if (!rc) + ops->put(buf + i, sizeof(buf) - i, p); + if (m->f & MP_BURN) + 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) < 8) + return (simple(m, radix, z, ops, p)); + + mp_div(&q, &m, m, pr[i]); + if (!MP_LEN(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); +} + +/* --- Main driver code --- */ + int mp_write(mp *m, int radix, const mptext_ops *ops, void *p) { - charv v = DA_INIT; - mpw b = radix; - mp bb; + int rc; /* --- Set various things up --- */ m = MP_COPY(m); - mp_build(&bb, &b, &b + 1); + MP_SPLIT(m); /* --- If the number is negative, sort that out --- */ if (m->f & MP_NEG) { if (ops->put("-", 1, p)) return (EOF); - MP_SPLIT(m); - m->f &= ~MP_NEG; } - /* --- Write digits to a temporary array --- */ - - do { - mp *q = MP_NEW, *r = MP_NEW; - int ch; - mpw x; + /* --- If the number is small, do it the easy way --- */ + + if (MP_LEN(m) < 8) + rc = simple(m, 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. + * + * The array size copes with the largest number possibly representable on + * the host machine. Such a large number shouldn't ever arise in real use. + */ + + else { + mp *pr[CHAR_BIT * sizeof(size_t)]; + size_t target = MP_LEN(m) / 2; + unsigned i = 0; + mp *z = mp_create(1); + + /* --- Set up the exponent table --- */ + + z->v[0] = radix; + z->f = 0; + for (;;) { + assert(((void)"Number is too unimaginably huge", + i < sizeof(pr) / sizeof(pr[0]))); + pr[i++] = z; + if (MP_LEN(z) > target) + break; + z = mp_sqr(MP_NEW, z); + } - mp_div(&q, &r, m, &bb); - x = r->v[0]; - if (x < 10) - ch = '0' + x; - else - ch = 'a' + x - 10; - DA_UNSHIFT(&v, ch); - MP_DROP(m); - MP_DROP(r); - m = q; - } while (MP_CMP(m, !=, MP_ZERO)); + /* --- Write out the answer --- */ - /* --- Finished that --- */ + rc = complicated(m, radix, pr, i - 1, 0, ops, p); - MP_DROP(m); + /* --- Tidy away the array --- */ - { - int rc = ops->put(DA(&v), DA_LEN(&v), p); - DA_DESTROY(&v); - return (rc ? EOF : 0); + while (i > 0) + mp_drop(pr[--i]); } + + /* --- Tidying up code --- */ + + MP_DROP(m); + return (rc); } /*----- Test rig ----------------------------------------------------------*/ -- 2.11.0