Objects/floatobject.c

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