3 * $Id: mptext.c,v 1.5 2000/06/17 11:46:19 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.5 2000/06/17 11:46:19 mdw
34 * New and much faster stack-based algorithm for reading integers. Support
35 * reading and writing binary integers in bases between 2 and 256.
37 * Revision 1.4 1999/12/22 15:56:56 mdw
38 * Use clever recursive algorithm for writing numbers out.
40 * Revision 1.3 1999/12/10 23:23:26 mdw
41 * Allocate slightly less memory.
43 * Revision 1.2 1999/11/20 22:24:15 mdw
44 * Use function versions of MPX_UMULN and MPX_UADDN.
46 * Revision 1.1 1999/11/17 18:02:16 mdw
47 * New multiprecision integer arithmetic suite.
51 /*----- Header files ------------------------------------------------------*/
61 /*----- Magical numbers ---------------------------------------------------*/
63 /* --- Maximum recursion depth --- *
65 * This is the number of bits in a @size_t@ object. Why?
67 * Just to convince yourself that this is correct: let @b = MPW_MAX + 1@.
68 * Then the largest possible @mp@ is %$M - 1$% where %$M = b^Z$%. Let %$r$%
69 * be a radix to read or write. Since the recursion squares the radix at
70 * each step, the highest number reached by the recursion is %$d$%, where:
74 * Solving gives that %$d = \lg \log_r b^Z$%. If %$r = 2$%, this is maximum,
75 * so choosing %$d = \lg \lg b^Z = \lg (Z \lg b) = \lg Z + \lg \lg b$%.
77 * Expressing %$\lg Z$% as @CHAR_BIT * sizeof(size_t)@ yields an
78 * overestimate, since a @size_t@ representation may contain `holes'.
79 * Choosing to represent %$\lg \lg b$% by 10 is almost certainly sufficient
80 * for `some time to come'.
83 #define DEPTH (CHAR_BIT * sizeof(size_t) + 10)
85 /*----- Main code ---------------------------------------------------------*/
87 /* --- @mp_read@ --- *
89 * Arguments: @mp *m@ = destination multiprecision number
90 * @int radix@ = base to assume for data (or zero to guess)
91 * @const mptext_ops *ops@ = pointer to operations block
92 * @void *p@ = data for the operations block
94 * Returns: The integer read, or zero if it didn't work.
96 * Use: Reads an integer from some source. If the @radix@ is
97 * specified, the number is assumed to be given in that radix,
98 * with the letters `a' (either upper- or lower-case) upwards
99 * standing for digits greater than 9. Otherwise, base 10 is
100 * assumed unless the number starts with `0' (octal), `0x' (hex)
101 * or `nnn_' (base `nnn'). An arbitrary amount of whitespace
102 * before the number is ignored.
105 /* --- About the algorithm --- *
107 * The algorithm here is rather aggressive. I maintain an array of
108 * successive squarings of the radix, and a stack of partial results, each
109 * with a counter attached indicating which radix square to multiply by.
110 * Once the item at the top of the stack reaches the same counter level as
111 * the next item down, they are combined together and the result is given a
112 * counter level one higher than either of the results.
114 * Gluing the results together at the end is slightly tricky. Pay attention
117 * This is more complicated because of the need to handle the slightly
121 mp
*mp_read(mp
*m
, int radix
, const mptext_ops
*ops
, void *p
)
123 int ch
; /* Current char being considered */
124 unsigned f
= 0; /* Flags about the current number */
125 int r
; /* Radix to switch over to */
126 mpw rd
; /* Radix as an @mp@ digit */
127 mp rr
; /* The @mp@ for the radix */
128 unsigned nf
= m ? m
->f
& MP_BURN
: 0; /* New @mp@ flags */
132 mp
*pow
[DEPTH
]; /* List of powers */
133 unsigned pows
; /* Next index to fill */
134 struct { unsigned i
; mp
*m
; } s
[DEPTH
]; /* Main stack */
135 unsigned sp
; /* Current stack pointer */
144 /* --- Initialize the stacks --- */
146 mp_build(&rr
, &rd
, &rd
+ 1);
152 /* --- Initialize the destination number --- */
157 /* --- Read an initial character --- */
163 /* --- Handle an initial sign --- */
172 /* --- If the radix is zero, look for leading zeros --- */
175 assert(((void)"ascii radix must be <= 36", radix
<= 36));
178 } else if (radix
< 0) {
180 assert(((void)"binary radix must fit in a byte ", rd
< UCHAR_MAX
));
182 } else if (ch
!= '0') {
197 /* --- Time to start --- */
199 for (;; ch
= ops
->get(p
)) {
202 /* --- An underscore indicates a numbered base --- */
204 if (ch
== '_' && r
> 0 && r
<= 36) {
207 /* --- Clear out the stacks --- */
209 for (i
= 1; i
< pows
; i
++)
212 for (i
= 0; i
< sp
; i
++)
216 /* --- Restart the search --- */
224 /* --- Check that the character is a digit and in range --- */
231 if (ch
>= '0' && ch
<= '9')
235 if (ch
>= 'a' && ch
<= 'z') /* ASCII dependent! */
242 /* --- Sort out what to do with the character --- */
244 if (x
>= 10 && r
>= 0)
252 /* --- Stick the character on the end of my integer --- */
254 assert(((void)"Number is too unimaginably huge", sp
< DEPTH
));
255 s
[sp
].m
= m
= mp_new(1, nf
);
259 /* --- Now grind through the stack --- */
261 while (sp
> 0 && s
[sp
- 1].i
== s
[sp
].i
) {
263 /* --- Combine the top two items --- */
267 m
= mp_mul(m
, m
, pow
[s
[sp
].i
]);
268 m
= mp_add(m
, m
, s
[sp
+ 1].m
);
270 MP_DROP(s
[sp
+ 1].m
);
273 /* --- Make a new radix power if necessary --- */
275 if (s
[sp
].i
>= pows
) {
276 assert(((void)"Number is too unimaginably huge", pows
< DEPTH
));
277 pow
[pows
] = mp_sqr(MP_NEW
, pow
[pows
- 1]);
287 /* --- If we're done, compute the rest of the number --- */
298 /* --- Combine the top two items --- */
302 z
= mp_mul(z
, z
, pow
[s
[sp
+ 1].i
]);
304 m
= mp_add(m
, m
, s
[sp
+ 1].m
);
306 MP_DROP(s
[sp
+ 1].m
);
308 /* --- Make a new radix power if necessary --- */
310 if (s
[sp
].i
>= pows
) {
311 assert(((void)"Number is too unimaginably huge", pows
< DEPTH
));
312 pow
[pows
] = mp_sqr(MP_NEW
, pow
[pows
- 1]);
321 for (i
= 0; i
< sp
; i
++)
325 /* --- Clear the radix power list --- */
329 for (i
= 1; i
< pows
; i
++)
333 /* --- Bail out if the number was bad --- */
338 /* --- Set the sign and return --- */
345 /* --- @mp_write@ --- *
347 * Arguments: @mp *m@ = pointer to a multi-precision integer
348 * @int radix@ = radix to use when writing the number out
349 * @const mptext_ops *ops@ = pointer to an operations block
350 * @void *p@ = data for the operations block
352 * Returns: Zero if it worked, nonzero otherwise.
354 * Use: Writes a large integer in textual form.
357 /* --- Simple case --- *
359 * Use a fixed-sized buffer and the simple single-precision division
360 * algorithm to pick off low-order digits. Put each digit in a buffer,
361 * working backwards from the end. If the buffer becomes full, recurse to
362 * get another one. Ensure that there are at least @z@ digits by writing
363 * leading zeroes if there aren't enough real digits.
366 static int simple(mp
*m
, int radix
, unsigned z
,
367 const mptext_ops
*ops
, void *p
)
371 unsigned i
= sizeof(buf
);
372 int rd
= radix
> 0 ? radix
: -radix
;
378 x
= mpx_udivn(m
->v
, m
->vl
, m
->v
, m
->vl
, rd
);
391 } while (i
&& MP_LEN(m
));
394 rc
= simple(m
, radix
, z
, ops
, p
);
396 static const char zero
[32] = "00000000000000000000000000000000";
397 while (!rc
&& z
>= sizeof(zero
)) {
398 rc
= ops
->put(zero
, sizeof(zero
), p
);
402 rc
= ops
->put(zero
, z
, p
);
405 ops
->put(buf
+ i
, sizeof(buf
) - i
, p
);
411 /* --- Complicated case --- *
413 * If the number is small, fall back to the simple case above. Otherwise
414 * divide and take remainder by current large power of the radix, and emit
415 * each separately. Don't emit a zero quotient. Be very careful about
416 * leading zeroes on the remainder part, because they're deeply significant.
419 static int complicated(mp
*m
, int radix
, mp
**pr
, unsigned i
, unsigned z
,
420 const mptext_ops
*ops
, void *p
)
427 return (simple(m
, radix
, z
, ops
, p
));
429 mp_div(&q
, &m
, m
, pr
[i
]);
437 rc
= complicated(q
, radix
, pr
, i
- 1, z
, ops
, p
);
440 rc
= complicated(m
, radix
, pr
, i
- 1, d
, ops
, p
);
445 /* --- Main driver code --- */
447 int mp_write(mp
*m
, int radix
, const mptext_ops
*ops
, void *p
)
451 /* --- Set various things up --- */
456 /* --- Check the radix for sensibleness --- */
459 assert(((void)"ascii radix must be <= 36", radix
<= 36));
461 assert(((void)"binary radix must fit in a byte", -radix
< UCHAR_MAX
));
463 assert(((void)"radix can't be zero in mp_write", 0));
465 /* --- If the number is negative, sort that out --- */
468 if (ops
->put("-", 1, p
))
473 /* --- If the number is small, do it the easy way --- */
476 rc
= simple(m
, radix
, 0, ops
, p
);
478 /* --- Use a clever algorithm --- *
480 * Square the radix repeatedly, remembering old results, until I get
481 * something more than half the size of the number @m@. Use this to divide
482 * the number: the quotient and remainder will be approximately the same
483 * size, and I'll have split them on a digit boundary, so I can just emit
484 * the quotient and remainder recursively, in order.
489 size_t target
= MP_LEN(m
) / 2;
491 mp
*z
= mp_new(1, 0);
493 /* --- Set up the exponent table --- */
495 z
->v
[0] = (radix
> 0 ? radix
: -radix
);
498 assert(((void)"Number is too unimaginably huge", i
< DEPTH
));
500 if (MP_LEN(z
) > target
)
502 z
= mp_sqr(MP_NEW
, z
);
505 /* --- Write out the answer --- */
507 rc
= complicated(m
, radix
, pr
, i
- 1, 0, ops
, p
);
509 /* --- Tidy away the array --- */
515 /* --- Tidying up code --- */
521 /*----- Test rig ----------------------------------------------------------*/
525 #include <mLib/testrig.h>
527 static int verify(dstr
*v
)
530 int ib
= *(int *)v
[0].buf
, ob
= *(int *)v
[2].buf
;
532 mp
*m
= mp_readdstr(MP_NEW
, &v
[1], 0, ib
);
535 fprintf(stderr
, "*** unexpected successful parse\n"
536 "*** input [%i] = ", ib
);
538 type_hex
.dump(&v
[1], stderr
);
540 fputs(v
[1].buf
, stderr
);
541 mp_writedstr(m
, &d
, 10);
542 fprintf(stderr
, "\n*** (value = %s)\n", d
.buf
);
545 mp_writedstr(m
, &d
, ob
);
546 if (d
.len
!= v
[3].len
|| memcmp(d
.buf
, v
[3].buf
, d
.len
) != 0) {
547 fprintf(stderr
, "*** failed read or write\n"
548 "*** input [%i] = ", ib
);
550 type_hex
.dump(&v
[1], stderr
);
552 fputs(v
[1].buf
, stderr
);
553 fprintf(stderr
, "\n*** output [%i] = ", ob
);
555 type_hex
.dump(&d
, stderr
);
557 fputs(d
.buf
, stderr
);
558 fprintf(stderr
, "\n*** expected [%i] = ", ob
);
560 type_hex
.dump(&v
[3], stderr
);
562 fputs(v
[3].buf
, stderr
);
570 fprintf(stderr
, "*** unexpected parse failure\n"
571 "*** input [%i] = ", ib
);
573 type_hex
.dump(&v
[1], stderr
);
575 fputs(v
[1].buf
, stderr
);
576 fprintf(stderr
, "\n*** expected [%i] = ", ob
);
578 type_hex
.dump(&v
[3], stderr
);
580 fputs(v
[3].buf
, stderr
);
587 assert(mparena_count(MPARENA_GLOBAL
) == 0);
591 static test_chunk tests
[] = {
592 { "mptext-ascii", verify
,
593 { &type_int
, &type_string
, &type_int
, &type_string
, 0 } },
594 { "mptext-bin-in", verify
,
595 { &type_int
, &type_hex
, &type_int
, &type_string
, 0 } },
596 { "mptext-bin-out", verify
,
597 { &type_int
, &type_string
, &type_int
, &type_hex
, 0 } },
601 int main(int argc
, char *argv
[])
604 test_run(argc
, argv
, tests
, SRCDIR
"/tests/mptext");
610 /*----- That's all, folks -------------------------------------------------*/