3 * Textual representation of multiprecision numbers
5 * (c) 1999 Straylight/Edgeware
8 /*----- Licensing notice --------------------------------------------------*
10 * This file is part of Catacomb.
12 * Catacomb is free software; you can redistribute it and/or modify
13 * it under the terms of the GNU Library General Public License as
14 * published by the Free Software Foundation; either version 2 of the
15 * License, or (at your option) any later version.
17 * Catacomb is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20 * GNU Library General Public License for more details.
22 * You should have received a copy of the GNU Library General Public
23 * License along with Catacomb; if not, write to the Free
24 * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
28 /*----- Header files ------------------------------------------------------*/
34 #include <mLib/macros.h>
40 /*----- Magical numbers ---------------------------------------------------*/
42 /* --- Maximum recursion depth --- *
44 * This is the number of bits in a @size_t@ object. Why?
46 * To see this, let %$b = \textit{MPW\_MAX} + 1$% and let %$Z$% be the
47 * largest @size_t@ value. Then the largest possible @mp@ is %$M - 1$% where
48 * %$M = b^Z$%. Let %$r$% be a radix to read or write. Since the recursion
49 * squares the radix at each step, the highest number reached by the
50 * recursion is %$d$%, where:
54 * Solving gives that %$d = \lg \log_r b^Z$%. If %$r = 2$%, this is maximum,
55 * so choosing %$d = \lg \lg b^Z = \lg (Z \lg b) = \lg Z + \lg \lg b$%.
57 * Expressing %$\lg Z$% as @CHAR_BIT * sizeof(size_t)@ yields an
58 * overestimate, since a @size_t@ representation may contain `holes'.
59 * Choosing to represent %$\lg \lg b$% by 10 is almost certainly sufficient
60 * for `some time to come'.
63 #define DEPTH (CHAR_BIT * sizeof(size_t) + 10)
65 /*----- Input -------------------------------------------------------------*/
67 /* --- @mp_read@ --- *
69 * Arguments: @mp *m@ = destination multiprecision number
70 * @int radix@ = base to assume for data (or zero to guess)
71 * @const mptext_ops *ops@ = pointer to operations block
72 * @void *p@ = data for the operations block
74 * Returns: The integer read, or zero if it didn't work.
76 * Use: Reads an integer from some source. If the @radix@ is
77 * specified, the number is assumed to be given in that radix,
78 * with the letters `a' (either upper- or lower-case) upwards
79 * standing for digits greater than 9. Otherwise, base 10 is
80 * assumed unless the number starts with `0' (octal), `0x' (hex)
81 * or `nnn_' (base `nnn'). An arbitrary amount of whitespace
82 * before the number is ignored.
85 /* --- About the algorithm --- *
87 * The algorithm here is rather aggressive. I maintain an array of
88 * successive squarings of the radix, and a stack of partial results, each
89 * with a counter attached indicating which radix square to multiply by.
90 * Once the item at the top of the stack reaches the same counter level as
91 * the next item down, they are combined together and the result is given a
92 * counter level one higher than either of the results.
94 * Gluing the results together at the end is slightly tricky. Pay attention
97 * This is more complicated because of the need to handle the slightly
101 static int char_digit(int ch
, int radix
)
103 int r
= radix
< 0 ?
-radix
: radix
;
106 if (ch
< 0) return (-1);
107 if (radix
< 0) d
= ch
;
108 else if ('0' <= ch
&& ch
<= '9') d
= ch
- '0';
109 else if ('a' <= ch
&& ch
<= 'z') d
= ch
- 'a' + 10;
110 else if ('A' <= ch
&& ch
<= 'Z') d
= ch
- 'A' + (radix
> 36 ?
36 : 10);
112 if (d
>= r
) return (-1);
116 static mp
*read_binary(int radix
, unsigned bit
, unsigned nf
,
117 const mptext_ops
*ops
, void *p
)
120 unsigned b
= MPW_BITS
;
127 /* --- The fast binary algorithm --- *
129 * We stack bits up starting at the top end of a word. When one word is
130 * full, we write it to the integer, and start another with the left-over
131 * bits. When the array in the integer is full, we resize using low-level
132 * calls and copy the current data to the top end. Finally, we do a single
133 * bit-shift when we know where the end of the number is.
136 m
= mp_dest(MP_NEW
, 1, nf
);
143 if ((d
= char_digit(ch
, radix
)) < 0) break;
145 /* --- Ignore leading zeroes, but notice that the number is valid --- */
148 if (!d
&& !nz
) continue;
151 /* --- Feed the digit into the accumulator --- */
157 a
|= MPW(d
) >> (bit
- b
);
162 v
= mpalloc(m
->a
, len
);
163 memcpy(v
+ n
, m
->v
, MPWS(n
));
165 m
->v
= v
; v
= m
->v
+ n
;
167 a
= (b
< MPW_BITS
) ?
MPW(d
) << b
: 0;
171 /* --- Finish up --- */
174 if (!any
) { mp_drop(m
); return (0); }
180 m
= mp_lsr(m
, m
, (unsigned long)n
* MPW_BITS
+ b
);
187 /* --- State for the general-base reader --- *
189 * There are two arrays. The @pow@ array is set so that @pow[i]@ contains
190 * %$R^{2^i}$% for @i < pows@. The stack @s@ contains partial results:
191 * each entry contains a value @m@ corresponding to %$2^i$% digits.
192 * Inductively, an empty stack represents zero; if a stack represents %$x$%
193 * then pushing a new entry on the top causes the stack to represent
196 * It is an invariant that each entry has a strictly smaller @i@ than the
197 * items beneath it. This is achieved by coaslescing entries at the top if
198 * they have equal %$i$% values: if the top items are %$(m, i)$%, and
199 * %$(M', i)$%, and the rest of the stack represents the integer %$x$%,
200 * then %$R^{2^i} (R^{2^i} x + M) + m = R^{2^{i+1}} x + (R^{2^i} M + m)$%,
201 * so we replace the top two items by %$((R^{2^i} M + m), i + 1)$%, and
202 * repeat if necessary.
206 struct { unsigned i
; mp
*m
; } s
[DEPTH
];
210 static void ensure_power(struct readstate
*rs
)
212 /* --- Make sure we have the necessary %$R^{2^i}$% computed --- */
214 if (rs
->s
[rs
->sp
].i
>= rs
->pows
) {
215 assert(rs
->pows
< DEPTH
);
216 rs
->pow
[rs
->pows
] = mp_sqr(MP_NEW
, rs
->pow
[rs
->pows
- 1]);
221 static void read_digit(struct readstate
*rs
, unsigned nf
, int d
)
223 mp
*m
= mp_new(1, nf
);
226 /* --- Put the new digit on top --- */
228 assert(rs
->sp
< DEPTH
);
232 /* --- Restore the stack invariant --- */
234 while (rs
->sp
&& rs
->s
[rs
->sp
- 1].i
<= rs
->s
[rs
->sp
].i
) {
240 m
= mp_mul(m
, m
, rs
->pow
[rs
->s
[rs
->sp
+ 1].i
]);
241 m
= mp_add(m
, m
, rs
->s
[rs
->sp
+ 1].m
);
242 MP_DROP(rs
->s
[rs
->sp
+ 1].m
);
247 /* --- Leave the stack pointer at an empty item --- */
252 static mp
*read_general(int radix
, unsigned t
, unsigned nf
,
253 const mptext_ops
*ops
, void *p
)
264 /* --- Prepare the stack --- */
266 r
= radix
< 0 ?
-radix
: radix
;
267 mp_build(&rr
, &r
, &r
+ 1);
272 /* --- If we've partially parsed some input then feed it in --- *
274 * Unfortunately, what we've got is backwards. Fortunately there's a
275 * fairly tight upper bound on how many digits @t@ might be, since we
276 * aborted that loop once it got too large.
281 while (t
) { assert(i
< sizeof(v
)); v
[i
++] = t
%r
; t
/= r
; }
282 while (i
) read_digit(&rs
, nf
, v
[--i
]);
286 /* --- Read more stuff --- */
290 if ((d
= char_digit(ch
, radix
)) < 0) break;
291 read_digit(&rs
, nf
, d
); any
= 1;
295 /* --- Stitch all of the numbers together --- *
297 * This is not the same code as @read_digit@. In particular, here we must
298 * cope with the partial result being some inconvenient power of %$R$%,
299 * rather than %$R^{2^i}$%.
302 if (!any
) return (0);
303 m
= MP_ZERO
; z
= MP_ONE
;
310 z
= mp_mul(z
, z
, rs
.pow
[rs
.s
[rs
.sp
].i
]);
313 for (i
= 0; i
< rs
.pows
; i
++) MP_DROP(rs
.pow
[i
]);
318 mp
*mp_read(mp
*m
, int radix
, const mptext_ops
*ops
, void *p
)
328 /* --- We don't actually need a destination so throw it away --- *
330 * But note the flags before we lose it entirely.
338 /* --- Maintain a lookahead character --- */
342 /* --- If we're reading text, skip leading space, and maybe a sign --- */
345 while (ISSPACE(ch
)) ch
= ops
->get(p
);
347 case '-': f
|= f_neg
; /* and on */
348 case '+': do ch
= ops
->get(p
); while (ISSPACE(ch
));
352 /* --- If we don't have a fixed radix, then parse one from the input --- *
354 * This is moderately easy if the input starts with `0x' or similar. If it
355 * starts with `0' and something else, then it might be octal, or just a
356 * plain old zero. Finally, it might start with a leading `NN_', in which
357 * case we carefully collect the decimal number until we're sure it's
358 * either a radix prefix (in which case we accept it and start over) or it
359 * isn't (in which case it's actually the start of a large number we need
367 case 'x': case 'X': radix
= 16; goto fetch
;
368 case 'o': case 'O': radix
= 8; goto fetch
;
369 case 'b': case 'B': radix
= 2; goto fetch
;
370 fetch
: ch
= ops
->get(p
); break;
371 default: radix
= 8; f
|= f_ok
; break;
374 if ((d
= char_digit(ch
, 10)) < 0) { ops
->unget(ch
, p
); return (0); }
379 if ((d
= char_digit(ch
, 10)) < 0) break;
381 if (ch
!= '_' || t
> 52) radix
= 10;
389 /* --- We're now ready to dispatch to the correct handler --- */
391 rd
= radix
< 0 ?
-radix
: radix
;
394 case 2: m
= read_binary(radix
, 1, nf
, ops
, p
); break;
395 case 4: m
= read_binary(radix
, 2, nf
, ops
, p
); break;
396 case 8: m
= read_binary(radix
, 3, nf
, ops
, p
); break;
397 case 16: m
= read_binary(radix
, 4, nf
, ops
, p
); break;
398 case 32: m
= read_binary(radix
, 5, nf
, ops
, p
); break;
399 case 64: m
= read_binary(radix
, 6, nf
, ops
, p
); break;
400 case 128: m
= read_binary(radix
, 7, nf
, ops
, p
); break;
401 default: m
= read_general(radix
, t
, nf
, ops
, p
); break;
404 /* --- That didn't work --- *
406 * If we've already read something then return that. Otherwise it's an
411 if (f
& f_ok
) return (MP_ZERO
);
415 /* --- Negate the result if we should do that --- */
417 if (f
& f_neg
) m
= mp_neg(m
, m
);
419 /* --- And we're all done --- */
427 /*----- Output ------------------------------------------------------------*/
429 /* --- @mp_write@ --- *
431 * Arguments: @mp *m@ = pointer to a multi-precision integer
432 * @int radix@ = radix to use when writing the number out
433 * @const mptext_ops *ops@ = pointer to an operations block
434 * @void *p@ = data for the operations block
436 * Returns: Zero if it worked, nonzero otherwise.
438 * Use: Writes a large integer in textual form.
441 static int digit_char(int d
, int radix
)
443 if (radix
< 0) return (d
);
444 else if (d
< 10) return (d
+ '0');
445 else if (d
< 26) return (d
- 10 + 'a');
446 else return (d
- 36 + 'A');
449 /* --- Simple case --- *
451 * Use a fixed-sized buffer and single-precision arithmetic to pick off
452 * low-order digits. Put each digit in a buffer, working backwards from the
453 * end. If the buffer becomes full, recurse to get another one. Ensure that
454 * there are at least @z@ digits by writing leading zeroes if there aren't
455 * enough real digits.
458 static int write_simple(mpw n
, int radix
, unsigned z
,
459 const mptext_ops
*ops
, void *p
)
463 unsigned i
= sizeof(buf
);
464 int rd
= radix
> 0 ? radix
: -radix
;
469 buf
[--i
] = digit_char(x
, radix
);
474 rc
= write_simple(n
, radix
, z
, ops
, p
);
477 memset(zbuf
, (radix
< 0) ?
0 : '0', sizeof(zbuf
));
478 while (!rc
&& z
>= sizeof(zbuf
)) {
479 rc
= ops
->put(zbuf
, sizeof(zbuf
), p
);
482 if (!rc
&& z
) rc
= ops
->put(zbuf
, z
, p
);
484 if (!rc
) rc
= ops
->put(buf
+ i
, sizeof(buf
) - i
, p
);
489 /* --- Complicated case --- *
491 * If the number is small, fall back to the simple case above. Otherwise
492 * divide and take remainder by current large power of the radix, and emit
493 * each separately. Don't emit a zero quotient. Be very careful about
494 * leading zeroes on the remainder part, because they're deeply significant.
497 static int write_complicated(mp
*m
, int radix
, mp
**pr
,
498 unsigned i
, unsigned z
,
499 const mptext_ops
*ops
, void *p
)
506 return (write_simple(MP_LEN(m
) ? m
->v
[0] : 0, radix
, z
, ops
, p
));
509 mp_div(&q
, &m
, m
, pr
[i
]);
510 if (MP_ZEROP(q
)) d
= z
;
514 rc
= write_complicated(q
, radix
, pr
, i
- 1, z
, ops
, p
);
516 if (!rc
) rc
= write_complicated(m
, radix
, pr
, i
- 1, d
, ops
, p
);
521 /* --- Binary case --- *
523 * Special case for binary output. Goes much faster.
526 static int write_binary(mp
*m
, int bit
, int radix
,
527 const mptext_ops
*ops
, void *p
)
541 /* --- Work out where to start --- */
544 if (n
% bit
) n
+= bit
- (n
% bit
);
548 if (n
>= MP_LEN(m
)) {
555 mask
= (1 << bit
) - 1;
558 /* --- Main code --- */
567 if (v
== m
->v
) break;
569 if (b
< MPW_BITS
) x
|= a
>> b
;
572 if (!x
&& !(f
& f_out
)) continue;
574 *q
++ = digit_char(x
, radix
);
575 if (q
>= buf
+ sizeof(buf
)) {
576 if ((rc
= ops
->put(buf
, sizeof(buf
), p
)) != 0) goto done
;
583 *q
++ = digit_char(x
, radix
);
584 rc
= ops
->put(buf
, q
- buf
, p
);
593 /* --- Main driver code --- */
595 int mp_write(mp
*m
, int radix
, const mptext_ops
*ops
, void *p
)
603 if (MP_EQ(m
, MP_ZERO
))
604 return (ops
->put(radix
> 0 ?
"0" : "\0", 1, p
));
606 /* --- Set various things up --- */
611 /* --- Check the radix for sensibleness --- */
614 assert(((void)"ascii radix must be <= 62", radix
<= 62));
616 assert(((void)"binary radix must fit in a byte", -radix
<= UCHAR_MAX
));
618 assert(((void)"radix can't be zero in mp_write", 0));
620 /* --- If the number is negative, sort that out --- */
624 if (ops
->put("-", 1, p
)) return (EOF
);
628 /* --- Handle binary radix --- */
631 case 2: case -2: return (write_binary(m
, 1, radix
, ops
, p
));
632 case 4: case -4: return (write_binary(m
, 2, radix
, ops
, p
));
633 case 8: case -8: return (write_binary(m
, 3, radix
, ops
, p
));
634 case 16: case -16: return (write_binary(m
, 4, radix
, ops
, p
));
635 case 32: case -32: return (write_binary(m
, 5, radix
, ops
, p
));
636 case -64: return (write_binary(m
, 6, radix
, ops
, p
));
637 case -128: return (write_binary(m
, 7, radix
, ops
, p
));
640 /* --- If the number is small, do it the easy way --- */
643 rc
= write_simple(MP_LEN(m
) ? m
->v
[0] : 0, radix
, 0, ops
, p
);
645 /* --- Use a clever algorithm --- *
647 * Square the radix repeatedly, remembering old results, until I get
648 * something more than half the size of the number @m@. Use this to divide
649 * the number: the quotient and remainder will be approximately the same
650 * size, and I'll have split them on a digit boundary, so I can just emit
651 * the quotient and remainder recursively, in order.
655 target
= (MP_LEN(m
) + 1) / 2;
658 /* --- Set up the exponent table --- */
660 z
->v
[0] = (radix
> 0 ? radix
: -radix
);
663 assert(((void)"Number is too unimaginably huge", i
< DEPTH
));
665 if (MP_LEN(z
) > target
) break;
666 z
= mp_sqr(MP_NEW
, z
);
669 /* --- Write out the answer --- */
671 rc
= write_complicated(m
, radix
, pr
, i
- 1, 0, ops
, p
);
673 /* --- Tidy away the array --- */
675 while (i
> 0) mp_drop(pr
[--i
]);
678 /* --- Tidying up code --- */
684 /*----- Test rig ----------------------------------------------------------*/
688 #include <mLib/testrig.h>
690 static int verify(dstr
*v
)
693 int ib
= *(int *)v
[0].buf
, ob
= *(int *)v
[2].buf
;
696 mp
*m
= mp_readdstr(MP_NEW
, &v
[1], &off
, ib
);
699 fprintf(stderr
, "*** unexpected successful parse\n"
700 "*** input [%2i] = ", ib
);
702 type_hex
.dump(&v
[1], stderr
);
704 fputs(v
[1].buf
, stderr
);
705 mp_writedstr(m
, &d
, 10);
706 fprintf(stderr
, "\n*** (value = %s)\n", d
.buf
);
709 mp_writedstr(m
, &d
, ob
);
710 if (d
.len
!= v
[3].len
|| MEMCMP(d
.buf
, !=, v
[3].buf
, d
.len
)) {
711 fprintf(stderr
, "*** failed read or write\n"
712 "*** input [%2i] = ", ib
);
714 type_hex
.dump(&v
[1], stderr
);
716 fputs(v
[1].buf
, stderr
);
717 fprintf(stderr
, "\n*** output [%2i] = ", ob
);
719 type_hex
.dump(&d
, stderr
);
721 fputs(d
.buf
, stderr
);
722 fprintf(stderr
, "\n*** expected [%2i] = ", ob
);
724 type_hex
.dump(&v
[3], stderr
);
726 fputs(v
[3].buf
, stderr
);
734 fprintf(stderr
, "*** unexpected parse failure\n"
735 "*** input [%2i] = ", ib
);
737 type_hex
.dump(&v
[1], stderr
);
739 fputs(v
[1].buf
, stderr
);
740 fprintf(stderr
, "\n*** expected [%2i] = ", ob
);
742 type_hex
.dump(&v
[3], stderr
);
744 fputs(v
[3].buf
, stderr
);
750 if (v
[1].len
- off
!= v
[4].len
||
751 MEMCMP(v
[1].buf
+ off
, !=, v
[4].buf
, v
[4].len
)) {
752 fprintf(stderr
, "*** leftovers incorrect\n"
753 "*** input [%2i] = ", ib
);
755 type_hex
.dump(&v
[1], stderr
);
757 fputs(v
[1].buf
, stderr
);
758 fprintf(stderr
, "\n*** expected `%s'\n"
760 v
[4].buf
, v
[1].buf
+ off
);
765 assert(mparena_count(MPARENA_GLOBAL
) == 0);
769 static test_chunk tests
[] = {
770 { "mptext-ascii", verify
,
771 { &type_int
, &type_string
, &type_int
, &type_string
, &type_string
, 0 } },
772 { "mptext-bin-in", verify
,
773 { &type_int
, &type_hex
, &type_int
, &type_string
, &type_string
, 0 } },
774 { "mptext-bin-out", verify
,
775 { &type_int
, &type_string
, &type_int
, &type_hex
, &type_string
, 0 } },
779 int main(int argc
, char *argv
[])
782 test_run(argc
, argv
, tests
, SRCDIR
"/t/mptext");
788 /*----- That's all, folks -------------------------------------------------*/