I finally got the time to update and merge Mark's and my trunk-math branch. The patch is collaborated work of Mark Dickinson and me. It was mostly done a few months ago. The patch fixes a lot of loose ends and edge cases related to operations with NaN, INF, very small values and complex math.

The patch also adds acosh, asinh, atanh, log1p and copysign to all platforms. Finally it fixes differences between platforms like different results or exceptions for edge cases. Have fun :)

5 files changed, 167 insertions, 729 deletions

diff --git a/Objects/complexobject.c b/Objects/complexobject.c
index 7f40ed6..6e0fbb2 100644
--- a/Objects/complexobject.c
+++ b/Objects/complexobject.c
@@ -187,6 +187,38 @@ c_powi(Py_complex x, long n)
 
 }
 
+double
+c_abs(Py_complex z)
+{
+	/* sets errno = ERANGE on overflow;  otherwise errno = 0 */
+	double result;
+
+	if (!Py_IS_FINITE(z.real) || !Py_IS_FINITE(z.imag)) {
+		/* C99 rules: if either the real or the imaginary part is an
+		   infinity, return infinity, even if the other part is a
+		   NaN. */
+		if (Py_IS_INFINITY(z.real)) {
+			result = fabs(z.real);
+			errno = 0;
+			return result;
+		}
+		if (Py_IS_INFINITY(z.imag)) {
+			result = fabs(z.imag);
+			errno = 0;
+			return result;
+		}
+		/* either the real or imaginary part is a NaN,
+		   and neither is infinite. Result should be NaN. */
+		return Py_NAN;
+	}
+	result = hypot(z.real, z.imag);
+	if (!Py_IS_FINITE(result))
+		errno = ERANGE;
+	else
+		errno = 0;
+	return result;
+}
+
 static PyObject *
 complex_subtype_from_c_complex(PyTypeObject *type, Py_complex cval)
 {
@@ -329,8 +361,7 @@ complex_to_buf(char *buf, int bufsz, PyComplexObject *v, int precision)
 		if (!Py_IS_FINITE(v->cval.imag)) {
 			if (Py_IS_NAN(v->cval.imag))
 				strncpy(buf, "nan*j", 6);
-			/* else if (copysign(1, v->cval.imag) == 1) */
-			else if (v->cval.imag > 0)
+			else if (copysign(1, v->cval.imag) == 1)
 				strncpy(buf, "inf*j", 6);
 			else
 				strncpy(buf, "-inf*j", 7);
@@ -492,6 +523,7 @@ static PyObject *
 complex_div(PyComplexObject *v, PyComplexObject *w)
 {
 	Py_complex quot;
+
 	PyFPE_START_PROTECT("complex_div", return 0)
 	errno = 0;
 	quot = c_quot(v->cval,w->cval);
@@ -655,9 +687,16 @@ static PyObject *
 complex_abs(PyComplexObject *v)
 {
 	double result;
+
 	PyFPE_START_PROTECT("complex_abs", return 0)
-	result = hypot(v->cval.real,v->cval.imag);
+	result = c_abs(v->cval);
 	PyFPE_END_PROTECT(result)
+
+	if (errno == ERANGE) {
+		PyErr_SetString(PyExc_OverflowError,
+				"absolute value too large");
+		return NULL;
+	}
 	return PyFloat_FromDouble(result);
 }
 
@@ -786,9 +825,29 @@ complex_getnewargs(PyComplexObject *v)
 	return Py_BuildValue("(D)", &v->cval);
 }
 
+#if 0
+static PyObject *
+complex_is_finite(PyObject *self)
+{
+	Py_complex c;
+	c = ((PyComplexObject *)self)->cval;
+	return PyBool_FromLong((long)(Py_IS_FINITE(c.real) &&
+				      Py_IS_FINITE(c.imag)));
+}
+
+PyDoc_STRVAR(complex_is_finite_doc,
+"complex.is_finite() -> bool\n"
+"\n"
+"Returns True if the real and the imaginary part is finite.");
+#endif
+
 static PyMethodDef complex_methods[] = {
 	{"conjugate",	(PyCFunction)complex_conjugate,	METH_NOARGS,
 	 complex_conjugate_doc},
+#if 0
+	{"is_finite",	(PyCFunction)complex_is_finite,	METH_NOARGS,
+	 complex_is_finite_doc},
+#endif
 	{"__getnewargs__",	(PyCFunction)complex_getnewargs,	METH_NOARGS},
 	{NULL,		NULL}		/* sentinel */
 };
diff --git a/Objects/doubledigits.c b/Objects/doubledigits.c
deleted file mode 100644
index 1f1c91c..0000000
--- a/Objects/doubledigits.c
+++ /dev/null
@@ -1,601 +0,0 @@
-/* Free-format floating point printer
- * 
- * Based on "Floating-Point Printer Sample Code", by Robert G. Burger,
- * http://www.cs.indiana.edu/~burger/fp/index.html
- */
-
-#include "Python.h"
-
-#if defined(__alpha) || defined(__i386) || defined(_M_IX86) || defined(_M_X64) || defined(_M_IA64)
-#define LITTLE_ENDIAN_IEEE_DOUBLE
-#elif !(defined(__ppc__) || defined(sparc) || defined(__sgi) || defined(_IBMR2) || defined(hpux))
-#error unknown machine type
-#endif
-
-#if defined(_M_IX86)
-#define UNSIGNED64 unsigned __int64
-#elif defined(__alpha)
-#define UNSIGNED64 unsigned long
-#else
-#define UNSIGNED64 unsigned long long
-#endif
-
-#ifndef U32
-#define U32 unsigned int
-#endif
-
-/* exponent stored + 1024, hidden bit to left of decimal point */
-#define bias 1023
-#define bitstoright 52
-#define m1mask 0xf
-#define hidden_bit 0x100000
-#ifdef LITTLE_ENDIAN_IEEE_DOUBLE
-struct dblflt {
-    unsigned int m4: 16;
-    unsigned int m3: 16;
-    unsigned int m2: 16;
-    unsigned int m1: 4;
-    unsigned int e: 11;
-    unsigned int s: 1;
-};
-#else
-/* Big Endian IEEE Double Floats */
-struct dblflt {
-    unsigned int s: 1;
-    unsigned int e: 11;
-    unsigned int m1: 4;
-    unsigned int m2: 16;
-    unsigned int m3: 16;
-    unsigned int m4: 16;
-};
-#endif
-#define float_radix 2.147483648e9
-
-
-typedef UNSIGNED64 Bigit;
-#define BIGSIZE 24
-#define MIN_E -1074
-#define MAX_FIVE 325
-#define B_P1 ((Bigit)1 << 52)
-
-typedef struct {
-   int l;
-   Bigit d[BIGSIZE];
-} Bignum;
-
-static Bignum R, S, MP, MM, five[MAX_FIVE];
-static Bignum S2, S3, S4, S5, S6, S7, S8, S9;
-static int ruf, k, s_n, use_mp, qr_shift, sl, slr;
-
-static void mul10(Bignum *x);
-static void big_short_mul(Bignum *x, Bigit y, Bignum *z);
-/*
-static void print_big(Bignum *x);
-*/
-static int estimate(int n);
-static void one_shift_left(int y, Bignum *z);
-static void short_shift_left(Bigit x, int y, Bignum *z);
-static void big_shift_left(Bignum *x, int y, Bignum *z);
-static int big_comp(Bignum *x, Bignum *y);
-static int sub_big(Bignum *x, Bignum *y, Bignum *z);
-static void add_big(Bignum *x, Bignum *y, Bignum *z);
-static int add_cmp(void);
-static int qr(void);
-
-/*static int _PyFloat_Digits(char *buf, double v, int *signum);*/
-/*static void _PyFloat_DigitsInit(void);*/
-
-#define ADD(x, y, z, k) {\
-   Bigit x_add, z_add;\
-   x_add = (x);\
-   if ((k))\
-     z_add = x_add + (y) + 1, (k) = (z_add <= x_add);\
-   else\
-     z_add = x_add + (y), (k) = (z_add < x_add);\
-   (z) = z_add;\
-}
-
-#define SUB(x, y, z, b) {\
-   Bigit x_sub, y_sub;\
-   x_sub = (x); y_sub = (y);\
-   if ((b))\
-     (z) = x_sub - y_sub - 1, b = (y_sub >= x_sub);\
-   else\
-     (z) = x_sub - y_sub, b = (y_sub > x_sub);\
-}
-
-#define MUL(x, y, z, k) {\
-   Bigit x_mul, low, high;\
-   x_mul = (x);\
-   low = (x_mul & 0xffffffff) * (y) + (k);\
-   high = (x_mul >> 32) * (y) + (low >> 32);\
-   (k) = high >> 32;\
-   (z) = (low & 0xffffffff) | (high << 32);\
-}
-
-#define SLL(x, y, z, k) {\
-   Bigit x_sll = (x);\
-   (z) = (x_sll << (y)) | (k);\
-   (k) = x_sll >> (64 - (y));\
-}
-
-static void
-mul10(Bignum *x)
-{
-   int i, l;
-   Bigit *p, k;
-
-   l = x->l;
-   for (i = l, p = &x->d[0], k = 0; i >= 0; i--)
-      MUL(*p, 10, *p++, k);
-   if (k != 0)
-      *p = k, x->l = l+1;
-}
-
-static void
-big_short_mul(Bignum *x, Bigit y, Bignum *z)
-{
-   int i, xl, zl;
-   Bigit *xp, *zp, k;
-   U32 high, low;
-
-   xl = x->l;
-   xp = &x->d[0];
-   zl = xl;
-   zp = &z->d[0];
-   high = y >> 32;
-   low = y & 0xffffffff;
-   for (i = xl, k = 0; i >= 0; i--, xp++, zp++) {
-      Bigit xlow, xhigh, z0, t, c, z1;
-      xlow = *xp & 0xffffffff;
-      xhigh = *xp >> 32;
-      z0 = (xlow * low) + k; /* Cout is (z0 < k) */
-      t = xhigh * low;
-      z1 = (xlow * high) + t;
-      c = (z1 < t);
-      t = z0 >> 32;
-      z1 += t;
-      c += (z1 < t);
-      *zp = (z1 << 32) | (z0 & 0xffffffff);
-      k = (xhigh * high) + (c << 32) + (z1 >> 32) + (z0 < k);
-   }
-   if (k != 0)
-      *zp = k, zl++;
-   z->l = zl;
-}
-
-/*
-static void
-print_big(Bignum *x)
-{
-   int i;
-   Bigit *p;
-
-   printf("#x");
-   i = x->l;
-   p = &x->d[i];
-   for (p = &x->d[i]; i >= 0; i--) {
-      Bigit b = *p--;
-      printf("%08x%08x", (int)(b >> 32), (int)(b & 0xffffffff));
-   }
-}
-*/
-
-static int
-estimate(int n)
-{
-   if (n < 0)
-      return (int)(n*0.3010299956639812);
-   else
-      return 1+(int)(n*0.3010299956639811);
-}
-
-static void
-one_shift_left(int y, Bignum *z)
-{
-   int n, m, i;
-   Bigit *zp;
-
-   n = y / 64;
-   m = y % 64;
-   zp = &z->d[0];
-   for (i = n; i > 0; i--) *zp++ = 0;
-   *zp = (Bigit)1 << m;
-   z->l = n;
-}
-
-static void
-short_shift_left(Bigit x, int y, Bignum *z)
-{
-   int n, m, i, zl;
-   Bigit *zp;
-   
-   n = y / 64;
-   m = y % 64;
-   zl = n;
-   zp = &(z->d[0]);
-   for (i = n; i > 0; i--) *zp++ = 0;
-   if (m == 0)
-      *zp = x;
-   else {
-      Bigit high = x >> (64 - m);
-      *zp = x << m;
-      if (high != 0)
-         *++zp = high, zl++;
-   }
-   z->l = zl;
-}
-
-static void
-big_shift_left(Bignum *x, int y, Bignum *z)
-{
-   int n, m, i, xl, zl;
-   Bigit *xp, *zp, k;
-   
-   n = y / 64;
-   m = y % 64;
-   xl = x->l;
-   xp = &(x->d[0]);
-   zl = xl + n;
-   zp = &(z->d[0]);
-   for (i = n; i > 0; i--) *zp++ = 0;
-   if (m == 0)
-      for (i = xl; i >= 0; i--) *zp++ = *xp++;
-   else {
-      for (i = xl, k = 0; i >= 0; i--)
-         SLL(*xp++, m, *zp++, k);
-      if (k != 0)
-         *zp = k, zl++;
-   }
-   z->l = zl;
-}
-
-
-static int
-big_comp(Bignum *x, Bignum *y)
-{
-   int i, xl, yl;
-   Bigit *xp, *yp;
-
-   xl = x->l;
-   yl = y->l;
-   if (xl > yl) return 1;
-   if (xl < yl) return -1;
-   xp = &x->d[xl];
-   yp = &y->d[xl];
-   for (i = xl; i >= 0; i--, xp--, yp--) {
-      Bigit a = *xp;
-      Bigit b = *yp;
-
-      if (a > b) return 1;
-      else if (a < b) return -1;
-   }
-   return 0;
-}
-
-static int
-sub_big(Bignum *x, Bignum *y, Bignum *z)
-{
-  int xl, yl, zl, b, i;
-  Bigit *xp, *yp, *zp;
-
-  xl = x->l;
-  yl = y->l;
-  if (yl > xl) return 1;
-  xp = &x->d[0];
-  yp = &y->d[0];
-  zp = &z->d[0];
-  
-  for (i = yl, b = 0; i >= 0; i--)
-    SUB(*xp++, *yp++, *zp++, b);
-  for (i = xl-yl; b && i > 0; i--) {
-    Bigit x_sub;
-    x_sub = *xp++;
-    *zp++ = x_sub - 1;
-    b = (x_sub == 0);
-  }
-  for (; i > 0; i--) *zp++ = *xp++;
-  if (b) return 1;
-  zl = xl;
-  while (*--zp == 0) zl--;
-  z->l = zl;
-  return 0;
-}
-
-static void
-add_big(Bignum *x, Bignum *y, Bignum *z)
-{
-  int xl, yl, k, i;
-  Bigit *xp, *yp, *zp;
-
-  xl = x->l;
-  yl = y->l;
-  if (yl > xl) {
-    int tl;
-    Bignum *tn;
-    tl = xl; xl = yl; yl = tl;
-    tn = x; x = y; y = tn;
-  }
-
-  xp = &x->d[0];
-  yp = &y->d[0];
-  zp = &z->d[0];
-  
-  for (i = yl, k = 0; i >= 0; i--)
-    ADD(*xp++, *yp++, *zp++, k);
-  for (i = xl-yl; k && i > 0; i--) {
-    Bigit z_add;
-    z_add = *xp++ + 1;
-    k = (z_add == 0);
-    *zp++ = z_add;
-  }
-  for (; i > 0; i--) *zp++ = *xp++;
-  if (k)
-    *zp = 1, z->l = xl+1;
-  else
-    z->l = xl;
-}
-
-static int
-add_cmp()
-{
-   int rl, ml, sl, suml;
-   static Bignum sum;
-
-   rl = R.l;
-   ml = (use_mp ? MP.l : MM.l);
-   sl = S.l;
-   
-   suml = rl >= ml ? rl : ml;
-   if ((sl > suml+1) || ((sl == suml+1) && (S.d[sl] > 1))) return -1;
-   if (sl < suml) return 1;
-
-   add_big(&R, (use_mp ? &MP : &MM), &sum);
-   return big_comp(&sum, &S);
-}
-
-static int
-qr()
-{
-  if (big_comp(&R, &S5) < 0)
-    if (big_comp(&R, &S2) < 0)
-      if (big_comp(&R, &S) < 0)
-        return 0;
-      else {
-        sub_big(&R, &S, &R);
-        return 1;
-      }
-    else if (big_comp(&R, &S3) < 0) {
-      sub_big(&R, &S2, &R);
-      return 2;
-    }
-    else if (big_comp(&R, &S4) < 0) {
-      sub_big(&R, &S3, &R);
-      return 3;
-    }
-    else {
-      sub_big(&R, &S4, &R);
-      return 4;
-    }
-  else if (big_comp(&R, &S7) < 0)
-    if (big_comp(&R, &S6) < 0) {
-      sub_big(&R, &S5, &R);
-      return 5;
-    }
-    else {
-      sub_big(&R, &S6, &R);
-      return 6;
-    }
-  else if (big_comp(&R, &S9) < 0)
-    if (big_comp(&R, &S8) < 0) {
-      sub_big(&R, &S7, &R);
-      return 7;
-    }
-    else {
-      sub_big(&R, &S8, &R);
-      return 8;
-    }
-  else {
-    sub_big(&R, &S9, &R);
-    return 9;
-  }
-}
-
-#define OUTDIG(d) { *buf++ = (d) + '0'; *buf = 0; return k; }
-
-int
-_PyFloat_Digits(char *buf, double v, int *signum)
-{
-   struct dblflt *x;
-   int sign, e, f_n, m_n, i, d, tc1, tc2;
-   Bigit f;
-
-   /* decompose float into sign, mantissa & exponent */
-   x = (struct dblflt *)&v;
-   sign = x->s;
-   e = x->e; 
-   f = (Bigit)(x->m1 << 16 | x->m2) << 32 | (U32)(x->m3 << 16 | x->m4);
-   if (e != 0) {
-      e = e - bias - bitstoright;
-      f |= (Bigit)hidden_bit << 32;
-   }
-   else if (f != 0)
-      /* denormalized */
-      e = 1 - bias - bitstoright;
-
-   *signum = sign;
-   if (f == 0) {
-     *buf++ = '0';
-     *buf = 0;
-     return 0;
-   }
-   
-   ruf = !(f & 1); /* ruf = (even? f) */
-
-   /* Compute the scaling factor estimate, k */
-   if (e > MIN_E)
-      k = estimate(e+52);
-   else {
-      int n;
-      Bigit y;
-      
-      for (n = e+52, y = (Bigit)1 << 52; f < y; n--) y >>= 1;
-      k = estimate(n);
-   }
-
-   if (e >= 0)
-      if (f != B_P1)
-         use_mp = 0, f_n = e+1, s_n = 1, m_n = e;
-      else
-         use_mp = 1, f_n = e+2, s_n = 2, m_n = e;
-   else
-      if ((e == MIN_E) || (f != B_P1))
-         use_mp = 0, f_n = 1, s_n = 1-e, m_n = 0;
-      else
-         use_mp = 1, f_n = 2, s_n = 2-e, m_n = 0;
-   
-   /* Scale it! */
-   if (k == 0) {
-      short_shift_left(f, f_n, &R);
-      one_shift_left(s_n, &S);
-      one_shift_left(m_n, &MM);
-      if (use_mp) one_shift_left(m_n+1, &MP);
-      qr_shift = 1;
-   }
-   else if (k > 0) {
-      s_n += k;
-      if (m_n >= s_n)
-         f_n -= s_n, m_n -= s_n, s_n = 0;
-      else 
-         f_n -= m_n, s_n -= m_n, m_n = 0;
-      short_shift_left(f, f_n, &R);
-      big_shift_left(&five[k-1], s_n, &S);
-      one_shift_left(m_n, &MM);
-      if (use_mp) one_shift_left(m_n+1, &MP);
-      qr_shift = 0;
-   }
-   else {
-      Bignum *power = &five[-k-1];
-
-      s_n += k;
-      big_short_mul(power, f, &S);
-      big_shift_left(&S, f_n, &R);
-      one_shift_left(s_n, &S);
-      big_shift_left(power, m_n, &MM);
-      if (use_mp) big_shift_left(power, m_n+1, &MP);
-      qr_shift = 1;
-   }
-
-   /* fixup */
-   if (add_cmp() <= -ruf) {
-     k--;
-     mul10(&R);
-     mul10(&MM);
-     if (use_mp) mul10(&MP);
-   }
-
-   /*
-   printf("\nk = %d\n", k);
-   printf("R = "); print_big(&R);
-   printf("\nS = "); print_big(&S);
-   printf("\nM- = "); print_big(&MM);
-   if (use_mp) printf("\nM+ = "), print_big(&MP);
-   putchar('\n');
-   fflush(0);
-   */
-   
-   if (qr_shift) {
-     sl = s_n / 64;
-     slr = s_n % 64;
-   }
-   else {
-     big_shift_left(&S, 1, &S2);
-     add_big(&S2, &S, &S3);
-     big_shift_left(&S2, 1, &S4);
-     add_big(&S4, &S, &S5);
-     add_big(&S4, &S2, &S6);
-     add_big(&S4, &S3, &S7);
-     big_shift_left(&S4, 1, &S8);
-     add_big(&S8, &S, &S9);
-   }
-
-again:
-   if (qr_shift) { /* Take advantage of the fact that S = (ash 1 s_n) */
-      if (R.l < sl)
-         d = 0;
-      else if (R.l == sl) {
-            Bigit *p;
-
-            p = &R.d[sl];
-            d = *p >> slr;
-            *p &= ((Bigit)1 << slr) - 1;
-            for (i = sl; (i > 0) && (*p == 0); i--) p--;
-            R.l = i;
-         }
-      else {
-         Bigit *p;
-         
-         p = &R.d[sl+1];
-         d = *p << (64 - slr) | *(p-1) >> slr;
-         p--;
-         *p &= ((Bigit)1 << slr) - 1;
-         for (i = sl; (i > 0) && (*p == 0); i--) p--;
-         R.l = i;
-      }
-   }
-   else /* We need to do quotient-remainder */
-     d = qr();
-
-   tc1 = big_comp(&R, &MM) < ruf;
-   tc2 = add_cmp() > -ruf;
-   if (!tc1)
-      if (!tc2) {
-         mul10(&R);
-         mul10(&MM);
-         if (use_mp) mul10(&MP);
-         *buf++ = d + '0';
-         goto again;
-      }
-      else
-        OUTDIG(d+1)
-   else
-      if (!tc2)
-         OUTDIG(d)
-      else {
-         big_shift_left(&R, 1, &MM);
-         if (big_comp(&MM, &S) == -1)
-            OUTDIG(d)
-         else
-            OUTDIG(d+1)
-      }
-}
-
-void
-_PyFloat_DigitsInit()
-{
-   int n, i, l;
-   Bignum *b;
-   Bigit *xp, *zp, k;
-
-   five[0].l = l = 0;
-   five[0].d[0] = 5;
-   for (n = MAX_FIVE-1, b = &five[0]; n > 0; n--) {
-      xp = &b->d[0];
-      b++;
-      zp = &b->d[0];
-      for (i = l, k = 0; i >= 0; i--)
-         MUL(*xp++, 5, *zp++, k);
-      if (k != 0)
-         *zp = k, l++;
-      b->l = l;
-   }
-
-   /*
-   for (n = 1, b = &five[0]; n <= MAX_FIVE; n++) {
-      big_shift_left(b++, n, &R);
-      print_big(&R);
-      putchar('\n');
-   }
-   fflush(0);
-   */
-}
diff --git a/Objects/floatobject.c b/Objects/floatobject.c
index b832f7a..f79995c 100644
--- a/Objects/floatobject.c
+++ b/Objects/floatobject.c
@@ -16,10 +16,6 @@
 
 #include "formatter_string.h"
 
