chiark / gitweb /
rand/rand-x86ish.S: Hoist argument register allocation outside.
[catacomb] / math / mptext-len.c
1 /* -*-c-*-
2  *
3  * Work out length of a number's string representation
4  *
5  * (c) 2002 Straylight/Edgeware
6  */
7
8 /*----- Licensing notice --------------------------------------------------*
9  *
10  * This file is part of Catacomb.
11  *
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.
16  *
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.
21  *
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,
25  * MA 02111-1307, USA.
26  */
27
28 /*----- Header files ------------------------------------------------------*/
29
30 #include "mp.h"
31 #include "mptext.h"
32
33 /*----- Main code ---------------------------------------------------------*/
34
35 /* --- @mptext_len@ --- *
36  *
37  * Arguments:   @mp *x@ = number to work on
38  *              @int r@ = radix the number will be expressed in
39  *
40  * Returns:     The number of digits needed to represent the number in the
41  *              given base.  This will not include space for a leading sign
42  *              (use @MP_NEGP@ to check that, or just add one on for luck);
43  *              neither will it add space for a terminating null.  In general
44  *              the answer will be an overestimate.
45  */
46
47 size_t mptext_len(mp *x, int r)
48 {
49   unsigned long b = mp_bits(x);
50   int s, ss = 2;
51   size_t n;
52   unsigned d = 0;
53
54   /* --- Huh? --- *
55    *
56    * The number of digits is at most %$\lceil b \log 2/\log r \rceil$%.  We
57    * produce an underestimate of %$\log_2 r = \log r/\log 2$% and divide by
58    * that.  How?  By linear interpolation between known points on the curve.
59    * The known points are precisely the powers of 2, so we can find a pair
60    * efficiently by doubling up.  The log curve is convex, so linear
61    * interpolation between points on the curve is always an underestimate.
62    *
63    * The integer maths here is a bit weird, so here's how it works.  If
64    * %$s = 2^d$% is the power of 2 below %$r$% then we want to compute
65    * %$\lceil b/(d + (r - s)/s) \rceil = \lceil (b s)/(s(d - 1) + r \rceil$%
66    * which is %$\lfloor (r + s (b + d - 1) - 1)/(r + s(d - 1)) \rfloor$%.
67    * Gluing the whole computation together like this makes the code hard to
68    * read, but means that there are fewer possibilities for rounding errors
69    * and thus we get a tighter bound.
70    */
71
72   /* --- Find the right pair of points --- */
73
74   if (r < 0) r = -r;
75   do {
76     s = ss;
77     d++;
78     if (r == s) {
79       n = (b + (d - 1))/d;
80       goto done;
81     }
82     ss = s << 1;
83   } while (ss <= r);
84
85   /* --- Do the interpolation --- */
86
87   n = (r + s*(b + d - 1) - 1)/(r + s*(d - 1));
88
89   /* --- Fixups --- */
90
91 done:
92   if (!n)
93     n = 1;
94   return (n);
95 }
96
97 /*----- That's all, folks -------------------------------------------------*/