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