[Python-checkins] r77237 - in python/branches/py3k: Include/longobject.h Modules/mathmodule.c Objects/longobject.c
mark.dickinson
python-checkins at python.org
Sat Jan 2 16:33:56 CET 2010
Author: mark.dickinson
Date: Sat Jan 2 16:33:56 2010
New Revision: 77237
Log:
Merged revisions 77234 via svnmerge from
svn+ssh://pythondev@svn.python.org/python/trunk
........
r77234 | mark.dickinson | 2010-01-02 14:45:40 +0000 (Sat, 02 Jan 2010) | 7 lines
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.
........
Modified:
python/branches/py3k/ (props changed)
python/branches/py3k/Include/longobject.h
python/branches/py3k/Modules/mathmodule.c
python/branches/py3k/Objects/longobject.c
Modified: python/branches/py3k/Include/longobject.h
==============================================================================
--- python/branches/py3k/Include/longobject.h (original)
+++ python/branches/py3k/Include/longobject.h Sat Jan 2 16:33:56 2010
@@ -44,13 +44,13 @@
/* For use by intobject.c only */
PyAPI_DATA(unsigned char) _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 *);
Modified: python/branches/py3k/Modules/mathmodule.c
==============================================================================
--- python/branches/py3k/Modules/mathmodule.c (original)
+++ python/branches/py3k/Modules/mathmodule.c Sat Jan 2 16:33:56 2010
@@ -54,7 +54,6 @@
#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 */
@@ -1342,11 +1341,12 @@
/* 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)
@@ -1354,18 +1354,21 @@
/* 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);
}
Modified: python/branches/py3k/Objects/longobject.c
==============================================================================
--- python/branches/py3k/Objects/longobject.c (original)
+++ python/branches/py3k/Objects/longobject.c Sat Jan 2 16:33:56 2010
@@ -98,9 +98,6 @@
#define SIGCHECK(PyTryBlock) \
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. */
@@ -911,224 +908,6 @@
}
-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,
- "Python int too large to convert to C double");
- return -1.0;
-}
-
/* Create a new long (or int) object from a C pointer */
PyObject *
@@ -2442,6 +2221,152 @@
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
More information about the Python-checkins
mailing list