3 * $Id: mptext.c,v 1.17 2002/10/19 11:59:04 mdw Exp $
5 * Textual representation of multiprecision numbers
7 * (c) 1999 Straylight/Edgeware
10 /*----- Licensing notice --------------------------------------------------*
12 * This file is part of Catacomb.
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.
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.
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,
30 /*----- Revision history --------------------------------------------------*
33 * Revision 1.17 2002/10/19 11:59:04 mdw
34 * Fix leftovers bug in reading.
36 * Revision 1.16 2002/10/15 22:57:43 mdw
37 * Bug fix: prevent negative zero.
39 * Revision 1.15 2002/10/15 19:18:15 mdw
40 * Fix fencepost bugs in binary radix writing.
42 * Revision 1.14 2002/10/09 00:33:44 mdw
43 * Allow `0o' and `0b' prefixes for octal and binary (from Haskell)
45 * Revision 1.13 2002/10/09 00:21:06 mdw
46 * Allow user-specified `r_xx' bases to be up to 62.
48 * Revision 1.12 2002/01/13 19:51:18 mdw
49 * Extend the textual format to bases up to 62 by distinguishing case.
51 * Revision 1.11 2001/06/16 23:42:17 mdw
54 * Revision 1.10 2001/06/16 13:22:39 mdw
55 * Added fast-track code for binary output bases, and tests.
57 * Revision 1.9 2001/02/03 16:05:17 mdw
58 * Make flags be unsigned. Improve the write algorithm: recurse until the
59 * parts are one word long and use single-precision arithmetic from there.
60 * Fix off-by-one bug when breaking the number apart.
62 * Revision 1.8 2000/12/06 20:32:42 mdw
63 * Reduce binary bytes (to allow marker bits to be ignored). Fix error
64 * message string a bit. Allow leading `+' signs.
66 * Revision 1.7 2000/07/15 10:01:08 mdw
67 * Bug fix in binary input.
69 * Revision 1.6 2000/06/25 12:58:23 mdw
70 * Fix the derivation of `depth' commentary.
72 * Revision 1.5 2000/06/17 11:46:19 mdw
73 * New and much faster stack-based algorithm for reading integers. Support
74 * reading and writing binary integers in bases between 2 and 256.
76 * Revision 1.4 1999/12/22 15:56:56 mdw
77 * Use clever recursive algorithm for writing numbers out.
79 * Revision 1.3 1999/12/10 23:23:26 mdw
80 * Allocate slightly less memory.
82 * Revision 1.2 1999/11/20 22:24:15 mdw
83 * Use function versions of MPX_UMULN and MPX_UADDN.
85 * Revision 1.1 1999/11/17 18:02:16 mdw
86 * New multiprecision integer arithmetic suite.
90 /*----- Header files ------------------------------------------------------*/
100 /*----- Magical numbers ---------------------------------------------------*/
102 /* --- Maximum recursion depth --- *
104 * This is the number of bits in a @size_t@ object. Why?
106 * To see this, let %$b = \textit{MPW\_MAX} + 1$% and let %$Z$% be the
107 * largest @size_t@ value. Then the largest possible @mp@ is %$M - 1$% where
108 * %$M = b^Z$%. Let %$r$% be a radix to read or write. Since the recursion
109 * squares the radix at each step, the highest number reached by the
110 * recursion is %$d$%, where:
114 * Solving gives that %$d = \lg \log_r b^Z$%. If %$r = 2$%, this is maximum,
115 * so choosing %$d = \lg \lg b^Z = \lg (Z \lg b) = \lg Z + \lg \lg b$%.
117 * Expressing %$\lg Z$% as @CHAR_BIT * sizeof(size_t)@ yields an
118 * overestimate, since a @size_t@ representation may contain `holes'.
119 * Choosing to represent %$\lg \lg b$% by 10 is almost certainly sufficient
120 * for `some time to come'.
123 #define DEPTH (CHAR_BIT * sizeof(size_t) + 10)
125 /*----- Main code ---------------------------------------------------------*/
127 /* --- @mp_read@ --- *
129 * Arguments: @mp *m@ = destination multiprecision number
130 * @int radix@ = base to assume for data (or zero to guess)
131 * @const mptext_ops *ops@ = pointer to operations block
132 * @void *p@ = data for the operations block
134 * Returns: The integer read, or zero if it didn't work.
136 * Use: Reads an integer from some source. If the @radix@ is
137 * specified, the number is assumed to be given in that radix,
138 * with the letters `a' (either upper- or lower-case) upwards
139 * standing for digits greater than 9. Otherwise, base 10 is
140 * assumed unless the number starts with `0' (octal), `0x' (hex)
141 * or `nnn_' (base `nnn'). An arbitrary amount of whitespace
142 * before the number is ignored.
145 /* --- About the algorithm --- *
147 * The algorithm here is rather aggressive. I maintain an array of
148 * successive squarings of the radix, and a stack of partial results, each
149 * with a counter attached indicating which radix square to multiply by.
150 * Once the item at the top of the stack reaches the same counter level as
151 * the next item down, they are combined together and the result is given a
152 * counter level one higher than either of the results.
154 * Gluing the results together at the end is slightly tricky. Pay attention
157 * This is more complicated because of the need to handle the slightly
161 mp
*mp_read(mp
*m
, int radix
, const mptext_ops
*ops
, void *p
)
163 int ch
; /* Current char being considered */
164 unsigned f
= 0; /* Flags about the current number */
165 int r
; /* Radix to switch over to */
166 mpw rd
; /* Radix as an @mp@ digit */
167 mp rr
; /* The @mp@ for the radix */
168 unsigned nf
= m ? m
->f
& MP_BURN
: 0; /* New @mp@ flags */
172 mp
*pow
[DEPTH
]; /* List of powers */
173 unsigned pows
; /* Next index to fill */
174 struct { unsigned i
; mp
*m
; } s
[DEPTH
]; /* Main stack */
175 unsigned sp
; /* Current stack pointer */
183 /* --- Initialize the stacks --- */
185 mp_build(&rr
, &rd
, &rd
+ 1);
191 /* --- Initialize the destination number --- */
196 /* --- Read an initial character --- */
202 /* --- Handle an initial sign --- */
204 if (radix
>= 0 && (ch
== '-' || ch
== '+')) {
207 do ch
= ops
->get(p
); while isspace(ch
);
210 /* --- If the radix is zero, look for leading zeros --- */
213 assert(((void)"ascii radix must be <= 62", radix
<= 62));
216 } else if (radix
< 0) {
218 assert(((void)"binary radix must fit in a byte", rd
< UCHAR_MAX
));
220 } else if (ch
!= '0') {
245 /* --- Use fast algorithm for binary radix --- *
247 * This is the restart point after having parsed a radix number from the
248 * input. We check whether the radix is binary, and if so use a fast
249 * algorithm which just stacks the bits up in the right order.
256 case 2: bit
= 1; goto bin
;
257 case 4: bit
= 2; goto bin
;
258 case 8: bit
= 3; goto bin
;
259 case 16: bit
= 4; goto bin
;
260 case 32: bit
= 5; goto bin
;
261 case 64: bit
= 6; goto bin
;
262 case 128: bit
= 7; goto bin
;
266 /* --- The fast binary algorithm --- *
268 * We stack bits up starting at the top end of a word. When one word is
269 * full, we write it to the integer, and start another with the left-over
270 * bits. When the array in the integer is full, we resize using low-level
271 * calls and copy the current data to the top end. Finally, we do a single
272 * bit-shift when we know where the end of the number is.
277 unsigned b
= MPW_BITS
;
281 m
= mp_dest(MP_NEW
, 1, nf
);
285 for (;; ch
= ops
->get(p
)) {
291 /* --- Check that the character is a digit and in range --- */
298 if (ch
>= '0' && ch
<= '9')
303 if (ch
>= 'a' && ch
<= 'z') /* ASCII dependent! */
305 else if (ch
>= 'A' && ch
<= 'Z')
314 /* --- Feed the digit into the accumulator --- */
317 if (!x
&& !(f
& f_start
))
324 a
|= MPW(x
) >> (bit
- b
);
331 v
= mpalloc(m
->a
, len
);
332 memcpy(v
+ n
, m
->v
, MPWS(n
));
337 a
= (b
< MPW_BITS
) ?
MPW(x
) << b
: 0;
341 /* --- Finish up --- */
352 m
= mp_lsr(m
, m
, (unsigned long)n
* MPW_BITS
+ b
);
358 /* --- Time to start --- */
360 for (;; ch
= ops
->get(p
)) {
366 /* --- An underscore indicates a numbered base --- */
368 if (ch
== '_' && r
> 0 && r
<= 62) {
371 /* --- Clear out the stacks --- */
373 for (i
= 1; i
< pows
; i
++)
376 for (i
= 0; i
< sp
; i
++)
380 /* --- Restart the search --- */
389 /* --- Check that the character is a digit and in range --- */
396 if (ch
>= '0' && ch
<= '9')
401 if (ch
>= 'a' && ch
<= 'z') /* ASCII dependent! */
403 else if (ch
>= 'A' && ch
<= 'Z')
410 /* --- Sort out what to do with the character --- */
412 if (x
>= 10 && r
>= 0)
420 /* --- Stick the character on the end of my integer --- */
422 assert(((void)"Number is too unimaginably huge", sp
< DEPTH
));
423 s
[sp
].m
= m
= mp_new(1, nf
);
427 /* --- Now grind through the stack --- */
429 while (sp
> 0 && s
[sp
- 1].i
== s
[sp
].i
) {
431 /* --- Combine the top two items --- */
435 m
= mp_mul(m
, m
, pow
[s
[sp
].i
]);
436 m
= mp_add(m
, m
, s
[sp
+ 1].m
);
438 MP_DROP(s
[sp
+ 1].m
);
441 /* --- Make a new radix power if necessary --- */
443 if (s
[sp
].i
>= pows
) {
444 assert(((void)"Number is too unimaginably huge", pows
< DEPTH
));
445 pow
[pows
] = mp_sqr(MP_NEW
, pow
[pows
- 1]);
455 /* --- If we're done, compute the rest of the number --- */
466 /* --- Combine the top two items --- */
470 z
= mp_mul(z
, z
, pow
[s
[sp
+ 1].i
]);
472 m
= mp_add(m
, m
, s
[sp
+ 1].m
);
474 MP_DROP(s
[sp
+ 1].m
);
476 /* --- Make a new radix power if necessary --- */
478 if (s
[sp
].i
>= pows
) {
479 assert(((void)"Number is too unimaginably huge", pows
< DEPTH
));
480 pow
[pows
] = mp_sqr(MP_NEW
, pow
[pows
- 1]);
489 for (i
= 0; i
< sp
; i
++)
493 /* --- Clear the radix power list --- */
497 for (i
= 1; i
< pows
; i
++)
501 /* --- Bail out if the number was bad --- */
507 /* --- Set the sign and return --- */
519 /* --- @mp_write@ --- *
521 * Arguments: @mp *m@ = pointer to a multi-precision integer
522 * @int radix@ = radix to use when writing the number out
523 * @const mptext_ops *ops@ = pointer to an operations block
524 * @void *p@ = data for the operations block
526 * Returns: Zero if it worked, nonzero otherwise.
528 * Use: Writes a large integer in textual form.
531 /* --- Simple case --- *
533 * Use a fixed-sized buffer and single-precision arithmetic to pick off
534 * low-order digits. Put each digit in a buffer, working backwards from the
535 * end. If the buffer becomes full, recurse to get another one. Ensure that
536 * there are at least @z@ digits by writing leading zeroes if there aren't
537 * enough real digits.
540 static int simple(mpw n
, int radix
, unsigned z
,
541 const mptext_ops
*ops
, void *p
)
545 unsigned i
= sizeof(buf
);
546 int rd
= radix
> 0 ? radix
: -radix
;
558 else if (x
< 36) /* Ascii specific */
568 rc
= simple(n
, radix
, z
, ops
, p
);
571 memset(zbuf
, (radix
< 0) ?
0 : '0', sizeof(zbuf
));
572 while (!rc
&& z
>= sizeof(zbuf
)) {
573 rc
= ops
->put(zbuf
, sizeof(zbuf
), p
);
577 rc
= ops
->put(zbuf
, z
, p
);
580 rc
= ops
->put(buf
+ i
, sizeof(buf
) - i
, p
);
585 /* --- Complicated case --- *
587 * If the number is small, fall back to the simple case above. Otherwise
588 * divide and take remainder by current large power of the radix, and emit
589 * each separately. Don't emit a zero quotient. Be very careful about
590 * leading zeroes on the remainder part, because they're deeply significant.
593 static int complicated(mp
*m
, int radix
, mp
**pr
, unsigned i
, unsigned z
,
594 const mptext_ops
*ops
, void *p
)
601 return (simple(MP_LEN(m
) ? m
->v
[0] : 0, radix
, z
, ops
, p
));
604 mp_div(&q
, &m
, m
, pr
[i
]);
612 rc
= complicated(q
, radix
, pr
, i
- 1, z
, ops
, p
);
615 rc
= complicated(m
, radix
, pr
, i
- 1, d
, ops
, p
);
620 /* --- Binary case --- *
622 * Special case for binary output. Goes much faster.
625 static int binary(mp
*m
, int bit
, int radix
, const mptext_ops
*ops
, void *p
)
640 /* --- Work out where to start --- */
644 n
+= bit
- (n
% bit
);
648 if (n
>= MP_LEN(m
)) {
655 mask
= (1 << bit
) - 1;
658 /* --- Main code --- */
674 if (!x
&& !(f
& f_out
))
682 ch
= 'a' + x
- 10; /* Ascii specific */
686 if (q
>= buf
+ sizeof(buf
)) {
687 if ((rc
= ops
->put(buf
, sizeof(buf
), p
)) != 0)
700 ch
= 'a' + x
- 10; /* Ascii specific */
704 rc
= ops
->put(buf
, q
- buf
, p
);
713 /* --- Main driver code --- */
715 int mp_write(mp
*m
, int radix
, const mptext_ops
*ops
, void *p
)
719 if (MP_EQ(m
, MP_ZERO
))
720 return (ops
->put("0", 1, p
));
722 /* --- Set various things up --- */
727 /* --- Check the radix for sensibleness --- */
730 assert(((void)"ascii radix must be <= 62", radix
<= 62));
732 assert(((void)"binary radix must fit in a byte", -radix
< UCHAR_MAX
));
734 assert(((void)"radix can't be zero in mp_write", 0));
736 /* --- If the number is negative, sort that out --- */
739 if (ops
->put("-", 1, p
))
744 /* --- Handle binary radix --- */
747 case 2: case -2: return (binary(m
, 1, radix
, ops
, p
));
748 case 4: case -4: return (binary(m
, 2, radix
, ops
, p
));
749 case 8: case -8: return (binary(m
, 3, radix
, ops
, p
));
750 case 16: case -16: return (binary(m
, 4, radix
, ops
, p
));
751 case 32: case -32: return (binary(m
, 5, radix
, ops
, p
));
752 case -64: return (binary(m
, 6, radix
, ops
, p
));
753 case -128: return (binary(m
, 7, radix
, ops
, p
));
756 /* --- If the number is small, do it the easy way --- */
759 rc
= simple(MP_LEN(m
) ? m
->v
[0] : 0, radix
, 0, ops
, p
);
761 /* --- Use a clever algorithm --- *
763 * Square the radix repeatedly, remembering old results, until I get
764 * something more than half the size of the number @m@. Use this to divide
765 * the number: the quotient and remainder will be approximately the same
766 * size, and I'll have split them on a digit boundary, so I can just emit
767 * the quotient and remainder recursively, in order.
772 size_t target
= (MP_LEN(m
) + 1) / 2;
774 mp
*z
= mp_new(1, 0);
776 /* --- Set up the exponent table --- */
778 z
->v
[0] = (radix
> 0 ? radix
: -radix
);
781 assert(((void)"Number is too unimaginably huge", i
< DEPTH
));
783 if (MP_LEN(z
) > target
)
785 z
= mp_sqr(MP_NEW
, z
);
788 /* --- Write out the answer --- */
790 rc
= complicated(m
, radix
, pr
, i
- 1, 0, ops
, p
);
792 /* --- Tidy away the array --- */
798 /* --- Tidying up code --- */
804 /*----- Test rig ----------------------------------------------------------*/
808 #include <mLib/testrig.h>
810 static int verify(dstr
*v
)
813 int ib
= *(int *)v
[0].buf
, ob
= *(int *)v
[2].buf
;
816 mp
*m
= mp_readdstr(MP_NEW
, &v
[1], &off
, ib
);
819 fprintf(stderr
, "*** unexpected successful parse\n"
820 "*** input [%2i] = ", ib
);
822 type_hex
.dump(&v
[1], stderr
);
824 fputs(v
[1].buf
, stderr
);
825 mp_writedstr(m
, &d
, 10);
826 fprintf(stderr
, "\n*** (value = %s)\n", d
.buf
);
829 mp_writedstr(m
, &d
, ob
);
830 if (d
.len
!= v
[3].len
|| memcmp(d
.buf
, v
[3].buf
, d
.len
) != 0) {
831 fprintf(stderr
, "*** failed read or write\n"
832 "*** input [%2i] = ", ib
);
834 type_hex
.dump(&v
[1], stderr
);
836 fputs(v
[1].buf
, stderr
);
837 fprintf(stderr
, "\n*** output [%2i] = ", ob
);
839 type_hex
.dump(&d
, stderr
);
841 fputs(d
.buf
, stderr
);
842 fprintf(stderr
, "\n*** expected [%2i] = ", ob
);
844 type_hex
.dump(&v
[3], stderr
);
846 fputs(v
[3].buf
, stderr
);
854 fprintf(stderr
, "*** unexpected parse failure\n"
855 "*** input [%2i] = ", ib
);
857 type_hex
.dump(&v
[1], stderr
);
859 fputs(v
[1].buf
, stderr
);
860 fprintf(stderr
, "\n*** expected [%2i] = ", ob
);
862 type_hex
.dump(&v
[3], stderr
);
864 fputs(v
[3].buf
, stderr
);
870 if (v
[1].len
- off
!= v
[4].len
||
871 memcmp(v
[1].buf
+ off
, v
[4].buf
, v
[4].len
) != 0) {
872 fprintf(stderr
, "*** leftovers incorrect\n"
873 "*** input [%2i] = ", ib
);
875 type_hex
.dump(&v
[1], stderr
);
877 fputs(v
[1].buf
, stderr
);
878 fprintf(stderr
, "\n*** expected `%s'\n"
880 v
[4].buf
, v
[1].buf
+ off
);
885 assert(mparena_count(MPARENA_GLOBAL
) == 0);
889 static test_chunk tests
[] = {
890 { "mptext-ascii", verify
,
891 { &type_int
, &type_string
, &type_int
, &type_string
, &type_string
, 0 } },
892 { "mptext-bin-in", verify
,
893 { &type_int
, &type_hex
, &type_int
, &type_string
, &type_string
, 0 } },
894 { "mptext-bin-out", verify
,
895 { &type_int
, &type_string
, &type_int
, &type_hex
, &type_string
, 0 } },
899 int main(int argc
, char *argv
[])
902 test_run(argc
, argv
, tests
, SRCDIR
"/tests/mptext");
908 /*----- That's all, folks -------------------------------------------------*/