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
   9 #include <ctype.h>
  10
  11 #if !defined(__STDC__)
  12 extern double fmod(double, double);
  13 extern double pow(double, double);
  14 #endif
  15
  16 /* Special free list -- see comments for same code in intobject.c. */
  17 #define BLOCK_SIZE      1000    /* 1K less typical malloc overhead */
  18 #define BHEAD_SIZE      8       /* Enough for a 64-bit pointer */
  19 #define N_FLOATOBJECTS  ((BLOCK_SIZE - BHEAD_SIZE) / sizeof(PyFloatObject))
  20
  21 struct _floatblock {
  22         struct _floatblock *next;
  23         PyFloatObject objects[N_FLOATOBJECTS];
  24 };
  25
  26 typedef struct _floatblock PyFloatBlock;
  27
  28 static PyFloatBlock *block_list = NULL;
  29 static PyFloatObject *free_list = NULL;
  30
  31 static PyFloatObject *
  32 fill_free_list(void)
  33 {
  34         PyFloatObject *p, *q;
  35         /* XXX Float blocks escape the object heap. Use PyObject_MALLOC ??? */
  36         p = (PyFloatObject *) PyMem_MALLOC(sizeof(PyFloatBlock));
  37         if (p == NULL)
  38                 return (PyFloatObject *) PyErr_NoMemory();
  39         ((PyFloatBlock *)p)->next = block_list;
  40         block_list = (PyFloatBlock *)p;
  41         p = &((PyFloatBlock *)p)->objects[0];
  42         q = p + N_FLOATOBJECTS;
  43         while (--q > p)
  44                 q->ob_type = (struct _typeobject *)(q-1);
  45         q->ob_type = NULL;
  46         return p + N_FLOATOBJECTS - 1;
  47 }
  48
  49 PyObject *
  50 PyFloat_FromDouble(double fval)
  51 {
  52         register PyFloatObject *op;
  53         if (free_list == NULL) {
  54                 if ((free_list = fill_free_list()) == NULL)
  55                         return NULL;
  56         }
  57         /* Inline PyObject_New */
  58         op = free_list;
  59         free_list = (PyFloatObject *)op->ob_type;
  60         PyObject_INIT(op, &PyFloat_Type);
  61         op->ob_fval = fval;
  62         return (PyObject *) op;
  63 }
  64
  65 /**************************************************************************
  66 RED_FLAG 22-Sep-2000 tim
  67 PyFloat_FromString's pend argument is braindead.  Prior to this RED_FLAG,
  68
  69 1.  If v was a regular string, *pend was set to point to its terminating
  70     null byte.  That's useless (the caller can find that without any
  71     help from this function!).
  72
  73 2.  If v was a Unicode string, or an object convertible to a character
  74     buffer, *pend was set to point into stack trash (the auto temp
  75     vector holding the character buffer).  That was downright dangerous.
  76
  77 Since we can't change the interface of a public API function, pend is
  78 still supported but now *officially* useless:  if pend is not NULL,
  79 *pend is set to NULL.
  80 **************************************************************************/
  81 PyObject *
  82 PyFloat_FromString(PyObject *v, char **pend)
  83 {
  84         const char *s, *last, *end;
  85         double x;
  86         char buffer[256]; /* for errors */
  87 #ifdef Py_USING_UNICODE
  88         char s_buffer[256]; /* for objects convertible to a char buffer */
  89 #endif
  90         int len;
  91
  92         if (pend)
  93                 *pend = NULL;
  94         if (PyString_Check(v)) {
  95                 s = PyString_AS_STRING(v);
  96                 len = PyString_GET_SIZE(v);
  97         }
  98 #ifdef Py_USING_UNICODE
  99         else if (PyUnicode_Check(v)) {
 100                 if (PyUnicode_GET_SIZE(v) >= sizeof(s_buffer)) {
 101                         PyErr_SetString(PyExc_ValueError,
 102                                 "Unicode float() literal too long to convert");
 103                         return NULL;
 104                 }
 105                 if (PyUnicode_EncodeDecimal(PyUnicode_AS_UNICODE(v),
 106                                             PyUnicode_GET_SIZE(v),
 107                                             s_buffer,
 108                                             NULL))
 109                         return NULL;
 110                 s = s_buffer;
 111                 len = (int)strlen(s);
 112         }
 113 #endif
 114         else if (PyObject_AsCharBuffer(v, &s, &len)) {
 115                 PyErr_SetString(PyExc_TypeError,
 116                                 "float() argument must be a string or a number");
 117                 return NULL;
 118         }
 119
 120         last = s + len;
 121         while (*s && isspace(Py_CHARMASK(*s)))
 122                 s++;
 123         if (*s == '\0') {
 124                 PyErr_SetString(PyExc_ValueError, "empty string for float()");
 125                 return NULL;
 126         }
 127         /* We don't care about overflow or underflow.  If the platform supports
 128          * them, infinities and signed zeroes (on underflow) are fine.
 129          * However, strtod can return 0 for denormalized numbers, where atof
 130          * does not.  So (alas!) we special-case a zero result.  Note that
 131          * whether strtod sets errno on underflow is not defined, so we can't
 132          * key off errno.
 133          */
 134         PyFPE_START_PROTECT("strtod", return NULL)
 135         x = PyOS_ascii_strtod(s, (char **)&end);
 136         PyFPE_END_PROTECT(x)
 137         errno = 0;
 138         /* Believe it or not, Solaris 2.6 can move end *beyond* the null
 139            byte at the end of the string, when the input is inf(inity). */
 140         if (end > last)
 141                 end = last;
 142         if (end == s) {
 143                 PyOS_snprintf(buffer, sizeof(buffer),
 144                               "invalid literal for float(): %.200s", s);
 145                 PyErr_SetString(PyExc_ValueError, buffer);
 146                 return NULL;
 147         }
 148         /* Since end != s, the platform made *some* kind of sense out
 149            of the input.  Trust it. */
 150         while (*end && isspace(Py_CHARMASK(*end)))
 151                 end++;
 152         if (*end != '\0') {
 153                 PyOS_snprintf(buffer, sizeof(buffer),
 154                               "invalid literal for float(): %.200s", s);
 155                 PyErr_SetString(PyExc_ValueError, buffer);
 156                 return NULL;
 157         }
 158         else if (end != last) {
 159                 PyErr_SetString(PyExc_ValueError,
 160                                 "null byte in argument for float()");
 161                 return NULL;
 162         }
 163         if (x == 0.0) {
 164                 /* See above -- may have been strtod being anal
 165                    about denorms. */
 166                 PyFPE_START_PROTECT("atof", return NULL)
 167                 x = PyOS_ascii_atof(s);
 168                 PyFPE_END_PROTECT(x)
 169                 errno = 0;    /* whether atof ever set errno is undefined */
 170         }
 171         return PyFloat_FromDouble(x);
 172 }
 173
 174 static void
 175 float_dealloc(PyFloatObject *op)
 176 {
 177         if (PyFloat_CheckExact(op)) {
 178                 op->ob_type = (struct _typeobject *)free_list;
 179                 free_list = op;
 180         }
 181         else
 182                 op->ob_type->tp_free((PyObject *)op);
 183 }
 184
 185 double
 186 PyFloat_AsDouble(PyObject *op)
 187 {
 188         PyNumberMethods *nb;
 189         PyFloatObject *fo;
 190         double val;
 191
 192         if (op && PyFloat_Check(op))
 193                 return PyFloat_AS_DOUBLE((PyFloatObject*) op);
 194
 195         if (op == NULL) {
 196                 PyErr_BadArgument();
 197                 return -1;
 198         }
 199
 200         if ((nb = op->ob_type->tp_as_number) == NULL || nb->nb_float == NULL) {
 201                 PyErr_SetString(PyExc_TypeError, "a float is required");
 202                 return -1;
 203         }
 204
 205         fo = (PyFloatObject*) (*nb->nb_float) (op);
 206         if (fo == NULL)
 207                 return -1;
 208         if (!PyFloat_Check(fo)) {
 209                 PyErr_SetString(PyExc_TypeError,
 210                                 "nb_float should return float object");
 211                 return -1;
 212         }
 213
 214         val = PyFloat_AS_DOUBLE(fo);
 215         Py_DECREF(fo);
 216
 217         return val;
 218 }
 219
 220 /* Methods */
 221
 222 static void
 223 format_float(char *buf, size_t buflen, PyFloatObject *v, int precision)
 224 {
 225         register char *cp;
 226         char format[32];
 227         /* Subroutine for float_repr and float_print.
 228            We want float numbers to be recognizable as such,
 229            i.e., they should contain a decimal point or an exponent.
 230            However, %g may print the number as an integer;
 231            in such cases, we append ".0" to the string. */
 232
 233         assert(PyFloat_Check(v));
 234         PyOS_snprintf(format, 32, "%%.%ig", precision);
 235         PyOS_ascii_formatd(buf, buflen, format, v->ob_fval);
 236         cp = buf;
 237         if (*cp == '-')
 238                 cp++;
 239         for (; *cp != '\0'; cp++) {
 240                 /* Any non-digit means it's not an integer;
 241                    this takes care of NAN and INF as well. */
 242                 if (!isdigit(Py_CHARMASK(*cp)))
 243                         break;
 244         }
 245         if (*cp == '\0') {
 246                 *cp++ = '.';
 247                 *cp++ = '0';
 248                 *cp++ = '\0';
 249         }
 250 }
 251
 252 /* XXX PyFloat_AsStringEx should not be a public API function (for one
 253    XXX thing, its signature passes a buffer without a length; for another,
 254    XXX it isn't useful outside this file).
 255 */
 256 void
 257 PyFloat_AsStringEx(char *buf, PyFloatObject *v, int precision)
 258 {
 259         format_float(buf, 100, v, precision);
 260 }
 261
 262 /* Macro and helper that convert PyObject obj to a C double and store
 263    the value in dbl; this replaces the functionality of the coercion
 264    slot function.  If conversion to double raises an exception, obj is
 265    set to NULL, and the function invoking this macro returns NULL.  If
 266    obj is not of float, int or long type, Py_NotImplemented is incref'ed,
 267    stored in obj, and returned from the function invoking this macro.
 268 */
 269 #define CONVERT_TO_DOUBLE(obj, dbl)                     \
 270         if (PyFloat_Check(obj))                         \
 271                 dbl = PyFloat_AS_DOUBLE(obj);           \
 272         else if (convert_to_double(&(obj), &(dbl)) < 0) \
 273                 return obj;
 274
 275 static int
 276 convert_to_double(PyObject **v, double *dbl)
 277 {
 278         register PyObject *obj = *v;
 279
 280         if (PyInt_Check(obj)) {
 281                 *dbl = (double)PyInt_AS_LONG(obj);
 282         }
 283         else if (PyLong_Check(obj)) {
 284                 *dbl = PyLong_AsDouble(obj);
 285                 if (*dbl == -1.0 && PyErr_Occurred()) {
 286                         *v = NULL;
 287                         return -1;
 288                 }
 289         }
 290         else {
 291                 Py_INCREF(Py_NotImplemented);
 292                 *v = Py_NotImplemented;
 293                 return -1;
 294         }
 295         return 0;
 296 }
 297
 298 /* Precisions used by repr() and str(), respectively.
 299
 300    The repr() precision (17 significant decimal digits) is the minimal number
 301    that is guaranteed to have enough precision so that if the number is read
 302    back in the exact same binary value is recreated.  This is true for IEEE
 303    floating point by design, and also happens to work for all other modern
 304    hardware.
 305
 306    The str() precision is chosen so that in most cases, the rounding noise
 307    created by various operations is suppressed, while giving plenty of
 308    precision for practical use.
 309
 310 */
 311
 312 #define PREC_REPR       17
 313 #define PREC_STR        12
 314
 315 /* XXX PyFloat_AsString and PyFloat_AsReprString should be deprecated:
 316    XXX they pass a char buffer without passing a length.
 317 */
 318 void
 319 PyFloat_AsString(char *buf, PyFloatObject *v)
 320 {
 321         format_float(buf, 100, v, PREC_STR);
 322 }
 323
 324 void
 325 PyFloat_AsReprString(char *buf, PyFloatObject *v)
 326 {
 327         format_float(buf, 100, v, PREC_REPR);
 328 }
 329
 330 /* ARGSUSED */
 331 static int
 332 float_print(PyFloatObject *v, FILE *fp, int flags)
 333 {
 334         char buf[100];
 335         format_float(buf, sizeof(buf), v,
 336                      (flags & Py_PRINT_RAW) ? PREC_STR : PREC_REPR);
 337         fputs(buf, fp);
 338         return 0;
 339 }
 340
 341 static PyObject *
 342 float_repr(PyFloatObject *v)
 343 {
 344         char buf[100];
 345         format_float(buf, sizeof(buf), v, PREC_REPR);
 346         return PyString_FromString(buf);
 347 }
 348
 349 static PyObject *
 350 float_str(PyFloatObject *v)
 351 {
 352         char buf[100];
 353         format_float(buf, sizeof(buf), v, PREC_STR);
 354         return PyString_FromString(buf);
 355 }
 356
 357 /* Comparison is pretty much a nightmare.  When comparing float to float,
 358  * we do it as straightforwardly (and long-windedly) as conceivable, so
 359  * that, e.g., Python x == y delivers the same result as the platform
 360  * C x == y when x and/or y is a NaN.
 361  * When mixing float with an integer type, there's no good *uniform* approach.
 362  * Converting the double to an integer obviously doesn't work, since we
 363  * may lose info from fractional bits.  Converting the integer to a double
 364  * also has two failure modes:  (1) a long int may trigger overflow (too
 365  * large to fit in the dynamic range of a C double); (2) even a C long may have
 366  * more bits than fit in a C double (e.g., on a a 64-bit box long may have
 367  * 63 bits of precision, but a C double probably has only 53), and then
 368  * we can falsely claim equality when low-order integer bits are lost by
 369  * coercion to double.  So this part is painful too.
 370  */
 371
 372 static PyObject*
 373 float_richcompare(PyObject *v, PyObject *w, int op)
 374 {
 375         double i, j;
 376         int r = 0;
 377
 378         assert(PyFloat_Check(v));
 379         i = PyFloat_AS_DOUBLE(v);
 380
 381         /* Switch on the type of w.  Set i and j to doubles to be compared,
 382          * and op to the richcomp to use.
 383          */
 384         if (PyFloat_Check(w))
 385                 j = PyFloat_AS_DOUBLE(w);
 386
 387         else if (Py_IS_INFINITY(i) || Py_IS_NAN(i)) {
 388                 if (PyInt_Check(w) || PyLong_Check(w))
 389                         /* If i is an infinity, its magnitude exceeds any
 390                          * finite integer, so it doesn't matter which int we
 391                          * compare i with.  If i is a NaN, similarly.
 392                          */
 393                         j = 0.0;
 394                 else
 395                         goto Unimplemented;
 396         }
 397
 398         else if (PyInt_Check(w)) {
 399                 long jj = PyInt_AS_LONG(w);
 400                 /* In the worst realistic case I can imagine, C double is a
 401                  * Cray single with 48 bits of precision, and long has 64
 402                  * bits.
 403                  */
 404 #if SIZEOF_LONG > 6
 405                 unsigned long abs = (unsigned long)(jj < 0 ? -jj : jj);
 406                 if (abs >> 48) {
 407                         /* Needs more than 48 bits.  Make it take the
 408                          * PyLong path.
 409                          */
 410                         PyObject *result;
 411                         PyObject *ww = PyLong_FromLong(jj);
 412
 413                         if (ww == NULL)
 414                                 return NULL;
 415                         result = float_richcompare(v, ww, op);
 416                         Py_DECREF(ww);
 417                         return result;
 418                 }
 419 #endif
 420                 j = (double)jj;
 421                 assert((long)j == jj);
 422         }
 423
 424         else if (PyLong_Check(w)) {
 425                 int vsign = i == 0.0 ? 0 : i < 0.0 ? -1 : 1;
 426                 int wsign = _PyLong_Sign(w);
 427                 size_t nbits;
 428                 int exponent;
 429
 430                 if (vsign != wsign) {
 431                         /* Magnitudes are irrelevant -- the signs alone
 432                          * determine the outcome.
 433                          */
 434                         i = (double)vsign;
 435                         j = (double)wsign;
 436                         goto Compare;
 437                 }
 438                 /* The signs are the same. */
 439                 /* Convert w to a double if it fits.  In particular, 0 fits. */
 440                 nbits = _PyLong_NumBits(w);
 441                 if (nbits == (size_t)-1 && PyErr_Occurred()) {
 442                         /* This long is so large that size_t isn't big enough
 443                          * to hold the # of bits.  Replace with little doubles
 444                          * that give the same outcome -- w is so large that
 445                          * its magnitude must exceed the magnitude of any
 446                          * finite float.
 447                          */
 448                         PyErr_Clear();
 449                         i = (double)vsign;
 450                         assert(wsign != 0);
 451                         j = wsign * 2.0;
 452                         goto Compare;
 453                 }
 454                 if (nbits <= 48) {
 455                         j = PyLong_AsDouble(w);
 456                         /* It's impossible that <= 48 bits overflowed. */
 457                         assert(j != -1.0 || ! PyErr_Occurred());
 458                         goto Compare;
 459                 }
 460                 assert(wsign != 0); /* else nbits was 0 */
 461                 assert(vsign != 0); /* if vsign were 0, then since wsign is
 462                                      * not 0, we would have taken the
 463                                      * vsign != wsign branch at the start */
 464                 /* We want to work with non-negative numbers. */
 465                 if (vsign < 0) {
 466                         /* "Multiply both sides" by -1; this also swaps the
 467                          * comparator.
 468                          */
 469                         i = -i;
 470                         op = _Py_SwappedOp[op];
 471                 }
 472                 assert(i > 0.0);
 473                 (void) frexp(i, &exponent);
 474                 /* exponent is the # of bits in v before the radix point;
 475                  * we know that nbits (the # of bits in w) > 48 at this point
 476                  */
 477                 if (exponent < 0 || (size_t)exponent < nbits) {
 478                         i = 1.0;
 479                         j = 2.0;
 480                         goto Compare;
 481                 }
 482                 if ((size_t)exponent > nbits) {
 483                         i = 2.0;
 484                         j = 1.0;
 485                         goto Compare;
 486                 }
 487                 /* v and w have the same number of bits before the radix
 488                  * point.  Construct two longs that have the same comparison
 489                  * outcome.
 490                  */
 491                 {
 492                         double fracpart;
 493                         double intpart;
 494                         PyObject *result = NULL;
 495                         PyObject *one = NULL;
 496                         PyObject *vv = NULL;
 497                         PyObject *ww = w;
 498
 499                         if (wsign < 0) {
 500                                 ww = PyNumber_Negative(w);
 501                                 if (ww == NULL)
 502                                         goto Error;
 503                         }
 504                         else
 505                                 Py_INCREF(ww);
 506
 507                         fracpart = modf(i, &intpart);
 508                         vv = PyLong_FromDouble(intpart);
 509                         if (vv == NULL)
 510                                 goto Error;
 511
 512                         if (fracpart != 0.0) {
 513                                 /* Shift left, and or a 1 bit into vv
 514                                  * to represent the lost fraction.
 515                                  */
 516                                 PyObject *temp;
 517
 518                                 one = PyInt_FromLong(1);
 519                                 if (one == NULL)
 520                                         goto Error;
 521
 522                                 temp = PyNumber_Lshift(ww, one);
 523                                 if (temp == NULL)
 524                                         goto Error;
 525                                 Py_DECREF(ww);
 526                                 ww = temp;
 527
 528                                 temp = PyNumber_Lshift(vv, one);
 529                                 if (temp == NULL)
 530                                         goto Error;
 531                                 Py_DECREF(vv);
 532                                 vv = temp;
 533
 534                                 temp = PyNumber_Or(vv, one);
 535                                 if (temp == NULL)
 536                                         goto Error;
 537                                 Py_DECREF(vv);
 538                                 vv = temp;
 539                         }
 540
 541                         r = PyObject_RichCompareBool(vv, ww, op);
 542                         if (r < 0)
 543                                 goto Error;
 544                         result = PyBool_FromLong(r);
 545                  Error:
 546                         Py_XDECREF(vv);
 547                         Py_XDECREF(ww);
 548                         Py_XDECREF(one);
 549                         return result;
 550                 }
 551         } /* else if (PyLong_Check(w)) */
 552
 553         else    /* w isn't float, int, or long */
 554                 goto Unimplemented;
 555
 556  Compare:
 557         PyFPE_START_PROTECT("richcompare", return NULL)
 558         switch (op) {
 559         case Py_EQ:
 560                 r = i == j;
 561                 break;
 562         case Py_NE:
 563                 r = i != j;
 564                 break;
 565         case Py_LE:
 566                 r = i <= j;
 567                 break;
 568         case Py_GE:
 569                 r = i >= j;
 570                 break;
 571         case Py_LT:
 572                 r = i < j;
 573                 break;
 574         case Py_GT:
 575                 r = i > j;
 576                 break;
 577         }
 578         PyFPE_END_PROTECT(r)
 579         return PyBool_FromLong(r);
 580
 581  Unimplemented:
 582         Py_INCREF(Py_NotImplemented);
 583         return Py_NotImplemented;
 584 }
 585
 586 static long
 587 float_hash(PyFloatObject *v)
 588 {
 589         return _Py_HashDouble(v->ob_fval);
 590 }
 591
 592 static PyObject *
 593 float_add(PyObject *v, PyObject *w)
 594 {
 595         double a,b;
 596         CONVERT_TO_DOUBLE(v, a);
 597         CONVERT_TO_DOUBLE(w, b);
 598         PyFPE_START_PROTECT("add", return 0)
 599         a = a + b;
 600         PyFPE_END_PROTECT(a)
 601         return PyFloat_FromDouble(a);
 602 }
 603
 604 static PyObject *
 605 float_sub(PyObject *v, PyObject *w)
 606 {
 607         double a,b;
 608         CONVERT_TO_DOUBLE(v, a);
 609         CONVERT_TO_DOUBLE(w, b);
 610         PyFPE_START_PROTECT("subtract", return 0)
 611         a = a - b;
 612         PyFPE_END_PROTECT(a)
 613         return PyFloat_FromDouble(a);
 614 }
 615
 616 static PyObject *
 617 float_mul(PyObject *v, PyObject *w)
 618 {
 619         double a,b;
 620         CONVERT_TO_DOUBLE(v, a);
 621         CONVERT_TO_DOUBLE(w, b);
 622         PyFPE_START_PROTECT("multiply", return 0)
 623         a = a * b;
 624         PyFPE_END_PROTECT(a)
 625         return PyFloat_FromDouble(a);
 626 }
 627
 628 static PyObject *
 629 float_div(PyObject *v, PyObject *w)
 630 {
 631         double a,b;
 632         CONVERT_TO_DOUBLE(v, a);
 633         CONVERT_TO_DOUBLE(w, b);
 634         if (b == 0.0) {
 635                 PyErr_SetString(PyExc_ZeroDivisionError, "float division");
 636                 return NULL;
 637         }
 638         PyFPE_START_PROTECT("divide", return 0)
 639         a = a / b;
 640         PyFPE_END_PROTECT(a)
 641         return PyFloat_FromDouble(a);
 642 }
 643
 644 static PyObject *
 645 float_classic_div(PyObject *v, PyObject *w)
 646 {
 647         double a,b;
 648         CONVERT_TO_DOUBLE(v, a);
 649         CONVERT_TO_DOUBLE(w, b);
 650         if (Py_DivisionWarningFlag >= 2 &&
 651             PyErr_Warn(PyExc_DeprecationWarning, "classic float division") < 0)
 652                 return NULL;
 653         if (b == 0.0) {
 654                 PyErr_SetString(PyExc_ZeroDivisionError, "float division");
 655                 return NULL;
 656         }
 657         PyFPE_START_PROTECT("divide", return 0)
 658         a = a / b;
 659         PyFPE_END_PROTECT(a)
 660         return PyFloat_FromDouble(a);
 661 }
 662
 663 static PyObject *
 664 float_rem(PyObject *v, PyObject *w)
 665 {
 666         double vx, wx;
 667         double mod;
 668         CONVERT_TO_DOUBLE(v, vx);
 669         CONVERT_TO_DOUBLE(w, wx);
 670         if (wx == 0.0) {
 671                 PyErr_SetString(PyExc_ZeroDivisionError, "float modulo");
 672                 return NULL;
 673         }
 674         PyFPE_START_PROTECT("modulo", return 0)
 675         mod = fmod(vx, wx);
 676         /* note: checking mod*wx < 0 is incorrect -- underflows to
 677            0 if wx < sqrt(smallest nonzero double) */
 678         if (mod && ((wx < 0) != (mod < 0))) {
 679                 mod += wx;
 680         }
 681         PyFPE_END_PROTECT(mod)
 682         return PyFloat_FromDouble(mod);
 683 }
 684
 685 static PyObject *
 686 float_divmod(PyObject *v, PyObject *w)
 687 {
 688         double vx, wx;
 689         double div, mod, floordiv;
 690         CONVERT_TO_DOUBLE(v, vx);
 691         CONVERT_TO_DOUBLE(w, wx);
 692         if (wx == 0.0) {
 693                 PyErr_SetString(PyExc_ZeroDivisionError, "float divmod()");
 694                 return NULL;
 695         }
 696         PyFPE_START_PROTECT("divmod", return 0)
 697         mod = fmod(vx, wx);
 698         /* fmod is typically exact, so vx-mod is *mathematically* an
 699            exact multiple of wx.  But this is fp arithmetic, and fp
 700            vx - mod is an approximation; the result is that div may
 701            not be an exact integral value after the division, although
 702            it will always be very close to one.
 703         */
 704         div = (vx - mod) / wx;
 705         if (mod) {
 706                 /* ensure the remainder has the same sign as the denominator */
 707                 if ((wx < 0) != (mod < 0)) {
 708                         mod += wx;
 709                         div -= 1.0;
 710                 }
 711         }
 712         else {
 713                 /* the remainder is zero, and in the presence of signed zeroes
 714                    fmod returns different results across platforms; ensure
 715                    it has the same sign as the denominator; we'd like to do
 716                    "mod = wx * 0.0", but that may get optimized away */
 717                 mod *= mod;  /* hide "mod = +0" from optimizer */
 718                 if (wx < 0.0)
 719                         mod = -mod;
 720         }
 721         /* snap quotient to nearest integral value */
 722         if (div) {
 723                 floordiv = floor(div);
 724                 if (div - floordiv > 0.5)
 725                         floordiv += 1.0;
 726         }
 727         else {
 728                 /* div is zero - get the same sign as the true quotient */
 729                 div *= div;     /* hide "div = +0" from optimizers */
 730                 floordiv = div * vx / wx; /* zero w/ sign of vx/wx */
 731         }
 732         PyFPE_END_PROTECT(floordiv)
 733         return Py_BuildValue("(dd)", floordiv, mod);
 734 }
 735
 736 static PyObject *
 737 float_floor_div(PyObject *v, PyObject *w)
 738 {
 739         PyObject *t, *r;
 740
 741         t = float_divmod(v, w);
 742         if (t == NULL || t == Py_NotImplemented)
 743                 return t;
 744         assert(PyTuple_CheckExact(t));
 745         r = PyTuple_GET_ITEM(t, 0);
 746         Py_INCREF(r);
 747         Py_DECREF(t);
 748         return r;
 749 }
 750
 751 static PyObject *
 752 float_pow(PyObject *v, PyObject *w, PyObject *z)
 753 {
 754         double iv, iw, ix;
 755
 756         if ((PyObject *)z != Py_None) {
 757                 PyErr_SetString(PyExc_TypeError, "pow() 3rd argument not "
 758                         "allowed unless all arguments are integers");
 759                 return NULL;
 760         }
 761
 762         CONVERT_TO_DOUBLE(v, iv);
 763         CONVERT_TO_DOUBLE(w, iw);
 764
 765         /* Sort out special cases here instead of relying on pow() */
 766         if (iw == 0) {          /* v**0 is 1, even 0**0 */
 767                 PyFPE_START_PROTECT("pow", return NULL)
 768                 if ((PyObject *)z != Py_None) {
 769                         double iz;
 770                         CONVERT_TO_DOUBLE(z, iz);
 771                         ix = fmod(1.0, iz);
 772                         if (ix != 0 && iz < 0)
 773                                 ix += iz;
 774                 }
 775                 else
 776                         ix = 1.0;
 777                 PyFPE_END_PROTECT(ix)
 778                 return PyFloat_FromDouble(ix);
 779         }
 780         if (iv == 0.0) {  /* 0**w is error if w<0, else 1 */
 781                 if (iw < 0.0) {
 782                         PyErr_SetString(PyExc_ZeroDivisionError,
 783                                         "0.0 cannot be raised to a negative power");
 784                         return NULL;
 785                 }
 786                 return PyFloat_FromDouble(0.0);
 787         }
 788         if (iv < 0.0) {
 789                 /* Whether this is an error is a mess, and bumps into libm
 790                  * bugs so we have to figure it out ourselves.
 791                  */
 792                 if (iw != floor(iw)) {
 793                         PyErr_SetString(PyExc_ValueError, "negative number "
 794                                 "cannot be raised to a fractional power");
 795                         return NULL;
 796                 }
 797                 /* iw is an exact integer, albeit perhaps a very large one.
 798                  * -1 raised to an exact integer should never be exceptional.
 799                  * Alas, some libms (chiefly glibc as of early 2003) return
 800                  * NaN and set EDOM on pow(-1, large_int) if the int doesn't
 801                  * happen to be representable in a *C* integer.  That's a
 802                  * bug; we let that slide in math.pow() (which currently
 803                  * reflects all platform accidents), but not for Python's **.
 804                  */
 805                  if (iv == -1.0 && !Py_IS_INFINITY(iw) && iw == iw) {
 806                         /* XXX the "iw == iw" was to weed out NaNs.  This
 807                          * XXX doesn't actually work on all platforms.
 808                          */
 809                         /* Return 1 if iw is even, -1 if iw is odd; there's
 810                          * no guarantee that any C integral type is big
 811                          * enough to hold iw, so we have to check this
 812                          * indirectly.
 813                          */
 814                         ix = floor(iw * 0.5) * 2.0;
 815                         return PyFloat_FromDouble(ix == iw ? 1.0 : -1.0);
 816                 }
 817                 /* Else iv != -1.0, and overflow or underflow are possible.
 818                  * Unless we're to write pow() ourselves, we have to trust
 819                  * the platform to do this correctly.
 820                  */
 821         }
 822         errno = 0;
 823         PyFPE_START_PROTECT("pow", return NULL)
 824         ix = pow(iv, iw);
 825         PyFPE_END_PROTECT(ix)
 826         Py_ADJUST_ERANGE1(ix);
 827         if (errno != 0) {
 828                 /* We don't expect any errno value other than ERANGE, but
 829                  * the range of libm bugs appears unbounded.
 830                  */
 831                 PyErr_SetFromErrno(errno == ERANGE ? PyExc_OverflowError :
 832                                                      PyExc_ValueError);
 833                 return NULL;
 834         }
 835         return PyFloat_FromDouble(ix);
 836 }
 837
 838 static PyObject *
 839 float_neg(PyFloatObject *v)
 840 {
 841         return PyFloat_FromDouble(-v->ob_fval);
 842 }
 843
 844 static PyObject *
 845 float_pos(PyFloatObject *v)
 846 {
 847         if (PyFloat_CheckExact(v)) {
 848                 Py_INCREF(v);
 849                 return (PyObject *)v;
 850         }
 851         else
 852                 return PyFloat_FromDouble(v->ob_fval);
 853 }
 854
 855 static PyObject *
 856 float_abs(PyFloatObject *v)
 857 {
 858         return PyFloat_FromDouble(fabs(v->ob_fval));
 859 }
 860
 861 static int
 862 float_nonzero(PyFloatObject *v)
 863 {
 864         return v->ob_fval != 0.0;
 865 }
 866
 867 static int
 868 float_coerce(PyObject **pv, PyObject **pw)
 869 {
 870         if (PyInt_Check(*pw)) {
 871                 long x = PyInt_AsLong(*pw);
 872                 *pw = PyFloat_FromDouble((double)x);
 873                 Py_INCREF(*pv);
 874                 return 0;
 875         }
 876         else if (PyLong_Check(*pw)) {
 877                 double x = PyLong_AsDouble(*pw);
 878                 if (x == -1.0 && PyErr_Occurred())
 879                         return -1;
 880                 *pw = PyFloat_FromDouble(x);
 881                 Py_INCREF(*pv);
 882                 return 0;
 883         }
 884         else if (PyFloat_Check(*pw)) {
 885                 Py_INCREF(*pv);
 886                 Py_INCREF(*pw);
 887                 return 0;
 888         }
 889         return 1; /* Can't do it */
 890 }
 891
 892 static PyObject *
 893 float_long(PyObject *v)
 894 {
 895         double x = PyFloat_AsDouble(v);
 896         return PyLong_FromDouble(x);
 897 }
 898
 899 static PyObject *
 900 float_int(PyObject *v)
 901 {
 902         double x = PyFloat_AsDouble(v);
 903         double wholepart;       /* integral portion of x, rounded toward 0 */
 904
 905         (void)modf(x, &wholepart);
 906         /* Try to get out cheap if this fits in a Python int.  The attempt
 907          * to cast to long must be protected, as C doesn't define what
 908          * happens if the double is too big to fit in a long.  Some rare
 909          * systems raise an exception then (RISCOS was mentioned as one,
 910          * and someone using a non-default option on Sun also bumped into
 911          * that).  Note that checking for >= and <= LONG_{MIN,MAX} would
 912          * still be vulnerable:  if a long has more bits of precision than
 913          * a double, casting MIN/MAX to double may yield an approximation,
 914          * and if that's rounded up, then, e.g., wholepart=LONG_MAX+1 would
 915          * yield true from the C expression wholepart<=LONG_MAX, despite
 916          * that wholepart is actually greater than LONG_MAX.
 917          */
 918         if (LONG_MIN < wholepart && wholepart < LONG_MAX) {
 919                 const long aslong = (long)wholepart;
 920                 return PyInt_FromLong(aslong);
 921         }
 922         return PyLong_FromDouble(wholepart);
 923 }
 924
 925 static PyObject *
 926 float_float(PyObject *v)
 927 {
 928         if (PyFloat_CheckExact(v))
 929                 Py_INCREF(v);
 930         else
 931                 v = PyFloat_FromDouble(((PyFloatObject *)v)->ob_fval);
 932         return v;
 933 }
 934
 935
 936 static PyObject *
 937 float_subtype_new(PyTypeObject *type, PyObject *args, PyObject *kwds);
 938
 939 static PyObject *
 940 float_new(PyTypeObject *type, PyObject *args, PyObject *kwds)
 941 {
 942         PyObject *x = Py_False; /* Integer zero */
 943         static const char *kwlist[] = {"x", 0};
 944
 945         if (type != &PyFloat_Type)
 946                 return float_subtype_new(type, args, kwds); /* Wimp out */
 947         if (!PyArg_ParseTupleAndKeywords(args, kwds, "|O:float", kwlist, &x))
 948                 return NULL;
 949         if (PyString_Check(x))
 950                 return PyFloat_FromString(x, NULL);
 951         return PyNumber_Float(x);
 952 }
 953
 954 /* Wimpy, slow approach to tp_new calls for subtypes of float:
 955    first create a regular float from whatever arguments we got,
 956    then allocate a subtype instance and initialize its ob_fval
 957    from the regular float.  The regular float is then thrown away.
 958 */
 959 static PyObject *
 960 float_subtype_new(PyTypeObject *type, PyObject *args, PyObject *kwds)
 961 {
 962         PyObject *tmp, *new;
 963
 964         assert(PyType_IsSubtype(type, &PyFloat_Type));
 965         tmp = float_new(&PyFloat_Type, args, kwds);
 966         if (tmp == NULL)
 967                 return NULL;
 968         assert(PyFloat_CheckExact(tmp));
 969         new = type->tp_alloc(type, 0);
 970         if (new == NULL) {
 971                 Py_DECREF(tmp);
 972                 return NULL;
 973         }
 974         ((PyFloatObject *)new)->ob_fval = ((PyFloatObject *)tmp)->ob_fval;
 975         Py_DECREF(tmp);
 976         return new;
 977 }
 978
 979 static PyObject *
 980 float_getnewargs(PyFloatObject *v)
 981 {
 982         return Py_BuildValue("(d)", v->ob_fval);
 983 }
 984
 985 /* this is for the benefit of the pack/unpack routines below */
 986
 987 typedef enum {
 988         unknown_format, ieee_big_endian_format, ieee_little_endian_format
 989 } float_format_type;
 990
 991 static float_format_type double_format, float_format;
 992 static float_format_type detected_double_format, detected_float_format;
 993
 994 static PyObject *
 995 float_getformat(PyTypeObject *v, PyObject* arg)
 996 {
 997         char* s;
 998         float_format_type r;
 999
1000         if (!PyString_Check(arg)) {
1001                 PyErr_Format(PyExc_TypeError,
1002              "__getformat__() argument must be string, not %.500s",
1003                              arg->ob_type->tp_name);
1004                 return NULL;
1005         }
1006         s = PyString_AS_STRING(arg);
1007         if (strcmp(s, "double") == 0) {
1008                 r = double_format;
1009         }
1010         else if (strcmp(s, "float") == 0) {
1011                 r = float_format;
1012         }
1013         else {
1014                 PyErr_SetString(PyExc_ValueError,
1015                                 "__getformat__() argument 1 must be "
1016                                 "'double' or 'float'");
1017                 return NULL;
1018         }
1019
1020         switch (r) {
1021         case unknown_format:
1022                 return PyString_FromString("unknown");
1023         case ieee_little_endian_format:
1024                 return PyString_FromString("IEEE, little-endian");
1025         case ieee_big_endian_format:
1026                 return PyString_FromString("IEEE, big-endian");
1027         default:
1028                 Py_FatalError("insane float_format or double_format");
1029                 return NULL;
1030         }
1031 }
1032
1033 PyDoc_STRVAR(float_getformat_doc,
1034 "float.__getformat__(typestr) -> string\n"
1035 "\n"
1036 "You probably don't want to use this function.  It exists mainly to be\n"
1037 "used in Python's test suite.\n"
1038 "\n"
1039 "typestr must be 'double' or 'float'.  This function returns whichever of\n"
1040 "'unknown', 'IEEE, big-endian' or 'IEEE, little-endian' best describes the\n"
1041 "format of floating point numbers used by the C type named by typestr.");
1042
1043 static PyObject *
1044 float_setformat(PyTypeObject *v, PyObject* args)
1045 {
1046         char* typestr;
1047         char* format;
1048         float_format_type f;
1049         float_format_type detected;
1050         float_format_type *p;
1051
1052         if (!PyArg_ParseTuple(args, "ss:__setformat__", &typestr, &format))
1053                 return NULL;
1054
1055         if (strcmp(typestr, "double") == 0) {
1056                 p = &double_format;
1057                 detected = detected_double_format;
1058         }
1059         else if (strcmp(typestr, "float") == 0) {
1060                 p = &float_format;
1061                 detected = detected_float_format;
1062         }
1063         else {
1064                 PyErr_SetString(PyExc_ValueError,
1065                                 "__setformat__() argument 1 must "
1066                                 "be 'double' or 'float'");
1067                 return NULL;
1068         }
1069
1070         if (strcmp(format, "unknown") == 0) {
1071                 f = unknown_format;
1072         }
1073         else if (strcmp(format, "IEEE, little-endian") == 0) {
1074                 f = ieee_little_endian_format;
1075         }
1076         else if (strcmp(format, "IEEE, big-endian") == 0) {
1077                 f = ieee_big_endian_format;
1078         }
1079         else {
1080                 PyErr_SetString(PyExc_ValueError,
1081                                 "__setformat__() argument 2 must be "
1082                                 "'unknown', 'IEEE, little-endian' or "
1083                                 "'IEEE, big-endian'");
1084                 return NULL;
1085
1086         }
1087
1088         if (f != unknown_format && f != detected) {
1089                 PyErr_Format(PyExc_ValueError,
1090                              "can only set %s format to 'unknown' or the "
1091                              "detected platform value", typestr);
1092                 return NULL;
1093         }
1094
1095         *p = f;
1096         Py_RETURN_NONE;
1097 }
1098
1099 PyDoc_STRVAR(float_setformat_doc,
1100 "float.__setformat__(typestr, fmt) -> None\n"
1101 "\n"
1102 "You probably don't want to use this function.  It exists mainly to be\n"
1103 "used in Python's test suite.\n"
1104 "\n"
1105 "typestr must be 'double' or 'float'.  fmt must be one of 'unknown',\n"
1106 "'IEEE, big-endian' or 'IEEE, little-endian', and in addition can only be\n"
1107 "one of the latter two if it appears to match the underlying C reality.\n"
1108 "\n"
1109 "Overrides the automatic determination of C-level floating point type.\n"
1110 "This affects how floats are converted to and from binary strings.");
1111
1112 static PyMethodDef float_methods[] = {
1113         {"__getnewargs__",      (PyCFunction)float_getnewargs,  METH_NOARGS},
1114         {"__getformat__",       (PyCFunction)float_getformat,
1115          METH_O|METH_CLASS,             float_getformat_doc},
1116         {"__setformat__",       (PyCFunction)float_setformat,
1117          METH_VARARGS|METH_CLASS,       float_setformat_doc},
1118         {NULL,          NULL}           /* sentinel */
1119 };
1120
1121 PyDoc_STRVAR(float_doc,
1122 "float(x) -> floating point number\n\
1123 \n\
1124 Convert a string or number to a floating point number, if possible.");
1125
1126
1127 static PyNumberMethods float_as_number = {
1128         (binaryfunc)float_add, /*nb_add*/
1129         (binaryfunc)float_sub, /*nb_subtract*/
1130         (binaryfunc)float_mul, /*nb_multiply*/
1131         (binaryfunc)float_classic_div, /*nb_divide*/
1132         (binaryfunc)float_rem, /*nb_remainder*/
1133         (binaryfunc)float_divmod, /*nb_divmod*/
1134         (ternaryfunc)float_pow, /*nb_power*/
1135         (unaryfunc)float_neg, /*nb_negative*/
1136         (unaryfunc)float_pos, /*nb_positive*/
1137         (unaryfunc)float_abs, /*nb_absolute*/
1138         (inquiry)float_nonzero, /*nb_nonzero*/
1139         0,              /*nb_invert*/
1140         0,              /*nb_lshift*/
1141         0,              /*nb_rshift*/
1142         0,              /*nb_and*/
1143         0,              /*nb_xor*/
1144         0,              /*nb_or*/
1145         (coercion)float_coerce, /*nb_coerce*/
1146         (unaryfunc)float_int, /*nb_int*/
1147         (unaryfunc)float_long, /*nb_long*/
1148         (unaryfunc)float_float, /*nb_float*/
1149         0,              /* nb_oct */
1150         0,              /* nb_hex */
1151         0,              /* nb_inplace_add */
1152         0,              /* nb_inplace_subtract */
1153         0,              /* nb_inplace_multiply */
1154         0,              /* nb_inplace_divide */
1155         0,              /* nb_inplace_remainder */
1156         0,              /* nb_inplace_power */
1157         0,              /* nb_inplace_lshift */
1158         0,              /* nb_inplace_rshift */
1159         0,              /* nb_inplace_and */
1160         0,              /* nb_inplace_xor */
1161         0,              /* nb_inplace_or */
1162         float_floor_div, /* nb_floor_divide */
1163         float_div,      /* nb_true_divide */
1164         0,              /* nb_inplace_floor_divide */
1165         0,              /* nb_inplace_true_divide */
1166 };
1167
1168 PyTypeObject PyFloat_Type = {
1169         PyObject_HEAD_INIT(&PyType_Type)
1170         0,
1171         "float",
1172         sizeof(PyFloatObject),
1173         0,
1174         (destructor)float_dealloc,              /* tp_dealloc */
1175         (printfunc)float_print,                 /* tp_print */
1176         0,                                      /* tp_getattr */
1177         0,                                      /* tp_setattr */
1178         0,                                      /* tp_compare */
1179         (reprfunc)float_repr,                   /* tp_repr */
1180         &float_as_number,                       /* tp_as_number */
1181         0,                                      /* tp_as_sequence */
1182         0,                                      /* tp_as_mapping */
1183         (hashfunc)float_hash,                   /* tp_hash */
1184         0,                                      /* tp_call */
1185         (reprfunc)float_str,                    /* tp_str */
1186         PyObject_GenericGetAttr,                /* tp_getattro */
1187         0,                                      /* tp_setattro */
1188         0,                                      /* tp_as_buffer */
1189         Py_TPFLAGS_DEFAULT | Py_TPFLAGS_CHECKTYPES |
1190                 Py_TPFLAGS_BASETYPE,            /* tp_flags */
1191         float_doc,                              /* tp_doc */
1192         0,                                      /* tp_traverse */
1193         0,                                      /* tp_clear */
1194         (richcmpfunc)float_richcompare,         /* tp_richcompare */
1195         0,                                      /* tp_weaklistoffset */
1196         0,                                      /* tp_iter */
1197         0,                                      /* tp_iternext */
1198         float_methods,                          /* tp_methods */
1199         0,                                      /* tp_members */
1200         0,                                      /* tp_getset */
1201         0,                                      /* tp_base */
1202         0,                                      /* tp_dict */
1203         0,                                      /* tp_descr_get */
1204         0,                                      /* tp_descr_set */
1205         0,                                      /* tp_dictoffset */
1206         0,                                      /* tp_init */
1207         0,                                      /* tp_alloc */
1208         float_new,                              /* tp_new */
1209 };
1210
1211 void
1212 _PyFloat_Init(void)
1213 {
1214         /* We attempt to determine if this machine is using IEEE
1215            floating point formats by peering at the bits of some
1216            carefully chosen values.  If it looks like we are on an
1217            IEEE platform, the float packing/unpacking routines can
1218            just copy bits, if not they resort to arithmetic & shifts
1219            and masks.  The shifts & masks approach works on all finite
1220            values, but what happens to infinities, NaNs and signed
1221            zeroes on packing is an accident, and attempting to unpack
1222            a NaN or an infinity will raise an exception.
1223
1224            Note that if we're on some whacked-out platform which uses
1225            IEEE formats but isn't strictly little-endian or big-
1226            endian, we will fall back to the portable shifts & masks
1227            method. */
1228
1229 #if SIZEOF_DOUBLE == 8
1230         {
1231                 double x = 9006104071832581.0;
1232                 if (memcmp(&x, "\x43\x3f\xff\x01\x02\x03\x04\x05", 8) == 0)
1233                         detected_double_format = ieee_big_endian_format;
1234                 else if (memcmp(&x, "\x05\x04\x03\x02\x01\xff\x3f\x43", 8) == 0)
1235                         detected_double_format = ieee_little_endian_format;
1236                 else
1237                         detected_double_format = unknown_format;
1238         }
1239 #else
1240         detected_double_format = unknown_format;
1241 #endif
1242
1243 #if SIZEOF_FLOAT == 4
1244         {
1245                 float y = 16711938.0;
1246                 if (memcmp(&y, "\x4b\x7f\x01\x02", 4) == 0)
1247                         detected_float_format = ieee_big_endian_format;
1248                 else if (memcmp(&y, "\x02\x01\x7f\x4b", 4) == 0)
1249                         detected_float_format = ieee_little_endian_format;
1250                 else
1251                         detected_float_format = unknown_format;
1252         }
1253 #else
1254         detected_float_format = unknown_format;
1255 #endif
1256
1257         double_format = detected_double_format;
1258         float_format = detected_float_format;
1259 }
1260
1261 void
1262 PyFloat_Fini(void)
1263 {
1264         PyFloatObject *p;
1265         PyFloatBlock *list, *next;
1266         unsigned i;
1267         int bc, bf;     /* block count, number of freed blocks */
1268         int frem, fsum; /* remaining unfreed floats per block, total */
1269
1270         bc = 0;
1271         bf = 0;
1272         fsum = 0;
1273         list = block_list;
1274         block_list = NULL;
1275         free_list = NULL;
1276         while (list != NULL) {
1277                 bc++;
1278                 frem = 0;
1279                 for (i = 0, p = &list->objects[0];
1280                      i < N_FLOATOBJECTS;
1281                      i++, p++) {
1282                         if (PyFloat_CheckExact(p) && p->ob_refcnt != 0)
1283                                 frem++;
1284                 }
1285                 next = list->next;
1286                 if (frem) {
1287                         list->next = block_list;
1288                         block_list = list;
1289                         for (i = 0, p = &list->objects[0];
1290                              i < N_FLOATOBJECTS;
1291                              i++, p++) {
1292                                 if (!PyFloat_CheckExact(p) ||
1293                                     p->ob_refcnt == 0) {
1294                                         p->ob_type = (struct _typeobject *)
1295                                                 free_list;
1296                                         free_list = p;
1297                                 }
1298                         }
1299                 }
1300                 else {
1301                         PyMem_FREE(list); /* XXX PyObject_FREE ??? */
1302                         bf++;
1303                 }
1304                 fsum += frem;
1305                 list = next;
1306         }
1307         if (!Py_VerboseFlag)
1308                 return;
1309         fprintf(stderr, "# cleanup floats");
1310         if (!fsum) {
1311                 fprintf(stderr, "\n");
1312         }
1313         else {
1314                 fprintf(stderr,
1315                         ": %d unfreed float%s in %d out of %d block%s\n",
1316                         fsum, fsum == 1 ? "" : "s",
1317                         bc - bf, bc, bc == 1 ? "" : "s");
1318         }
1319         if (Py_VerboseFlag > 1) {
1320                 list = block_list;
1321                 while (list != NULL) {
1322                         for (i = 0, p = &list->objects[0];
1323                              i < N_FLOATOBJECTS;
1324                              i++, p++) {
1325                                 if (PyFloat_CheckExact(p) &&
1326                                     p->ob_refcnt != 0) {
1327                                         char buf[100];
1328                                         PyFloat_AsString(buf, p);
1329                                         fprintf(stderr,
1330                              "#   <float at %p, refcnt=%d, val=%s>\n",
1331                                                 p, p->ob_refcnt, buf);
1332                                 }
1333                         }
1334                         list = list->next;
1335                 }
1336         }
1337 }
1338
1339 /*----------------------------------------------------------------------------
1340  * _PyFloat_{Pack,Unpack}{4,8}.  See floatobject.h.
1341  *
1342  * TODO:  On platforms that use the standard IEEE-754 single and double
1343  * formats natively, these routines could simply copy the bytes.
1344  */
1345 int
1346 _PyFloat_Pack4(double x, unsigned char *p, int le)
1347 {
1348         if (float_format == unknown_format) {
1349                 unsigned char sign;
1350                 int e;
1351                 double f;
1352                 unsigned int fbits;
1353                 int incr = 1;
1354
1355                 if (le) {
1356                         p += 3;
1357                         incr = -1;
1358                 }
1359
1360                 if (x < 0) {
1361                         sign = 1;
1362                         x = -x;
1363                 }
1364                 else
1365                         sign = 0;
1366
1367                 f = frexp(x, &e);
1368
1369                 /* Normalize f to be in the range [1.0, 2.0) */
1370                 if (0.5 <= f && f < 1.0) {
1371                         f *= 2.0;
1372                         e--;
1373                 }
1374                 else if (f == 0.0)
1375                         e = 0;
1376                 else {
1377                         PyErr_SetString(PyExc_SystemError,
1378                                         "frexp() result out of range");
1379                         return -1;
1380                 }
1381
1382                 if (e >= 128)
1383                         goto Overflow;
1384                 else if (e < -126) {
1385                         /* Gradual underflow */
1386                         f = ldexp(f, 126 + e);
1387                         e = 0;
1388                 }
1389                 else if (!(e == 0 && f == 0.0)) {
1390                         e += 127;
1391                         f -= 1.0; /* Get rid of leading 1 */
1392                 }
1393
1394                 f *= 8388608.0; /* 2**23 */
1395                 fbits = (unsigned int)(f + 0.5); /* Round */
1396                 assert(fbits <= 8388608);
1397                 if (fbits >> 23) {
1398                         /* The carry propagated out of a string of 23 1 bits. */
1399                         fbits = 0;
1400                         ++e;
1401                         if (e >= 255)
1402                                 goto Overflow;
1403                 }
1404
1405                 /* First byte */
1406                 *p = (sign << 7) | (e >> 1);
1407                 p += incr;
1408
1409                 /* Second byte */
1410                 *p = (char) (((e & 1) << 7) | (fbits >> 16));
1411                 p += incr;
1412
1413                 /* Third byte */
1414                 *p = (fbits >> 8) & 0xFF;
1415                 p += incr;
1416
1417                 /* Fourth byte */
1418                 *p = fbits & 0xFF;
1419
1420                 /* Done */
1421                 return 0;
1422
1423           Overflow:
1424                 PyErr_SetString(PyExc_OverflowError,
1425                                 "float too large to pack with f format");
1426                 return -1;
1427         }
1428         else {
1429                 float y = (float)x;
1430                 const char *s = (char*)&y;
1431                 int i, incr = 1;
1432
1433                 if ((float_format == ieee_little_endian_format && !le)
1434                     || (float_format == ieee_big_endian_format && le)) {
1435                         p += 3;
1436                         incr = -1;
1437                 }
1438
1439                 for (i = 0; i < 4; i++) {
1440                         *p = *s++;
1441                         p += incr;
1442                 }
1443                 return 0;
1444         }
1445 }
1446
1447 int
1448 _PyFloat_Pack8(double x, unsigned char *p, int le)
1449 {
1450         if (double_format == unknown_format) {
1451                 unsigned char sign;
1452                 int e;
1453                 double f;
1454                 unsigned int fhi, flo;
1455                 int incr = 1;
1456
1457                 if (le) {
1458                         p += 7;
1459                         incr = -1;
1460                 }
1461
1462                 if (x < 0) {
1463                         sign = 1;
1464                         x = -x;
1465                 }
1466                 else
1467                         sign = 0;
1468
1469                 f = frexp(x, &e);
1470
1471                 /* Normalize f to be in the range [1.0, 2.0) */
1472                 if (0.5 <= f && f < 1.0) {
1473                         f *= 2.0;
1474                         e--;
1475                 }
1476                 else if (f == 0.0)
1477                         e = 0;
1478                 else {
1479                         PyErr_SetString(PyExc_SystemError,
1480                                         "frexp() result out of range");
1481                         return -1;
1482                 }
1483
1484                 if (e >= 1024)
1485                         goto Overflow;
1486                 else if (e < -1022) {
1487                         /* Gradual underflow */
1488                         f = ldexp(f, 1022 + e);
1489                         e = 0;
1490                 }
1491                 else if (!(e == 0 && f == 0.0)) {
1492                         e += 1023;
1493                         f -= 1.0; /* Get rid of leading 1 */
1494                 }
1495
1496                 /* fhi receives the high 28 bits; flo the low 24 bits (== 52 bits) */
1497                 f *= 268435456.0; /* 2**28 */
1498                 fhi = (unsigned int)f; /* Truncate */
1499                 assert(fhi < 268435456);
1500
1501                 f -= (double)fhi;
1502                 f *= 16777216.0; /* 2**24 */
1503                 flo = (unsigned int)(f + 0.5); /* Round */
1504                 assert(flo <= 16777216);
1505                 if (flo >> 24) {
1506                         /* The carry propagated out of a string of 24 1 bits. */
1507                         flo = 0;
1508                         ++fhi;
1509                         if (fhi >> 28) {
1510                                 /* And it also progagated out of the next 28 bits. */
1511                                 fhi = 0;
1512                                 ++e;
1513                                 if (e >= 2047)
1514                                         goto Overflow;
1515                         }
1516                 }
1517
1518                 /* First byte */
1519                 *p = (sign << 7) | (e >> 4);
1520                 p += incr;
1521
1522                 /* Second byte */
1523                 *p = (unsigned char) (((e & 0xF) << 4) | (fhi >> 24));
1524                 p += incr;
1525
1526                 /* Third byte */
1527                 *p = (fhi >> 16) & 0xFF;
1528                 p += incr;
1529
1530                 /* Fourth byte */
1531                 *p = (fhi >> 8) & 0xFF;
1532                 p += incr;
1533
1534                 /* Fifth byte */
1535                 *p = fhi & 0xFF;
1536                 p += incr;
1537
1538                 /* Sixth byte */
1539                 *p = (flo >> 16) & 0xFF;
1540                 p += incr;
1541
1542                 /* Seventh byte */
1543                 *p = (flo >> 8) & 0xFF;
1544                 p += incr;
1545
1546                 /* Eighth byte */
1547                 *p = flo & 0xFF;
1548                 p += incr;
1549
1550                 /* Done */
1551                 return 0;
1552
1553           Overflow:
1554                 PyErr_SetString(PyExc_OverflowError,
1555                                 "float too large to pack with d format");
1556                 return -1;
1557         }
1558         else {
1559                 const char *s = (char*)&x;
1560                 int i, incr = 1;
1561
1562                 if ((double_format == ieee_little_endian_format && !le)
1563                     || (double_format == ieee_big_endian_format && le)) {
1564                         p += 7;
1565                         incr = -1;
1566                 }
1567
1568                 for (i = 0; i < 8; i++) {
1569                         *p = *s++;
1570                         p += incr;
1571                 }
1572                 return 0;
1573         }
1574 }
1575
1576 double
1577 _PyFloat_Unpack4(const unsigned char *p, int le)
1578 {
1579         if (float_format == unknown_format) {
1580                 unsigned char sign;
1581                 int e;
1582                 unsigned int f;
1583                 double x;
1584                 int incr = 1;
1585
1586                 if (le) {
1587                         p += 3;
1588                         incr = -1;
1589                 }
1590
1591                 /* First byte */
1592                 sign = (*p >> 7) & 1;
1593                 e = (*p & 0x7F) << 1;
1594                 p += incr;
1595
1596                 /* Second byte */
1597                 e |= (*p >> 7) & 1;
1598                 f = (*p & 0x7F) << 16;
1599                 p += incr;
1600
1601                 if (e == 255) {
1602                         PyErr_SetString(
1603                                 PyExc_ValueError,
1604                                 "can't unpack IEEE 754 special value "
1605                                 "on non-IEEE platform");
1606                         return -1;
1607                 }
1608
1609                 /* Third byte */
1610                 f |= *p << 8;
1611                 p += incr;
1612
1613                 /* Fourth byte */
1614                 f |= *p;
1615
1616                 x = (double)f / 8388608.0;
1617
1618                 /* XXX This sadly ignores Inf/NaN issues */
1619                 if (e == 0)
1620                         e = -126;
1621                 else {
1622                         x += 1.0;
1623                         e -= 127;
1624                 }
1625                 x = ldexp(x, e);
1626
1627                 if (sign)
1628                         x = -x;
1629
1630                 return x;
1631         }
1632         else {
1633                 float x;
1634
1635                 if ((float_format == ieee_little_endian_format && !le)
1636                     || (float_format == ieee_big_endian_format && le)) {
1637                         char buf[4];
1638                         char *d = &buf[3];
1639                         int i;
1640
1641                         for (i = 0; i < 4; i++) {
1642                                 *d-- = *p++;
1643                         }
1644                         memcpy(&x, buf, 4);
1645                 }
1646                 else {
1647                         memcpy(&x, p, 4);
1648                 }
1649
1650                 return x;
1651         }
1652 }
1653
1654 double
1655 _PyFloat_Unpack8(const unsigned char *p, int le)
1656 {
1657         if (double_format == unknown_format) {
1658                 unsigned char sign;
1659                 int e;
1660                 unsigned int fhi, flo;
1661                 double x;
1662                 int incr = 1;
1663
1664                 if (le) {
1665                         p += 7;
1666                         incr = -1;
1667                 }
1668
1669                 /* First byte */
1670                 sign = (*p >> 7) & 1;
1671                 e = (*p & 0x7F) << 4;
1672
1673                 p += incr;
1674
1675                 /* Second byte */
1676                 e |= (*p >> 4) & 0xF;
1677                 fhi = (*p & 0xF) << 24;
1678                 p += incr;
1679
1680                 if (e == 2047) {
1681                         PyErr_SetString(
1682                                 PyExc_ValueError,
1683                                 "can't unpack IEEE 754 special value "
1684                                 "on non-IEEE platform");
1685                         return -1.0;
1686                 }
1687
1688                 /* Third byte */
1689                 fhi |= *p << 16;
1690                 p += incr;
1691
1692                 /* Fourth byte */
1693                 fhi |= *p  << 8;
1694                 p += incr;
1695
1696                 /* Fifth byte */
1697                 fhi |= *p;
1698                 p += incr;
1699
1700                 /* Sixth byte */
1701                 flo = *p << 16;
1702                 p += incr;
1703
1704                 /* Seventh byte */
1705                 flo |= *p << 8;
1706                 p += incr;
1707
1708                 /* Eighth byte */
1709                 flo |= *p;
1710
1711                 x = (double)fhi + (double)flo / 16777216.0; /* 2**24 */
1712                 x /= 268435456.0; /* 2**28 */
1713
1714                 if (e == 0)
1715                         e = -1022;
1716                 else {
1717                         x += 1.0;
1718                         e -= 1023;
1719                 }
1720                 x = ldexp(x, e);
1721
1722                 if (sign)
1723                         x = -x;
1724
1725                 return x;
1726         }
1727         else {
1728                 double x;
1729
1730                 if ((double_format == ieee_little_endian_format && !le)
1731                     || (double_format == ieee_big_endian_format && le)) {
1732                         char buf[8];
1733                         char *d = &buf[7];
1734                         int i;
1735
1736                         for (i = 0; i < 8; i++) {
1737                                 *d-- = *p++;
1738                         }
1739                         memcpy(&x, buf, 8);
1740                 }
1741                 else {
1742                         memcpy(&x, p, 8);
1743                 }
1744
1745                 return x;
1746         }
1747 }