-#if !defined(__STDC__)
-extern double fmod(double, double);
-extern double pow(double, double);
-#endif
 
 #ifdef _OSF_SOURCE
 /* OSF1 5.1 doesn't make this available with XOPEN_SOURCE_EXTENDED defined */
@@ -251,11 +247,11 @@ PyFloat_FromString(PyObject *v, char **pend)
 			p++;
 		}
 		if (PyOS_strnicmp(p, "inf", 4) == 0) {
-			return PyFloat_FromDouble(sign * Py_HUGE_VAL);
+			Py_RETURN_INF(sign);
 		}
 #ifdef Py_NAN
 		if(PyOS_strnicmp(p, "nan", 4) == 0) {
-			return PyFloat_FromDouble(Py_NAN);
+			Py_RETURN_NAN;
 		}
 #endif
 		PyOS_snprintf(buffer, sizeof(buffer),
@@ -405,110 +401,6 @@ PyFloat_AsStringEx(char *buf, PyFloatObject *v, int precision)
 	format_float(buf, 100, v, precision);
 }
 
-#ifdef Py_BROKEN_REPR
-/* The following function is based on Tcl_PrintDouble,
- * from tclUtil.c.
- */
-
-#define is_infinite(d)	( (d) > DBL_MAX || (d) < -DBL_MAX )
-#define is_nan(d)		((d) != (d))
-
-static void
-format_double_repr(char *dst, double value)
-{
-    char *p, c;
-    int exp;
-    int signum;
-    char buffer[30];
-
-	/*
-	 * Handle NaN.
-	 */
-
-	if (is_nan(value)) {
-	    strcpy(dst, "nan");
-	    return;
-	}
-
-	/*
-	 * Handle infinities.
-	 */
-
-	if (is_infinite(value)) {
-	    if (value < 0) {
-		strcpy(dst, "-inf");
-	    } else {
-		strcpy(dst, "inf");
-	    }
-	    return;
-	}
-
-	/*
-	 * Ordinary (normal and denormal) values.
-	 */
-
-	exp = _PyFloat_Digits(buffer, value, &signum)+1;
-	if (signum) {
-	    *dst++ = '-';
-	}
-	p = buffer;
-	if (exp < -3 || exp > 17) {
-	    /*
-	     * E format for numbers < 1e-3 or >= 1e17.
-	     */
-
-	    *dst++ = *p++;
-	    c = *p;
-	    if (c != '\0') {
-		*dst++ = '.';
-		while (c != '\0') {
-		    *dst++ = c;
-		    c = *++p;
-		}
-	    }
-	    sprintf(dst, "e%+d", exp-1);
-	} else {
-	    /*
-	     * F format for others.
-	     */
-
-	    if (exp <= 0) {
-		*dst++ = '0';
-	    }
-	    c = *p;
-	    while (exp-- > 0) {
-		if (c != '\0') {
-		    *dst++ = c;
-		    c = *++p;
-		} else {
-		    *dst++ = '0';
-		}
-	    }
-	    *dst++ = '.';
-	    if (c == '\0') {
-		*dst++ = '0';
-	    } else {
-		while (++exp < 0) {
-		    *dst++ = '0';
-		}
-		while (c != '\0') {
-		    *dst++ = c;
-		    c = *++p;
-		}
-	    }
-	    *dst++ = '\0';
-	}
-}
-
-static void
-format_float_repr(char *buf, PyFloatObject *v)
-{
-	assert(PyFloat_Check(v));
-	format_double_repr(buf, PyFloat_AS_DOUBLE(v));
-}
-
-#endif /* Py_BROKEN_REPR */
-
 /* Macro and helper that convert PyObject obj to a C double and store
    the value in dbl; this replaces the functionality of the coercion
    slot function.  If conversion to double raises an exception, obj is
@@ -593,13 +485,8 @@ float_print(PyFloatObject *v, FILE *fp, int flags)
 static PyObject *
 float_repr(PyFloatObject *v)
 {
-#ifdef Py_BROKEN_REPR
-	char buf[30];
-	format_float_repr(buf, v);
-#else
 	char buf[100];
 	format_float(buf, sizeof(buf), v, PREC_REPR);
-#endif
 
 	return PyString_FromString(buf);
 }
@@ -889,10 +776,13 @@ float_div(PyObject *v, PyObject *w)
 	double a,b;
 	CONVERT_TO_DOUBLE(v, a);
 	CONVERT_TO_DOUBLE(w, b);
+#ifdef Py_NAN
 	if (b == 0.0) {
-		PyErr_SetString(PyExc_ZeroDivisionError, "float division");
+		PyErr_SetString(PyExc_ZeroDivisionError,
+				"float division");
 		return NULL;
 	}
+#endif
 	PyFPE_START_PROTECT("divide", return 0)
 	a = a / b;
 	PyFPE_END_PROTECT(a)
@@ -908,10 +798,13 @@ float_classic_div(PyObject *v, PyObject *w)
 	if (Py_DivisionWarningFlag >= 2 &&
 	    PyErr_Warn(PyExc_DeprecationWarning, "classic float division") < 0)
 		return NULL;
+#ifdef Py_NAN
 	if (b == 0.0) {
-		PyErr_SetString(PyExc_ZeroDivisionError, "float division");
+		PyErr_SetString(PyExc_ZeroDivisionError,
+				"float division");
 		return NULL;
 	}
+#endif
 	PyFPE_START_PROTECT("divide", return 0)
 	a = a / b;
 	PyFPE_END_PROTECT(a)
@@ -923,12 +816,15 @@ float_rem(PyObject *v, PyObject *w)
 {
 	double vx, wx;
 	double mod;
- 	CONVERT_TO_DOUBLE(v, vx);
- 	CONVERT_TO_DOUBLE(w, wx);
+	CONVERT_TO_DOUBLE(v, vx);
+	CONVERT_TO_DOUBLE(w, wx);
+#ifdef Py_NAN
 	if (wx == 0.0) {
-		PyErr_SetString(PyExc_ZeroDivisionError, "float modulo");
+		PyErr_SetString(PyExc_ZeroDivisionError,
+				"float modulo");
 		return NULL;
 	}
+#endif
 	PyFPE_START_PROTECT("modulo", return 0)
 	mod = fmod(vx, wx);
 	/* note: checking mod*wx < 0 is incorrect -- underflows to
@@ -1032,6 +928,9 @@ float_pow(PyObject *v, PyObject *w, PyObject *z)
 		}
 		return PyFloat_FromDouble(0.0);
 	}
+	if (iv == 1.0) { /* 1**w is 1, even 1**inf and 1**nan */
+		return PyFloat_FromDouble(1.0);
+	}
 	if (iv < 0.0) {
 		/* Whether this is an error is a mess, and bumps into libm
 		 * bugs so we have to figure it out ourselves.
@@ -1123,6 +1022,57 @@ float_coerce(PyObject **pv, PyObject **pw)
 }
 
 static PyObject *
+float_is_integer(PyObject *v)
+{
+	double x = PyFloat_AsDouble(v);
+	PyObject *o;
+	
+	if (x == -1.0 && PyErr_Occurred())
+		return NULL;
+	if (!Py_IS_FINITE(x))
+		Py_RETURN_FALSE;
+	PyFPE_START_PROTECT("is_integer", return NULL)
+	o = (floor(x) == x) ? Py_True : Py_False;
+	PyFPE_END_PROTECT(x)
+	if (errno != 0) {
+		PyErr_SetFromErrno(errno == ERANGE ? PyExc_OverflowError :
+						     PyExc_ValueError);
+		return NULL;
+	}
+	Py_INCREF(o);
+	return o;
+}
+
+#if 0
+static PyObject *
+float_is_inf(PyObject *v)
+{
+	double x = PyFloat_AsDouble(v);
+	if (x == -1.0 && PyErr_Occurred())
+		return NULL;
+	return PyBool_FromLong((long)Py_IS_INFINITY(x));
+}
+
+static PyObject *
+float_is_nan(PyObject *v)
+{
+	double x = PyFloat_AsDouble(v);
+	if (x == -1.0 && PyErr_Occurred())
+		return NULL;
+	return PyBool_FromLong((long)Py_IS_NAN(x));
+}
+
+static PyObject *
+float_is_finite(PyObject *v)
+{
+	double x = PyFloat_AsDouble(v);
+	if (x == -1.0 && PyErr_Occurred())
+		return NULL;
+	return PyBool_FromLong((long)Py_IS_FINITE(x));
+}
+#endif
+
+static PyObject *
 float_trunc(PyObject *v)
 {
 	double x = PyFloat_AsDouble(v);
@@ -1480,12 +1430,22 @@ PyDoc_STRVAR(float__format__doc,
 
 
 static PyMethodDef float_methods[] = {
-  	{"conjugate",	(PyCFunction)float_float,	METH_NOARGS,
+	{"conjugate",	(PyCFunction)float_float,	METH_NOARGS,
 	 "Returns self, the complex conjugate of any float."},
 	{"__trunc__",	(PyCFunction)float_trunc, METH_NOARGS,
          "Returns the Integral closest to x between 0 and x."},
 	{"as_integer_ratio", (PyCFunction)float_as_integer_ratio, METH_NOARGS,
 	 float_as_integer_ratio_doc},
+	{"is_integer",	(PyCFunction)float_is_integer,	METH_NOARGS,
+	 "Returns True if the float is an integer."},
+#if 0
+	{"is_inf",	(PyCFunction)float_is_inf,	METH_NOARGS,
+	 "Returns True if the float is positive or negative infinite."},
+	{"is_finite",	(PyCFunction)float_is_finite,	METH_NOARGS,
+	 "Returns True if the float is finite, neither infinite nor NaN."},
+	{"is_nan",	(PyCFunction)float_is_nan,	METH_NOARGS,
+	 "Returns True if the float is not a number (NaN)."},
+#endif
 	{"__getnewargs__",	(PyCFunction)float_getnewargs,	METH_NOARGS},
 	{"__getformat__",	(PyCFunction)float_getformat,	
 	 METH_O|METH_CLASS,		float_getformat_doc},
@@ -1646,10 +1606,6 @@ _PyFloat_Init(void)
 	double_format = detected_double_format;
 	float_format = detected_float_format;
 
-#ifdef Py_BROKEN_REPR	
-	/* Initialize floating point repr */
-	_PyFloat_DigitsInit();
-#endif
 	/* Init float info */
 	if (FloatInfoType.tp_name == 0)
 		PyStructSequence_InitType(&FloatInfoType, &floatinfo_desc);
diff --git a/Objects/intobject.c b/Objects/intobject.c
index e313954..eacad9d 100644
--- a/Objects/intobject.c
+++ b/Objects/intobject.c
@@ -1143,9 +1143,21 @@ int__format__(PyObject *self, PyObject *args)
 	return NULL;
 }
 
