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         Py_ssize_t 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) >= (Py_ssize_t)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 = 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_FINITE(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_FINITE(iw)) {
 806                         /* Return 1 if iw is even, -1 if iw is odd; there's
 807                          * no guarantee that any C integral type is big
 808                          * enough to hold iw, so we have to check this
 809                          * indirectly.
 810                          */
 811                         ix = floor(iw * 0.5) * 2.0;
 812                         return PyFloat_FromDouble(ix == iw ? 1.0 : -1.0);
 813                 }
 814                 /* Else iv != -1.0, and overflow or underflow are possible.
 815                  * Unless we're to write pow() ourselves, we have to trust
 816                  * the platform to do this correctly.
 817                  */
 818         }
 819         errno = 0;
 820         PyFPE_START_PROTECT("pow", return NULL)
 821         ix = pow(iv, iw);
 822         PyFPE_END_PROTECT(ix)
 823         Py_ADJUST_ERANGE1(ix);
 824         if (errno != 0) {
 825                 /* We don't expect any errno value other than ERANGE, but
 826                  * the range of libm bugs appears unbounded.
 827                  */
 828                 PyErr_SetFromErrno(errno == ERANGE ? PyExc_OverflowError :
 829                                                      PyExc_ValueError);
 830                 return NULL;
 831         }
 832         return PyFloat_FromDouble(ix);
 833 }
 834
 835 static PyObject *
 836 float_neg(PyFloatObject *v)
 837 {
 838         return PyFloat_FromDouble(-v->ob_fval);
 839 }
 840
 841 static PyObject *
 842 float_pos(PyFloatObject *v)
 843 {
 844         if (PyFloat_CheckExact(v)) {
 845                 Py_INCREF(v);
 846                 return (PyObject *)v;
 847         }
 848         else
 849                 return PyFloat_FromDouble(v->ob_fval);
 850 }
 851
 852 static PyObject *
 853 float_abs(PyFloatObject *v)
 854 {
 855         return PyFloat_FromDouble(fabs(v->ob_fval));
 856 }
 857
 858 static int
 859 float_nonzero(PyFloatObject *v)
 860 {
 861         return v->ob_fval != 0.0;
 862 }
 863
 864 static int
 865 float_coerce(PyObject **pv, PyObject **pw)
 866 {
 867         if (PyInt_Check(*pw)) {
 868                 long x = PyInt_AsLong(*pw);
 869                 *pw = PyFloat_FromDouble((double)x);
 870                 Py_INCREF(*pv);
 871                 return 0;
 872         }
 873         else if (PyLong_Check(*pw)) {
 874                 double x = PyLong_AsDouble(*pw);
 875                 if (x == -1.0 && PyErr_Occurred())
 876                         return -1;
 877                 *pw = PyFloat_FromDouble(x);
 878                 Py_INCREF(*pv);
 879                 return 0;
 880         }
 881         else if (PyFloat_Check(*pw)) {
 882                 Py_INCREF(*pv);
 883                 Py_INCREF(*pw);
 884                 return 0;
 885         }
 886         return 1; /* Can't do it */
 887 }
 888
 889 static PyObject *
 890 float_long(PyObject *v)
 891 {
 892         double x = PyFloat_AsDouble(v);
 893         return PyLong_FromDouble(x);
 894 }
 895
 896 static PyObject *
 897 float_int(PyObject *v)
 898 {
 899         double x = PyFloat_AsDouble(v);
 900         double wholepart;       /* integral portion of x, rounded toward 0 */
 901
 902         (void)modf(x, &wholepart);
 903         /* Try to get out cheap if this fits in a Python int.  The attempt
 904          * to cast to long must be protected, as C doesn't define what
 905          * happens if the double is too big to fit in a long.  Some rare
 906          * systems raise an exception then (RISCOS was mentioned as one,
 907          * and someone using a non-default option on Sun also bumped into
 908          * that).  Note that checking for >= and <= LONG_{MIN,MAX} would
 909          * still be vulnerable:  if a long has more bits of precision than
 910          * a double, casting MIN/MAX to double may yield an approximation,
 911          * and if that's rounded up, then, e.g., wholepart=LONG_MAX+1 would
 912          * yield true from the C expression wholepart<=LONG_MAX, despite
 913          * that wholepart is actually greater than LONG_MAX.
 914          */
 915         if (LONG_MIN < wholepart && wholepart < LONG_MAX) {
 916                 const long aslong = (long)wholepart;
 917                 return PyInt_FromLong(aslong);
 918         }
 919         return PyLong_FromDouble(wholepart);
 920 }
 921
 922 static PyObject *
 923 float_float(PyObject *v)
 924 {
 925         if (PyFloat_CheckExact(v))
 926                 Py_INCREF(v);
 927         else
 928                 v = PyFloat_FromDouble(((PyFloatObject *)v)->ob_fval);
 929         return v;
 930 }
 931
 932
 933 static PyObject *
 934 float_subtype_new(PyTypeObject *type, PyObject *args, PyObject *kwds);
 935
 936 static PyObject *
 937 float_new(PyTypeObject *type, PyObject *args, PyObject *kwds)
 938 {
 939         PyObject *x = Py_False; /* Integer zero */
 940         static char *kwlist[] = {"x", 0};
 941
 942         if (type != &PyFloat_Type)
 943                 return float_subtype_new(type, args, kwds); /* Wimp out */
 944         if (!PyArg_ParseTupleAndKeywords(args, kwds, "|O:float", kwlist, &x))
 945                 return NULL;
 946         if (PyString_Check(x))
 947                 return PyFloat_FromString(x, NULL);
 948         return PyNumber_Float(x);
 949 }
 950
 951 /* Wimpy, slow approach to tp_new calls for subtypes of float:
 952    first create a regular float from whatever arguments we got,
 953    then allocate a subtype instance and initialize its ob_fval
 954    from the regular float.  The regular float is then thrown away.
 955 */
 956 static PyObject *
 957 float_subtype_new(PyTypeObject *type, PyObject *args, PyObject *kwds)
 958 {
 959         PyObject *tmp, *newobj;
 960
 961         assert(PyType_IsSubtype(type, &PyFloat_Type));
 962         tmp = float_new(&PyFloat_Type, args, kwds);
 963         if (tmp == NULL)
 964                 return NULL;
 965         assert(PyFloat_CheckExact(tmp));
 966         newobj = type->tp_alloc(type, 0);
 967         if (newobj == NULL) {
 968                 Py_DECREF(tmp);
 969                 return NULL;
 970         }
 971         ((PyFloatObject *)newobj)->ob_fval = ((PyFloatObject *)tmp)->ob_fval;
 972         Py_DECREF(tmp);
 973         return newobj;
 974 }
 975
 976 static PyObject *
 977 float_getnewargs(PyFloatObject *v)
 978 {
 979         return Py_BuildValue("(d)", v->ob_fval);
 980 }
 981
 982 /* this is for the benefit of the pack/unpack routines below */
 983
 984 typedef enum {
 985         unknown_format, ieee_big_endian_format, ieee_little_endian_format
 986 } float_format_type;
 987
 988 static float_format_type double_format, float_format;
 989 static float_format_type detected_double_format, detected_float_format;
 990
 991 static PyObject *
 992 float_getformat(PyTypeObject *v, PyObject* arg)
 993 {
 994         char* s;
 995         float_format_type r;
 996
 997         if (!PyString_Check(arg)) {
 998                 PyErr_Format(PyExc_TypeError,
 999              "__getformat__() argument must be string, not %.500s",
1000                              arg->ob_type->tp_name);
1001                 return NULL;
1002         }
1003         s = PyString_AS_STRING(arg);
1004         if (strcmp(s, "double") == 0) {
1005                 r = double_format;
1006         }
1007         else if (strcmp(s, "float") == 0) {
1008                 r = float_format;
1009         }
1010         else {
1011                 PyErr_SetString(PyExc_ValueError,
1012                                 "__getformat__() argument 1 must be "
1013                                 "'double' or 'float'");
1014                 return NULL;
1015         }
1016
1017         switch (r) {
1018         case unknown_format:
1019                 return PyString_FromString("unknown");
1020         case ieee_little_endian_format:
1021                 return PyString_FromString("IEEE, little-endian");
1022         case ieee_big_endian_format:
1023                 return PyString_FromString("IEEE, big-endian");
1024         default:
1025                 Py_FatalError("insane float_format or double_format");
1026                 return NULL;
1027         }
1028 }
1029
1030 PyDoc_STRVAR(float_getformat_doc,
1031 "float.__getformat__(typestr) -> string\n"
1032 "\n"
1033 "You probably don't want to use this function.  It exists mainly to be\n"
1034 "used in Python's test suite.\n"
1035 "\n"
1036 "typestr must be 'double' or 'float'.  This function returns whichever of\n"
1037 "'unknown', 'IEEE, big-endian' or 'IEEE, little-endian' best describes the\n"
1038 "format of floating point numbers used by the C type named by typestr.");
1039
1040 static PyObject *
1041 float_setformat(PyTypeObject *v, PyObject* args)
1042 {
1043         char* typestr;
1044         char* format;
1045         float_format_type f;
1046         float_format_type detected;
1047         float_format_type *p;
1048
1049         if (!PyArg_ParseTuple(args, "ss:__setformat__", &typestr, &format))
1050                 return NULL;
1051
1052         if (strcmp(typestr, "double") == 0) {
1053                 p = &double_format;
1054                 detected = detected_double_format;
1055         }
1056         else if (strcmp(typestr, "float") == 0) {
1057                 p = &float_format;
1058                 detected = detected_float_format;
1059         }
1060         else {
1061                 PyErr_SetString(PyExc_ValueError,
1062                                 "__setformat__() argument 1 must "
1063                                 "be 'double' or 'float'");
1064                 return NULL;
1065         }
1066
1067         if (strcmp(format, "unknown") == 0) {
1068                 f = unknown_format;
1069         }
1070         else if (strcmp(format, "IEEE, little-endian") == 0) {
1071                 f = ieee_little_endian_format;
1072         }
1073         else if (strcmp(format, "IEEE, big-endian") == 0) {
1074                 f = ieee_big_endian_format;
1075         }
1076         else {
1077                 PyErr_SetString(PyExc_ValueError,
1078                                 "__setformat__() argument 2 must be "
1079                                 "'unknown', 'IEEE, little-endian' or "
1080                                 "'IEEE, big-endian'");
1081                 return NULL;
1082
1083         }
1084
1085         if (f != unknown_format && f != detected) {
1086                 PyErr_Format(PyExc_ValueError,
1087                              "can only set %s format to 'unknown' or the "
1088                              "detected platform value", typestr);
1089                 return NULL;
1090         }
1091
1092         *p = f;
1093         Py_RETURN_NONE;
1094 }
1095
1096 PyDoc_STRVAR(float_setformat_doc,
1097 "float.__setformat__(typestr, fmt) -> None\n"
1098 "\n"
1099 "You probably don't want to use this function.  It exists mainly to be\n"
1100 "used in Python's test suite.\n"
1101 "\n"
1102 "typestr must be 'double' or 'float'.  fmt must be one of 'unknown',\n"
1103 "'IEEE, big-endian' or 'IEEE, little-endian', and in addition can only be\n"
1104 "one of the latter two if it appears to match the underlying C reality.\n"
1105 "\n"
1106 "Overrides the automatic determination of C-level floating point type.\n"
1107 "This affects how floats are converted to and from binary strings.");
1108
1109 static PyMethodDef float_methods[] = {
1110         {"__getnewargs__",      (PyCFunction)float_getnewargs,  METH_NOARGS},
1111         {"__getformat__",       (PyCFunction)float_getformat,
1112          METH_O|METH_CLASS,             float_getformat_doc},
1113         {"__setformat__",       (PyCFunction)float_setformat,
1114          METH_VARARGS|METH_CLASS,       float_setformat_doc},
1115         {NULL,          NULL}           /* sentinel */
1116 };
1117
1118 PyDoc_STRVAR(float_doc,
1119 "float(x) -> floating point number\n\
1120 \n\
1121 Convert a string or number to a floating point number, if possible.");
1122
1123
1124 static PyNumberMethods float_as_number = {
1125         float_add,      /*nb_add*/
1126         float_sub,      /*nb_subtract*/
1127         float_mul,      /*nb_multiply*/
1128         float_classic_div, /*nb_divide*/
1129         float_rem,      /*nb_remainder*/
1130         float_divmod,   /*nb_divmod*/
1131         float_pow,      /*nb_power*/
1132         (unaryfunc)float_neg, /*nb_negative*/
1133         (unaryfunc)float_pos, /*nb_positive*/
1134         (unaryfunc)float_abs, /*nb_absolute*/
1135         (inquiry)float_nonzero, /*nb_nonzero*/
1136         0,              /*nb_invert*/
1137         0,              /*nb_lshift*/
1138         0,              /*nb_rshift*/
1139         0,              /*nb_and*/
1140         0,              /*nb_xor*/
1141         0,              /*nb_or*/
1142         float_coerce,   /*nb_coerce*/
1143         float_int,      /*nb_int*/
1144         float_long,     /*nb_long*/
1145         float_float,    /*nb_float*/
1146         0,              /* nb_oct */
1147         0,              /* nb_hex */
1148         0,              /* nb_inplace_add */
1149         0,              /* nb_inplace_subtract */
1150         0,              /* nb_inplace_multiply */
1151         0,              /* nb_inplace_divide */
1152         0,              /* nb_inplace_remainder */
1153         0,              /* nb_inplace_power */
1154         0,              /* nb_inplace_lshift */
1155         0,              /* nb_inplace_rshift */
1156         0,              /* nb_inplace_and */
1157         0,              /* nb_inplace_xor */
1158         0,              /* nb_inplace_or */
1159         float_floor_div, /* nb_floor_divide */
1160         float_div,      /* nb_true_divide */
1161         0,              /* nb_inplace_floor_divide */
1162         0,              /* nb_inplace_true_divide */
1163 };
1164
1165 PyTypeObject PyFloat_Type = {
1166         PyObject_HEAD_INIT(&PyType_Type)
1167         0,
1168         "float",
1169         sizeof(PyFloatObject),
1170         0,
1171         (destructor)float_dealloc,              /* tp_dealloc */
1172         (printfunc)float_print,                 /* tp_print */
1173         0,                                      /* tp_getattr */
1174         0,                                      /* tp_setattr */
1175         0,                                      /* tp_compare */
1176         (reprfunc)float_repr,                   /* tp_repr */
1177         &float_as_number,                       /* tp_as_number */
1178         0,                                      /* tp_as_sequence */
1179         0,                                      /* tp_as_mapping */
1180         (hashfunc)float_hash,                   /* tp_hash */
1181         0,                                      /* tp_call */
1182         (reprfunc)float_str,                    /* tp_str */
1183         PyObject_GenericGetAttr,                /* tp_getattro */
1184         0,                                      /* tp_setattro */
1185         0,                                      /* tp_as_buffer */
1186         Py_TPFLAGS_DEFAULT | Py_TPFLAGS_CHECKTYPES |
1187                 Py_TPFLAGS_BASETYPE,            /* tp_flags */
1188         float_doc,                              /* tp_doc */
1189         0,                                      /* tp_traverse */
1190         0,                                      /* tp_clear */
1191         float_richcompare,                      /* tp_richcompare */
1192         0,                                      /* tp_weaklistoffset */
1193         0,                                      /* tp_iter */
1194         0,                                      /* tp_iternext */
1195         float_methods,                          /* tp_methods */
1196         0,                                      /* tp_members */
1197         0,                                      /* tp_getset */
1198         0,                                      /* tp_base */
1199         0,                                      /* tp_dict */
1200         0,                                      /* tp_descr_get */
1201         0,                                      /* tp_descr_set */
1202         0,                                      /* tp_dictoffset */
1203         0,                                      /* tp_init */
1204         0,                                      /* tp_alloc */
1205         float_new,                              /* tp_new */
1206 };
1207
1208 void
1209 _PyFloat_Init(void)
1210 {
1211         /* We attempt to determine if this machine is using IEEE
1212            floating point formats by peering at the bits of some
1213            carefully chosen values.  If it looks like we are on an
1214            IEEE platform, the float packing/unpacking routines can
1215            just copy bits, if not they resort to arithmetic & shifts
1216            and masks.  The shifts & masks approach works on all finite
1217            values, but what happens to infinities, NaNs and signed
1218            zeroes on packing is an accident, and attempting to unpack
1219            a NaN or an infinity will raise an exception.
1220
1221            Note that if we're on some whacked-out platform which uses
1222            IEEE formats but isn't strictly little-endian or big-
1223            endian, we will fall back to the portable shifts & masks
1224            method. */
1225
1226 #if SIZEOF_DOUBLE == 8
1227         {
1228                 double x = 9006104071832581.0;
1229                 if (memcmp(&x, "\x43\x3f\xff\x01\x02\x03\x04\x05", 8) == 0)
1230                         detected_double_format = ieee_big_endian_format;
1231                 else if (memcmp(&x, "\x05\x04\x03\x02\x01\xff\x3f\x43", 8) == 0)
1232                         detected_double_format = ieee_little_endian_format;
1233                 else
1234                         detected_double_format = unknown_format;
1235         }
1236 #else
1237         detected_double_format = unknown_format;
1238 #endif
1239
1240 #if SIZEOF_FLOAT == 4
1241         {
1242                 float y = 16711938.0;
1243                 if (memcmp(&y, "\x4b\x7f\x01\x02", 4) == 0)
1244                         detected_float_format = ieee_big_endian_format;
1245                 else if (memcmp(&y, "\x02\x01\x7f\x4b", 4) == 0)
1246                         detected_float_format = ieee_little_endian_format;
1247                 else
1248                         detected_float_format = unknown_format;
1249         }
1250 #else
1251         detected_float_format = unknown_format;
1252 #endif
1253
1254         double_format = detected_double_format;
1255         float_format = detected_float_format;
1256 }
1257
1258 void
1259 PyFloat_Fini(void)
1260 {
1261         PyFloatObject *p;
1262         PyFloatBlock *list, *next;
1263         unsigned i;
1264         int bc, bf;     /* block count, number of freed blocks */
1265         int frem, fsum; /* remaining unfreed floats per block, total */
1266
1267         bc = 0;
1268         bf = 0;
1269         fsum = 0;
1270         list = block_list;
1271         block_list = NULL;
1272         free_list = NULL;
1273         while (list != NULL) {
1274                 bc++;
1275                 frem = 0;
1276                 for (i = 0, p = &list->objects[0];
1277                      i < N_FLOATOBJECTS;
1278                      i++, p++) {
1279                         if (PyFloat_CheckExact(p) && p->ob_refcnt != 0)
1280                                 frem++;
1281                 }
1282                 next = list->next;
1283                 if (frem) {
1284                         list->next = block_list;
1285                         block_list = list;
1286                         for (i = 0, p = &list->objects[0];
1287                              i < N_FLOATOBJECTS;
1288                              i++, p++) {
1289                                 if (!PyFloat_CheckExact(p) ||
1290                                     p->ob_refcnt == 0) {
1291                                         p->ob_type = (struct _typeobject *)
1292                                                 free_list;
1293                                         free_list = p;
1294                                 }
1295                         }
1296                 }
1297                 else {
1298                         PyMem_FREE(list); /* XXX PyObject_FREE ??? */
1299                         bf++;
1300                 }
1301                 fsum += frem;
1302                 list = next;
1303         }
1304         if (!Py_VerboseFlag)
1305                 return;
1306         fprintf(stderr, "# cleanup floats");
1307         if (!fsum) {
1308                 fprintf(stderr, "\n");
1309         }
1310         else {
1311                 fprintf(stderr,
1312                         ": %d unfreed float%s in %d out of %d block%s\n",
1313                         fsum, fsum == 1 ? "" : "s",
1314                         bc - bf, bc, bc == 1 ? "" : "s");
1315         }
1316         if (Py_VerboseFlag > 1) {
1317                 list = block_list;
1318                 while (list != NULL) {
1319                         for (i = 0, p = &list->objects[0];
1320                              i < N_FLOATOBJECTS;
1321                              i++, p++) {
1322                                 if (PyFloat_CheckExact(p) &&
1323                                     p->ob_refcnt != 0) {
1324                                         char buf[100];
1325                                         PyFloat_AsString(buf, p);
1326                                         /* XXX(twouters) cast refcount to
1327                                            long until %zd is universally
1328                                            available
1329                                          */
1330                                         fprintf(stderr,
1331                              "#   <float at %p, refcnt=%ld, val=%s>\n",
1332                                                 p, (long)p->ob_refcnt, buf);
1333                                 }
1334                         }
1335                         list = list->next;
1336                 }
1337         }
1338 }
1339
1340 /*----------------------------------------------------------------------------
1341  * _PyFloat_{Pack,Unpack}{4,8}.  See floatobject.h.
1342  *
1343  * TODO:  On platforms that use the standard IEEE-754 single and double
1344  * formats natively, these routines could simply copy the bytes.
1345  */
1346 int
1347 _PyFloat_Pack4(double x, unsigned char *p, int le)
1348 {
1349         if (float_format == unknown_format) {
1350                 unsigned char sign;
1351                 int e;
1352                 double f;
1353                 unsigned int fbits;
1354                 int incr = 1;
1355
1356                 if (le) {
1357                         p += 3;
1358                         incr = -1;
1359                 }
1360
1361                 if (x < 0) {
1362                         sign = 1;
1363                         x = -x;
1364                 }
1365                 else
1366                         sign = 0;
1367
1368                 f = frexp(x, &e);
1369
1370                 /* Normalize f to be in the range [1.0, 2.0) */
1371                 if (0.5 <= f && f < 1.0) {
1372                         f *= 2.0;
1373                         e--;
1374                 }
1375                 else if (f == 0.0)
1376                         e = 0;
1377                 else {
1378                         PyErr_SetString(PyExc_SystemError,
1379                                         "frexp() result out of range");
1380                         return -1;
1381                 }
1382
1383                 if (e >= 128)
1384                         goto Overflow;
1385                 else if (e < -126) {
1386                         /* Gradual underflow */
1387                         f = ldexp(f, 126 + e);
1388                         e = 0;
1389                 }
1390                 else if (!(e == 0 && f == 0.0)) {
1391                         e += 127;
1392                         f -= 1.0; /* Get rid of leading 1 */
1393                 }
1394
1395                 f *= 8388608.0; /* 2**23 */
1396                 fbits = (unsigned int)(f + 0.5); /* Round */
1397                 assert(fbits <= 8388608);
1398                 if (fbits >> 23) {
1399                         /* The carry propagated out of a string of 23 1 bits. */
1400                         fbits = 0;
1401                         ++e;
1402                         if (e >= 255)
1403                                 goto Overflow;
1404                 }
1405
1406                 /* First byte */
1407                 *p = (sign << 7) | (e >> 1);
1408                 p += incr;
1409
1410                 /* Second byte */
1411                 *p = (char) (((e & 1) << 7) | (fbits >> 16));
1412                 p += incr;
1413
1414                 /* Third byte */
1415                 *p = (fbits >> 8) & 0xFF;
1416                 p += incr;
1417
1418                 /* Fourth byte */
1419                 *p = fbits & 0xFF;
1420
1421                 /* Done */
1422                 return 0;
1423
1424           Overflow:
1425                 PyErr_SetString(PyExc_OverflowError,
1426                                 "float too large to pack with f format");
1427                 return -1;
1428         }
1429         else {
1430                 float y = (float)x;
1431                 const char *s = (char*)&y;
1432                 int i, incr = 1;
1433
1434                 if ((float_format == ieee_little_endian_format && !le)
1435                     || (float_format == ieee_big_endian_format && le)) {
1436                         p += 3;
1437                         incr = -1;
1438                 }
1439
1440                 for (i = 0; i < 4; i++) {
1441                         *p = *s++;
1442                         p += incr;
1443                 }
1444                 return 0;
1445         }
1446 }
1447
1448 int
1449 _PyFloat_Pack8(double x, unsigned char *p, int le)
1450 {
1451         if (double_format == unknown_format) {
1452                 unsigned char sign;
1453                 int e;
1454                 double f;
1455                 unsigned int fhi, flo;
1456                 int incr = 1;
1457
1458                 if (le) {
1459                         p += 7;
1460                         incr = -1;
1461                 }
1462
1463                 if (x < 0) {
1464                         sign = 1;
1465                         x = -x;
1466                 }
1467                 else
1468                         sign = 0;
1469
1470                 f = frexp(x, &e);
1471
1472                 /* Normalize f to be in the range [1.0, 2.0) */
1473                 if (0.5 <= f && f < 1.0) {
1474                         f *= 2.0;
1475                         e--;
1476                 }
1477                 else if (f == 0.0)
1478                         e = 0;
1479                 else {
1480                         PyErr_SetString(PyExc_SystemError,
1481                                         "frexp() result out of range");
1482                         return -1;
1483                 }
1484
1485                 if (e >= 1024)
1486                         goto Overflow;
1487                 else if (e < -1022) {
1488                         /* Gradual underflow */
1489                         f = ldexp(f, 1022 + e);
1490                         e = 0;
1491                 }
1492                 else if (!(e == 0 && f == 0.0)) {
1493                         e += 1023;
1494                         f -= 1.0; /* Get rid of leading 1 */
1495                 }
1496
1497                 /* fhi receives the high 28 bits; flo the low 24 bits (== 52 bits) */
1498                 f *= 268435456.0; /* 2**28 */
1499                 fhi = (unsigned int)f; /* Truncate */
1500                 assert(fhi < 268435456);
1501
1502                 f -= (double)fhi;
1503                 f *= 16777216.0; /* 2**24 */
1504                 flo = (unsigned int)(f + 0.5); /* Round */
1505                 assert(flo <= 16777216);
1506                 if (flo >> 24) {
1507                         /* The carry propagated out of a string of 24 1 bits. */
1508                         flo = 0;
1509                         ++fhi;
1510                         if (fhi >> 28) {
1511                                 /* And it also progagated out of the next 28 bits. */
1512                                 fhi = 0;
1513                                 ++e;
1514                                 if (e >= 2047)
1515                                         goto Overflow;
1516                         }
1517                 }
1518
1519                 /* First byte */
1520                 *p = (sign << 7) | (e >> 4);
1521                 p += incr;
1522
1523                 /* Second byte */
1524                 *p = (unsigned char) (((e & 0xF) << 4) | (fhi >> 24));
1525                 p += incr;
1526
1527                 /* Third byte */
1528                 *p = (fhi >> 16) & 0xFF;
1529                 p += incr;
1530
1531                 /* Fourth byte */
1532                 *p = (fhi >> 8) & 0xFF;
1533                 p += incr;
1534
1535                 /* Fifth byte */
1536                 *p = fhi & 0xFF;
1537                 p += incr;
1538
1539                 /* Sixth byte */
1540                 *p = (flo >> 16) & 0xFF;
1541                 p += incr;
1542
1543                 /* Seventh byte */
1544                 *p = (flo >> 8) & 0xFF;
1545                 p += incr;
1546
1547                 /* Eighth byte */
1548                 *p = flo & 0xFF;
1549                 p += incr;
1550
1551                 /* Done */
1552                 return 0;
1553
1554           Overflow:
1555                 PyErr_SetString(PyExc_OverflowError,
1556                                 "float too large to pack with d format");
1557                 return -1;
1558         }
1559         else {
1560                 const char *s = (char*)&x;
1561                 int i, incr = 1;
1562
1563                 if ((double_format == ieee_little_endian_format && !le)
1564                     || (double_format == ieee_big_endian_format && le)) {
1565                         p += 7;
1566                         incr = -1;
1567                 }
1568
1569                 for (i = 0; i < 8; i++) {
1570                         *p = *s++;
1571                         p += incr;
1572                 }
1573                 return 0;
1574         }
1575 }
1576
1577 double
1578 _PyFloat_Unpack4(const unsigned char *p, int le)
1579 {
1580         if (float_format == unknown_format) {
1581                 unsigned char sign;
1582                 int e;
1583                 unsigned int f;
1584                 double x;
1585                 int incr = 1;
1586
1587                 if (le) {
1588                         p += 3;
1589                         incr = -1;
1590                 }
1591
1592                 /* First byte */
1593                 sign = (*p >> 7) & 1;
1594                 e = (*p & 0x7F) << 1;
1595                 p += incr;
1596
1597                 /* Second byte */
1598                 e |= (*p >> 7) & 1;
1599                 f = (*p & 0x7F) << 16;
1600                 p += incr;
1601
1602                 if (e == 255) {
1603                         PyErr_SetString(
1604                                 PyExc_ValueError,
1605                                 "can't unpack IEEE 754 special value "
1606                                 "on non-IEEE platform");
1607                         return -1;
1608                 }
1609
1610                 /* Third byte */
1611                 f |= *p << 8;
1612                 p += incr;
1613
1614                 /* Fourth byte */
1615                 f |= *p;
1616
1617                 x = (double)f / 8388608.0;
1618
1619                 /* XXX This sadly ignores Inf/NaN issues */
1620                 if (e == 0)
1621                         e = -126;
1622                 else {
1623                         x += 1.0;
1624                         e -= 127;
1625                 }
1626                 x = ldexp(x, e);
1627
1628                 if (sign)
1629                         x = -x;
1630
1631                 return x;
1632         }
1633         else {
1634                 float x;
1635
1636                 if ((float_format == ieee_little_endian_format && !le)
1637                     || (float_format == ieee_big_endian_format && le)) {
1638                         char buf[4];
1639                         char *d = &buf[3];
1640                         int i;
1641
1642                         for (i = 0; i < 4; i++) {
1643                                 *d-- = *p++;
1644                         }
1645                         memcpy(&x, buf, 4);
1646                 }
1647                 else {
1648                         memcpy(&x, p, 4);
1649                 }
1650
1651                 return x;
1652         }
1653 }
1654
1655 double
1656 _PyFloat_Unpack8(const unsigned char *p, int le)
1657 {
1658         if (double_format == unknown_format) {
1659                 unsigned char sign;
1660                 int e;
1661                 unsigned int fhi, flo;
1662                 double x;
1663                 int incr = 1;
1664
1665                 if (le) {
1666                         p += 7;
1667                         incr = -1;
1668                 }
1669
1670                 /* First byte */
1671                 sign = (*p >> 7) & 1;
1672                 e = (*p & 0x7F) << 4;
1673
1674                 p += incr;
1675
1676                 /* Second byte */
1677                 e |= (*p >> 4) & 0xF;
1678                 fhi = (*p & 0xF) << 24;
1679                 p += incr;
1680
1681                 if (e == 2047) {
1682                         PyErr_SetString(
1683                                 PyExc_ValueError,
1684                                 "can't unpack IEEE 754 special value "
1685                                 "on non-IEEE platform");
1686                         return -1.0;
1687                 }
1688
1689                 /* Third byte */
1690                 fhi |= *p << 16;
1691                 p += incr;
1692
1693                 /* Fourth byte */
1694                 fhi |= *p  << 8;
1695                 p += incr;
1696
1697                 /* Fifth byte */
1698                 fhi |= *p;
1699                 p += incr;
1700
1701                 /* Sixth byte */
1702                 flo = *p << 16;
1703                 p += incr;
1704
1705                 /* Seventh byte */
1706                 flo |= *p << 8;
1707                 p += incr;
1708
1709                 /* Eighth byte */
1710                 flo |= *p;
1711
1712                 x = (double)fhi + (double)flo / 16777216.0; /* 2**24 */
1713                 x /= 268435456.0; /* 2**28 */
1714
1715                 if (e == 0)
1716                         e = -1022;
1717                 else {
1718                         x += 1.0;
1719                         e -= 1023;
1720                 }
1721                 x = ldexp(x, e);
1722
1723                 if (sign)
1724                         x = -x;
1725
1726                 return x;
1727         }
1728         else {
1729                 double x;
1730
1731                 if ((double_format == ieee_little_endian_format && !le)
1732                     || (double_format == ieee_big_endian_format && le)) {
1733                         char buf[8];
1734                         char *d = &buf[7];
1735                         int i;
1736
1737                         for (i = 0; i < 8; i++) {
1738                                 *d-- = *p++;
1739                         }
1740                         memcpy(&x, buf, 8);
1741                 }
1742                 else {
1743                         memcpy(&x, p, 8);
1744                 }
1745
1746                 return x;
1747         }
1748 }