Refactor some longobject internals: PyLong_AsDouble and _PyLong_AsScaledDouble

(the latter renamed to _PyLong_Frexp) now use the same core code.  The
exponent produced by _PyLong_Frexp now has type Py_ssize_t instead of the
previously used int, and no longer needs scaling by PyLong_SHIFT.  This
frees the math module from having to know anything about the PyLong
implementation.  This closes issue #5576.

3 files changed, 168 insertions, 240 deletions

diff --git a/Include/longobject.h b/Include/longobject.h
index 0cb8e2f..71815e5 100644
--- a/Include/longobject.h
+++ b/Include/longobject.h
@@ -33,13 +33,13 @@ PyAPI_FUNC(PyObject *) PyLong_GetInfo(void);
 #define _PyLong_FromSsize_t PyLong_FromSsize_t
 PyAPI_DATA(int) _PyLong_DigitValue[256];
 
-/* _PyLong_AsScaledDouble returns a double x and an exponent e such that
-   the true value is approximately equal to x * 2**(SHIFT*e).  e is >= 0.
-   x is 0.0 if and only if the input is 0 (in which case, e and x are both
-   zeroes).  Overflow is impossible.  Note that the exponent returned must
-   be multiplied by SHIFT!  There may not be enough room in an int to store
-   e*SHIFT directly. */
-PyAPI_FUNC(double) _PyLong_AsScaledDouble(PyObject *vv, int *e);
+/* _PyLong_Frexp returns a double x and an exponent e such that the
+   true value is approximately equal to x * 2**e.  e is >= 0.  x is
+   0.0 if and only if the input is 0 (in which case, e and x are both
+   zeroes); otherwise, 0.5 <= abs(x) < 1.0.  On overflow, which is
+   possible if the number of bits doesn't fit into a Py_ssize_t, sets
+   OverflowError and returns -1.0 for x, 0 for e. */
+PyAPI_FUNC(double) _PyLong_Frexp(PyLongObject *a, Py_ssize_t *e);
 
 PyAPI_FUNC(double) PyLong_AsDouble(PyObject *);
 PyAPI_FUNC(PyObject *) PyLong_FromVoidPtr(void *);
diff --git a/Modules/mathmodule.c b/Modules/mathmodule.c
index 6eb9e91..849f7f0 100644
--- a/Modules/mathmodule.c
+++ b/Modules/mathmodule.c
@@ -54,7 +54,6 @@ raised for division by zero and mod by zero.
 
 #include "Python.h"
 #include "_math.h"
-#include "longintrepr.h" /* just for SHIFT */
 
 #ifdef _OSF_SOURCE
 /* OSF1 5.1 doesn't make this available with XOPEN_SOURCE_EXTENDED defined */
@@ -1269,11 +1268,12 @@ PyDoc_STRVAR(math_modf_doc,
 
 /* A decent logarithm is easy to compute even for huge longs, but libm can't
    do that by itself -- loghelper can.  func is log or log10, and name is
-   "log" or "log10".  Note that overflow isn't possible:  a long can contain
-   no more than INT_MAX * SHIFT bits, so has value certainly less than
-   2**(2**64 * 2**16) == 2**2**80, and log2 of that is 2**80, which is
+   "log" or "log10".  Note that overflow of the result isn't possible: a long
+   can contain no more than INT_MAX * SHIFT bits, so has value certainly less
+   than 2**(2**64 * 2**16) == 2**2**80, and log2 of that is 2**80, which is
    small enough to fit in an IEEE single.  log and log10 are even smaller.
-*/
+   However, intermediate overflow is possible for a long if the number of bits
+   in that long is larger than PY_SSIZE_T_MAX. */
 
 static PyObject*
 loghelper(PyObject* arg, double (*func)(double), char *funcname)
@@ -1281,18 +1281,21 @@ loghelper(PyObject* arg, double (*func)(double), char *funcname)
 	/* If it is long, do it ourselves. */
 	if (PyLong_Check(arg)) {
 		double x;
-		int e;
-		x = _PyLong_AsScaledDouble(arg, &e);
+		Py_ssize_t e;
+		x = _PyLong_Frexp((PyLongObject *)arg, &e);
+		if (x == -1.0 && PyErr_Occurred())
+			return NULL;
 		if (x <= 0.0) {
 			PyErr_SetString(PyExc_ValueError,
 					"math domain error");
 			return NULL;
 		}
-		/* Value is ~= x * 2**(e*PyLong_SHIFT), so the log ~=
-		   log(x) + log(2) * e * PyLong_SHIFT.
-		   CAUTION:  e*PyLong_SHIFT may overflow using int arithmetic,
-		   so force use of double. */
-		x = func(x) + (e * (double)PyLong_SHIFT) * func(2.0);
+		/* Special case for log(1), to make sure we get an
+		   exact result there. */
+		if (e == 1 && x == 0.5)
+			return PyFloat_FromDouble(0.0);
+		/* Value is ~= x * 2**e, so the log ~= log(x) + log(2) * e. */
+		x = func(x) + func(2.0) * e;
 		return PyFloat_FromDouble(x);
 	}
 
diff --git a/Objects/longobject.c b/Objects/longobject.c
index 606b902..51442d3 100644
--- a/Objects/longobject.c
+++ b/Objects/longobject.c
@@ -39,9 +39,6 @@
 		if (PyErr_CheckSignals()) PyTryBlock \
 	}
 