+#if 0
+static PyObject *
+int_is_finite(PyObject *v)
+{
+	Py_RETURN_TRUE;
+}
+#endif
+
 static PyMethodDef int_methods[] = {
 	{"conjugate",	(PyCFunction)int_int,	METH_NOARGS,
 	 "Returns self, the complex conjugate of any int."},
+#if 0
+	{"is_finite",	(PyCFunction)int_is_finite,	METH_NOARGS,
+	 "Returns always True."},
+#endif
 	{"__trunc__",	(PyCFunction)int_int,	METH_NOARGS,
          "Truncating an Integral returns itself."},
 	{"__getnewargs__",	(PyCFunction)int_getnewargs,	METH_NOARGS},
diff --git a/Objects/longobject.c b/Objects/longobject.c
index 4ab8f1e..7fc4576 100644
--- a/Objects/longobject.c
+++ b/Objects/longobject.c
@@ -3441,9 +3441,21 @@ long__format__(PyObject *self, PyObject *args)
 	return NULL;
 }
 
+#if 0
+static PyObject *
+long_is_finite(PyObject *v)
+{
+	Py_RETURN_TRUE;
+}
+#endif
+
 static PyMethodDef long_methods[] = {
 	{"conjugate",	(PyCFunction)long_long,	METH_NOARGS,
 	 "Returns self, the complex conjugate of any long."},
+#if 0
+	{"is_finite",	(PyCFunction)long_is_finite,	METH_NOARGS,
+	 "Returns always True."},
+#endif
 	{"__trunc__",	(PyCFunction)long_long,	METH_NOARGS,
          "Truncating an Integral returns itself."},
 	{"__getnewargs__",	(PyCFunction)long_getnewargs,	METH_NOARGS},