Objects/floatobject.c

   1
   2 /* Float object implementation */
   3
   4 /* XXX There should be overflow checks here, but it's hard to check
   5    for any kind of float exception without losing portability. */
   6
   7 #include "Python.h"
   8 #include "structseq.h"
   9
  10 #include <ctype.h>
  11 #include <float.h>
  12
  13 #undef MAX
  14 #undef MIN
  15 #define MAX(x, y) ((x) < (y) ? (y) : (x))
  16 #define MIN(x, y) ((x) < (y) ? (x) : (y))
  17
  18 #ifdef _OSF_SOURCE
  19 /* OSF1 5.1 doesn't make this available with XOPEN_SOURCE_EXTENDED defined */
  20 extern int finite(double);
  21 #endif
  22
  23 /* Special free list -- see comments for same code in intobject.c. */
  24 #define BLOCK_SIZE      1000    /* 1K less typical malloc overhead */
  25 #define BHEAD_SIZE      8       /* Enough for a 64-bit pointer */
  26 #define N_FLOATOBJECTS  ((BLOCK_SIZE - BHEAD_SIZE) / sizeof(PyFloatObject))
  27
  28 struct _floatblock {
  29     struct _floatblock *next;
  30     PyFloatObject objects[N_FLOATOBJECTS];
  31 };
  32
  33 typedef struct _floatblock PyFloatBlock;
  34
  35 static PyFloatBlock *block_list = NULL;
  36 static PyFloatObject *free_list = NULL;
  37
  38 static PyFloatObject *
  39 fill_free_list(void)
  40 {
  41     PyFloatObject *p, *q;
  42     /* XXX Float blocks escape the object heap. Use PyObject_MALLOC ??? */
  43     p = (PyFloatObject *) PyMem_MALLOC(sizeof(PyFloatBlock));
  44     if (p == NULL)
  45         return (PyFloatObject *) PyErr_NoMemory();
  46     ((PyFloatBlock *)p)->next = block_list;
  47     block_list = (PyFloatBlock *)p;
  48     p = &((PyFloatBlock *)p)->objects[0];
  49     q = p + N_FLOATOBJECTS;
  50     while (--q > p)
  51         Py_TYPE(q) = (struct _typeobject *)(q-1);
  52     Py_TYPE(q) = NULL;
  53     return p + N_FLOATOBJECTS - 1;
  54 }
  55
  56 double
  57 PyFloat_GetMax(void)
  58 {
  59     return DBL_MAX;
  60 }
  61
  62 double
  63 PyFloat_GetMin(void)
  64 {
  65     return DBL_MIN;
  66 }
  67
  68 static PyTypeObject FloatInfoType = {0, 0, 0, 0, 0, 0};
  69
  70 PyDoc_STRVAR(floatinfo__doc__,
  71 "sys.float_info\n\
  72 \n\
  73 A structseq holding information about the float type. It contains low level\n\
  74 information about the precision and internal representation. Please study\n\
  75 your system's :file:`float.h` for more information.");
  76
  77 static PyStructSequence_Field floatinfo_fields[] = {
  78     {"max",             "DBL_MAX -- maximum representable finite float"},
  79     {"max_exp",         "DBL_MAX_EXP -- maximum int e such that radix**(e-1) "
  80                     "is representable"},
  81     {"max_10_exp",      "DBL_MAX_10_EXP -- maximum int e such that 10**e "
  82                     "is representable"},
  83     {"min",             "DBL_MIN -- Minimum positive normalizer float"},
  84     {"min_exp",         "DBL_MIN_EXP -- minimum int e such that radix**(e-1) "
  85                     "is a normalized float"},
  86     {"min_10_exp",      "DBL_MIN_10_EXP -- minimum int e such that 10**e is "
  87                     "a normalized"},
  88     {"dig",             "DBL_DIG -- digits"},
  89     {"mant_dig",        "DBL_MANT_DIG -- mantissa digits"},
  90     {"epsilon",         "DBL_EPSILON -- Difference between 1 and the next "
  91                     "representable float"},
  92     {"radix",           "FLT_RADIX -- radix of exponent"},
  93     {"rounds",          "FLT_ROUNDS -- addition rounds"},
  94     {0}
  95 };
  96
  97 static PyStructSequence_Desc floatinfo_desc = {
  98     "sys.float_info",           /* name */
  99     floatinfo__doc__,           /* doc */
 100     floatinfo_fields,           /* fields */
 101     11
 102 };
 103
 104 PyObject *
 105 PyFloat_GetInfo(void)
 106 {
 107     PyObject* floatinfo;
 108     int pos = 0;
 109
 110     floatinfo = PyStructSequence_New(&FloatInfoType);
 111     if (floatinfo == NULL) {
 112         return NULL;
 113     }
 114
 115 #define SetIntFlag(flag) \
 116     PyStructSequence_SET_ITEM(floatinfo, pos++, PyInt_FromLong(flag))
 117 #define SetDblFlag(flag) \
 118     PyStructSequence_SET_ITEM(floatinfo, pos++, PyFloat_FromDouble(flag))
 119
 120     SetDblFlag(DBL_MAX);
 121     SetIntFlag(DBL_MAX_EXP);
 122     SetIntFlag(DBL_MAX_10_EXP);
 123     SetDblFlag(DBL_MIN);
 124     SetIntFlag(DBL_MIN_EXP);
 125     SetIntFlag(DBL_MIN_10_EXP);
 126     SetIntFlag(DBL_DIG);
 127     SetIntFlag(DBL_MANT_DIG);
 128     SetDblFlag(DBL_EPSILON);
 129     SetIntFlag(FLT_RADIX);
 130     SetIntFlag(FLT_ROUNDS);
 131 #undef SetIntFlag
 132 #undef SetDblFlag
 133
 134     if (PyErr_Occurred()) {
 135         Py_CLEAR(floatinfo);
 136         return NULL;
 137     }
 138     return floatinfo;
 139 }
 140
 141 PyObject *
 142 PyFloat_FromDouble(double fval)
 143 {
 144     register PyFloatObject *op;
 145     if (free_list == NULL) {
 146         if ((free_list = fill_free_list()) == NULL)
 147             return NULL;
 148     }
 149     /* Inline PyObject_New */
 150     op = free_list;
 151     free_list = (PyFloatObject *)Py_TYPE(op);
 152     PyObject_INIT(op, &PyFloat_Type);
 153     op->ob_fval = fval;
 154     return (PyObject *) op;
 155 }
 156
 157 /**************************************************************************
 158 RED_FLAG 22-Sep-2000 tim
 159 PyFloat_FromString's pend argument is braindead.  Prior to this RED_FLAG,
 160
 161 1.  If v was a regular string, *pend was set to point to its terminating
 162     null byte.  That's useless (the caller can find that without any
 163     help from this function!).
 164
 165 2.  If v was a Unicode string, or an object convertible to a character
 166     buffer, *pend was set to point into stack trash (the auto temp
 167     vector holding the character buffer).  That was downright dangerous.
 168
 169 Since we can't change the interface of a public API function, pend is
 170 still supported but now *officially* useless:  if pend is not NULL,
 171 *pend is set to NULL.
 172 **************************************************************************/
 173 PyObject *
 174 PyFloat_FromString(PyObject *v, char **pend)
 175 {
 176     const char *s, *last, *end;
 177     double x;
 178     char buffer[256]; /* for errors */
 179 #ifdef Py_USING_UNICODE
 180     char *s_buffer = NULL;
 181 #endif
 182     Py_ssize_t len;
 183     PyObject *result = NULL;
 184
 185     if (pend)
 186         *pend = NULL;
 187     if (PyString_Check(v)) {
 188         s = PyString_AS_STRING(v);
 189         len = PyString_GET_SIZE(v);
 190     }
 191 #ifdef Py_USING_UNICODE
 192     else if (PyUnicode_Check(v)) {
 193         s_buffer = (char *)PyMem_MALLOC(PyUnicode_GET_SIZE(v)+1);
 194         if (s_buffer == NULL)
 195             return PyErr_NoMemory();
 196         if (PyUnicode_EncodeDecimal(PyUnicode_AS_UNICODE(v),
 197                                     PyUnicode_GET_SIZE(v),
 198                                     s_buffer,
 199                                     NULL))
 200             goto error;
 201         s = s_buffer;
 202         len = strlen(s);
 203     }
 204 #endif
 205     else if (PyObject_AsCharBuffer(v, &s, &len)) {
 206         PyErr_SetString(PyExc_TypeError,
 207             "float() argument must be a string or a number");
 208         return NULL;
 209     }
 210     last = s + len;
 211
 212     while (Py_ISSPACE(*s))
 213         s++;
 214     /* We don't care about overflow or underflow.  If the platform
 215      * supports them, infinities and signed zeroes (on underflow) are
 216      * fine. */
 217     x = PyOS_string_to_double(s, (char **)&end, NULL);
 218     if (x == -1.0 && PyErr_Occurred())
 219         goto error;
 220     while (Py_ISSPACE(*end))
 221         end++;
 222     if (end == last)
 223         result = PyFloat_FromDouble(x);
 224     else {
 225         PyOS_snprintf(buffer, sizeof(buffer),
 226                       "invalid literal for float(): %.200s", s);
 227         PyErr_SetString(PyExc_ValueError, buffer);
 228         result = NULL;
 229     }
 230
 231   error:
 232 #ifdef Py_USING_UNICODE
 233     if (s_buffer)
 234         PyMem_FREE(s_buffer);
 235 #endif
 236     return result;
 237 }
 238
 239 static void
 240 float_dealloc(PyFloatObject *op)
 241 {
 242     if (PyFloat_CheckExact(op)) {
 243         Py_TYPE(op) = (struct _typeobject *)free_list;
 244         free_list = op;
 245     }
 246     else
 247         Py_TYPE(op)->tp_free((PyObject *)op);
 248 }
 249
 250 double
 251 PyFloat_AsDouble(PyObject *op)
 252 {
 253     PyNumberMethods *nb;
 254     PyFloatObject *fo;
 255     double val;
 256
 257     if (op && PyFloat_Check(op))
 258         return PyFloat_AS_DOUBLE((PyFloatObject*) op);
 259
 260     if (op == NULL) {
 261         PyErr_BadArgument();
 262         return -1;
 263     }
 264
 265     if ((nb = Py_TYPE(op)->tp_as_number) == NULL || nb->nb_float == NULL) {
 266         PyErr_SetString(PyExc_TypeError, "a float is required");
 267         return -1;
 268     }
 269
 270     fo = (PyFloatObject*) (*nb->nb_float) (op);
 271     if (fo == NULL)
 272         return -1;
 273     if (!PyFloat_Check(fo)) {
 274         PyErr_SetString(PyExc_TypeError,
 275                         "nb_float should return float object");
 276         return -1;
 277     }
 278
 279     val = PyFloat_AS_DOUBLE(fo);
 280     Py_DECREF(fo);
 281
 282     return val;
 283 }
 284
 285 /* Methods */
 286
 287 /* Macro and helper that convert PyObject obj to a C double and store
 288    the value in dbl; this replaces the functionality of the coercion
 289    slot function.  If conversion to double raises an exception, obj is
 290    set to NULL, and the function invoking this macro returns NULL.  If
 291    obj is not of float, int or long type, Py_NotImplemented is incref'ed,
 292    stored in obj, and returned from the function invoking this macro.
 293 */
 294 #define CONVERT_TO_DOUBLE(obj, dbl)                     \
 295     if (PyFloat_Check(obj))                             \
 296         dbl = PyFloat_AS_DOUBLE(obj);                   \
 297     else if (convert_to_double(&(obj), &(dbl)) < 0)     \
 298         return obj;
 299
 300 static int
 301 convert_to_double(PyObject **v, double *dbl)
 302 {
 303     register PyObject *obj = *v;
 304
 305     if (PyInt_Check(obj)) {
 306         *dbl = (double)PyInt_AS_LONG(obj);
 307     }
 308     else if (PyLong_Check(obj)) {
 309         *dbl = PyLong_AsDouble(obj);
 310         if (*dbl == -1.0 && PyErr_Occurred()) {
 311             *v = NULL;
 312             return -1;
 313         }
 314     }
 315     else {
 316         Py_INCREF(Py_NotImplemented);
 317         *v = Py_NotImplemented;
 318         return -1;
 319     }
 320     return 0;
 321 }
 322
 323 /* XXX PyFloat_AsString and PyFloat_AsReprString are deprecated:
 324    XXX they pass a char buffer without passing a length.
 325 */
 326 void
 327 PyFloat_AsString(char *buf, PyFloatObject *v)
 328 {
 329     char *tmp = PyOS_double_to_string(v->ob_fval, 'g',
 330                     PyFloat_STR_PRECISION,
 331                     Py_DTSF_ADD_DOT_0, NULL);
 332     strcpy(buf, tmp);
 333     PyMem_Free(tmp);
 334 }
 335
 336 void
 337 PyFloat_AsReprString(char *buf, PyFloatObject *v)
 338 {
 339     char * tmp = PyOS_double_to_string(v->ob_fval, 'r', 0,
 340                     Py_DTSF_ADD_DOT_0, NULL);
 341     strcpy(buf, tmp);
 342     PyMem_Free(tmp);
 343 }
 344
 345 /* ARGSUSED */
 346 static int
 347 float_print(PyFloatObject *v, FILE *fp, int flags)
 348 {
 349     char *buf;
 350     if (flags & Py_PRINT_RAW)
 351         buf = PyOS_double_to_string(v->ob_fval,
 352                                     'g', PyFloat_STR_PRECISION,
 353                                     Py_DTSF_ADD_DOT_0, NULL);
 354     else
 355         buf = PyOS_double_to_string(v->ob_fval,
 356                             'r', 0, Py_DTSF_ADD_DOT_0, NULL);
 357     Py_BEGIN_ALLOW_THREADS
 358     fputs(buf, fp);
 359     Py_END_ALLOW_THREADS
 360     PyMem_Free(buf);
 361     return 0;
 362 }
 363
 364 static PyObject *
 365 float_str_or_repr(PyFloatObject *v, int precision, char format_code)
 366 {
 367     PyObject *result;
 368     char *buf = PyOS_double_to_string(PyFloat_AS_DOUBLE(v),
 369                                   format_code, precision,
 370                                   Py_DTSF_ADD_DOT_0,
 371                                   NULL);
 372     if (!buf)
 373         return PyErr_NoMemory();
 374     result = PyString_FromString(buf);
 375     PyMem_Free(buf);
 376     return result;
 377 }
 378
 379 static PyObject *
 380 float_repr(PyFloatObject *v)
 381 {
 382     return float_str_or_repr(v, 0, 'r');
 383 }
 384
 385 static PyObject *
 386 float_str(PyFloatObject *v)
 387 {
 388     return float_str_or_repr(v, PyFloat_STR_PRECISION, 'g');
 389 }
 390
 391 /* Comparison is pretty much a nightmare.  When comparing float to float,
 392  * we do it as straightforwardly (and long-windedly) as conceivable, so
 393  * that, e.g., Python x == y delivers the same result as the platform
 394  * C x == y when x and/or y is a NaN.
 395  * When mixing float with an integer type, there's no good *uniform* approach.
 396  * Converting the double to an integer obviously doesn't work, since we
 397  * may lose info from fractional bits.  Converting the integer to a double
 398  * also has two failure modes:  (1) a long int may trigger overflow (too
 399  * large to fit in the dynamic range of a C double); (2) even a C long may have
 400  * more bits than fit in a C double (e.g., on a a 64-bit box long may have
 401  * 63 bits of precision, but a C double probably has only 53), and then
 402  * we can falsely claim equality when low-order integer bits are lost by
 403  * coercion to double.  So this part is painful too.
 404  */
 405
 406 static PyObject*
 407 float_richcompare(PyObject *v, PyObject *w, int op)
 408 {
 409     double i, j;
 410     int r = 0;
 411
 412     assert(PyFloat_Check(v));
 413     i = PyFloat_AS_DOUBLE(v);
 414
 415     /* Switch on the type of w.  Set i and j to doubles to be compared,
 416      * and op to the richcomp to use.
 417      */
 418     if (PyFloat_Check(w))
 419         j = PyFloat_AS_DOUBLE(w);
 420
 421     else if (!Py_IS_FINITE(i)) {
 422         if (PyInt_Check(w) || PyLong_Check(w))
 423             /* If i is an infinity, its magnitude exceeds any
 424              * finite integer, so it doesn't matter which int we
 425              * compare i with.  If i is a NaN, similarly.
 426              */
 427             j = 0.0;
 428         else
 429             goto Unimplemented;
 430     }
 431
 432     else if (PyInt_Check(w)) {
 433         long jj = PyInt_AS_LONG(w);
 434         /* In the worst realistic case I can imagine, C double is a
 435          * Cray single with 48 bits of precision, and long has 64
 436          * bits.
 437          */
 438 #if SIZEOF_LONG > 6
 439         unsigned long abs = (unsigned long)(jj < 0 ? -jj : jj);
 440         if (abs >> 48) {
 441             /* Needs more than 48 bits.  Make it take the
 442              * PyLong path.
 443              */
 444             PyObject *result;
 445             PyObject *ww = PyLong_FromLong(jj);
 446
 447             if (ww == NULL)
 448                 return NULL;
 449             result = float_richcompare(v, ww, op);
 450             Py_DECREF(ww);
 451             return result;
 452         }
 453 #endif
 454         j = (double)jj;
 455         assert((long)j == jj);
 456     }
 457
 458     else if (PyLong_Check(w)) {
 459         int vsign = i == 0.0 ? 0 : i < 0.0 ? -1 : 1;
 460         int wsign = _PyLong_Sign(w);
 461         size_t nbits;
 462         int exponent;
 463
 464         if (vsign != wsign) {
 465             /* Magnitudes are irrelevant -- the signs alone
 466              * determine the outcome.
 467              */
 468             i = (double)vsign;
 469             j = (double)wsign;
 470             goto Compare;
 471         }
 472         /* The signs are the same. */
 473         /* Convert w to a double if it fits.  In particular, 0 fits. */
 474         nbits = _PyLong_NumBits(w);
 475         if (nbits == (size_t)-1 && PyErr_Occurred()) {
 476             /* This long is so large that size_t isn't big enough
 477              * to hold the # of bits.  Replace with little doubles
 478              * that give the same outcome -- w is so large that
 479              * its magnitude must exceed the magnitude of any
 480              * finite float.
 481              */
 482             PyErr_Clear();
 483             i = (double)vsign;
 484             assert(wsign != 0);
 485             j = wsign * 2.0;
 486             goto Compare;
 487         }
 488         if (nbits <= 48) {
 489             j = PyLong_AsDouble(w);
 490             /* It's impossible that <= 48 bits overflowed. */
 491             assert(j != -1.0 || ! PyErr_Occurred());
 492             goto Compare;
 493         }
 494         assert(wsign != 0); /* else nbits was 0 */
 495         assert(vsign != 0); /* if vsign were 0, then since wsign is
 496                              * not 0, we would have taken the
 497                              * vsign != wsign branch at the start */
 498         /* We want to work with non-negative numbers. */
 499         if (vsign < 0) {
 500             /* "Multiply both sides" by -1; this also swaps the
 501              * comparator.
 502              */
 503             i = -i;
 504             op = _Py_SwappedOp[op];
 505         }
 506         assert(i > 0.0);
 507         (void) frexp(i, &exponent);
 508         /* exponent is the # of bits in v before the radix point;
 509          * we know that nbits (the # of bits in w) > 48 at this point
 510          */
 511         if (exponent < 0 || (size_t)exponent < nbits) {
 512             i = 1.0;
 513             j = 2.0;
 514             goto Compare;
 515         }
 516         if ((size_t)exponent > nbits) {
 517             i = 2.0;
 518             j = 1.0;
 519             goto Compare;
 520         }
 521         /* v and w have the same number of bits before the radix
 522          * point.  Construct two longs that have the same comparison
 523          * outcome.
 524          */
 525         {
 526             double fracpart;
 527             double intpart;
 528             PyObject *result = NULL;
 529             PyObject *one = NULL;
 530             PyObject *vv = NULL;
 531             PyObject *ww = w;
 532
 533             if (wsign < 0) {
 534                 ww = PyNumber_Negative(w);
 535                 if (ww == NULL)
 536                     goto Error;
 537             }
 538             else
 539                 Py_INCREF(ww);
 540
 541             fracpart = modf(i, &intpart);
 542             vv = PyLong_FromDouble(intpart);
 543             if (vv == NULL)
 544                 goto Error;
 545
 546             if (fracpart != 0.0) {
 547                 /* Shift left, and or a 1 bit into vv
 548                  * to represent the lost fraction.
 549                  */
 550                 PyObject *temp;
 551
 552                 one = PyInt_FromLong(1);
 553                 if (one == NULL)
 554                     goto Error;
 555
 556                 temp = PyNumber_Lshift(ww, one);
 557                 if (temp == NULL)
 558                     goto Error;
 559                 Py_DECREF(ww);
 560                 ww = temp;
 561
 562                 temp = PyNumber_Lshift(vv, one);
 563                 if (temp == NULL)
 564                     goto Error;
 565                 Py_DECREF(vv);
 566                 vv = temp;
 567
 568                 temp = PyNumber_Or(vv, one);
 569                 if (temp == NULL)
 570                     goto Error;
 571                 Py_DECREF(vv);
 572                 vv = temp;
 573             }
 574
 575             r = PyObject_RichCompareBool(vv, ww, op);
 576             if (r < 0)
 577                 goto Error;
 578             result = PyBool_FromLong(r);
 579          Error:
 580             Py_XDECREF(vv);
 581             Py_XDECREF(ww);
 582             Py_XDECREF(one);
 583             return result;
 584         }
 585     } /* else if (PyLong_Check(w)) */
 586
 587     else        /* w isn't float, int, or long */
 588         goto Unimplemented;
 589
 590  Compare:
 591     PyFPE_START_PROTECT("richcompare", return NULL)
 592     switch (op) {
 593     case Py_EQ:
 594         r = i == j;
 595         break;
 596     case Py_NE:
 597         r = i != j;
 598         break;
 599     case Py_LE:
 600         r = i <= j;
 601         break;
 602     case Py_GE:
 603         r = i >= j;
 604         break;
 605     case Py_LT:
 606         r = i < j;
 607         break;
 608     case Py_GT:
 609         r = i > j;
 610         break;
 611     }
 612     PyFPE_END_PROTECT(r)
 613     return PyBool_FromLong(r);
 614
 615  Unimplemented:
 616     Py_INCREF(Py_NotImplemented);
 617     return Py_NotImplemented;
 618 }
 619
 620 static long
 621 float_hash(PyFloatObject *v)
 622 {
 623     return _Py_HashDouble(v->ob_fval);
 624 }
 625
 626 static PyObject *
 627 float_add(PyObject *v, PyObject *w)
 628 {
 629     double a,b;
 630     CONVERT_TO_DOUBLE(v, a);
 631     CONVERT_TO_DOUBLE(w, b);
 632     PyFPE_START_PROTECT("add", return 0)
 633     a = a + b;
 634     PyFPE_END_PROTECT(a)
 635     return PyFloat_FromDouble(a);
 636 }
 637
 638 static PyObject *
 639 float_sub(PyObject *v, PyObject *w)
 640 {
 641     double a,b;
 642     CONVERT_TO_DOUBLE(v, a);
 643     CONVERT_TO_DOUBLE(w, b);
 644     PyFPE_START_PROTECT("subtract", return 0)
 645     a = a - b;
 646     PyFPE_END_PROTECT(a)
 647     return PyFloat_FromDouble(a);
 648 }
 649
 650 static PyObject *
 651 float_mul(PyObject *v, PyObject *w)
 652 {
 653     double a,b;
 654     CONVERT_TO_DOUBLE(v, a);
 655     CONVERT_TO_DOUBLE(w, b);
 656     PyFPE_START_PROTECT("multiply", return 0)
 657     a = a * b;
 658     PyFPE_END_PROTECT(a)
 659     return PyFloat_FromDouble(a);
 660 }
 661
 662 static PyObject *
 663 float_div(PyObject *v, PyObject *w)
 664 {
 665     double a,b;
 666     CONVERT_TO_DOUBLE(v, a);
 667     CONVERT_TO_DOUBLE(w, b);
 668 #ifdef Py_NAN
 669     if (b == 0.0) {
 670         PyErr_SetString(PyExc_ZeroDivisionError,
 671                         "float division by zero");
 672         return NULL;
 673     }
 674 #endif
 675     PyFPE_START_PROTECT("divide", return 0)
 676     a = a / b;
 677     PyFPE_END_PROTECT(a)
 678     return PyFloat_FromDouble(a);
 679 }
 680
 681 static PyObject *
 682 float_classic_div(PyObject *v, PyObject *w)
 683 {
 684     double a,b;
 685     CONVERT_TO_DOUBLE(v, a);
 686     CONVERT_TO_DOUBLE(w, b);
 687     if (Py_DivisionWarningFlag >= 2 &&
 688         PyErr_Warn(PyExc_DeprecationWarning, "classic float division") < 0)
 689         return NULL;
 690 #ifdef Py_NAN
 691     if (b == 0.0) {
 692         PyErr_SetString(PyExc_ZeroDivisionError,
 693                         "float division by zero");
 694         return NULL;
 695     }
 696 #endif
 697     PyFPE_START_PROTECT("divide", return 0)
 698     a = a / b;
 699     PyFPE_END_PROTECT(a)
 700     return PyFloat_FromDouble(a);
 701 }
 702
 703 static PyObject *
 704 float_rem(PyObject *v, PyObject *w)
 705 {
 706     double vx, wx;
 707     double mod;
 708     CONVERT_TO_DOUBLE(v, vx);
 709     CONVERT_TO_DOUBLE(w, wx);
 710 #ifdef Py_NAN
 711     if (wx == 0.0) {
 712         PyErr_SetString(PyExc_ZeroDivisionError,
 713                         "float modulo");
 714         return NULL;
 715     }
 716 #endif
 717     PyFPE_START_PROTECT("modulo", return 0)
 718     mod = fmod(vx, wx);
 719     /* note: checking mod*wx < 0 is incorrect -- underflows to
 720        0 if wx < sqrt(smallest nonzero double) */
 721     if (mod && ((wx < 0) != (mod < 0))) {
 722         mod += wx;
 723     }
 724     PyFPE_END_PROTECT(mod)
 725     return PyFloat_FromDouble(mod);
 726 }
 727
 728 static PyObject *
 729 float_divmod(PyObject *v, PyObject *w)
 730 {
 731     double vx, wx;
 732     double div, mod, floordiv;
 733     CONVERT_TO_DOUBLE(v, vx);
 734     CONVERT_TO_DOUBLE(w, wx);
 735     if (wx == 0.0) {
 736         PyErr_SetString(PyExc_ZeroDivisionError, "float divmod()");
 737         return NULL;
 738     }
 739     PyFPE_START_PROTECT("divmod", return 0)
 740     mod = fmod(vx, wx);
 741     /* fmod is typically exact, so vx-mod is *mathematically* an
 742        exact multiple of wx.  But this is fp arithmetic, and fp
 743        vx - mod is an approximation; the result is that div may
 744        not be an exact integral value after the division, although
 745        it will always be very close to one.
 746     */
 747     div = (vx - mod) / wx;
 748     if (mod) {
 749         /* ensure the remainder has the same sign as the denominator */
 750         if ((wx < 0) != (mod < 0)) {
 751             mod += wx;
 752             div -= 1.0;
 753         }
 754     }
 755     else {
 756         /* the remainder is zero, and in the presence of signed zeroes
 757            fmod returns different results across platforms; ensure
 758            it has the same sign as the denominator; we'd like to do
 759            "mod = wx * 0.0", but that may get optimized away */
 760         mod *= mod;  /* hide "mod = +0" from optimizer */
 761         if (wx < 0.0)
 762             mod = -mod;
 763     }
 764     /* snap quotient to nearest integral value */
 765     if (div) {
 766         floordiv = floor(div);
 767         if (div - floordiv > 0.5)
 768             floordiv += 1.0;
 769     }
 770     else {
 771         /* div is zero - get the same sign as the true quotient */
 772         div *= div;             /* hide "div = +0" from optimizers */
 773         floordiv = div * vx / wx; /* zero w/ sign of vx/wx */
 774     }
 775     PyFPE_END_PROTECT(floordiv)
 776     return Py_BuildValue("(dd)", floordiv, mod);
 777 }
 778
 779 static PyObject *
 780 float_floor_div(PyObject *v, PyObject *w)
 781 {
 782     PyObject *t, *r;
 783
 784     t = float_divmod(v, w);
 785     if (t == NULL || t == Py_NotImplemented)
 786         return t;
 787     assert(PyTuple_CheckExact(t));
 788     r = PyTuple_GET_ITEM(t, 0);
 789     Py_INCREF(r);
 790     Py_DECREF(t);
 791     return r;
 792 }
 793
 794 /* determine whether x is an odd integer or not;  assumes that
 795    x is not an infinity or nan. */
 796 #define DOUBLE_IS_ODD_INTEGER(x) (fmod(fabs(x), 2.0) == 1.0)
 797
 798 static PyObject *
 799 float_pow(PyObject *v, PyObject *w, PyObject *z)
 800 {
 801     double iv, iw, ix;
 802     int negate_result = 0;
 803
 804     if ((PyObject *)z != Py_None) {
 805         PyErr_SetString(PyExc_TypeError, "pow() 3rd argument not "
 806             "allowed unless all arguments are integers");
 807         return NULL;
 808     }
 809
 810     CONVERT_TO_DOUBLE(v, iv);
 811     CONVERT_TO_DOUBLE(w, iw);
 812
 813     /* Sort out special cases here instead of relying on pow() */
 814     if (iw == 0) {              /* v**0 is 1, even 0**0 */
 815         return PyFloat_FromDouble(1.0);
 816     }
 817     if (Py_IS_NAN(iv)) {        /* nan**w = nan, unless w == 0 */
 818         return PyFloat_FromDouble(iv);
 819     }
 820     if (Py_IS_NAN(iw)) {        /* v**nan = nan, unless v == 1; 1**nan = 1 */
 821         return PyFloat_FromDouble(iv == 1.0 ? 1.0 : iw);
 822     }
 823     if (Py_IS_INFINITY(iw)) {
 824         /* v**inf is: 0.0 if abs(v) < 1; 1.0 if abs(v) == 1; inf if
 825          *     abs(v) > 1 (including case where v infinite)
 826          *
 827          * v**-inf is: inf if abs(v) < 1; 1.0 if abs(v) == 1; 0.0 if
 828          *     abs(v) > 1 (including case where v infinite)
 829          */
 830         iv = fabs(iv);
 831         if (iv == 1.0)
 832             return PyFloat_FromDouble(1.0);
 833         else if ((iw > 0.0) == (iv > 1.0))
 834             return PyFloat_FromDouble(fabs(iw)); /* return inf */
 835         else
 836             return PyFloat_FromDouble(0.0);
 837     }
 838     if (Py_IS_INFINITY(iv)) {
 839         /* (+-inf)**w is: inf for w positive, 0 for w negative; in
 840          *     both cases, we need to add the appropriate sign if w is
 841          *     an odd integer.
 842          */
 843         int iw_is_odd = DOUBLE_IS_ODD_INTEGER(iw);
 844         if (iw > 0.0)
 845             return PyFloat_FromDouble(iw_is_odd ? iv : fabs(iv));
 846         else
 847             return PyFloat_FromDouble(iw_is_odd ?
 848                                       copysign(0.0, iv) : 0.0);
 849     }
 850     if (iv == 0.0) {  /* 0**w is: 0 for w positive, 1 for w zero
 851                          (already dealt with above), and an error
 852                          if w is negative. */
 853         int iw_is_odd = DOUBLE_IS_ODD_INTEGER(iw);
 854         if (iw < 0.0) {
 855             PyErr_SetString(PyExc_ZeroDivisionError,
 856                             "0.0 cannot be raised to a "
 857                             "negative power");
 858             return NULL;
 859         }
 860         /* use correct sign if iw is odd */
 861         return PyFloat_FromDouble(iw_is_odd ? iv : 0.0);
 862     }
 863
 864     if (iv < 0.0) {
 865         /* Whether this is an error is a mess, and bumps into libm
 866          * bugs so we have to figure it out ourselves.
 867          */
 868         if (iw != floor(iw)) {
 869             PyErr_SetString(PyExc_ValueError, "negative number "
 870                 "cannot be raised to a fractional power");
 871             return NULL;
 872         }
 873         /* iw is an exact integer, albeit perhaps a very large
 874          * one.  Replace iv by its absolute value and remember
 875          * to negate the pow result if iw is odd.
 876          */
 877         iv = -iv;
 878         negate_result = DOUBLE_IS_ODD_INTEGER(iw);
 879     }
 880
 881     if (iv == 1.0) { /* 1**w is 1, even 1**inf and 1**nan */
 882         /* (-1) ** large_integer also ends up here.  Here's an
 883          * extract from the comments for the previous
 884          * implementation explaining why this special case is
 885          * necessary:
 886          *
 887          * -1 raised to an exact integer should never be exceptional.
 888          * Alas, some libms (chiefly glibc as of early 2003) return
 889          * NaN and set EDOM on pow(-1, large_int) if the int doesn't
 890          * happen to be representable in a *C* integer.  That's a
 891          * bug.
 892          */
 893         return PyFloat_FromDouble(negate_result ? -1.0 : 1.0);
 894     }
 895
 896     /* Now iv and iw are finite, iw is nonzero, and iv is
 897      * positive and not equal to 1.0.  We finally allow
 898      * the platform pow to step in and do the rest.
 899      */
 900     errno = 0;
 901     PyFPE_START_PROTECT("pow", return NULL)
 902     ix = pow(iv, iw);
 903     PyFPE_END_PROTECT(ix)
 904     Py_ADJUST_ERANGE1(ix);
 905     if (negate_result)
 906         ix = -ix;
 907
 908     if (errno != 0) {
 909         /* We don't expect any errno value other than ERANGE, but
 910          * the range of libm bugs appears unbounded.
 911          */
 912         PyErr_SetFromErrno(errno == ERANGE ? PyExc_OverflowError :
 913                              PyExc_ValueError);
 914         return NULL;
 915     }
 916     return PyFloat_FromDouble(ix);
 917 }
 918
 919 #undef DOUBLE_IS_ODD_INTEGER
 920
 921 static PyObject *
 922 float_neg(PyFloatObject *v)
 923 {
 924     return PyFloat_FromDouble(-v->ob_fval);
 925 }
 926
 927 static PyObject *
 928 float_abs(PyFloatObject *v)
 929 {
 930     return PyFloat_FromDouble(fabs(v->ob_fval));
 931 }
 932
 933 static int
 934 float_nonzero(PyFloatObject *v)
 935 {
 936     return v->ob_fval != 0.0;
 937 }
 938
 939 static int
 940 float_coerce(PyObject **pv, PyObject **pw)
 941 {
 942     if (PyInt_Check(*pw)) {
 943         long x = PyInt_AsLong(*pw);
 944         *pw = PyFloat_FromDouble((double)x);
 945         Py_INCREF(*pv);
 946         return 0;
 947     }
 948     else if (PyLong_Check(*pw)) {
 949         double x = PyLong_AsDouble(*pw);
 950         if (x == -1.0 && PyErr_Occurred())
 951             return -1;
 952         *pw = PyFloat_FromDouble(x);
 953         Py_INCREF(*pv);
 954         return 0;
 955     }
 956     else if (PyFloat_Check(*pw)) {
 957         Py_INCREF(*pv);
 958         Py_INCREF(*pw);
 959         return 0;
 960     }
 961     return 1; /* Can't do it */
 962 }
 963
 964 static PyObject *
 965 float_is_integer(PyObject *v)
 966 {
 967     double x = PyFloat_AsDouble(v);
 968     PyObject *o;
 969
 970     if (x == -1.0 && PyErr_Occurred())
 971         return NULL;
 972     if (!Py_IS_FINITE(x))
 973         Py_RETURN_FALSE;
 974     errno = 0;
 975     PyFPE_START_PROTECT("is_integer", return NULL)
 976     o = (floor(x) == x) ? Py_True : Py_False;
 977     PyFPE_END_PROTECT(x)
 978     if (errno != 0) {
 979         PyErr_SetFromErrno(errno == ERANGE ? PyExc_OverflowError :
 980                              PyExc_ValueError);
 981         return NULL;
 982     }
 983     Py_INCREF(o);
 984     return o;
 985 }
 986
 987 #if 0
 988 static PyObject *
 989 float_is_inf(PyObject *v)
 990 {
 991     double x = PyFloat_AsDouble(v);
 992     if (x == -1.0 && PyErr_Occurred())
 993         return NULL;
 994     return PyBool_FromLong((long)Py_IS_INFINITY(x));
 995 }
 996
 997 static PyObject *
 998 float_is_nan(PyObject *v)
 999 {
1000     double x = PyFloat_AsDouble(v);
1001     if (x == -1.0 && PyErr_Occurred())
1002         return NULL;
1003     return PyBool_FromLong((long)Py_IS_NAN(x));
1004 }
1005
1006 static PyObject *
1007 float_is_finite(PyObject *v)
1008 {
1009     double x = PyFloat_AsDouble(v);
1010     if (x == -1.0 && PyErr_Occurred())
1011         return NULL;
1012     return PyBool_FromLong((long)Py_IS_FINITE(x));
1013 }
1014 #endif
1015
1016 static PyObject *
1017 float_trunc(PyObject *v)
1018 {
1019     double x = PyFloat_AsDouble(v);
1020     double wholepart;           /* integral portion of x, rounded toward 0 */
1021
1022     (void)modf(x, &wholepart);
1023     /* Try to get out cheap if this fits in a Python int.  The attempt
1024      * to cast to long must be protected, as C doesn't define what
1025      * happens if the double is too big to fit in a long.  Some rare
1026      * systems raise an exception then (RISCOS was mentioned as one,
1027      * and someone using a non-default option on Sun also bumped into
1028      * that).  Note that checking for >= and <= LONG_{MIN,MAX} would
1029      * still be vulnerable:  if a long has more bits of precision than
1030      * a double, casting MIN/MAX to double may yield an approximation,
1031      * and if that's rounded up, then, e.g., wholepart=LONG_MAX+1 would
1032      * yield true from the C expression wholepart<=LONG_MAX, despite
1033      * that wholepart is actually greater than LONG_MAX.
1034      */
1035     if (LONG_MIN < wholepart && wholepart < LONG_MAX) {
1036         const long aslong = (long)wholepart;
1037         return PyInt_FromLong(aslong);
1038     }
1039     return PyLong_FromDouble(wholepart);
1040 }
1041
1042 static PyObject *
1043 float_long(PyObject *v)
1044 {
1045     double x = PyFloat_AsDouble(v);
1046     return PyLong_FromDouble(x);
1047 }
1048
1049 /* _Py_double_round: rounds a finite nonzero double to the closest multiple of
1050    10**-ndigits; here ndigits is within reasonable bounds (typically, -308 <=
1051    ndigits <= 323).  Returns a Python float, or sets a Python error and
1052    returns NULL on failure (OverflowError and memory errors are possible). */
1053
1054 #ifndef PY_NO_SHORT_FLOAT_REPR
1055 /* version of _Py_double_round that uses the correctly-rounded string<->double
1056    conversions from Python/dtoa.c */
1057
1058 /* FIVE_POW_LIMIT is the largest k such that 5**k is exactly representable as
1059    a double.  Since we're using the code in Python/dtoa.c, it should be safe
1060    to assume that C doubles are IEEE 754 binary64 format.  To be on the safe
1061    side, we check this. */
1062 #if DBL_MANT_DIG == 53
1063 #define FIVE_POW_LIMIT 22
1064 #else
1065 #error "C doubles do not appear to be IEEE 754 binary64 format"
1066 #endif
1067
1068 PyObject *
1069 _Py_double_round(double x, int ndigits) {
1070
1071     double rounded, m;
1072     Py_ssize_t buflen, mybuflen=100;
1073     char *buf, *buf_end, shortbuf[100], *mybuf=shortbuf;
1074     int decpt, sign, val, halfway_case;
1075     PyObject *result = NULL;
1076
1077     /* The basic idea is very simple: convert and round the double to a
1078        decimal string using _Py_dg_dtoa, then convert that decimal string
1079        back to a double with _Py_dg_strtod.  There's one minor difficulty:
1080        Python 2.x expects round to do round-half-away-from-zero, while
1081        _Py_dg_dtoa does round-half-to-even.  So we need some way to detect
1082        and correct the halfway cases.
1083
1084        Detection: a halfway value has the form k * 0.5 * 10**-ndigits for
1085        some odd integer k.  Or in other words, a rational number x is
1086        exactly halfway between two multiples of 10**-ndigits if its
1087        2-valuation is exactly -ndigits-1 and its 5-valuation is at least
1088        -ndigits.  For ndigits >= 0 the latter condition is automatically
1089        satisfied for a binary float x, since any such float has
1090        nonnegative 5-valuation.  For 0 > ndigits >= -22, x needs to be an
1091        integral multiple of 5**-ndigits; we can check this using fmod.
1092        For -22 > ndigits, there are no halfway cases: 5**23 takes 54 bits
1093        to represent exactly, so any odd multiple of 0.5 * 10**n for n >=
1094        23 takes at least 54 bits of precision to represent exactly.
1095
1096        Correction: a simple strategy for dealing with halfway cases is to
1097        (for the halfway cases only) call _Py_dg_dtoa with an argument of
1098        ndigits+1 instead of ndigits (thus doing an exact conversion to
1099        decimal), round the resulting string manually, and then convert
1100        back using _Py_dg_strtod.
1101     */
1102
1103     /* nans, infinities and zeros should have already been dealt
1104        with by the caller (in this case, builtin_round) */
1105     assert(Py_IS_FINITE(x) && x != 0.0);
1106
1107     /* find 2-valuation val of x */
1108     m = frexp(x, &val);
1109     while (m != floor(m)) {
1110         m *= 2.0;
1111         val--;
1112     }
1113
1114     /* determine whether this is a halfway case */
1115     if (val == -ndigits-1) {
1116         if (ndigits >= 0)
1117             halfway_case = 1;
1118         else if (ndigits >= -FIVE_POW_LIMIT) {
1119             double five_pow = 1.0;
1120             int i;
1121             for (i=0; i < -ndigits; i++)
1122                 five_pow *= 5.0;
1123             halfway_case = fmod(x, five_pow) == 0.0;
1124         }
1125         else
1126             halfway_case = 0;
1127     }
1128     else
1129         halfway_case = 0;
1130
1131     /* round to a decimal string; use an extra place for halfway case */
1132     buf = _Py_dg_dtoa(x, 3, ndigits+halfway_case, &decpt, &sign, &buf_end);
1133     if (buf == NULL) {
1134         PyErr_NoMemory();
1135         return NULL;
1136     }
1137     buflen = buf_end - buf;
1138
1139     /* in halfway case, do the round-half-away-from-zero manually */
1140     if (halfway_case) {
1141         int i, carry;
1142         /* sanity check: _Py_dg_dtoa should not have stripped
1143            any zeros from the result: there should be exactly
1144            ndigits+1 places following the decimal point, and
1145            the last digit in the buffer should be a '5'.*/
1146         assert(buflen - decpt == ndigits+1);
1147         assert(buf[buflen-1] == '5');
1148
1149         /* increment and shift right at the same time. */
1150         decpt += 1;
1151         carry = 1;
1152         for (i=buflen-1; i-- > 0;) {
1153             carry += buf[i] - '0';
1154             buf[i+1] = carry % 10 + '0';
1155             carry /= 10;
1156         }
1157         buf[0] = carry + '0';
1158     }
1159
1160     /* Get new buffer if shortbuf is too small.  Space needed <= buf_end -
1161        buf + 8: (1 extra for '0', 1 for sign, 5 for exp, 1 for '\0'). */
1162     if (buflen + 8 > mybuflen) {
1163         mybuflen = buflen+8;
1164         mybuf = (char *)PyMem_Malloc(mybuflen);
1165         if (mybuf == NULL) {
1166             PyErr_NoMemory();
1167             goto exit;
1168         }
1169     }
1170     /* copy buf to mybuf, adding exponent, sign and leading 0 */
1171     PyOS_snprintf(mybuf, mybuflen, "%s0%se%d", (sign ? "-" : ""),
1172                   buf, decpt - (int)buflen);
1173
1174     /* and convert the resulting string back to a double */
1175     errno = 0;
1176     rounded = _Py_dg_strtod(mybuf, NULL);
1177     if (errno == ERANGE && fabs(rounded) >= 1.)
1178         PyErr_SetString(PyExc_OverflowError,
1179                         "rounded value too large to represent");
1180     else
1181         result = PyFloat_FromDouble(rounded);
1182
1183     /* done computing value;  now clean up */
1184     if (mybuf != shortbuf)
1185         PyMem_Free(mybuf);
1186   exit:
1187     _Py_dg_freedtoa(buf);
1188     return result;
1189 }
1190
1191 #undef FIVE_POW_LIMIT
1192
1193 #else /* PY_NO_SHORT_FLOAT_REPR */
1194
1195 /* fallback version, to be used when correctly rounded binary<->decimal
1196    conversions aren't available */
1197
1198 PyObject *
1199 _Py_double_round(double x, int ndigits) {
1200     double pow1, pow2, y, z;
1201     if (ndigits >= 0) {
1202         if (ndigits > 22) {
1203             /* pow1 and pow2 are each safe from overflow, but
1204                pow1*pow2 ~= pow(10.0, ndigits) might overflow */
1205             pow1 = pow(10.0, (double)(ndigits-22));
1206             pow2 = 1e22;
1207         }
1208         else {
1209             pow1 = pow(10.0, (double)ndigits);
1210             pow2 = 1.0;
1211         }
1212         y = (x*pow1)*pow2;
1213         /* if y overflows, then rounded value is exactly x */
1214         if (!Py_IS_FINITE(y))
1215             return PyFloat_FromDouble(x);
1216     }
1217     else {
1218         pow1 = pow(10.0, (double)-ndigits);
1219         pow2 = 1.0; /* unused; silences a gcc compiler warning */
1220         y = x / pow1;
1221     }
1222
1223     z = round(y);
1224     if (fabs(y-z) == 0.5)
1225         /* halfway between two integers; use round-away-from-zero */
1226         z = y + copysign(0.5, y);
1227
1228     if (ndigits >= 0)
1229         z = (z / pow2) / pow1;
1230     else
1231         z *= pow1;
1232
1233     /* if computation resulted in overflow, raise OverflowError */
1234     if (!Py_IS_FINITE(z)) {
1235         PyErr_SetString(PyExc_OverflowError,
1236                         "overflow occurred during round");
1237         return NULL;
1238     }
1239
1240     return PyFloat_FromDouble(z);
1241 }
1242
1243 #endif /* PY_NO_SHORT_FLOAT_REPR */
1244
1245 static PyObject *
1246 float_float(PyObject *v)
1247 {
1248     if (PyFloat_CheckExact(v))
1249         Py_INCREF(v);
1250     else
1251         v = PyFloat_FromDouble(((PyFloatObject *)v)->ob_fval);
1252     return v;
1253 }
1254
1255 /* turn ASCII hex characters into integer values and vice versa */
1256
1257 static char
1258 char_from_hex(int x)
1259 {
1260     assert(0 <= x && x < 16);
1261     return "0123456789abcdef"[x];
1262 }
1263
1264 static int
1265 hex_from_char(char c) {
1266     int x;
1267     switch(c) {
1268     case '0':
1269         x = 0;
1270         break;
1271     case '1':
1272         x = 1;
1273         break;
1274     case '2':
1275         x = 2;
1276         break;
1277     case '3':
1278         x = 3;
1279         break;
1280     case '4':
1281         x = 4;
1282         break;
1283     case '5':
1284         x = 5;
1285         break;
1286     case '6':
1287         x = 6;
1288         break;
1289     case '7':
1290         x = 7;
1291         break;
1292     case '8':
1293         x = 8;
1294         break;
1295     case '9':
1296         x = 9;
1297         break;
1298     case 'a':
1299     case 'A':
1300         x = 10;
1301         break;
1302     case 'b':
1303     case 'B':
1304         x = 11;
1305         break;
1306     case 'c':
1307     case 'C':
1308         x = 12;
1309         break;
1310     case 'd':
1311     case 'D':
1312         x = 13;
1313         break;
1314     case 'e':
1315     case 'E':
1316         x = 14;
1317         break;
1318     case 'f':
1319     case 'F':
1320         x = 15;
1321         break;
1322     default:
1323         x = -1;
1324         break;
1325     }
1326     return x;
1327 }
1328
1329 /* convert a float to a hexadecimal string */
1330
1331 /* TOHEX_NBITS is DBL_MANT_DIG rounded up to the next integer
1332    of the form 4k+1. */
1333 #define TOHEX_NBITS DBL_MANT_DIG + 3 - (DBL_MANT_DIG+2)%4
1334
1335 static PyObject *
1336 float_hex(PyObject *v)
1337 {
1338     double x, m;
1339     int e, shift, i, si, esign;
1340     /* Space for 1+(TOHEX_NBITS-1)/4 digits, a decimal point, and the
1341        trailing NUL byte. */
1342     char s[(TOHEX_NBITS-1)/4+3];
1343
1344     CONVERT_TO_DOUBLE(v, x);
1345
1346     if (Py_IS_NAN(x) || Py_IS_INFINITY(x))
1347         return float_str((PyFloatObject *)v);
1348
1349     if (x == 0.0) {
1350         if (copysign(1.0, x) == -1.0)
1351             return PyString_FromString("-0x0.0p+0");
1352         else
1353             return PyString_FromString("0x0.0p+0");
1354     }
1355
1356     m = frexp(fabs(x), &e);
1357     shift = 1 - MAX(DBL_MIN_EXP - e, 0);
1358     m = ldexp(m, shift);
1359     e -= shift;
1360
1361     si = 0;
1362     s[si] = char_from_hex((int)m);
1363     si++;
1364     m -= (int)m;
1365     s[si] = '.';
1366     si++;
1367     for (i=0; i < (TOHEX_NBITS-1)/4; i++) {
1368         m *= 16.0;
1369         s[si] = char_from_hex((int)m);
1370         si++;
1371         m -= (int)m;
1372     }
1373     s[si] = '\0';
1374
1375     if (e < 0) {
1376         esign = (int)'-';
1377         e = -e;
1378     }
1379     else
1380         esign = (int)'+';
1381
1382     if (x < 0.0)
1383         return PyString_FromFormat("-0x%sp%c%d", s, esign, e);
1384     else
1385         return PyString_FromFormat("0x%sp%c%d", s, esign, e);
1386 }
1387
1388 PyDoc_STRVAR(float_hex_doc,
1389 "float.hex() -> string\n\
1390 \n\
1391 Return a hexadecimal representation of a floating-point number.\n\
1392 >>> (-0.1).hex()\n\
1393 '-0x1.999999999999ap-4'\n\
1394 >>> 3.14159.hex()\n\
1395 '0x1.921f9f01b866ep+1'");
1396
1397 /* Case-insensitive locale-independent string match used for nan and inf
1398    detection. t should be lower-case and null-terminated.  Return a nonzero
1399    result if the first strlen(t) characters of s match t and 0 otherwise. */
1400
1401 static int
1402 case_insensitive_match(const char *s, const char *t)
1403 {
1404     while(*t && Py_TOLOWER(*s) == *t) {
1405         s++;
1406         t++;
1407     }
1408     return *t ? 0 : 1;
1409 }
1410
1411 /* Convert a hexadecimal string to a float. */
1412
1413 static PyObject *
1414 float_fromhex(PyObject *cls, PyObject *arg)
1415 {
1416     PyObject *result_as_float, *result;
1417     double x;
1418     long exp, top_exp, lsb, key_digit;
1419     char *s, *coeff_start, *s_store, *coeff_end, *exp_start, *s_end;
1420     int half_eps, digit, round_up, sign=1;
1421     Py_ssize_t length, ndigits, fdigits, i;
1422
1423     /*
1424      * For the sake of simplicity and correctness, we impose an artificial
1425      * limit on ndigits, the total number of hex digits in the coefficient
1426      * The limit is chosen to ensure that, writing exp for the exponent,
1427      *
1428      *   (1) if exp > LONG_MAX/2 then the value of the hex string is
1429      *   guaranteed to overflow (provided it's nonzero)
1430      *
1431      *   (2) if exp < LONG_MIN/2 then the value of the hex string is
1432      *   guaranteed to underflow to 0.
1433      *
1434      *   (3) if LONG_MIN/2 <= exp <= LONG_MAX/2 then there's no danger of
1435      *   overflow in the calculation of exp and top_exp below.
1436      *
1437      * More specifically, ndigits is assumed to satisfy the following
1438      * inequalities:
1439      *
1440      *   4*ndigits <= DBL_MIN_EXP - DBL_MANT_DIG - LONG_MIN/2
1441      *   4*ndigits <= LONG_MAX/2 + 1 - DBL_MAX_EXP
1442      *
1443      * If either of these inequalities is not satisfied, a ValueError is
1444      * raised.  Otherwise, write x for the value of the hex string, and
1445      * assume x is nonzero.  Then
1446      *
1447      *   2**(exp-4*ndigits) <= |x| < 2**(exp+4*ndigits).
1448      *
1449      * Now if exp > LONG_MAX/2 then:
1450      *
1451      *   exp - 4*ndigits >= LONG_MAX/2 + 1 - (LONG_MAX/2 + 1 - DBL_MAX_EXP)
1452      *                    = DBL_MAX_EXP
1453      *
1454      * so |x| >= 2**DBL_MAX_EXP, which is too large to be stored in C
1455      * double, so overflows.  If exp < LONG_MIN/2, then
1456      *
1457      *   exp + 4*ndigits <= LONG_MIN/2 - 1 + (
1458      *                      DBL_MIN_EXP - DBL_MANT_DIG - LONG_MIN/2)
1459      *                    = DBL_MIN_EXP - DBL_MANT_DIG - 1
1460      *
1461      * and so |x| < 2**(DBL_MIN_EXP-DBL_MANT_DIG-1), hence underflows to 0
1462      * when converted to a C double.
1463      *
1464      * It's easy to show that if LONG_MIN/2 <= exp <= LONG_MAX/2 then both
1465      * exp+4*ndigits and exp-4*ndigits are within the range of a long.
1466      */
1467
1468     if (PyString_AsStringAndSize(arg, &s, &length))
1469         return NULL;
1470     s_end = s + length;
1471
1472     /********************
1473      * Parse the string *
1474      ********************/
1475
1476     /* leading whitespace and optional sign */
1477     while (Py_ISSPACE(*s))
1478         s++;
1479     if (*s == '-') {
1480         s++;
1481         sign = -1;
1482     }
1483     else if (*s == '+')
1484         s++;
1485
1486     /* infinities and nans */
1487     if (*s == 'i' || *s == 'I') {
1488         if (!case_insensitive_match(s+1, "nf"))
1489             goto parse_error;
1490         s += 3;
1491         x = Py_HUGE_VAL;
1492         if (case_insensitive_match(s, "inity"))
1493             s += 5;
1494         goto finished;
1495     }
1496     if (*s == 'n' || *s == 'N') {
1497         if (!case_insensitive_match(s+1, "an"))
1498             goto parse_error;
1499         s += 3;
1500         x = Py_NAN;
1501         goto finished;
1502     }
1503
1504     /* [0x] */
1505     s_store = s;
1506     if (*s == '0') {
1507         s++;
1508         if (*s == 'x' || *s == 'X')
1509             s++;
1510         else
1511             s = s_store;
1512     }
1513
1514     /* coefficient: <integer> [. <fraction>] */
1515     coeff_start = s;
1516     while (hex_from_char(*s) >= 0)
1517         s++;
1518     s_store = s;
1519     if (*s == '.') {
1520         s++;
1521         while (hex_from_char(*s) >= 0)
1522             s++;
1523         coeff_end = s-1;
1524     }
1525     else
1526         coeff_end = s;
1527
1528     /* ndigits = total # of hex digits; fdigits = # after point */
1529     ndigits = coeff_end - coeff_start;
1530     fdigits = coeff_end - s_store;
1531     if (ndigits == 0)
1532         goto parse_error;
1533     if (ndigits > MIN(DBL_MIN_EXP - DBL_MANT_DIG - LONG_MIN/2,
1534                       LONG_MAX/2 + 1 - DBL_MAX_EXP)/4)
1535         goto insane_length_error;
1536
1537     /* [p <exponent>] */
1538     if (*s == 'p' || *s == 'P') {
1539         s++;
1540         exp_start = s;
1541         if (*s == '-' || *s == '+')
1542             s++;
1543         if (!('0' <= *s && *s <= '9'))
1544             goto parse_error;
1545         s++;
1546         while ('0' <= *s && *s <= '9')
1547             s++;
1548         exp = strtol(exp_start, NULL, 10);
1549     }
1550     else
1551         exp = 0;
1552
1553 /* for 0 <= j < ndigits, HEX_DIGIT(j) gives the jth most significant digit */
1554 #define HEX_DIGIT(j) hex_from_char(*((j) < fdigits ?            \
1555                      coeff_end-(j) :                                    \
1556                      coeff_end-1-(j)))
1557
1558     /*******************************************
1559      * Compute rounded value of the hex string *
1560      *******************************************/
1561
1562     /* Discard leading zeros, and catch extreme overflow and underflow */
1563     while (ndigits > 0 && HEX_DIGIT(ndigits-1) == 0)
1564         ndigits--;
1565     if (ndigits == 0 || exp < LONG_MIN/2) {
1566         x = 0.0;
1567         goto finished;
1568     }
1569     if (exp > LONG_MAX/2)
1570         goto overflow_error;
1571
1572     /* Adjust exponent for fractional part. */
1573     exp = exp - 4*((long)fdigits);
1574
1575     /* top_exp = 1 more than exponent of most sig. bit of coefficient */
1576     top_exp = exp + 4*((long)ndigits - 1);
1577     for (digit = HEX_DIGIT(ndigits-1); digit != 0; digit /= 2)
1578         top_exp++;
1579
1580     /* catch almost all nonextreme cases of overflow and underflow here */
1581     if (top_exp < DBL_MIN_EXP - DBL_MANT_DIG) {
1582         x = 0.0;
1583         goto finished;
1584     }
1585     if (top_exp > DBL_MAX_EXP)
1586         goto overflow_error;
1587
1588     /* lsb = exponent of least significant bit of the *rounded* value.
1589        This is top_exp - DBL_MANT_DIG unless result is subnormal. */
1590     lsb = MAX(top_exp, (long)DBL_MIN_EXP) - DBL_MANT_DIG;
1591
1592     x = 0.0;
1593     if (exp >= lsb) {
1594         /* no rounding required */
1595         for (i = ndigits-1; i >= 0; i--)
1596             x = 16.0*x + HEX_DIGIT(i);
1597         x = ldexp(x, (int)(exp));
1598         goto finished;
1599     }
1600     /* rounding required.  key_digit is the index of the hex digit
1601        containing the first bit to be rounded away. */
1602     half_eps = 1 << (int)((lsb - exp - 1) % 4);
1603     key_digit = (lsb - exp - 1) / 4;
1604     for (i = ndigits-1; i > key_digit; i--)
1605         x = 16.0*x + HEX_DIGIT(i);
1606     digit = HEX_DIGIT(key_digit);
1607     x = 16.0*x + (double)(digit & (16-2*half_eps));
1608
1609     /* round-half-even: round up if bit lsb-1 is 1 and at least one of
1610        bits lsb, lsb-2, lsb-3, lsb-4, ... is 1. */
1611     if ((digit & half_eps) != 0) {
1612         round_up = 0;
1613         if ((digit & (3*half_eps-1)) != 0 ||
1614             (half_eps == 8 && (HEX_DIGIT(key_digit+1) & 1) != 0))
1615             round_up = 1;
1616         else
1617             for (i = key_digit-1; i >= 0; i--)
1618                 if (HEX_DIGIT(i) != 0) {
1619                     round_up = 1;
1620                     break;
1621                 }
1622         if (round_up == 1) {
1623             x += 2*half_eps;
1624             if (top_exp == DBL_MAX_EXP &&
1625                 x == ldexp((double)(2*half_eps), DBL_MANT_DIG))
1626                 /* overflow corner case: pre-rounded value <
1627                    2**DBL_MAX_EXP; rounded=2**DBL_MAX_EXP. */
1628                 goto overflow_error;
1629         }
1630     }
1631     x = ldexp(x, (int)(exp+4*key_digit));
1632
1633   finished:
1634     /* optional trailing whitespace leading to the end of the string */
1635     while (Py_ISSPACE(*s))
1636         s++;
1637     if (s != s_end)
1638         goto parse_error;
1639     result_as_float = Py_BuildValue("(d)", sign * x);
1640     if (result_as_float == NULL)
1641         return NULL;
1642     result = PyObject_CallObject(cls, result_as_float);
1643     Py_DECREF(result_as_float);
1644     return result;
1645
1646   overflow_error:
1647     PyErr_SetString(PyExc_OverflowError,
1648                     "hexadecimal value too large to represent as a float");
1649     return NULL;
1650
1651   parse_error:
1652     PyErr_SetString(PyExc_ValueError,
1653                     "invalid hexadecimal floating-point string");
1654     return NULL;
1655
1656   insane_length_error:
1657     PyErr_SetString(PyExc_ValueError,
1658                     "hexadecimal string too long to convert");
1659     return NULL;
1660 }
1661
1662 PyDoc_STRVAR(float_fromhex_doc,
1663 "float.fromhex(string) -> float\n\
1664 \n\
1665 Create a floating-point number from a hexadecimal string.\n\
1666 >>> float.fromhex('0x1.ffffp10')\n\
1667 2047.984375\n\
1668 >>> float.fromhex('-0x1p-1074')\n\
1669 -4.9406564584124654e-324");
1670
1671
1672 static PyObject *
1673 float_as_integer_ratio(PyObject *v, PyObject *unused)
1674 {
1675     double self;
1676     double float_part;
1677     int exponent;
1678     int i;
1679
1680     PyObject *prev;
1681     PyObject *py_exponent = NULL;
1682     PyObject *numerator = NULL;
1683     PyObject *denominator = NULL;
1684     PyObject *result_pair = NULL;
1685     PyNumberMethods *long_methods = PyLong_Type.tp_as_number;
1686
1687 #define INPLACE_UPDATE(obj, call) \
1688     prev = obj; \
1689     obj = call; \
1690     Py_DECREF(prev); \
1691
1692     CONVERT_TO_DOUBLE(v, self);
1693
1694     if (Py_IS_INFINITY(self)) {
1695       PyErr_SetString(PyExc_OverflowError,
1696                       "Cannot pass infinity to float.as_integer_ratio.");
1697       return NULL;
1698     }
1699 #ifdef Py_NAN
1700     if (Py_IS_NAN(self)) {
1701       PyErr_SetString(PyExc_ValueError,
1702                       "Cannot pass NaN to float.as_integer_ratio.");
1703       return NULL;
1704     }
1705 #endif
1706
1707     PyFPE_START_PROTECT("as_integer_ratio", goto error);
1708     float_part = frexp(self, &exponent);        /* self == float_part * 2**exponent exactly */
1709     PyFPE_END_PROTECT(float_part);
1710
1711     for (i=0; i<300 && float_part != floor(float_part) ; i++) {
1712         float_part *= 2.0;
1713         exponent--;
1714     }
1715     /* self == float_part * 2**exponent exactly and float_part is integral.
1716        If FLT_RADIX != 2, the 300 steps may leave a tiny fractional part
1717        to be truncated by PyLong_FromDouble(). */
1718
1719     numerator = PyLong_FromDouble(float_part);
1720     if (numerator == NULL) goto error;
1721
1722     /* fold in 2**exponent */
1723     denominator = PyLong_FromLong(1);
1724     py_exponent = PyLong_FromLong(labs((long)exponent));
1725     if (py_exponent == NULL) goto error;
1726     INPLACE_UPDATE(py_exponent,
1727                    long_methods->nb_lshift(denominator, py_exponent));
1728     if (py_exponent == NULL) goto error;
1729     if (exponent > 0) {
1730         INPLACE_UPDATE(numerator,
1731                        long_methods->nb_multiply(numerator, py_exponent));
1732         if (numerator == NULL) goto error;
1733     }
1734     else {
1735         Py_DECREF(denominator);
1736         denominator = py_exponent;
1737         py_exponent = NULL;
1738     }
1739
1740     /* Returns ints instead of longs where possible */
1741     INPLACE_UPDATE(numerator, PyNumber_Int(numerator));
1742     if (numerator == NULL) goto error;
1743     INPLACE_UPDATE(denominator, PyNumber_Int(denominator));
1744     if (denominator == NULL) goto error;
1745
1746     result_pair = PyTuple_Pack(2, numerator, denominator);
1747
1748 #undef INPLACE_UPDATE
1749 error:
1750     Py_XDECREF(py_exponent);
1751     Py_XDECREF(denominator);
1752     Py_XDECREF(numerator);
1753     return result_pair;
1754 }
1755
1756 PyDoc_STRVAR(float_as_integer_ratio_doc,
1757 "float.as_integer_ratio() -> (int, int)\n"
1758 "\n"
1759 "Returns a pair of integers, whose ratio is exactly equal to the original\n"
1760 "float and with a positive denominator.\n"
1761 "Raises OverflowError on infinities and a ValueError on NaNs.\n"
1762 "\n"
1763 ">>> (10.0).as_integer_ratio()\n"
1764 "(10, 1)\n"
1765 ">>> (0.0).as_integer_ratio()\n"
1766 "(0, 1)\n"
1767 ">>> (-.25).as_integer_ratio()\n"
1768 "(-1, 4)");
1769
1770
1771 static PyObject *
1772 float_subtype_new(PyTypeObject *type, PyObject *args, PyObject *kwds);
1773
1774 static PyObject *
1775 float_new(PyTypeObject *type, PyObject *args, PyObject *kwds)
1776 {
1777     PyObject *x = Py_False; /* Integer zero */
1778     static char *kwlist[] = {"x", 0};
1779
1780     if (type != &PyFloat_Type)
1781         return float_subtype_new(type, args, kwds); /* Wimp out */
1782     if (!PyArg_ParseTupleAndKeywords(args, kwds, "|O:float", kwlist, &x))
1783         return NULL;
1784     /* If it's a string, but not a string subclass, use
1785        PyFloat_FromString. */
1786     if (PyString_CheckExact(x))
1787         return PyFloat_FromString(x, NULL);
1788     return PyNumber_Float(x);
1789 }
1790
1791 /* Wimpy, slow approach to tp_new calls for subtypes of float:
1792    first create a regular float from whatever arguments we got,
1793    then allocate a subtype instance and initialize its ob_fval
1794    from the regular float.  The regular float is then thrown away.
1795 */
1796 static PyObject *
1797 float_subtype_new(PyTypeObject *type, PyObject *args, PyObject *kwds)
1798 {
1799     PyObject *tmp, *newobj;
1800
1801     assert(PyType_IsSubtype(type, &PyFloat_Type));
1802     tmp = float_new(&PyFloat_Type, args, kwds);
1803     if (tmp == NULL)
1804         return NULL;
1805     assert(PyFloat_CheckExact(tmp));
1806     newobj = type->tp_alloc(type, 0);
1807     if (newobj == NULL) {
1808         Py_DECREF(tmp);
1809         return NULL;
1810     }
1811     ((PyFloatObject *)newobj)->ob_fval = ((PyFloatObject *)tmp)->ob_fval;
1812     Py_DECREF(tmp);
1813     return newobj;
1814 }
1815
1816 static PyObject *
1817 float_getnewargs(PyFloatObject *v)
1818 {
1819     return Py_BuildValue("(d)", v->ob_fval);
1820 }
1821
1822 /* this is for the benefit of the pack/unpack routines below */
1823
1824 typedef enum {
1825     unknown_format, ieee_big_endian_format, ieee_little_endian_format
1826 } float_format_type;
1827
1828 static float_format_type double_format, float_format;
1829 static float_format_type detected_double_format, detected_float_format;
1830
1831 static PyObject *
1832 float_getformat(PyTypeObject *v, PyObject* arg)
1833 {
1834     char* s;
1835     float_format_type r;
1836
1837     if (!PyString_Check(arg)) {
1838         PyErr_Format(PyExc_TypeError,
1839          "__getformat__() argument must be string, not %.500s",
1840                          Py_TYPE(arg)->tp_name);
1841         return NULL;
1842     }
1843     s = PyString_AS_STRING(arg);
1844     if (strcmp(s, "double") == 0) {
1845         r = double_format;
1846     }
1847     else if (strcmp(s, "float") == 0) {
1848         r = float_format;
1849     }
1850     else {
1851         PyErr_SetString(PyExc_ValueError,
1852                         "__getformat__() argument 1 must be "
1853                         "'double' or 'float'");
1854         return NULL;
1855     }
1856
1857     switch (r) {
1858     case unknown_format:
1859         return PyString_FromString("unknown");
1860     case ieee_little_endian_format:
1861         return PyString_FromString("IEEE, little-endian");
1862     case ieee_big_endian_format:
1863         return PyString_FromString("IEEE, big-endian");
1864     default:
1865         Py_FatalError("insane float_format or double_format");
1866         return NULL;
1867     }
1868 }
1869
1870 PyDoc_STRVAR(float_getformat_doc,
1871 "float.__getformat__(typestr) -> string\n"
1872 "\n"
1873 "You probably don't want to use this function.  It exists mainly to be\n"
1874 "used in Python's test suite.\n"
1875 "\n"
1876 "typestr must be 'double' or 'float'.  This function returns whichever of\n"
1877 "'unknown', 'IEEE, big-endian' or 'IEEE, little-endian' best describes the\n"
1878 "format of floating point numbers used by the C type named by typestr.");
1879
1880 static PyObject *
1881 float_setformat(PyTypeObject *v, PyObject* args)
1882 {
1883     char* typestr;
1884     char* format;
1885     float_format_type f;
1886     float_format_type detected;
1887     float_format_type *p;
1888
1889     if (!PyArg_ParseTuple(args, "ss:__setformat__", &typestr, &format))
1890         return NULL;
1891
1892     if (strcmp(typestr, "double") == 0) {
1893         p = &double_format;
1894         detected = detected_double_format;
1895     }
1896     else if (strcmp(typestr, "float") == 0) {
1897         p = &float_format;
1898         detected = detected_float_format;
1899     }
1900     else {
1901         PyErr_SetString(PyExc_ValueError,
1902                         "__setformat__() argument 1 must "
1903                         "be 'double' or 'float'");
1904         return NULL;
1905     }
1906
1907     if (strcmp(format, "unknown") == 0) {
1908         f = unknown_format;
1909     }
1910     else if (strcmp(format, "IEEE, little-endian") == 0) {
1911         f = ieee_little_endian_format;
1912     }
1913     else if (strcmp(format, "IEEE, big-endian") == 0) {
1914         f = ieee_big_endian_format;
1915     }
1916     else {
1917         PyErr_SetString(PyExc_ValueError,
1918                         "__setformat__() argument 2 must be "
1919                         "'unknown', 'IEEE, little-endian' or "
1920                         "'IEEE, big-endian'");
1921         return NULL;
1922
1923     }
1924
1925     if (f != unknown_format && f != detected) {
1926         PyErr_Format(PyExc_ValueError,
1927                      "can only set %s format to 'unknown' or the "
1928                      "detected platform value", typestr);
1929         return NULL;
1930     }
1931
1932     *p = f;
1933     Py_RETURN_NONE;
1934 }
1935
1936 PyDoc_STRVAR(float_setformat_doc,
1937 "float.__setformat__(typestr, fmt) -> None\n"
1938 "\n"
1939 "You probably don't want to use this function.  It exists mainly to be\n"
1940 "used in Python's test suite.\n"
1941 "\n"
1942 "typestr must be 'double' or 'float'.  fmt must be one of 'unknown',\n"
1943 "'IEEE, big-endian' or 'IEEE, little-endian', and in addition can only be\n"
1944 "one of the latter two if it appears to match the underlying C reality.\n"
1945 "\n"
1946 "Overrides the automatic determination of C-level floating point type.\n"
1947 "This affects how floats are converted to and from binary strings.");
1948
1949 static PyObject *
1950 float_getzero(PyObject *v, void *closure)
1951 {
1952     return PyFloat_FromDouble(0.0);
1953 }
1954
1955 static PyObject *
1956 float__format__(PyObject *self, PyObject *args)
1957 {
1958     PyObject *format_spec;
1959
1960     if (!PyArg_ParseTuple(args, "O:__format__", &format_spec))
1961         return NULL;
1962     if (PyBytes_Check(format_spec))
1963         return _PyFloat_FormatAdvanced(self,
1964                                        PyBytes_AS_STRING(format_spec),
1965                                        PyBytes_GET_SIZE(format_spec));
1966     if (PyUnicode_Check(format_spec)) {
1967         /* Convert format_spec to a str */
1968         PyObject *result;
1969         PyObject *str_spec = PyObject_Str(format_spec);
1970
1971         if (str_spec == NULL)
1972             return NULL;
1973
1974         result = _PyFloat_FormatAdvanced(self,
1975                                          PyBytes_AS_STRING(str_spec),
1976                                          PyBytes_GET_SIZE(str_spec));
1977
1978         Py_DECREF(str_spec);
1979         return result;
1980     }
1981     PyErr_SetString(PyExc_TypeError, "__format__ requires str or unicode");
1982     return NULL;
1983 }
1984
1985 PyDoc_STRVAR(float__format__doc,
1986 "float.__format__(format_spec) -> string\n"
1987 "\n"
1988 "Formats the float according to format_spec.");
1989
1990
1991 static PyMethodDef float_methods[] = {
1992     {"conjugate",       (PyCFunction)float_float,       METH_NOARGS,
1993      "Returns self, the complex conjugate of any float."},
1994     {"__trunc__",       (PyCFunction)float_trunc, METH_NOARGS,
1995      "Returns the Integral closest to x between 0 and x."},
1996     {"as_integer_ratio", (PyCFunction)float_as_integer_ratio, METH_NOARGS,
1997      float_as_integer_ratio_doc},
1998     {"fromhex", (PyCFunction)float_fromhex,
1999      METH_O|METH_CLASS, float_fromhex_doc},
2000     {"hex", (PyCFunction)float_hex,
2001      METH_NOARGS, float_hex_doc},
2002     {"is_integer",      (PyCFunction)float_is_integer,  METH_NOARGS,
2003      "Returns True if the float is an integer."},
2004 #if 0
2005     {"is_inf",          (PyCFunction)float_is_inf,      METH_NOARGS,
2006      "Returns True if the float is positive or negative infinite."},
2007     {"is_finite",       (PyCFunction)float_is_finite,   METH_NOARGS,
2008      "Returns True if the float is finite, neither infinite nor NaN."},
2009     {"is_nan",          (PyCFunction)float_is_nan,      METH_NOARGS,
2010      "Returns True if the float is not a number (NaN)."},
2011 #endif
2012     {"__getnewargs__",          (PyCFunction)float_getnewargs,  METH_NOARGS},
2013     {"__getformat__",           (PyCFunction)float_getformat,
2014      METH_O|METH_CLASS,                 float_getformat_doc},
2015     {"__setformat__",           (PyCFunction)float_setformat,
2016      METH_VARARGS|METH_CLASS,           float_setformat_doc},
2017     {"__format__",          (PyCFunction)float__format__,
2018      METH_VARARGS,                  float__format__doc},
2019     {NULL,              NULL}           /* sentinel */
2020 };
2021
2022 static PyGetSetDef float_getset[] = {
2023     {"real",
2024      (getter)float_float, (setter)NULL,
2025      "the real part of a complex number",
2026      NULL},
2027     {"imag",
2028      (getter)float_getzero, (setter)NULL,
2029      "the imaginary part of a complex number",
2030      NULL},
2031     {NULL}  /* Sentinel */
2032 };
2033
2034 PyDoc_STRVAR(float_doc,
2035 "float(x) -> floating point number\n\
2036 \n\
2037 Convert a string or number to a floating point number, if possible.");
2038
2039
2040 static PyNumberMethods float_as_number = {
2041     float_add,          /*nb_add*/
2042     float_sub,          /*nb_subtract*/
2043     float_mul,          /*nb_multiply*/
2044     float_classic_div, /*nb_divide*/
2045     float_rem,          /*nb_remainder*/
2046     float_divmod,       /*nb_divmod*/
2047     float_pow,          /*nb_power*/
2048     (unaryfunc)float_neg, /*nb_negative*/
2049     (unaryfunc)float_float, /*nb_positive*/
2050     (unaryfunc)float_abs, /*nb_absolute*/
2051     (inquiry)float_nonzero, /*nb_nonzero*/
2052     0,                  /*nb_invert*/
2053     0,                  /*nb_lshift*/
2054     0,                  /*nb_rshift*/
2055     0,                  /*nb_and*/
2056     0,                  /*nb_xor*/
2057     0,                  /*nb_or*/
2058     float_coerce,       /*nb_coerce*/
2059     float_trunc,        /*nb_int*/
2060     float_long,         /*nb_long*/
2061     float_float,        /*nb_float*/
2062     0,                  /* nb_oct */
2063     0,                  /* nb_hex */
2064     0,                  /* nb_inplace_add */
2065     0,                  /* nb_inplace_subtract */
2066     0,                  /* nb_inplace_multiply */
2067     0,                  /* nb_inplace_divide */
2068     0,                  /* nb_inplace_remainder */
2069     0,                  /* nb_inplace_power */
2070     0,                  /* nb_inplace_lshift */
2071     0,                  /* nb_inplace_rshift */
2072     0,                  /* nb_inplace_and */
2073     0,                  /* nb_inplace_xor */
2074     0,                  /* nb_inplace_or */
2075     float_floor_div, /* nb_floor_divide */
2076     float_div,          /* nb_true_divide */
2077     0,                  /* nb_inplace_floor_divide */
2078     0,                  /* nb_inplace_true_divide */
2079 };
2080
2081 PyTypeObject PyFloat_Type = {
2082     PyVarObject_HEAD_INIT(&PyType_Type, 0)
2083     "float",
2084     sizeof(PyFloatObject),
2085     0,
2086     (destructor)float_dealloc,                  /* tp_dealloc */
2087     (printfunc)float_print,                     /* tp_print */
2088     0,                                          /* tp_getattr */
2089     0,                                          /* tp_setattr */
2090     0,                                          /* tp_compare */
2091     (reprfunc)float_repr,                       /* tp_repr */
2092     &float_as_number,                           /* tp_as_number */
2093     0,                                          /* tp_as_sequence */
2094     0,                                          /* tp_as_mapping */
2095     (hashfunc)float_hash,                       /* tp_hash */
2096     0,                                          /* tp_call */
2097     (reprfunc)float_str,                        /* tp_str */
2098     PyObject_GenericGetAttr,                    /* tp_getattro */
2099     0,                                          /* tp_setattro */
2100     0,                                          /* tp_as_buffer */
2101     Py_TPFLAGS_DEFAULT | Py_TPFLAGS_CHECKTYPES |
2102         Py_TPFLAGS_BASETYPE,                    /* tp_flags */
2103     float_doc,                                  /* tp_doc */
2104     0,                                          /* tp_traverse */
2105     0,                                          /* tp_clear */
2106     float_richcompare,                          /* tp_richcompare */
2107     0,                                          /* tp_weaklistoffset */
2108     0,                                          /* tp_iter */
2109     0,                                          /* tp_iternext */
2110     float_methods,                              /* tp_methods */
2111     0,                                          /* tp_members */
2112     float_getset,                               /* tp_getset */
2113     0,                                          /* tp_base */
2114     0,                                          /* tp_dict */
2115     0,                                          /* tp_descr_get */
2116     0,                                          /* tp_descr_set */
2117     0,                                          /* tp_dictoffset */
2118     0,                                          /* tp_init */
2119     0,                                          /* tp_alloc */
2120     float_new,                                  /* tp_new */
2121 };
2122
2123 void
2124 _PyFloat_Init(void)
2125 {
2126     /* We attempt to determine if this machine is using IEEE
2127        floating point formats by peering at the bits of some
2128        carefully chosen values.  If it looks like we are on an
2129        IEEE platform, the float packing/unpacking routines can
2130        just copy bits, if not they resort to arithmetic & shifts
2131        and masks.  The shifts & masks approach works on all finite
2132        values, but what happens to infinities, NaNs and signed
2133        zeroes on packing is an accident, and attempting to unpack
2134        a NaN or an infinity will raise an exception.
2135
2136        Note that if we're on some whacked-out platform which uses
2137        IEEE formats but isn't strictly little-endian or big-
2138        endian, we will fall back to the portable shifts & masks
2139        method. */
2140
2141 #if SIZEOF_DOUBLE == 8
2142     {
2143         double x = 9006104071832581.0;
2144         if (memcmp(&x, "\x43\x3f\xff\x01\x02\x03\x04\x05", 8) == 0)
2145             detected_double_format = ieee_big_endian_format;
2146         else if (memcmp(&x, "\x05\x04\x03\x02\x01\xff\x3f\x43", 8) == 0)
2147             detected_double_format = ieee_little_endian_format;
2148         else
2149             detected_double_format = unknown_format;
2150     }
2151 #else
2152     detected_double_format = unknown_format;
2153 #endif
2154
2155 #if SIZEOF_FLOAT == 4
2156     {
2157         float y = 16711938.0;
2158         if (memcmp(&y, "\x4b\x7f\x01\x02", 4) == 0)
2159             detected_float_format = ieee_big_endian_format;
2160         else if (memcmp(&y, "\x02\x01\x7f\x4b", 4) == 0)
2161             detected_float_format = ieee_little_endian_format;
2162         else
2163             detected_float_format = unknown_format;
2164     }
2165 #else
2166     detected_float_format = unknown_format;
2167 #endif
2168
2169     double_format = detected_double_format;
2170     float_format = detected_float_format;
2171
2172     /* Init float info */
2173     if (FloatInfoType.tp_name == 0)
2174         PyStructSequence_InitType(&FloatInfoType, &floatinfo_desc);
2175 }
2176
2177 int
2178 PyFloat_ClearFreeList(void)
2179 {
2180     PyFloatObject *p;
2181     PyFloatBlock *list, *next;
2182     int i;
2183     int u;                      /* remaining unfreed ints per block */
2184     int freelist_size = 0;
2185
2186     list = block_list;
2187     block_list = NULL;
2188     free_list = NULL;
2189     while (list != NULL) {
2190         u = 0;
2191         for (i = 0, p = &list->objects[0];
2192              i < N_FLOATOBJECTS;
2193              i++, p++) {
2194             if (PyFloat_CheckExact(p) && Py_REFCNT(p) != 0)
2195                 u++;
2196         }
2197         next = list->next;
2198         if (u) {
2199             list->next = block_list;
2200             block_list = list;
2201             for (i = 0, p = &list->objects[0];
2202                  i < N_FLOATOBJECTS;
2203                  i++, p++) {
2204                 if (!PyFloat_CheckExact(p) ||
2205                     Py_REFCNT(p) == 0) {
2206                     Py_TYPE(p) = (struct _typeobject *)
2207                         free_list;
2208                     free_list = p;
2209                 }
2210             }
2211         }
2212         else {
2213             PyMem_FREE(list);
2214         }
2215         freelist_size += u;
2216         list = next;
2217     }
2218     return freelist_size;
2219 }
2220
2221 void
2222 PyFloat_Fini(void)
2223 {
2224     PyFloatObject *p;
2225     PyFloatBlock *list;
2226     int i;
2227     int u;                      /* total unfreed floats per block */
2228
2229     u = PyFloat_ClearFreeList();
2230
2231     if (!Py_VerboseFlag)
2232         return;
2233     fprintf(stderr, "# cleanup floats");
2234     if (!u) {
2235         fprintf(stderr, "\n");
2236     }
2237     else {
2238         fprintf(stderr,
2239             ": %d unfreed float%s\n",
2240             u, u == 1 ? "" : "s");
2241     }
2242     if (Py_VerboseFlag > 1) {
2243         list = block_list;
2244         while (list != NULL) {
2245             for (i = 0, p = &list->objects[0];
2246                  i < N_FLOATOBJECTS;
2247                  i++, p++) {
2248                 if (PyFloat_CheckExact(p) &&
2249                     Py_REFCNT(p) != 0) {
2250                     char *buf = PyOS_double_to_string(
2251                         PyFloat_AS_DOUBLE(p), 'r',
2252                         0, 0, NULL);
2253                     if (buf) {
2254                         /* XXX(twouters) cast
2255                            refcount to long
2256                            until %zd is
2257                            universally
2258                            available
2259                         */
2260                         fprintf(stderr,
2261                  "#   <float at %p, refcnt=%ld, val=%s>\n",
2262                                     p, (long)Py_REFCNT(p), buf);
2263                                     PyMem_Free(buf);
2264                             }
2265                 }
2266             }
2267             list = list->next;
2268         }
2269     }
2270 }
2271
2272 /*----------------------------------------------------------------------------
2273  * _PyFloat_{Pack,Unpack}{4,8}.  See floatobject.h.
2274  */
2275 int
2276 _PyFloat_Pack4(double x, unsigned char *p, int le)
2277 {
2278     if (float_format == unknown_format) {
2279         unsigned char sign;
2280         int e;
2281         double f;
2282         unsigned int fbits;
2283         int incr = 1;
2284
2285         if (le) {
2286             p += 3;
2287             incr = -1;
2288         }
2289
2290         if (x < 0) {
2291             sign = 1;
2292             x = -x;
2293         }
2294         else
2295             sign = 0;
2296
2297         f = frexp(x, &e);
2298
2299         /* Normalize f to be in the range [1.0, 2.0) */
2300         if (0.5 <= f && f < 1.0) {
2301             f *= 2.0;
2302             e--;
2303         }
2304         else if (f == 0.0)
2305             e = 0;
2306         else {
2307             PyErr_SetString(PyExc_SystemError,
2308                             "frexp() result out of range");
2309             return -1;
2310         }
2311
2312         if (e >= 128)
2313             goto Overflow;
2314         else if (e < -126) {
2315             /* Gradual underflow */
2316             f = ldexp(f, 126 + e);
2317             e = 0;
2318         }
2319         else if (!(e == 0 && f == 0.0)) {
2320             e += 127;
2321             f -= 1.0; /* Get rid of leading 1 */
2322         }
2323
2324         f *= 8388608.0; /* 2**23 */
2325         fbits = (unsigned int)(f + 0.5); /* Round */
2326         assert(fbits <= 8388608);
2327         if (fbits >> 23) {
2328             /* The carry propagated out of a string of 23 1 bits. */
2329             fbits = 0;
2330             ++e;
2331             if (e >= 255)
2332                 goto Overflow;
2333         }
2334
2335         /* First byte */
2336         *p = (sign << 7) | (e >> 1);
2337         p += incr;
2338
2339         /* Second byte */
2340         *p = (char) (((e & 1) << 7) | (fbits >> 16));
2341         p += incr;
2342
2343         /* Third byte */
2344         *p = (fbits >> 8) & 0xFF;
2345         p += incr;
2346
2347         /* Fourth byte */
2348         *p = fbits & 0xFF;
2349
2350         /* Done */
2351         return 0;
2352
2353     }
2354     else {
2355         float y = (float)x;
2356         const char *s = (char*)&y;
2357         int i, incr = 1;
2358
2359         if (Py_IS_INFINITY(y) && !Py_IS_INFINITY(x))
2360             goto Overflow;
2361
2362         if ((float_format == ieee_little_endian_format && !le)
2363             || (float_format == ieee_big_endian_format && le)) {
2364             p += 3;
2365             incr = -1;
2366         }
2367
2368         for (i = 0; i < 4; i++) {
2369             *p = *s++;
2370             p += incr;
2371         }
2372         return 0;
2373     }
2374   Overflow:
2375     PyErr_SetString(PyExc_OverflowError,
2376                     "float too large to pack with f format");
2377     return -1;
2378 }
2379
2380 int
2381 _PyFloat_Pack8(double x, unsigned char *p, int le)
2382 {
2383     if (double_format == unknown_format) {
2384         unsigned char sign;
2385         int e;
2386         double f;
2387         unsigned int fhi, flo;
2388         int incr = 1;
2389
2390         if (le) {
2391             p += 7;
2392             incr = -1;
2393         }
2394
2395         if (x < 0) {
2396             sign = 1;
2397             x = -x;
2398         }
2399         else
2400             sign = 0;
2401
2402         f = frexp(x, &e);
2403
2404         /* Normalize f to be in the range [1.0, 2.0) */
2405         if (0.5 <= f && f < 1.0) {
2406             f *= 2.0;
2407             e--;
2408         }
2409         else if (f == 0.0)
2410             e = 0;
2411         else {
2412             PyErr_SetString(PyExc_SystemError,
2413                             "frexp() result out of range");
2414             return -1;
2415         }
2416
2417         if (e >= 1024)
2418             goto Overflow;
2419         else if (e < -1022) {
2420             /* Gradual underflow */
2421             f = ldexp(f, 1022 + e);
2422             e = 0;
2423         }
2424         else if (!(e == 0 && f == 0.0)) {
2425             e += 1023;
2426             f -= 1.0; /* Get rid of leading 1 */
2427         }
2428
2429         /* fhi receives the high 28 bits; flo the low 24 bits (== 52 bits) */
2430         f *= 268435456.0; /* 2**28 */
2431         fhi = (unsigned int)f; /* Truncate */
2432         assert(fhi < 268435456);
2433
2434         f -= (double)fhi;
2435         f *= 16777216.0; /* 2**24 */
2436         flo = (unsigned int)(f + 0.5); /* Round */
2437         assert(flo <= 16777216);
2438         if (flo >> 24) {
2439             /* The carry propagated out of a string of 24 1 bits. */
2440             flo = 0;
2441             ++fhi;
2442             if (fhi >> 28) {
2443                 /* And it also progagated out of the next 28 bits. */
2444                 fhi = 0;
2445                 ++e;
2446                 if (e >= 2047)
2447                     goto Overflow;
2448             }
2449         }
2450
2451         /* First byte */
2452         *p = (sign << 7) | (e >> 4);
2453         p += incr;
2454
2455         /* Second byte */
2456         *p = (unsigned char) (((e & 0xF) << 4) | (fhi >> 24));
2457         p += incr;
2458
2459         /* Third byte */
2460         *p = (fhi >> 16) & 0xFF;
2461         p += incr;
2462
2463         /* Fourth byte */
2464         *p = (fhi >> 8) & 0xFF;
2465         p += incr;
2466
2467         /* Fifth byte */
2468         *p = fhi & 0xFF;
2469         p += incr;
2470
2471         /* Sixth byte */
2472         *p = (flo >> 16) & 0xFF;
2473         p += incr;
2474
2475         /* Seventh byte */
2476         *p = (flo >> 8) & 0xFF;
2477         p += incr;
2478
2479         /* Eighth byte */
2480         *p = flo & 0xFF;
2481         /* p += incr; Unneeded (for now) */
2482
2483         /* Done */
2484         return 0;
2485
2486       Overflow:
2487         PyErr_SetString(PyExc_OverflowError,
2488                         "float too large to pack with d format");
2489         return -1;
2490     }
2491     else {
2492         const char *s = (char*)&x;
2493         int i, incr = 1;
2494
2495         if ((double_format == ieee_little_endian_format && !le)
2496             || (double_format == ieee_big_endian_format && le)) {
2497             p += 7;
2498             incr = -1;
2499         }
2500
2501         for (i = 0; i < 8; i++) {
2502             *p = *s++;
2503             p += incr;
2504         }
2505         return 0;
2506     }
2507 }
2508
2509 double
2510 _PyFloat_Unpack4(const unsigned char *p, int le)
2511 {
2512     if (float_format == unknown_format) {
2513         unsigned char sign;
2514         int e;
2515         unsigned int f;
2516         double x;
2517         int incr = 1;
2518
2519         if (le) {
2520             p += 3;
2521             incr = -1;
2522         }
2523
2524         /* First byte */
2525         sign = (*p >> 7) & 1;
2526         e = (*p & 0x7F) << 1;
2527         p += incr;
2528
2529         /* Second byte */
2530         e |= (*p >> 7) & 1;
2531         f = (*p & 0x7F) << 16;
2532         p += incr;
2533
2534         if (e == 255) {
2535             PyErr_SetString(
2536                 PyExc_ValueError,
2537                 "can't unpack IEEE 754 special value "
2538                 "on non-IEEE platform");
2539             return -1;
2540         }
2541
2542         /* Third byte */
2543         f |= *p << 8;
2544         p += incr;
2545
2546         /* Fourth byte */
2547         f |= *p;
2548
2549         x = (double)f / 8388608.0;
2550
2551         /* XXX This sadly ignores Inf/NaN issues */
2552         if (e == 0)
2553             e = -126;
2554         else {
2555             x += 1.0;
2556             e -= 127;
2557         }
2558         x = ldexp(x, e);
2559
2560         if (sign)
2561             x = -x;
2562
2563         return x;
2564     }
2565     else {
2566         float x;
2567
2568         if ((float_format == ieee_little_endian_format && !le)
2569             || (float_format == ieee_big_endian_format && le)) {
2570             char buf[4];
2571             char *d = &buf[3];
2572             int i;
2573
2574             for (i = 0; i < 4; i++) {
2575                 *d-- = *p++;
2576             }
2577             memcpy(&x, buf, 4);
2578         }
2579         else {
2580             memcpy(&x, p, 4);
2581         }
2582
2583         return x;
2584     }
2585 }
2586
2587 double
2588 _PyFloat_Unpack8(const unsigned char *p, int le)
2589 {
2590     if (double_format == unknown_format) {
2591         unsigned char sign;
2592         int e;
2593         unsigned int fhi, flo;
2594         double x;
2595         int incr = 1;
2596
2597         if (le) {
2598             p += 7;
2599             incr = -1;
2600         }
2601
2602         /* First byte */
2603         sign = (*p >> 7) & 1;
2604         e = (*p & 0x7F) << 4;
2605
2606         p += incr;
2607
2608         /* Second byte */
2609         e |= (*p >> 4) & 0xF;
2610         fhi = (*p & 0xF) << 24;
2611         p += incr;
2612
2613         if (e == 2047) {
2614             PyErr_SetString(
2615                 PyExc_ValueError,
2616                 "can't unpack IEEE 754 special value "
2617                 "on non-IEEE platform");
2618             return -1.0;
2619         }
2620
2621         /* Third byte */
2622         fhi |= *p << 16;
2623         p += incr;
2624
2625         /* Fourth byte */
2626         fhi |= *p  << 8;
2627         p += incr;
2628
2629         /* Fifth byte */
2630         fhi |= *p;
2631         p += incr;
2632
2633         /* Sixth byte */
2634         flo = *p << 16;
2635         p += incr;
2636
2637         /* Seventh byte */
2638         flo |= *p << 8;
2639         p += incr;
2640
2641         /* Eighth byte */
2642         flo |= *p;
2643
2644         x = (double)fhi + (double)flo / 16777216.0; /* 2**24 */
2645         x /= 268435456.0; /* 2**28 */
2646
2647         if (e == 0)
2648             e = -1022;
2649         else {
2650             x += 1.0;
2651             e -= 1023;
2652         }
2653         x = ldexp(x, e);
2654
2655         if (sign)
2656             x = -x;
2657
2658         return x;
2659     }
2660     else {
2661         double x;
2662
2663         if ((double_format == ieee_little_endian_format && !le)
2664             || (double_format == ieee_big_endian_format && le)) {
2665             char buf[8];
2666             char *d = &buf[7];
2667             int i;
2668
2669             for (i = 0; i < 8; i++) {
2670                 *d-- = *p++;
2671             }
2672             memcpy(&x, buf, 8);
2673         }
2674         else {
2675             memcpy(&x, p, 8);
2676         }
2677
2678         return x;
2679     }
2680 }