-/* forward declaration */
-static int bits_in_digit(digit d);
-
 /* Normalize (remove leading zeros from) a long int object.
    Doesn't attempt to free the storage--in most cases, due to the nature
    of the algorithms used, this could save at most be one word anyway. */
@@ -750,224 +747,6 @@ Overflow:
 
 }
 
-double
-_PyLong_AsScaledDouble(PyObject *vv, int *exponent)
-{
-/* NBITS_WANTED should be > the number of bits in a double's precision,
-   but small enough so that 2**NBITS_WANTED is within the normal double
-   range.  nbitsneeded is set to 1 less than that because the most-significant
-   Python digit contains at least 1 significant bit, but we don't want to
-   bother counting them (catering to the worst case cheaply).
-
-   57 is one more than VAX-D double precision; I (Tim) don't know of a double
-   format with more precision than that; it's 1 larger so that we add in at
-   least one round bit to stand in for the ignored least-significant bits.
-*/
-#define NBITS_WANTED 57
-	PyLongObject *v;
-	double x;
-	const double multiplier = (double)(1L << PyLong_SHIFT);
-	Py_ssize_t i;
-	int sign;
-	int nbitsneeded;
-
-	if (vv == NULL || !PyLong_Check(vv)) {
-		PyErr_BadInternalCall();
-		return -1;
-	}
-	v = (PyLongObject *)vv;
-	i = Py_SIZE(v);
-	sign = 1;
-	if (i < 0) {
-		sign = -1;
-		i = -(i);
-	}
-	else if (i == 0) {
-		*exponent = 0;
-		return 0.0;
-	}
-	--i;
-	x = (double)v->ob_digit[i];
-	nbitsneeded = NBITS_WANTED - 1;
-	/* Invariant:  i Python digits remain unaccounted for. */
-	while (i > 0 && nbitsneeded > 0) {
-		--i;
-		x = x * multiplier + (double)v->ob_digit[i];
-		nbitsneeded -= PyLong_SHIFT;
-	}
-	/* There are i digits we didn't shift in.  Pretending they're all
-	   zeroes, the true value is x * 2**(i*PyLong_SHIFT). */
-	*exponent = i;
-	assert(x > 0.0);
-	return x * sign;
-#undef NBITS_WANTED
-}
-
-/* Get a C double from a long int object.  Rounds to the nearest double,
-   using the round-half-to-even rule in the case of a tie. */
-
-double
-PyLong_AsDouble(PyObject *vv)
-{
-	PyLongObject *v = (PyLongObject *)vv;
-	Py_ssize_t rnd_digit, rnd_bit, m, n;
-	digit lsb, *d;
-	int round_up = 0;
-	double x;
-
-	if (vv == NULL || !PyLong_Check(vv)) {
-		PyErr_BadInternalCall();
-		return -1.0;
-	}
-
-	/* Notes on the method: for simplicity, assume v is positive and >=
-	   2**DBL_MANT_DIG. (For negative v we just ignore the sign until the
-	   end; for small v no rounding is necessary.)  Write n for the number
-	   of bits in v, so that 2**(n-1) <= v < 2**n, and n > DBL_MANT_DIG.
-
-	   Some terminology: the *rounding bit* of v is the 1st bit of v that
-	   will be rounded away (bit n - DBL_MANT_DIG - 1); the *parity bit*
-	   is the bit immediately above.  The round-half-to-even rule says
-	   that we round up if the rounding bit is set, unless v is exactly
-	   halfway between two floats and the parity bit is zero.
-
-	   Write d[0] ... d[m] for the digits of v, least to most significant.
-	   Let rnd_bit be the index of the rounding bit, and rnd_digit the
-	   index of the PyLong digit containing the rounding bit.  Then the
-	   bits of the digit d[rnd_digit] look something like:
-
-	              rounding bit
-	                  |
-	                  v
-	      msb -> sssssrttttttttt <- lsb
-	                 ^
-	                 |
-	              parity bit
-
-	   where 's' represents a 'significant bit' that will be included in
-	   the mantissa of the result, 'r' is the rounding bit, and 't'
-	   represents a 'trailing bit' following the rounding bit.  Note that
-	   if the rounding bit is at the top of d[rnd_digit] then the parity
-	   bit will be the lsb of d[rnd_digit+1].  If we set
-
-	      lsb = 1 << (rnd_bit % PyLong_SHIFT)
-
-	   then d[rnd_digit] & (PyLong_BASE - 2*lsb) selects just the
-	   significant bits of d[rnd_digit], d[rnd_digit] & (lsb-1) gets the
-	   trailing bits, and d[rnd_digit] & lsb gives the rounding bit.
-
-	   We initialize the double x to the integer given by digits
-	   d[rnd_digit:m-1], but with the rounding bit and trailing bits of
-	   d[rnd_digit] masked out.  So the value of x comes from the top
-	   DBL_MANT_DIG bits of v, multiplied by 2*lsb.  Note that in the loop
-	   that produces x, all floating-point operations are exact (assuming
-	   that FLT_RADIX==2).  Now if we're rounding down, the value we want
-	   to return is simply
-
-	      x * 2**(PyLong_SHIFT * rnd_digit).
-
-	   and if we're rounding up, it's
-
-	      (x + 2*lsb) * 2**(PyLong_SHIFT * rnd_digit).
-
-	   Under the round-half-to-even rule, we round up if, and only
-	   if, the rounding bit is set *and* at least one of the
-	   following three conditions is satisfied:
-
-	      (1) the parity bit is set, or
-	      (2) at least one of the trailing bits of d[rnd_digit] is set, or
-	      (3) at least one of the digits d[i], 0 <= i < rnd_digit
-	         is nonzero.
-
-	   Finally, we have to worry about overflow.  If v >= 2**DBL_MAX_EXP,
-	   or equivalently n > DBL_MAX_EXP, then overflow occurs.  If v <
-	   2**DBL_MAX_EXP then we're usually safe, but there's a corner case
-	   to consider: if v is very close to 2**DBL_MAX_EXP then it's
-	   possible that v is rounded up to exactly 2**DBL_MAX_EXP, and then
-	   again overflow occurs.
-	*/
-
-	if (Py_SIZE(v) == 0)
-		return 0.0;
-	m = ABS(Py_SIZE(v)) - 1;
-	d = v->ob_digit;
-	assert(d[m]);  /* v should be normalized */
-
-	/* fast path for case where 0 < abs(v) < 2**DBL_MANT_DIG */
-	if (m < DBL_MANT_DIG / PyLong_SHIFT ||
-	    (m == DBL_MANT_DIG / PyLong_SHIFT &&
-	     d[m] < (digit)1 << DBL_MANT_DIG%PyLong_SHIFT)) {
-		x = d[m];
-		while (--m >= 0)
-			x = x*PyLong_BASE + d[m];
-		return Py_SIZE(v) < 0 ? -x : x;
-	}
-
-	/* if m is huge then overflow immediately; otherwise, compute the
-	   number of bits n in v.  The condition below implies n (= #bits) >=
-	   m * PyLong_SHIFT + 1 > DBL_MAX_EXP, hence v >= 2**DBL_MAX_EXP. */
-	if (m > (DBL_MAX_EXP-1)/PyLong_SHIFT)
-		goto overflow;
-	n = m * PyLong_SHIFT + bits_in_digit(d[m]);
-	if (n > DBL_MAX_EXP)
-		goto overflow;
-
-	/* find location of rounding bit */
-	assert(n > DBL_MANT_DIG); /* dealt with |v| < 2**DBL_MANT_DIG above */
-	rnd_bit = n - DBL_MANT_DIG - 1;
-	rnd_digit = rnd_bit/PyLong_SHIFT;
-	lsb = (digit)1 << (rnd_bit%PyLong_SHIFT);
-
-	/* Get top DBL_MANT_DIG bits of v.  Assumes PyLong_SHIFT <
-	   DBL_MANT_DIG, so we'll need bits from at least 2 digits of v. */
-	x = d[m];
-	assert(m > rnd_digit);
-	while (--m > rnd_digit)
-		x = x*PyLong_BASE + d[m];
-	x = x*PyLong_BASE + (d[m] & (PyLong_BASE-2*lsb));
-
-	/* decide whether to round up, using round-half-to-even */
-	assert(m == rnd_digit);
-	if (d[m] & lsb) { /* if (rounding bit is set) */
-		digit parity_bit;
-		if (lsb == PyLong_BASE/2)
-			parity_bit = d[m+1] & 1;
-		else
-			parity_bit = d[m] & 2*lsb;
-		if (parity_bit)
-			round_up = 1;
-		else if (d[m] & (lsb-1))
-			round_up = 1;
-		else {
-			while (--m >= 0) {
-				if (d[m]) {
-					round_up = 1;
-					break;
-				}
-			}
-		}
-	}
-
-	/* and round up if necessary */
-	if (round_up) {
-		x += 2*lsb;
-		if (n == DBL_MAX_EXP &&
-		    x == ldexp((double)(2*lsb), DBL_MANT_DIG)) {
-			/* overflow corner case */
-			goto overflow;
-		}
-	}
-
-	/* shift, adjust for sign, and return */
-	x = ldexp(x, rnd_digit*PyLong_SHIFT);
-	return Py_SIZE(v) < 0 ? -x : x;
-
-  overflow:
-	PyErr_SetString(PyExc_OverflowError,
-		"long int too large to convert to float");
-	return -1.0;
-}
-
 /* Create a new long (or int) object from a C pointer */
 
 PyObject *
