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

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