Fast estimation of number representation lengths.
[u/mdw/catacomb] / mptext.c
CommitLineData
d3409d5e 1/* -*-c-*-
2 *
9bca44cb 3 * $Id: mptext.c,v 1.16 2002/10/15 22:57:43 mdw Exp $
d3409d5e 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 $
9bca44cb 33 * Revision 1.16 2002/10/15 22:57:43 mdw
34 * Bug fix: prevent negative zero.
35 *
afd054c1 36 * Revision 1.15 2002/10/15 19:18:15 mdw
37 * Fix fencepost bugs in binary radix writing.
38 *
3db58009 39 * Revision 1.14 2002/10/09 00:33:44 mdw
40 * Allow `0o' and `0b' prefixes for octal and binary (from Haskell)
41 *
6ea6fe51 42 * Revision 1.13 2002/10/09 00:21:06 mdw
43 * Allow user-specified `r_xx' bases to be up to 62.
44 *
631673a1 45 * Revision 1.12 2002/01/13 19:51:18 mdw
46 * Extend the textual format to bases up to 62 by distinguishing case.
47 *
eaa515d8 48 * Revision 1.11 2001/06/16 23:42:17 mdw
49 * Typesetting fixes.
50 *
a951033d 51 * Revision 1.10 2001/06/16 13:22:39 mdw
52 * Added fast-track code for binary output bases, and tests.
53 *
3bc9cb53 54 * Revision 1.9 2001/02/03 16:05:17 mdw
55 * Make flags be unsigned. Improve the write algorithm: recurse until the
56 * parts are one word long and use single-precision arithmetic from there.
57 * Fix off-by-one bug when breaking the number apart.
58 *
9d3838a0 59 * Revision 1.8 2000/12/06 20:32:42 mdw
60 * Reduce binary bytes (to allow marker bits to be ignored). Fix error
61 * message string a bit. Allow leading `+' signs.
62 *
7d45ed6c 63 * Revision 1.7 2000/07/15 10:01:08 mdw
64 * Bug fix in binary input.
65 *
dd9199f0 66 * Revision 1.6 2000/06/25 12:58:23 mdw
67 * Fix the derivation of `depth' commentary.
68 *
2b26f2d7 69 * Revision 1.5 2000/06/17 11:46:19 mdw
70 * New and much faster stack-based algorithm for reading integers. Support
71 * reading and writing binary integers in bases between 2 and 256.
72 *
e360a4f2 73 * Revision 1.4 1999/12/22 15:56:56 mdw
74 * Use clever recursive algorithm for writing numbers out.
75 *
9c3df6c0 76 * Revision 1.3 1999/12/10 23:23:26 mdw
77 * Allocate slightly less memory.
78 *
90b6f0be 79 * Revision 1.2 1999/11/20 22:24:15 mdw
80 * Use function versions of MPX_UMULN and MPX_UADDN.
81 *
d3409d5e 82 * Revision 1.1 1999/11/17 18:02:16 mdw
83 * New multiprecision integer arithmetic suite.
84 *
85 */
86
87/*----- Header files ------------------------------------------------------*/
88
89#include <ctype.h>
2b26f2d7 90#include <limits.h>
d3409d5e 91#include <stdio.h>
92
d3409d5e 93#include "mp.h"
94#include "mptext.h"
e360a4f2 95#include "paranoia.h"
d3409d5e 96
2b26f2d7 97/*----- Magical numbers ---------------------------------------------------*/
98
99/* --- Maximum recursion depth --- *
100 *
101 * This is the number of bits in a @size_t@ object. Why?
102 *
eaa515d8 103 * To see this, let %$b = \textit{MPW\_MAX} + 1$% and let %$Z$% be the
dd9199f0 104 * largest @size_t@ value. Then the largest possible @mp@ is %$M - 1$% where
105 * %$M = b^Z$%. Let %$r$% be a radix to read or write. Since the recursion
106 * squares the radix at each step, the highest number reached by the
107 * recursion is %$d$%, where:
2b26f2d7 108 *
dd9199f0 109 * %$r^{2^d} = b^Z$%.
2b26f2d7 110 *
111 * Solving gives that %$d = \lg \log_r b^Z$%. If %$r = 2$%, this is maximum,
112 * so choosing %$d = \lg \lg b^Z = \lg (Z \lg b) = \lg Z + \lg \lg b$%.
113 *
114 * Expressing %$\lg Z$% as @CHAR_BIT * sizeof(size_t)@ yields an
115 * overestimate, since a @size_t@ representation may contain `holes'.
116 * Choosing to represent %$\lg \lg b$% by 10 is almost certainly sufficient
117 * for `some time to come'.
118 */
119
120#define DEPTH (CHAR_BIT * sizeof(size_t) + 10)
121
d3409d5e 122/*----- Main code ---------------------------------------------------------*/
123
124/* --- @mp_read@ --- *
125 *
126 * Arguments: @mp *m@ = destination multiprecision number
127 * @int radix@ = base to assume for data (or zero to guess)
128 * @const mptext_ops *ops@ = pointer to operations block
129 * @void *p@ = data for the operations block
130 *
131 * Returns: The integer read, or zero if it didn't work.
132 *
133 * Use: Reads an integer from some source. If the @radix@ is
134 * specified, the number is assumed to be given in that radix,
135 * with the letters `a' (either upper- or lower-case) upwards
136 * standing for digits greater than 9. Otherwise, base 10 is
137 * assumed unless the number starts with `0' (octal), `0x' (hex)
138 * or `nnn_' (base `nnn'). An arbitrary amount of whitespace
139 * before the number is ignored.
140 */
141
2b26f2d7 142/* --- About the algorithm --- *
143 *
144 * The algorithm here is rather aggressive. I maintain an array of
145 * successive squarings of the radix, and a stack of partial results, each
146 * with a counter attached indicating which radix square to multiply by.
147 * Once the item at the top of the stack reaches the same counter level as
148 * the next item down, they are combined together and the result is given a
149 * counter level one higher than either of the results.
150 *
151 * Gluing the results together at the end is slightly tricky. Pay attention
152 * to the code.
153 *
154 * This is more complicated because of the need to handle the slightly
155 * bizarre syntax.
156 */
157
d3409d5e 158mp *mp_read(mp *m, int radix, const mptext_ops *ops, void *p)
159{
2b26f2d7 160 int ch; /* Current char being considered */
161 unsigned f = 0; /* Flags about the current number */
162 int r; /* Radix to switch over to */
163 mpw rd; /* Radix as an @mp@ digit */
164 mp rr; /* The @mp@ for the radix */
165 unsigned nf = m ? m->f & MP_BURN : 0; /* New @mp@ flags */
166
167 /* --- Stacks --- */
168
169 mp *pow[DEPTH]; /* List of powers */
170 unsigned pows; /* Next index to fill */
171 struct { unsigned i; mp *m; } s[DEPTH]; /* Main stack */
172 unsigned sp; /* Current stack pointer */
173
174 /* --- Flags --- */
d3409d5e 175
3bc9cb53 176#define f_neg 1u
177#define f_ok 2u
a951033d 178#define f_start 4u
d3409d5e 179
2b26f2d7 180 /* --- Initialize the stacks --- */
181
182 mp_build(&rr, &rd, &rd + 1);
183 pow[0] = &rr;
184 pows = 1;
185
186 sp = 0;
187
d3409d5e 188 /* --- Initialize the destination number --- */
189
2b26f2d7 190 if (m)
191 MP_DROP(m);
d3409d5e 192
193 /* --- Read an initial character --- */
194
195 ch = ops->get(p);
196 while (isspace(ch))
197 ch = ops->get(p);
198
199 /* --- Handle an initial sign --- */
200
9d3838a0 201 if (radix >= 0 && (ch == '-' || ch == '+')) {
202 if (ch == '-')
203 f |= f_neg;
204 do ch = ops->get(p); while isspace(ch);
d3409d5e 205 }
206
207 /* --- If the radix is zero, look for leading zeros --- */
208
2b26f2d7 209 if (radix > 0) {
631673a1 210 assert(((void)"ascii radix must be <= 62", radix <= 62));
2b26f2d7 211 rd = radix;
212 r = -1;
213 } else if (radix < 0) {
214 rd = -radix;
9d3838a0 215 assert(((void)"binary radix must fit in a byte", rd < UCHAR_MAX));
d3409d5e 216 r = -1;
2b26f2d7 217 } else if (ch != '0') {
218 rd = 10;
d3409d5e 219 r = 0;
220 } else {
221 ch = ops->get(p);
3db58009 222 switch (ch) {
223 case 'x':
224 rd = 16;
225 goto prefix;
226 case 'o':
227 rd = 8;
228 goto prefix;
229 case 'b':
230 rd = 2;
231 goto prefix;
232 prefix:
233 ch = ops->get(p);
234 break;
235 default:
236 rd = 8;
237 f |= f_ok;
d3409d5e 238 }
239 r = -1;
240 }
241
a951033d 242 /* --- Use fast algorithm for binary radix --- *
243 *
244 * This is the restart point after having parsed a radix number from the
245 * input. We check whether the radix is binary, and if so use a fast
246 * algorithm which just stacks the bits up in the right order.
247 */
248
249restart:
250 switch (rd) {
251 unsigned bit;
252
253 case 2: bit = 1; goto bin;
254 case 4: bit = 2; goto bin;
255 case 8: bit = 3; goto bin;
256 case 16: bit = 4; goto bin;
257 case 32: bit = 5; goto bin;
258 case 64: bit = 6; goto bin;
259 case 128: bit = 7; goto bin;
260 default:
261 break;
262
263 /* --- The fast binary algorithm --- *
264 *
265 * We stack bits up starting at the top end of a word. When one word is
266 * full, we write it to the integer, and start another with the left-over
267 * bits. When the array in the integer is full, we resize using low-level
268 * calls and copy the current data to the top end. Finally, we do a single
269 * bit-shift when we know where the end of the number is.
270 */
271
272 bin: {
273 mpw a = 0;
274 unsigned b = MPW_BITS;
275 size_t len, n;
276 mpw *v;
277
278 m = mp_dest(MP_NEW, 1, nf);
279 len = n = m->sz;
280 n = len;
281 v = m->v + n;
282 for (;; ch = ops->get(p)) {
283 unsigned x;
284
285 if (ch < 0)
286 break;
287
288 /* --- Check that the character is a digit and in range --- */
289
290 if (radix < 0)
291 x = ch % rd;
292 else {
293 if (!isalnum(ch))
294 break;
295 if (ch >= '0' && ch <= '9')
296 x = ch - '0';
297 else {
631673a1 298 if (rd <= 36)
299 ch = tolower(ch);
a951033d 300 if (ch >= 'a' && ch <= 'z') /* ASCII dependent! */
301 x = ch - 'a' + 10;
631673a1 302 else if (ch >= 'A' && ch <= 'Z')
303 x = ch - 'A' + 36;
a951033d 304 else
305 break;
306 }
307 }
308 if (x >= rd)
309 break;
310
311 /* --- Feed the digit into the accumulator --- */
312
313 f |= f_ok;
314 if (!x && !(f & f_start))
315 continue;
316 f |= f_start;
317 if (b > bit) {
318 b -= bit;
319 a |= MPW(x) << b;
320 } else {
321 a |= MPW(x) >> (bit - b);
322 b += MPW_BITS - bit;
323 *--v = MPW(a);
324 n--;
325 if (!n) {
326 n = len;
327 len <<= 1;
328 v = mpalloc(m->a, len);
329 memcpy(v + n, m->v, MPWS(n));
330 mpfree(m->a, m->v);
331 m->v = v;
332 v = m->v + n;
333 }
334 a = (b < MPW_BITS) ? MPW(x) << b : 0;
335 }
336 }
337
338 /* --- Finish up --- */
339
340 if (!(f & f_ok)) {
341 mp_drop(m);
342 m = 0;
343 } else {
344 *--v = MPW(a);
345 n--;
346 m->sz = len;
347 m->vl = m->v + len;
348 m->f &= ~MP_UNDEF;
349 m = mp_lsr(m, m, (unsigned long)n * MPW_BITS + b);
350 }
351 goto done;
352 }}
353
d3409d5e 354 /* --- Time to start --- */
355
356 for (;; ch = ops->get(p)) {
a951033d 357 unsigned x;
d3409d5e 358
7d45ed6c 359 if (ch < 0)
360 break;
361
d3409d5e 362 /* --- An underscore indicates a numbered base --- */
363
6ea6fe51 364 if (ch == '_' && r > 0 && r <= 62) {
2b26f2d7 365 unsigned i;
366
367 /* --- Clear out the stacks --- */
368
369 for (i = 1; i < pows; i++)
370 MP_DROP(pow[i]);
371 pows = 1;
372 for (i = 0; i < sp; i++)
373 MP_DROP(s[i].m);
374 sp = 0;
375
376 /* --- Restart the search --- */
377
378 rd = r;
d3409d5e 379 r = -1;
380 f &= ~f_ok;
a951033d 381 ch = ops->get(p);
382 goto restart;
d3409d5e 383 }
384
385 /* --- Check that the character is a digit and in range --- */
386
2b26f2d7 387 if (radix < 0)
9d3838a0 388 x = ch % rd;
d3409d5e 389 else {
2b26f2d7 390 if (!isalnum(ch))
d3409d5e 391 break;
2b26f2d7 392 if (ch >= '0' && ch <= '9')
393 x = ch - '0';
394 else {
631673a1 395 if (rd <= 36)
396 ch = tolower(ch);
2b26f2d7 397 if (ch >= 'a' && ch <= 'z') /* ASCII dependent! */
398 x = ch - 'a' + 10;
631673a1 399 else if (ch >= 'A' && ch <= 'Z')
400 x = ch - 'A' + 36;
2b26f2d7 401 else
402 break;
403 }
d3409d5e 404 }
405
406 /* --- Sort out what to do with the character --- */
407
408 if (x >= 10 && r >= 0)
409 r = -1;
2b26f2d7 410 if (x >= rd)
d3409d5e 411 break;
412
413 if (r >= 0)
414 r = r * 10 + x;
415
416 /* --- Stick the character on the end of my integer --- */
417
2b26f2d7 418 assert(((void)"Number is too unimaginably huge", sp < DEPTH));
419 s[sp].m = m = mp_new(1, nf);
420 m->v[0] = x;
421 s[sp].i = 0;
422
423 /* --- Now grind through the stack --- */
424
425 while (sp > 0 && s[sp - 1].i == s[sp].i) {
426
427 /* --- Combine the top two items --- */
428
429 sp--;
430 m = s[sp].m;
431 m = mp_mul(m, m, pow[s[sp].i]);
432 m = mp_add(m, m, s[sp + 1].m);
433 s[sp].m = m;
434 MP_DROP(s[sp + 1].m);
435 s[sp].i++;
436
437 /* --- Make a new radix power if necessary --- */
438
439 if (s[sp].i >= pows) {
440 assert(((void)"Number is too unimaginably huge", pows < DEPTH));
441 pow[pows] = mp_sqr(MP_NEW, pow[pows - 1]);
442 pows++;
443 }
444 }
d3409d5e 445 f |= f_ok;
2b26f2d7 446 sp++;
d3409d5e 447 }
448
449 ops->unget(ch, p);
450
2b26f2d7 451 /* --- If we're done, compute the rest of the number --- */
452
453 if (f & f_ok) {
454 if (!sp)
455 return (MP_ZERO);
456 else {
457 mp *z = MP_ONE;
458 sp--;
459
460 while (sp > 0) {
461
462 /* --- Combine the top two items --- */
463
464 sp--;
465 m = s[sp].m;
466 z = mp_mul(z, z, pow[s[sp + 1].i]);
467 m = mp_mul(m, m, z);
468 m = mp_add(m, m, s[sp + 1].m);
469 s[sp].m = m;
470 MP_DROP(s[sp + 1].m);
471
472 /* --- Make a new radix power if necessary --- */
473
474 if (s[sp].i >= pows) {
475 assert(((void)"Number is too unimaginably huge", pows < DEPTH));
476 pow[pows] = mp_sqr(MP_NEW, pow[pows - 1]);
477 pows++;
478 }
479 }
480 MP_DROP(z);
481 m = s[0].m;
482 }
483 } else {
484 unsigned i;
485 for (i = 0; i < sp; i++)
486 MP_DROP(s[i].m);
487 }
488
489 /* --- Clear the radix power list --- */
490
491 {
492 unsigned i;
493 for (i = 1; i < pows; i++)
494 MP_DROP(pow[i]);
495 }
496
d3409d5e 497 /* --- Bail out if the number was bad --- */
498
a951033d 499done:
2b26f2d7 500 if (!(f & f_ok))
d3409d5e 501 return (0);
d3409d5e 502
503 /* --- Set the sign and return --- */
504
d3409d5e 505 if (f & f_neg)
506 m->f |= MP_NEG;
9bca44cb 507 MP_SHRINK(m);
d3409d5e 508 return (m);
3bc9cb53 509
a951033d 510#undef f_start
3bc9cb53 511#undef f_neg
512#undef f_ok
d3409d5e 513}
514
515/* --- @mp_write@ --- *
516 *
517 * Arguments: @mp *m@ = pointer to a multi-precision integer
518 * @int radix@ = radix to use when writing the number out
519 * @const mptext_ops *ops@ = pointer to an operations block
520 * @void *p@ = data for the operations block
521 *
522 * Returns: Zero if it worked, nonzero otherwise.
523 *
524 * Use: Writes a large integer in textual form.
525 */
526
e360a4f2 527/* --- Simple case --- *
528 *
3bc9cb53 529 * Use a fixed-sized buffer and single-precision arithmetic to pick off
530 * low-order digits. Put each digit in a buffer, working backwards from the
531 * end. If the buffer becomes full, recurse to get another one. Ensure that
532 * there are at least @z@ digits by writing leading zeroes if there aren't
533 * enough real digits.
e360a4f2 534 */
535
3bc9cb53 536static int simple(mpw n, int radix, unsigned z,
e360a4f2 537 const mptext_ops *ops, void *p)
538{
539 int rc = 0;
540 char buf[64];
541 unsigned i = sizeof(buf);
2b26f2d7 542 int rd = radix > 0 ? radix : -radix;
e360a4f2 543
544 do {
545 int ch;
546 mpw x;
547
3bc9cb53 548 x = n % rd;
549 n /= rd;
2b26f2d7 550 if (radix < 0)
551 ch = x;
3bc9cb53 552 else if (x < 10)
553 ch = '0' + x;
631673a1 554 else if (x < 36) /* Ascii specific */
3bc9cb53 555 ch = 'a' + x - 10;
631673a1 556 else
557 ch = 'A' + x - 36;
e360a4f2 558 buf[--i] = ch;
559 if (z)
560 z--;
3bc9cb53 561 } while (i && n);
e360a4f2 562
3bc9cb53 563 if (n)
564 rc = simple(n, radix, z, ops, p);
e360a4f2 565 else {
a951033d 566 char zbuf[32];
567 memset(zbuf, (radix < 0) ? 0 : '0', sizeof(zbuf));
568 while (!rc && z >= sizeof(zbuf)) {
569 rc = ops->put(zbuf, sizeof(zbuf), p);
570 z -= sizeof(zbuf);
e360a4f2 571 }
572 if (!rc && z)
a951033d 573 rc = ops->put(zbuf, z, p);
e360a4f2 574 }
575 if (!rc)
3bc9cb53 576 rc = ops->put(buf + i, sizeof(buf) - i, p);
577 BURN(buf);
e360a4f2 578 return (rc);
579}
580
581/* --- Complicated case --- *
582 *
583 * If the number is small, fall back to the simple case above. Otherwise
584 * divide and take remainder by current large power of the radix, and emit
585 * each separately. Don't emit a zero quotient. Be very careful about
586 * leading zeroes on the remainder part, because they're deeply significant.
587 */
588
589static int complicated(mp *m, int radix, mp **pr, unsigned i, unsigned z,
590 const mptext_ops *ops, void *p)
591{
592 int rc = 0;
593 mp *q = MP_NEW;
594 unsigned d = 1 << i;
595
3bc9cb53 596 if (MP_LEN(m) < 2)
597 return (simple(MP_LEN(m) ? m->v[0] : 0, radix, z, ops, p));
e360a4f2 598
3bc9cb53 599 assert(i);
e360a4f2 600 mp_div(&q, &m, m, pr[i]);
601 if (!MP_LEN(q))
602 d = z;
603 else {
604 if (z > d)
605 z -= d;
606 else
607 z = 0;
608 rc = complicated(q, radix, pr, i - 1, z, ops, p);
609 }
610 if (!rc)
611 rc = complicated(m, radix, pr, i - 1, d, ops, p);
612 mp_drop(q);
613 return (rc);
614}
615
a951033d 616/* --- Binary case --- *
617 *
618 * Special case for binary output. Goes much faster.
619 */
620
621static int binary(mp *m, int bit, int radix, const mptext_ops *ops, void *p)
622{
623 mpw *v;
624 mpw a;
625 int rc = 0;
626 unsigned b;
627 unsigned mask;
628 unsigned long n;
629 unsigned f = 0;
630 char buf[8], *q;
631 unsigned x;
632 int ch;
633
634#define f_out 1u
635
636 /* --- Work out where to start --- */
637
638 n = mp_bits(m);
afd054c1 639 if (n % bit)
640 n += bit - (n % bit);
a951033d 641 b = n % MPW_BITS;
642 n /= MPW_BITS;
afd054c1 643
644 if (n >= MP_LEN(m)) {
a951033d 645 n--;
646 b += MPW_BITS;
647 }
648
649 v = m->v + n;
650 a = *v;
651 mask = (1 << bit) - 1;
652 q = buf;
653
654 /* --- Main code --- */
655
656 for (;;) {
657 if (b > bit) {
658 b -= bit;
659 x = a >> b;
660 } else {
661 x = a << (bit - b);
662 b += MPW_BITS - bit;
663 if (v == m->v)
664 break;
665 a = *--v;
666 if (b < MPW_BITS)
667 x |= a >> b;
668 }
669 x &= mask;
670 if (!x && !(f & f_out))
671 continue;
672
673 if (radix < 0)
674 ch = x;
675 else if (x < 10)
676 ch = '0' + x;
631673a1 677 else if (x < 36)
678 ch = 'a' + x - 10; /* Ascii specific */
a951033d 679 else
631673a1 680 ch = 'A' + x - 36;
a951033d 681 *q++ = ch;
682 if (q >= buf + sizeof(buf)) {
683 if ((rc = ops->put(buf, sizeof(buf), p)) != 0)
684 goto done;
685 q = buf;
686 }
687 f |= f_out;
688 }
689
690 x &= mask;
691 if (radix < 0)
692 ch = x;
693 else if (x < 10)
694 ch = '0' + x;
631673a1 695 else if (x < 36)
696 ch = 'a' + x - 10; /* Ascii specific */
a951033d 697 else
631673a1 698 ch = 'A' + x - 36;
a951033d 699 *q++ = ch;
700 rc = ops->put(buf, q - buf, p);
701
702done:
703 mp_drop(m);
704 return (rc);
705
706#undef f_out
707}
708
e360a4f2 709/* --- Main driver code --- */
710
d3409d5e 711int mp_write(mp *m, int radix, const mptext_ops *ops, void *p)
712{
e360a4f2 713 int rc;
d3409d5e 714
afd054c1 715 if (MP_EQ(m, MP_ZERO))
716 return (ops->put("0", 1, p));
717
d3409d5e 718 /* --- Set various things up --- */
719
720 m = MP_COPY(m);
e360a4f2 721 MP_SPLIT(m);
d3409d5e 722
2b26f2d7 723 /* --- Check the radix for sensibleness --- */
724
725 if (radix > 0)
631673a1 726 assert(((void)"ascii radix must be <= 62", radix <= 62));
2b26f2d7 727 else if (radix < 0)
728 assert(((void)"binary radix must fit in a byte", -radix < UCHAR_MAX));
729 else
730 assert(((void)"radix can't be zero in mp_write", 0));
731
d3409d5e 732 /* --- If the number is negative, sort that out --- */
733
734 if (m->f & MP_NEG) {
735 if (ops->put("-", 1, p))
736 return (EOF);
2b26f2d7 737 m->f &= ~MP_NEG;
d3409d5e 738 }
739
a951033d 740 /* --- Handle binary radix --- */
741
742 switch (radix) {
743 case 2: case -2: return (binary(m, 1, radix, ops, p));
744 case 4: case -4: return (binary(m, 2, radix, ops, p));
745 case 8: case -8: return (binary(m, 3, radix, ops, p));
746 case 16: case -16: return (binary(m, 4, radix, ops, p));
747 case 32: case -32: return (binary(m, 5, radix, ops, p));
748 case -64: return (binary(m, 6, radix, ops, p));
749 case -128: return (binary(m, 7, radix, ops, p));
750 }
751
e360a4f2 752 /* --- If the number is small, do it the easy way --- */
753
3bc9cb53 754 if (MP_LEN(m) < 2)
755 rc = simple(MP_LEN(m) ? m->v[0] : 0, radix, 0, ops, p);
e360a4f2 756
757 /* --- Use a clever algorithm --- *
758 *
759 * Square the radix repeatedly, remembering old results, until I get
760 * something more than half the size of the number @m@. Use this to divide
761 * the number: the quotient and remainder will be approximately the same
762 * size, and I'll have split them on a digit boundary, so I can just emit
763 * the quotient and remainder recursively, in order.
e360a4f2 764 */
765
766 else {
2b26f2d7 767 mp *pr[DEPTH];
3bc9cb53 768 size_t target = (MP_LEN(m) + 1) / 2;
e360a4f2 769 unsigned i = 0;
2b26f2d7 770 mp *z = mp_new(1, 0);
e360a4f2 771
772 /* --- Set up the exponent table --- */
773
2b26f2d7 774 z->v[0] = (radix > 0 ? radix : -radix);
e360a4f2 775 z->f = 0;
776 for (;;) {
2b26f2d7 777 assert(((void)"Number is too unimaginably huge", i < DEPTH));
e360a4f2 778 pr[i++] = z;
779 if (MP_LEN(z) > target)
780 break;
781 z = mp_sqr(MP_NEW, z);
782 }
d3409d5e 783
e360a4f2 784 /* --- Write out the answer --- */
d3409d5e 785
e360a4f2 786 rc = complicated(m, radix, pr, i - 1, 0, ops, p);
d3409d5e 787
e360a4f2 788 /* --- Tidy away the array --- */
d3409d5e 789
e360a4f2 790 while (i > 0)
791 mp_drop(pr[--i]);
d3409d5e 792 }
e360a4f2 793
794 /* --- Tidying up code --- */
795
796 MP_DROP(m);
797 return (rc);
d3409d5e 798}
799
800/*----- Test rig ----------------------------------------------------------*/
801
802#ifdef TEST_RIG
803
804#include <mLib/testrig.h>
805
806static int verify(dstr *v)
807{
808 int ok = 1;
809 int ib = *(int *)v[0].buf, ob = *(int *)v[2].buf;
810 dstr d = DSTR_INIT;
811 mp *m = mp_readdstr(MP_NEW, &v[1], 0, ib);
812 if (m) {
813 if (!ob) {
814 fprintf(stderr, "*** unexpected successful parse\n"
a951033d 815 "*** input [%2i] = ", ib);
2b26f2d7 816 if (ib < 0)
817 type_hex.dump(&v[1], stderr);
818 else
819 fputs(v[1].buf, stderr);
d3409d5e 820 mp_writedstr(m, &d, 10);
2b26f2d7 821 fprintf(stderr, "\n*** (value = %s)\n", d.buf);
d3409d5e 822 ok = 0;
823 } else {
824 mp_writedstr(m, &d, ob);
825 if (d.len != v[3].len || memcmp(d.buf, v[3].buf, d.len) != 0) {
826 fprintf(stderr, "*** failed read or write\n"
a951033d 827 "*** input [%2i] = ", ib);
2b26f2d7 828 if (ib < 0)
829 type_hex.dump(&v[1], stderr);
830 else
831 fputs(v[1].buf, stderr);
a951033d 832 fprintf(stderr, "\n*** output [%2i] = ", ob);
2b26f2d7 833 if (ob < 0)
834 type_hex.dump(&d, stderr);
835 else
836 fputs(d.buf, stderr);
a951033d 837 fprintf(stderr, "\n*** expected [%2i] = ", ob);
2b26f2d7 838 if (ob < 0)
839 type_hex.dump(&v[3], stderr);
840 else
841 fputs(v[3].buf, stderr);
842 fputc('\n', stderr);
d3409d5e 843 ok = 0;
844 }
845 }
846 mp_drop(m);
847 } else {
848 if (ob) {
849 fprintf(stderr, "*** unexpected parse failure\n"
2b26f2d7 850 "*** input [%i] = ", ib);
851 if (ib < 0)
852 type_hex.dump(&v[1], stderr);
853 else
854 fputs(v[1].buf, stderr);
855 fprintf(stderr, "\n*** expected [%i] = ", ob);
856 if (ob < 0)
857 type_hex.dump(&v[3], stderr);
858 else
859 fputs(v[3].buf, stderr);
860 fputc('\n', stderr);
d3409d5e 861 ok = 0;
862 }
863 }
864
865 dstr_destroy(&d);
9c3df6c0 866 assert(mparena_count(MPARENA_GLOBAL) == 0);
d3409d5e 867 return (ok);
868}
869
870static test_chunk tests[] = {
2b26f2d7 871 { "mptext-ascii", verify,
d3409d5e 872 { &type_int, &type_string, &type_int, &type_string, 0 } },
2b26f2d7 873 { "mptext-bin-in", verify,
874 { &type_int, &type_hex, &type_int, &type_string, 0 } },
875 { "mptext-bin-out", verify,
876 { &type_int, &type_string, &type_int, &type_hex, 0 } },
d3409d5e 877 { 0, 0, { 0 } }
878};
879
880int main(int argc, char *argv[])
881{
882 sub_init();
883 test_run(argc, argv, tests, SRCDIR "/tests/mptext");
884 return (0);
885}
886
887#endif
888
889/*----- That's all, folks -------------------------------------------------*/