@@ -2299,6 +2078,152 @@ x_divrem(PyLongObject *v1, PyLongObject *w1, PyLongObject **prem)
 	return long_normalize(a);
 }
 
+/* For a nonzero PyLong a, express a in the form x * 2**e, with 0.5 <=
+   abs(x) < 1.0 and e >= 0; return x and put e in *e.  Here x is
+   rounded to DBL_MANT_DIG significant bits using round-half-to-even.
+   If a == 0, return 0.0 and set *e = 0.  If the resulting exponent
+   e is larger than PY_SSIZE_T_MAX, raise OverflowError and return
+   -1.0. */
+
+/* attempt to define 2.0**DBL_MANT_DIG as a compile-time constant */
+#if DBL_MANT_DIG == 53
+#define EXP2_DBL_MANT_DIG 9007199254740992.0
+#else
+#define EXP2_DBL_MANT_DIG (ldexp(1.0, DBL_MANT_DIG))
+#endif
+
+double
+_PyLong_Frexp(PyLongObject *a, Py_ssize_t *e)
+{
+	Py_ssize_t a_size, a_bits, shift_digits, shift_bits, x_size;
+	/* See below for why x_digits is always large enough. */
+	digit rem, x_digits[2 + (DBL_MANT_DIG + 1) / PyLong_SHIFT];
+	double dx;
+	/* Correction term for round-half-to-even rounding.  For a digit x,
+	   "x + half_even_correction[x & 7]" gives x rounded to the nearest
+	   multiple of 4, rounding ties to a multiple of 8. */
+	static const int half_even_correction[8] = {0, -1, -2, 1, 0, -1, 2, 1};
+
+	a_size = ABS(Py_SIZE(a));
+	if (a_size == 0) {
+		/* Special case for 0: significand 0.0, exponent 0. */
+		*e = 0;
+		return 0.0;
+	}
+	a_bits = bits_in_digit(a->ob_digit[a_size-1]);
+	/* The following is an overflow-free version of the check
+	   "if ((a_size - 1) * PyLong_SHIFT + a_bits > PY_SSIZE_T_MAX) ..." */
+	if (a_size >= (PY_SSIZE_T_MAX - 1) / PyLong_SHIFT + 1 &&
+	    (a_size > (PY_SSIZE_T_MAX - 1) / PyLong_SHIFT + 1 ||
+	     a_bits > (PY_SSIZE_T_MAX - 1) % PyLong_SHIFT + 1))
+		 goto overflow;
+	a_bits = (a_size - 1) * PyLong_SHIFT + a_bits;
+
+	/* Shift the first DBL_MANT_DIG + 2 bits of a into x_digits[0:x_size]
+	   (shifting left if a_bits <= DBL_MANT_DIG + 2).
+
+	   Number of digits needed for result: write // for floor division.
+	   Then if shifting left, we end up using
+
+	     1 + a_size + (DBL_MANT_DIG + 2 - a_bits) // PyLong_SHIFT
+
+	   digits.  If shifting right, we use
+
+	     a_size - (a_bits - DBL_MANT_DIG - 2) // PyLong_SHIFT
+
+	   digits.  Using a_size = 1 + (a_bits - 1) // PyLong_SHIFT along with
+	   the inequalities
+
+	     m // PyLong_SHIFT + n // PyLong_SHIFT <= (m + n) // PyLong_SHIFT
+	     m // PyLong_SHIFT - n // PyLong_SHIFT <=
+	                                      1 + (m - n - 1) // PyLong_SHIFT,
+
+	   valid for any integers m and n, we find that x_size satisfies
+
+	     x_size <= 2 + (DBL_MANT_DIG + 1) // PyLong_SHIFT
+
+	   in both cases.
+	*/
+	if (a_bits <= DBL_MANT_DIG + 2) {
+		shift_digits = (DBL_MANT_DIG + 2 - a_bits) / PyLong_SHIFT;
+		shift_bits = (DBL_MANT_DIG + 2 - a_bits) % PyLong_SHIFT;
+		x_size = 0;
+		while (x_size < shift_digits)
+			x_digits[x_size++] = 0;
+		rem = v_lshift(x_digits + x_size, a->ob_digit, a_size,
+			       shift_bits);
+		x_size += a_size;
+		x_digits[x_size++] = rem;
+	}
+	else {
+		shift_digits = (a_bits - DBL_MANT_DIG - 2) / PyLong_SHIFT;
+		shift_bits = (a_bits - DBL_MANT_DIG - 2) % PyLong_SHIFT;
+		rem = v_rshift(x_digits, a->ob_digit + shift_digits,
+			       a_size - shift_digits, shift_bits);
+		x_size = a_size - shift_digits;
+		/* For correct rounding below, we need the least significant
+		   bit of x to be 'sticky' for this shift: if any of the bits
+		   shifted out was nonzero, we set the least significant bit
+		   of x. */
+		if (rem)
+			x_digits[0] |= 1;
+		else
+			while (shift_digits > 0)
+				if (a->ob_digit[--shift_digits]) {
+					x_digits[0] |= 1;
+					break;
+				}
+	}
+	assert(1 <= x_size && x_size <= sizeof(x_digits)/sizeof(digit));
+
+	/* Round, and convert to double. */
+	x_digits[0] += half_even_correction[x_digits[0] & 7];
+	dx = x_digits[--x_size];
+	while (x_size > 0)
+		dx = dx * PyLong_BASE + x_digits[--x_size];
+
+	/* Rescale;  make correction if result is 1.0. */
+	dx /= 4.0 * EXP2_DBL_MANT_DIG;
+	if (dx == 1.0) {
+		if (a_bits == PY_SSIZE_T_MAX)
+			goto overflow;
+		dx = 0.5;
+		a_bits += 1;
+	}
+
+	*e = a_bits;
+	return Py_SIZE(a) < 0 ? -dx : dx;
+
+  overflow:
+	/* exponent > PY_SSIZE_T_MAX */
+	PyErr_SetString(PyExc_OverflowError,
+			"huge integer: number of bits overflows a Py_ssize_t");
+	*e = 0;
+	return -1.0;
+}
+
+/* Get a C double from a long int object.  Rounds to the nearest double,
+   using the round-half-to-even rule in the case of a tie. */
+
+double
+PyLong_AsDouble(PyObject *v)
+{
+	Py_ssize_t exponent;
+	double x;
+
+	if (v == NULL || !PyLong_Check(v)) {
+		PyErr_BadInternalCall();
+		return -1.0;
+	}
+	x = _PyLong_Frexp((PyLongObject *)v, &exponent);
+	if ((x == -1.0 && PyErr_Occurred()) || exponent > DBL_MAX_EXP) {
+		PyErr_SetString(PyExc_OverflowError,
+				"long int too large to convert to float");
+		return -1.0;
+	}
+	return ldexp(x, exponent);
+}
+
 /* Methods */
 
 static void