chiark - git - mdw - catacomb/blob - mptext.c

   1 /* -*-c-*-
   2  *
   3  * $Id: mptext.c,v 1.6 2000/06/25 12:58:23 mdw Exp $
   4  *
   5  * Textual representation of multiprecision numbers
   6  *
   7  * (c) 1999 Straylight/Edgeware
   8  */
   9
  10 /*----- Licensing notice --------------------------------------------------*
  11  *
  12  * This file is part of Catacomb.
  13  *
  14  * Catacomb is free software; you can redistribute it and/or modify
  15  * it under the terms of the GNU Library General Public License as
  16  * published by the Free Software Foundation; either version 2 of the
  17  * License, or (at your option) any later version.
  18  *
  19  * Catacomb is distributed in the hope that it will be useful,
  20  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  21  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  22  * GNU Library General Public License for more details.
  23  *
  24  * You should have received a copy of the GNU Library General Public
  25  * License along with Catacomb; if not, write to the Free
  26  * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
  27  * MA 02111-1307, USA.
  28  */
  29
  30 /*----- Revision history --------------------------------------------------*
  31  *
  32  * $Log: mptext.c,v $
  33  * Revision 1.6  2000/06/25 12:58:23  mdw
  34  * Fix the derivation of `depth' commentary.
  35  *
  36  * Revision 1.5  2000/06/17 11:46:19  mdw
  37  * New and much faster stack-based algorithm for reading integers.  Support
  38  * reading and writing binary integers in bases between 2 and 256.
  39  *
  40  * Revision 1.4  1999/12/22 15:56:56  mdw
  41  * Use clever recursive algorithm for writing numbers out.
  42  *
  43  * Revision 1.3  1999/12/10 23:23:26  mdw
  44  * Allocate slightly less memory.
  45  *
  46  * Revision 1.2  1999/11/20 22:24:15  mdw
  47  * Use function versions of MPX_UMULN and MPX_UADDN.
  48  *
  49  * Revision 1.1  1999/11/17 18:02:16  mdw
  50  * New multiprecision integer arithmetic suite.
  51  *
  52  */
  53
  54 /*----- Header files ------------------------------------------------------*/
  55
  56 #include <ctype.h>
  57 #include <limits.h>
  58 #include <stdio.h>
  59
  60 #include "mp.h"
  61 #include "mptext.h"
  62 #include "paranoia.h"
  63
  64 /*----- Magical numbers ---------------------------------------------------*/
  65
  66 /* --- Maximum recursion depth --- *
  67  *
  68  * This is the number of bits in a @size_t@ object.  Why?
  69  *
  70  * To see this, let %$b = \mathit{MPW\_MAX} + 1$% and let %$Z$% be the
  71  * largest @size_t@ value.  Then the largest possible @mp@ is %$M - 1$% where
  72  * %$M = b^Z$%.  Let %$r$% be a radix to read or write.  Since the recursion
  73  * squares the radix at each step, the highest number reached by the
  74  * recursion is %$d$%, where:
  75  *
  76  *   %$r^{2^d} = b^Z$%.
  77  *
  78  * Solving gives that %$d = \lg \log_r b^Z$%.  If %$r = 2$%, this is maximum,
  79  * so choosing %$d = \lg \lg b^Z = \lg (Z \lg b) = \lg Z + \lg \lg b$%.
  80  *
  81  * Expressing %$\lg Z$% as @CHAR_BIT * sizeof(size_t)@ yields an
  82  * overestimate, since a @size_t@ representation may contain `holes'.
  83  * Choosing to represent %$\lg \lg b$% by 10 is almost certainly sufficient
  84  * for `some time to come'.
  85  */
  86
  87 #define DEPTH (CHAR_BIT * sizeof(size_t) + 10)
  88
  89 /*----- Main code ---------------------------------------------------------*/
  90
  91 /* --- @mp_read@ --- *
  92  *
  93  * Arguments:   @mp *m@ = destination multiprecision number
  94  *              @int radix@ = base to assume for data (or zero to guess)
  95  *              @const mptext_ops *ops@ = pointer to operations block
  96  *              @void *p@ = data for the operations block
  97  *
  98  * Returns:     The integer read, or zero if it didn't work.
  99  *
 100  * Use:         Reads an integer from some source.  If the @radix@ is
 101  *              specified, the number is assumed to be given in that radix,
 102  *              with the letters `a' (either upper- or lower-case) upwards
 103  *              standing for digits greater than 9.  Otherwise, base 10 is
 104  *              assumed unless the number starts with `0' (octal), `0x' (hex)
 105  *              or `nnn_' (base `nnn').  An arbitrary amount of whitespace
 106  *              before the number is ignored.
 107  */
 108
 109 /* --- About the algorithm --- *
 110  *
 111  * The algorithm here is rather aggressive.  I maintain an array of
 112  * successive squarings of the radix, and a stack of partial results, each
 113  * with a counter attached indicating which radix square to multiply by.
 114  * Once the item at the top of the stack reaches the same counter level as
 115  * the next item down, they are combined together and the result is given a
 116  * counter level one higher than either of the results.
 117  *
 118  * Gluing the results together at the end is slightly tricky.  Pay attention
 119  * to the code.
 120  *
 121  * This is more complicated because of the need to handle the slightly
 122  * bizarre syntax.
 123  */
 124
 125 mp *mp_read(mp *m, int radix, const mptext_ops *ops, void *p)
 126 {
 127   int ch;                               /* Current char being considered */
 128   unsigned f = 0;                       /* Flags about the current number */
 129   int r;                                /* Radix to switch over to */
 130   mpw rd;                               /* Radix as an @mp@ digit */
 131   mp rr;                                /* The @mp@ for the radix */
 132   unsigned nf = m ? m->f & MP_BURN : 0; /* New @mp@ flags */
 133
 134   /* --- Stacks --- */
 135
 136   mp *pow[DEPTH];                       /* List of powers */
 137   unsigned pows;                        /* Next index to fill */
 138   struct { unsigned i; mp *m; } s[DEPTH]; /* Main stack */
 139   unsigned sp;                          /* Current stack pointer */
 140
 141   /* --- Flags --- */
 142
 143   enum {
 144     f_neg = 1u,
 145     f_ok = 2u
 146   };
 147
 148   /* --- Initialize the stacks --- */
 149
 150   mp_build(&rr, &rd, &rd + 1);
 151   pow[0] = &rr;
 152   pows = 1;
 153
 154   sp = 0;
 155
 156   /* --- Initialize the destination number --- */
 157
 158   if (m)
 159     MP_DROP(m);
 160
 161   /* --- Read an initial character --- */
 162
 163   ch = ops->get(p);
 164   while (isspace(ch))
 165     ch = ops->get(p);
 166
 167   /* --- Handle an initial sign --- */
 168
 169   if (ch == '-') {
 170     f |= f_neg;
 171     ch = ops->get(p);
 172     while (isspace(ch))
 173       ch = ops->get(p);
 174   }
 175
 176   /* --- If the radix is zero, look for leading zeros --- */
 177
 178   if (radix > 0) {
 179     assert(((void)"ascii radix must be <= 36", radix <= 36));
 180     rd = radix;
 181     r = -1;
 182   } else if (radix < 0) {
 183     rd = -radix;
 184     assert(((void)"binary radix must fit in a byte ", rd < UCHAR_MAX));
 185     r = -1;
 186   } else if (ch != '0') {
 187     rd = 10;
 188     r = 0;
 189   } else {
 190     ch = ops->get(p);
 191     if (ch == 'x') {
 192       ch = ops->get(p);
 193       rd = 16;
 194     } else {
 195       rd = 8;
 196       f |= f_ok;
 197     }
 198     r = -1;
 199   }
 200
 201   /* --- Time to start --- */
 202
 203   for (;; ch = ops->get(p)) {
 204     int x;
 205
 206     /* --- An underscore indicates a numbered base --- */
 207
 208     if (ch == '_' && r > 0 && r <= 36) {
 209       unsigned i;
 210
 211       /* --- Clear out the stacks --- */
 212
 213       for (i = 1; i < pows; i++)
 214         MP_DROP(pow[i]);
 215       pows = 1;
 216       for (i = 0; i < sp; i++)
 217         MP_DROP(s[i].m);
 218       sp = 0;
 219
 220       /* --- Restart the search --- */
 221
 222       rd = r;
 223       r = -1;
 224       f &= ~f_ok;
 225       continue;
 226     }
 227
 228     /* --- Check that the character is a digit and in range --- */
 229
 230     if (radix < 0)
 231       x = ch;
 232     else {
 233       if (!isalnum(ch))
 234         break;
 235       if (ch >= '0' && ch <= '9')
 236         x = ch - '0';
 237       else {
 238         ch = tolower(ch);
 239         if (ch >= 'a' && ch <= 'z')     /* ASCII dependent! */
 240           x = ch - 'a' + 10;
 241         else
 242           break;
 243       }
 244     }
 245
 246     /* --- Sort out what to do with the character --- */
 247
 248     if (x >= 10 && r >= 0)
 249       r = -1;
 250     if (x >= rd)
 251       break;
 252
 253     if (r >= 0)
 254       r = r * 10 + x;
 255
 256     /* --- Stick the character on the end of my integer --- */
 257
 258     assert(((void)"Number is too unimaginably huge", sp < DEPTH));
 259     s[sp].m = m = mp_new(1, nf);
 260     m->v[0] = x;
 261     s[sp].i = 0;
 262
 263     /* --- Now grind through the stack --- */
 264
 265     while (sp > 0 && s[sp - 1].i == s[sp].i) {
 266
 267       /* --- Combine the top two items --- */
 268
 269       sp--;
 270       m = s[sp].m;
 271       m = mp_mul(m, m, pow[s[sp].i]);
 272       m = mp_add(m, m, s[sp + 1].m);
 273       s[sp].m = m;
 274       MP_DROP(s[sp + 1].m);
 275       s[sp].i++;
 276
 277       /* --- Make a new radix power if necessary --- */
 278
 279       if (s[sp].i >= pows) {
 280         assert(((void)"Number is too unimaginably huge", pows < DEPTH));
 281         pow[pows] = mp_sqr(MP_NEW, pow[pows - 1]);
 282         pows++;
 283       }
 284     }
 285     f |= f_ok;
 286     sp++;
 287   }
 288
 289   ops->unget(ch, p);
 290
 291   /* --- If we're done, compute the rest of the number --- */
 292
 293   if (f & f_ok) {
 294     if (!sp)
 295       return (MP_ZERO);
 296     else {
 297       mp *z = MP_ONE;
 298       sp--;
 299
 300       while (sp > 0) {
 301
 302         /* --- Combine the top two items --- */
 303
 304         sp--;
 305         m = s[sp].m;
 306         z = mp_mul(z, z, pow[s[sp + 1].i]);
 307         m = mp_mul(m, m, z);
 308         m = mp_add(m, m, s[sp + 1].m);
 309         s[sp].m = m;
 310         MP_DROP(s[sp + 1].m);
 311
 312         /* --- Make a new radix power if necessary --- */
 313
 314         if (s[sp].i >= pows) {
 315           assert(((void)"Number is too unimaginably huge", pows < DEPTH));
 316           pow[pows] = mp_sqr(MP_NEW, pow[pows - 1]);
 317           pows++;
 318         }
 319       }
 320       MP_DROP(z);
 321       m = s[0].m;
 322     }
 323   } else {
 324     unsigned i;
 325     for (i = 0; i < sp; i++)
 326       MP_DROP(s[i].m);
 327   }
 328
 329   /* --- Clear the radix power list --- */
 330
 331   {
 332     unsigned i;
 333     for (i = 1; i < pows; i++)
 334       MP_DROP(pow[i]);
 335   }
 336
 337   /* --- Bail out if the number was bad --- */
 338
 339   if (!(f & f_ok))
 340     return (0);
 341
 342   /* --- Set the sign and return --- */
 343
 344   if (f & f_neg)
 345     m->f |= MP_NEG;
 346   return (m);
 347 }
 348
 349 /* --- @mp_write@ --- *
 350  *
 351  * Arguments:   @mp *m@ = pointer to a multi-precision integer
 352  *              @int radix@ = radix to use when writing the number out
 353  *              @const mptext_ops *ops@ = pointer to an operations block
 354  *              @void *p@ = data for the operations block
 355  *
 356  * Returns:     Zero if it worked, nonzero otherwise.
 357  *
 358  * Use:         Writes a large integer in textual form.
 359  */
 360
 361 /* --- Simple case --- *
 362  *
 363  * Use a fixed-sized buffer and the simple single-precision division
 364  * algorithm to pick off low-order digits.  Put each digit in a buffer,
 365  * working backwards from the end.  If the buffer becomes full, recurse to
 366  * get another one.  Ensure that there are at least @z@ digits by writing
 367  * leading zeroes if there aren't enough real digits.
 368  */
 369
 370 static int simple(mp *m, int radix, unsigned z,
 371                   const mptext_ops *ops, void *p)
 372 {
 373   int rc = 0;
 374   char buf[64];
 375   unsigned i = sizeof(buf);
 376   int rd = radix > 0 ? radix : -radix;
 377
 378   do {
 379     int ch;
 380     mpw x;
 381
 382     x = mpx_udivn(m->v, m->vl, m->v, m->vl, rd);
 383     MP_SHRINK(m);
 384     if (radix < 0)
 385       ch = x;
 386     else {
 387       if (x < 10)
 388         ch = '0' + x;
 389       else
 390         ch = 'a' + x - 10;
 391     }
 392     buf[--i] = ch;
 393     if (z)
 394       z--;
 395   } while (i && MP_LEN(m));
 396
 397   if (MP_LEN(m))
 398     rc = simple(m, radix, z, ops, p);
 399   else {
 400     static const char zero[32] = "00000000000000000000000000000000";
 401     while (!rc && z >= sizeof(zero)) {
 402       rc = ops->put(zero, sizeof(zero), p);
 403       z -= sizeof(zero);
 404     }
 405     if (!rc && z)
 406       rc = ops->put(zero, z, p);
 407   }
 408   if (!rc)
 409     ops->put(buf + i, sizeof(buf) - i, p);
 410   if (m->f & MP_BURN)
 411     BURN(buf);
 412   return (rc);
 413 }
 414
 415 /* --- Complicated case --- *
 416  *
 417  * If the number is small, fall back to the simple case above.  Otherwise
 418  * divide and take remainder by current large power of the radix, and emit
 419  * each separately.  Don't emit a zero quotient.  Be very careful about
 420  * leading zeroes on the remainder part, because they're deeply significant.
 421  */
 422
 423 static int complicated(mp *m, int radix, mp **pr, unsigned i, unsigned z,
 424                        const mptext_ops *ops, void *p)
 425 {
 426   int rc = 0;
 427   mp *q = MP_NEW;
 428   unsigned d = 1 << i;
 429
 430   if (MP_LEN(m) < 8)
 431     return (simple(m, radix, z, ops, p));
 432
 433   mp_div(&q, &m, m, pr[i]);
 434   if (!MP_LEN(q))
 435     d = z;
 436   else {
 437     if (z > d)
 438       z -= d;
 439     else
 440       z = 0;
 441     rc = complicated(q, radix, pr, i - 1, z, ops, p);
 442   }
 443   if (!rc)
 444     rc = complicated(m, radix, pr, i - 1, d, ops, p);
 445   mp_drop(q);
 446   return (rc);
 447 }
 448
 449 /* --- Main driver code --- */
 450
 451 int mp_write(mp *m, int radix, const mptext_ops *ops, void *p)
 452 {
 453   int rc;
 454
 455   /* --- Set various things up --- */
 456
 457   m = MP_COPY(m);
 458   MP_SPLIT(m);
 459
 460   /* --- Check the radix for sensibleness --- */
 461
 462   if (radix > 0)
 463     assert(((void)"ascii radix must be <= 36", radix <= 36));
 464   else if (radix < 0)
 465     assert(((void)"binary radix must fit in a byte", -radix < UCHAR_MAX));
 466   else
 467     assert(((void)"radix can't be zero in mp_write", 0));
 468
 469   /* --- If the number is negative, sort that out --- */
 470
 471   if (m->f & MP_NEG) {
 472     if (ops->put("-", 1, p))
 473       return (EOF);
 474     m->f &= ~MP_NEG;
 475   }
 476
 477   /* --- If the number is small, do it the easy way --- */
 478
 479   if (MP_LEN(m) < 8)
 480     rc = simple(m, radix, 0, ops, p);
 481
 482   /* --- Use a clever algorithm --- *
 483    *
 484    * Square the radix repeatedly, remembering old results, until I get
 485    * something more than half the size of the number @m@.  Use this to divide
 486    * the number: the quotient and remainder will be approximately the same
 487    * size, and I'll have split them on a digit boundary, so I can just emit
 488    * the quotient and remainder recursively, in order.
 489    */
 490
 491   else {
 492     mp *pr[DEPTH];
 493     size_t target = MP_LEN(m) / 2;
 494     unsigned i = 0;
 495     mp *z = mp_new(1, 0);
 496
 497     /* --- Set up the exponent table --- */
 498
 499     z->v[0] = (radix > 0 ? radix : -radix);
 500     z->f = 0;
 501     for (;;) {
 502       assert(((void)"Number is too unimaginably huge", i < DEPTH));
 503       pr[i++] = z;
 504       if (MP_LEN(z) > target)
 505         break;
 506       z = mp_sqr(MP_NEW, z);
 507     }
 508
 509     /* --- Write out the answer --- */
 510
 511     rc = complicated(m, radix, pr, i - 1, 0, ops, p);
 512
 513     /* --- Tidy away the array --- */
 514
 515     while (i > 0)
 516       mp_drop(pr[--i]);
 517   }
 518
 519   /* --- Tidying up code --- */
 520
 521   MP_DROP(m);
 522   return (rc);
 523 }
 524
 525 /*----- Test rig ----------------------------------------------------------*/
 526
 527 #ifdef TEST_RIG
 528
 529 #include <mLib/testrig.h>
 530
 531 static int verify(dstr *v)
 532 {
 533   int ok = 1;
 534   int ib = *(int *)v[0].buf, ob = *(int *)v[2].buf;
 535   dstr d = DSTR_INIT;
 536   mp *m = mp_readdstr(MP_NEW, &v[1], 0, ib);
 537   if (m) {
 538     if (!ob) {
 539       fprintf(stderr, "*** unexpected successful parse\n"
 540                       "*** input [%i] = ", ib);
 541       if (ib < 0)
 542         type_hex.dump(&v[1], stderr);
 543       else
 544         fputs(v[1].buf, stderr);
 545       mp_writedstr(m, &d, 10);
 546       fprintf(stderr, "\n*** (value = %s)\n", d.buf);
 547       ok = 0;
 548     } else {
 549       mp_writedstr(m, &d, ob);
 550       if (d.len != v[3].len || memcmp(d.buf, v[3].buf, d.len) != 0) {
 551         fprintf(stderr, "*** failed read or write\n"
 552                         "*** input [%i]    = ", ib);
 553         if (ib < 0)
 554           type_hex.dump(&v[1], stderr);
 555         else
 556           fputs(v[1].buf, stderr);
 557         fprintf(stderr, "\n*** output [%i]   = ", ob);
 558         if (ob < 0)
 559           type_hex.dump(&d, stderr);
 560         else
 561           fputs(d.buf, stderr);
 562         fprintf(stderr, "\n*** expected [%i]   = ", ob);
 563         if (ob < 0)
 564           type_hex.dump(&v[3], stderr);
 565         else
 566           fputs(v[3].buf, stderr);
 567         fputc('\n', stderr);
 568         ok = 0;
 569       }
 570     }
 571     mp_drop(m);
 572   } else {
 573     if (ob) {
 574       fprintf(stderr, "*** unexpected parse failure\n"
 575                       "*** input [%i]    = ", ib);
 576       if (ib < 0)
 577         type_hex.dump(&v[1], stderr);
 578       else
 579         fputs(v[1].buf, stderr);
 580       fprintf(stderr, "\n*** expected [%i]   = ", ob);
 581       if (ob < 0)
 582         type_hex.dump(&v[3], stderr);
 583       else
 584         fputs(v[3].buf, stderr);
 585       fputc('\n', stderr);
 586       ok = 0;
 587     }
 588   }
 589
 590   dstr_destroy(&d);
 591   assert(mparena_count(MPARENA_GLOBAL) == 0);
 592   return (ok);
 593 }
 594
 595 static test_chunk tests[] = {
 596   { "mptext-ascii", verify,
 597     { &type_int, &type_string, &type_int, &type_string, 0 } },
 598   { "mptext-bin-in", verify,
 599     { &type_int, &type_hex, &type_int, &type_string, 0 } },
 600   { "mptext-bin-out", verify,
 601     { &type_int, &type_string, &type_int, &type_hex, 0 } },
 602   { 0, 0, { 0 } }
 603 };
 604
 605 int main(int argc, char *argv[])
 606 {
 607   sub_init();
 608   test_run(argc, argv, tests, SRCDIR "/tests/mptext");
 609   return (0);
 610 }
 611
 612 #endif
 613
 614 /*----- That's all, folks -------------------------------------------------*/