[u/mdw/catacomb] / mptext-len.c

/* -*-c-*-
 *
 * $Id$
 *
 * Work out length of a number's string representation
 *
 * (c) 2002 Straylight/Edgeware
 */

/*----- Licensing notice --------------------------------------------------*
 *
 * This file is part of Catacomb.
 *
 * Catacomb is free software; you can redistribute it and/or modify
 * it under the terms of the GNU Library General Public License as
 * published by the Free Software Foundation; either version 2 of the
 * License, or (at your option) any later version.
 *
 * Catacomb is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU Library General Public License for more details.
 *
 * You should have received a copy of the GNU Library General Public
 * License along with Catacomb; if not, write to the Free
 * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
 * MA 02111-1307, USA.
 */

/*----- Header files ------------------------------------------------------*/

#include "mp.h"
#include "mptext.h"

/*----- Main code ---------------------------------------------------------*/

/* --- @mptext_len@ --- *
 *
 * Arguments:	@mp *x@ = number to work on
 *		@int r@ = radix the number will be expressed in
 *
 * Returns:	The number of digits needed to represent the number in the
 *		given base.  This will not include space for a leading sign
 *		(use @MP_NEGP@ to check that, or just add one on for luck);
 *		neither will it add space for a terminating null.  In general
 *		the answer will be an overestimate.
 */

size_t mptext_len(mp *x, int r)
{
  unsigned long b = mp_bits(x);
  int s, ss = 2;
  size_t n;
  unsigned d = 0;

  /* --- Huh? --- *
   *
   * The number of digits is at most %$\lceil b \log 2/\log r \rceil$%.  We
   * produce an underestimate of %$\log_2 r = \log r/\log 2$% and divide by
   * that.  How?  By linear interpolation between known points on the curve.
   * The known points are precisely the powers of 2, so we can find a pair
   * efficiently by doubling up.  The log curve is convex, so linear
   * interpolation between points on the curve is always an underestimate.
   *
   * The integer maths here is a bit weird, so here's how it works.  If
   * %$s = 2^d$% is the power of 2 below %$r$% then we want to compute
   * %$\lceil b/(d + (r - s)/s) \rceil = \lceil (b s)/(s(d - 1) + r \rceil$%
   * which is %$\lfloor (r + s (b + d - 1) - 1)/(r + s(d - 1)) \rfloor$%.
   * Gluing the whole computation together like this makes the code hard to
   * read, but means that there are fewer possibilities for rounding errors
   * and thus we get a tighter bound.
   */

  /* --- Find the right pair of points --- */

  if (r < 0) r = -r;
  do {
    s = ss;
    d++;
    if (r == s) {
      n = (b + (d - 1))/d;
      goto done;
    }
    ss = s << 1;
  } while (ss <= r);

  /* --- Do the interpolation --- */

  n = (r + s*(b + d - 1) - 1)/(r + s*(d - 1));

  /* --- Fixups --- */

done:
  if (!n)
    n = 1;
  return (n);
}

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